Heavy-quark form factors in the large \(\beta _0\) limit

Eur. Phys. J. C Heavy-quark form factors in the large β0 limit Andrey G. Grozin 0 1 2 3 0 THEP, Institut für Physik, Universität Mainz , Mainz , Germany 1 Institut für Theoretische Teilchenphysik, Karlsruher Institut für Technologie , Karlsruhe , Germany 2 Budker Institute of Nuclear Physics , Novosibirsk , Russia 3 Novosibirsk State University , Novosibirsk , Russia Heavy-quark form factors are calculated at β0αs ∼ 1 to all orders in αs at the first order in 1/β0. Using the inversion relation generalized to vertex functions, we reduce the massive on-shell Feynman integral to the HQET one. This HQET vertex integral can be expressed via a 2 F1 function; the nth term of its ε expansion is explicitly known. We confirm existing results for nlL−1αsL terms in the form factors (up to L = 3), and we present results for higher L. 1 Introduction Quark form factors are building blocks for various production cross sections and decay widths in QCD. Massive quark form factors are known up to two loops [ 1 ]; recently they have been calculated at three loops in the large Nc limit [ 2 ]. We shall consider heavy-quark form factors in the large β0 limit, where β0αs ∼ 1, and 1/β0 is an expansion parameter (see the reviews [ 3–5 ]). A bare form factor can be written as F = 1 + ∞ L−1 L=1 n=0 aLn β0n g02 (4π )d/2 L Keeping terms with the highest degree of β0 in each order of perturbation theory, we get 1 F = 1 + β0 f β0g02 (4π )d/2 + O The leading coefficients aL,L−1 can easily be obtained from n L−1 terms (Fig. 1). We shall consider only the first 1/β0 f order.1 1 In some cases it is possible to obtain results for 1/β02 corrections; see, e.g. [ 6–8 ]. (1) (2) k1,2 10]: J (μ) = 2 Heavy-quark bilinear currents We consider the QCD currents J0 = Q¯ 0 Q0 = Z (αs(n f )(μ)) J (μ), = γ [μ1 · · · γ μn], where Q0 is a bare heavy-quark field. The antisymmetrized product of n γ matrices has the property γ μ γμ = η(d − 2n) , η = (−1)n. All results for form factors of this current will explicitly depend on n and η. In situations when the initial heavy-quark momentum p1 and the final one p2 can be written as p1,2 = mv1,2 + k1,2 (m 2 is the on-shell mass, v1,2 = 1) with small residual momenta m, these currents can be expanded in HQET ones [9, 2 i=0 Gi (μ, μ )O˜ i (μ ) + O where the leading HQET currents are J˜i0 = h¯v20 i hv10 = Z˜ (αs(nl )(μ)) J˜i (μ), i = , /v1 + /v2, /v1 /v2, and the O˜ i are local and bilocal dimension-4 HQET operators with appropriate quantum numbers. Here hv1,20 are two (unrelated) bare fields describing HQET quarks with the velocities v1,2 having small (variable) residual momenta; the HQET Lagrangian explicitly contains v1,2. These reference velocities can be changed by arbitrary small vectors of order ki /m (reparametrization invariance). The HQET (3) (4) (5) (6) mv k1 k2 a k3 k4 mv ω k1 k2 b k3 k4 ω Hi (μ, μ ) = Hˆi current renormalization constant Z˜ does not depend on the Dirac structure and is a function of the Minkowski angle ϑ : v1 · v2 = cosh ϑ = w. For our purpose it is convenient to choose v1,2 = p1,2/m, i. e., both residual momenta k1,2 = 0. Then the matrix elements of O˜ i vanish: non-zero expressions for these matrix elements (having dimensionality of energy) cannot be constructed, because we have no non-zero dimensionful parameters. The coefficients Hi in (5) can be obtained by matching the on-shell matrix elements (k1,2 = 0) in QCD and HQET: 2 i=0 Q( p2 = mv2)| J0|Q( p1 = mv1) = Fi u¯2 i u1, Q(k2 = 0)| J˜i0|Q(k1 = 0) = F˜i u¯2 i u1, F˜i = 1, (7) where u1,2 are the Dirac spinors of the initial quark and the final one (all loop corrections to F˜i vanish because they contain no scale). Therefore the bare matching coefficients (in the relation similar to (5) but for the bare currents) are H 0 i = Fi /F˜i = Fi . The renormalized matching coefficients are Hi (μ, μ ) = H 0 Z˜ (αs(nl )(μ )) Fi Z˜ i Z (αs(n f )(μ)) = F˜i Z UV divergences cancel in the ratio Fi /Z as well as in the ratio F˜i /Z˜ . Both Fi and F˜i contain IR divergences which cancel in the ratio Fi /F˜i because HQET is constructed to reproduce the IR behaviour of QCD (F˜i have no loop corrections because their UV and IR divergences cancel each other). The dependence of Hi (μ, μ ) on μ and μ is determined by the RG equations. Their solution can be written as × αs(n f )(μ) αs(n f )(μ0) αs(nl )(μ ) αs(nl )(μ0) γn0/(2β0(n f )) −γ˜0/(2β0(nl )) Kγ(nn f )(αs(n f )(μ)) K −(nγ˜l)(αs(nl )(μ )), where for any anomalous dimension γ (αs ) = γ0αs /(4π ) + γ1(αs /(4π ))2 + · · · we define Kγ (αs ) = exp Matrix elements of the currents with n = 0, 1 can be written via smaller numbers of form factors: Q(mv2)| J |Q(mv1) = F Su¯2u1, F S = F0 + 2F1 + (2w − 1)F2, where Fi with n = 0, η = 1 are used, and Q(mv2)| J μ|Q(mv1) = (F1V + F2V )u¯2γ μu1 (v1 + v2)μ −F2V u¯2u1 2 , F V 1 = F0 + 2F1 − (2w − 3)F2, F V 2 = −4(F1 + F2), where Fi with n = 1, η = −1 are used. 3 Inversion relations On-shell massive self-energy integrals with one massive line and any number of massless ones in some cases can be expressed via similar off-shell HQET integrals. Suppose all massless lines can be drawn at one side of the massive one and the resulting graph is planar (e.g., the diagram in Fig. 2a). Lines of such a diagram subdivide the plane into a number of polygonal cells (plus the exterior); with each cell we can associate a loop momentum (flowing counterclockwise). Then outer massless edges of the diagram correspond to the denominators −ki2 − i 0; inner massless edges to −(ki − k j )2 − i 0; and massive edges to m2 − (ki + mv)2 − i 0 (Table 1). The corresponding HQET diagram (Fig. 2b) has HQET denominators −2ki · v − 2ω − i 0 instead of massive ones. First we perform a Wick rotation of all loop momenta ki0 → i ki0 (in the v rest frame). Then, in Euclidean momentum space, we invert each loop momentum [ 11 ]: (10) (11) (12) (13) (14) (9) ki → kk2i . i see Table 1. As a result, a massive on-shell diagram (Fig. 2a) becomes m− ni (the sum runs over all massive line segments, ni are their indices, i. e. the powers of the denominators) times the off-shell HQET diagram (Fig. 2b) with ω = −(2m)−1 (15). The indices of all inner massless edges, as well as of all massive edges (which become HQET ones), remain intact (see Table 1). From the same table it is clear that the index of an outer massless edge becomes d − ni , where the sum runs over all edges of the cell to which this outer edge belongs (they can be all massless, or one of them can be massive). If there is a cell ki bounded only by inner massless edges, and maybe one massive one, then the denominator (ki2)d− n j will appear (Fig. 3). This denominator does not correspond to any line, and hence the resulting integral is not a Feynman integral at all; in this case, the discussed relation becomes rather useless (though formally correct). The inversion relations [ 11 ] were used, e.g., in [ 12–14 ]). The inversion relations can be generalized to similar vertex integrals; the masses of the initial particle and the final one may differ. At one loop (Fig. 4), the integrals I (n1, n2, n; ϑ; ω1, ω2) = M(n1, n2, n; ϑ; m1, m2) = dd k iπd/2 1 × [−k2 − 2m1v1 · k − i0]n1 [−k2 − 2m2v2 · k − i0]n2 (−k2 − i0)n dd k iπd/2 (16) Euclidean ki2 (ki − k j )2 The integrals I (17) have been investigated in [ 15 ]. Here we need only the integrals M (16) with m1 = m2; they reduce to the integrals I (17) with ω1 = ω2, which are especially simple [ 15 ]: I (n1, n2, n; ϑ ; ω, ω) = (−2ω)d−n1−n2−2n I (n1 + n2, n) × 3 F2 nn1 +12,nn2,2,n1d2+−n22+n1 1 − c2osh ϑ , where I (n1, n) = (−d + n1 + 2n) (d/2 − n) (n1) (n) (17) (19) (20) k + m2v2 m2v2 − q is the one-loop HQET self-energy integral. We only need integer n1,2; in this case all I reduce by IBP to 2 master integrals [ 15 ]: I (1, 0, n) (trivial) and I (1, 1, n) (given by (19)). Inversion relations can be generalized to diagrams with more external legs. For example, the one-loop massive box diagram with two on-shell legs and the corresponding offshell HQET one (Fig. 5) M(n1, n2, n3, n4; ϑ; m1, m2; q2, q · v1, q · v2) = 1 × (−k2 − 2m1v1 · k)n1 (−k2 − 2m2v2 · k)n2 (−(k + q)2)n3 (−k2)n4 , (21) I (n1, n2, n3, n4; ϑ; ω1, ω2; q2, q · v1, q · v2) = 1 × (−2k · v1 − 2ω1)n1 (−2k · v2 − 2ω2)n2 (−(k + q)2)n3 (−k2)n4 dd k iπd/2 dd k iπd/2 We need only terms with the highest degree of n f ; therefore, there is no need to distinguish between n f and nl = n f − 1, or any n f + const. The gluon propagator can be written as Dμν (k) = k2(1 − 1 (k2)) gμν − kμkν k2 , where the gluon self-energy is (k2) = β0 (4πg)02d/2 e−γ ε D(ε) (−k2)−ε, ε m2v2 Fig. 5 Box diagrams D(ε) = eγ ε (1 − ε) (1 + ε) 2(1 − ε) 5 (1 − 2ε)(1 − 23 ε) (1 − 2ε) = 1 + 3 ε + · · · (the same function appears also in the one-loop self-energy integral with arbitrary masses m1,2 and arbitrary p2, where At this leading large β0 order, the coupling constant renormalization is simple: β0 (4πg)02d/2 e−γ ε = b Zα(b)μ2ε, αs (μ) 1 b = β0 4π , Zα = 1 + b/ε . The bare QCD matrix elements can be written in the form [ 6,16 ] 1 ∞ fi (ε, Lε) Fi = δi0 + β0 L=1 L (−m2)L + O . (27) It is convenient to write the functions fi (ε, u) in the form usual for on-shell massive QCD problems (see [ 5 ]) eγ ε fi (ε, u) = CF D(ε) (1 − 2u) (1 + u) (3 − u − ε) Ni (ε, u). (28) We calculate the vertex function (Fig. 1) and multiply it by Z oQs with the 1/β0 accuracy (see [ 5 ]). Reducing on-shell massive QCD integrals to off-shell HQET ones by the inversion relation (18) and then to the master integrals by IBP [ 15 ], we obtain N0(ε, u) = −ηu n − 2 + ε − 2(w + 1)u(n − 2)2 w − 1 −u ηu + 4(w + 1)ε (n − 2) +2(2 − u) w + (w + 1)u −(6w + 2u + ηu2)ε +2 w − (w + 1)u ε2 F +ηu n − 2 + ε + 2(n − 2)2 + 4ε(n − 2) w − 1 −6(1 − u2) + 2(1 − u)(5 + 2u)ε −2(1 − 2u)ε2, n − 2 + ε N1(ε, u) = u ηw w − 1 − ηu(n − 2) − 2 + u + ε −ηuε F − ηu n − 2 + ε , w − 1 n − 2 + ε N2(ε, u) = ηu w − 1 ×[1 − (1 + (w − 1)u)F ], (26) (29) (30) (22) (23) (24) (25) where F = 2 F1 1, 1 + u 1 − w 3/2 2 both indices are equal to 1 [ 17 ]). At ϑ = 0 this result agrees with the result of [ 18 ] at m1 = m2; see also [ 5 ].2 Re-expressing the bare form factors (27) via the renormalized coupling we obtain 1 ∞ fi (ε, Lε) Fi = δi0 + β0 L=1 L D(ε) is regular at the origin; expanding (b/(ε+b))L in b, we obtain a quadruple sum. In the coefficient of ε−1 all fnm except fn0 cancel; differentiating this coefficient in log b (and using the fact that F (30) at u = 0 is ϑ/ sinh ϑ ) we obtain the anomalous dimension corresponding to Z /Z˜ [ 6,16 ]: These anomalous dimensions at the 1/β0 order are [ 19,20 ] γn = 4CF b (1 + 23 b) (2 + 2b) β0 (1 + b)2(2 + b) 3(1 + b) (1 − b) × (n − 1)(3 − n + 2b) + O γ˜ = 4CF b (1 + 23 b) (2 + 2b) β0 (1 + b) 3(1 + b) (1 − b) (37) (38) (39) (40) In the coefficient of ε0 all fnm except fn0 and f0m cancel. The coefficients fn0 form Kγn−γ˜ (αs (μ)), see (9); we have [ 6 ] where where the Borel images of the perturbative series for Hˆi are 1 Si (u) = u The integral (37) is not well defined because of poles at the integration contour. The leading renormalon ambiguities are given by the residues at u = 1/2 [ 21 ] (see also [ 5 ]). It is easy to calculate these residues because F (30) at u = 1/2 is just 2/(w + 1): As demonstrated in [ 21 ], matrix elements of the QCD currents between ground-state mesons (pseudoscalar or vector) are unambiguous: the IR renormalon ambiguities of the leading matching coefficients Hi are compensated by the UV renormalon ambiguities in the matrix elements of the 1/m suppressed HQET operators O˜ i in (5) (see also [ 5 ]). The hypergeometric function F (30) has been expanded in u to all orders [ 17 ], the coefficients are expressed via Nielsen polylogarithms Snm (x ). The result [ 17 ] is written for the case of an Euclidean angle3; its analytical continuation to Minkowski angles is 1 F = sinh ϑ (2 cosh(ϑ/2))2u n ∞ sinh(ϑ u) u un (−2)n−m Sm,n−m+1(−eϑ ) + O Our results satisfy this requirement ( f1,2(−b, 0) = 0 because the QCD current J does not mix with currents with other Dirac structures). 2 There are a few typos in Sect. 8.8 of [ 5 ]. The unnumbered formula below (8.93) should read R0 = cosh(Lu), R1 = sinh[(1 − 2u)L/2] sinh(L/2) In the second formula in (8.95), the coefficient of R0 should contain an extra factor 3. In both formulae in (8.96), their right-hand sides should be 1 + αs correction. −e−ϑu +eϑu n=1 ∞ un n=1 m=1 m=1 n (−2)n−m Sm,n−m+1(−e−ϑ ) . (41) It is possible to re-express this expansion in terms of Nielsen polylogarithms of just one argument, see [ 23 ], but then the symmetry ϑ → −ϑ will not be explicit. 3 M. Yu. Kalmykov has informed me that there is a typo: the power of cos ϑ in (2.7) should be 1 + 2ε. This typo has been corrected in [ 22 ]. Acknowledgements I am grateful to M. Steinhauser for useful comments and hospitality in Karlsruhe, where the major part of this work was done; to J. M. Henn for useful discussions and hospitality in Mainz; and to M. Yu. Kalmykov for bringing Ref. [ 17 ] to my attention and discussions related to it. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP3. Appendix A: Anticommuting γ5 and ’t Hooft–Veltman γ5 For flavour-nonsinglet currents one may use the anticommuting γ5 without encountering contradictions; they are related to the currents with the ’t Hooft–Veltman γ5 by a finite renormalization [ 24–26 ]: (q¯ γ5AC nτ q)μ = Z2−n(α(n f )(μ))(q¯ γ5HV nτ q)μ, (A.1) where τ is a flavour matrix with Tr τ = 0. The currents with γ5AC n have anomalous dimensions γn, because they can be obtained from the case of massless quarks; γ5HV n is just 4−n with reshuffled components. Equating the derivatives in d log μ we obtain Z2−n(αs ) = Kγ(nn −f)γ4−n (αs ), where the anomalous dimensions γn and γ4−n differ starting from two loops. In particular, Z0(αs ) = 1. In HQET currents with γ5AC and with γ5HV have the same anomalous dimension γ˜ , and the finite renormalization factor similar to (A.2) is 1. In the large β0 limit (see (35)) 8n Zn(αs ) = exp − β0 (1 + 23 b) (2 + 2b) (1 + b)2(2 + b) 3(1 + b) (1 − b) 0 b × +O db 1 β02 . At the leading 1/β0 order we may use these formulae for flavour singlet currents, too. The matrix γ5AC n has the same property (4) but with η = −(−1)n. From our results (27)– (29) we see that, indeed, Hˆγ5AC n = Hˆ n η→−η = Hˆγ5HV n = Hˆ 4−n . Matrix elements of the currents with γ5AC and n = 0, 1 can be written via smaller numbers of form factors: Q(mv2)| J |Q(mv1) = F P u¯2γ5ACu1, F P = F0 − 2F1 − (2w + 1)F2, (A.2) (A.3) (A.4) (A.5) 5 Appendix B: Expansion of the hypergeometric function F We can also find several terms of this expansion using the Mathematica package HypExp [ 27,28 ] (which uses HPL [ 29,30 ]). This results in 1 F = sinh ϑ u2 ϑ − H−+(τ )u − (H−+−(τ ) − 2H−+(τ )l) 2 3 −(H−+−−(τ ) − 2H−+−(τ )l + 2H−+(τ )l2) u 3 H−+−−−(τ ) − 2H−+−−(τ )l + 2H−+−(τ )l2 − where Fi with n = 0, η = −1 are used, and , where Fi with n = 1, η = 1 are used. The divergence of the axial current is i ∂μ(Q¯ 0γ5ACγ μ Q0) = 2m0 Q¯ 0γ5AC Q0, where the bare mass m0 = Z mosm. Taking the matrix element of this equation we obtain F A 1 + w − 1 F A 2 = Z mos F P . 2 The on-shell mass renormalization constant Z mos at the first 1/β0 order is given by the formula similar to (27), (28) with Nm (ε, u) = −2(3 − 2ε)(1 − u); see, e.g., [ 5 ]. And indeed, from (29), (A.5)–(A.6) we obtain N A 1 + w −2 1 N2A = N P + Nm . (A.6) (A.7) (A.8) (A.9) (B.10) where τ = tanh ϑ2 , l = 21 H−(τ ) = log cosh ϑ2 H+(τ ) = ϑ, , and H···(τ ) are harmonic polylogarithms (see [ 29–31 ]). Only one new polylogarithm appears at each order. In order to compare the expansion coefficients in (41) and in (B.10), we need to transform them to harmonic polylogarithms of the same argument, which we choose as x = e−ϑ . In (41), we first rewrite Snm (−x −1) via Snm (−x ) using the formula from [ 23 ]; then we rewrite Snm (−x ) via H···(−x ) and then via H···(x ); we rewrite log cosh(ϑ/2) (B.11) via H···(x ); and finally we re-express products of harmonic polylogarithms via their linear combinations. In (B.10) we rewrite harmonic polylogarithms with ± indices [ 30 ] via normal ones with indices 0, ±1; substitute τ = (1 − x )/(1 + x ) and re-express via H···(x ); and finally convert products of harmonic polylogarithms to sums. All these steps are done in Mathematica using HPL [ 29,30 ]. We have checked that all the coefficients presented in (B.10) agree with (41). (B.11) Appendix C: Vector form factors The vector form factors F1V,2 (13) can be written in the form (27), (28); from (29), (13) we obtain N1V (ε, u) = 2[2w + u − 3u2 − 3wε + 2wuε − (w − 3)u2ε +wε2 − (w + 1)uε2]F −2 2 + u − 3u2 − 3ε + 2uε +2u2ε + ε2 − 2uε2 , N2V (ε, u) = 4u(1 + u − 2uε)F. All loop corrections to F1V vanish at ϑ = 0, and hence N V 1 = 0 at w = 1. The form factor F1V = H1V /Z˜ , where Z˜ at the 1/β0 order is determined by the anomalous dimension (36), and H1V contains only non-negative powers of ε. We choose μ = μ = μ0 = m. H1V at ε = 0 is given by the coefficients fn0 (which produce K−γ˜ (10)) and f0n (which produce Hˆ1V (37)); εn terms (n > 0) require all fnm . Writing the expansion (B.10) as F = f0 − f1u − f2u2/2 − f3u3/3 − · · · we obtain up to four loops (C.12) (C.13) −2w f1 + (3w + 1) f0 − 4 H V 1 = 1 + CF β0 b − w f2 + (3w + 1) f1 − 2 − 3 w f3 + 3w + 1 2 f2 + 2 π 2 π 2 3 ζ3w − 4 w − 12 −16w 2 π 2 π 2 3 ζ3w − 4 w − 12 − 16w f1 π 4w + 92 nlL−1αsL terms with L = 1, 2, 3 in F1V,2 from [ 2 ]. Using HPL [ 29, 30 ] we have successfully reproduced all π 4 + 269ζ3 + 206 3 + · · · . (C.15) 1. W. Bernreuther , R. Bonciani , T. Gehrmann , R. Heinesch , T. Leineweber , P. Mastrolia , E. Remiddi, Two-loop QCD corrections to the heavy quark form-factors: the vector contributions . Nucl. Phys. B 706 , 245 ( 2005 ). arXiv:hep-ph/0406046 2. J.M. Henn , A.V. Smirnov , V.A. Smirnov , M. Steinhauser , Massive three-loop form factor in the planar limit . JHEP 01 , 074 ( 2017 ). arXiv: 1611 .07535 [hep-ph] 3. M. Beneke , Renormalons. Phys. Rep . 317 , 1 ( 1999 ). arXiv:hep-ph/9807443 4. M. Beneke , V.M. Braun , Renormalons and power corrections , in At the Frontier of Particle Physics: Handbook of QCD, ed. by M.A. Shifman (World Scientific , Singapore, 2001 ), vol 3 , p. 1719 . arXiv:hep-ph/0010208 5. A.G. Grozin , Heavy quark effective theory . in Springer Tracts in Modern Physics, vol 201 , Chap 8 (Springer, Berlin, 2004 ). arXiv:hep-ph/0311050 6. D.J. Broadhurst , Large N expansion of QED: asymptotic photon propagator and contributions to the muon anomaly, for any number of loops . Z. Phys. C 58 , 339 ( 1993 ) 7. A.G. Grozin , J.M. Henn , G.P. Korchemsky , P. Marquard , The threeloop cusp anomalous dimension in QCD and its supersymmetric extensions . JHEP 01 , 140 ( 2016 ). arXiv: 1510 .07803 [hep-ph] 8. A.G. Grozin , Leading and next-to-leading large-n f terms in the cusp anomalous dimension and quark-antiquark potential . PoS LL 2016 , 053 ( 2016 ). arXiv: 1605 .03886 [hep-ph] 9. A.F. Falk , B. Grinstein , Heavy meson pair production in e+eannihilation from the static quark effective theory . Phys. Lett. B 249 , 314 ( 1990 ) 10. M. Neubert , Short distance expansion of heavy quark currents . Phys. Rev. D 46 , 2212 ( 1992 ) 11. D.J. Broadhurst , A.G. Grozin , Multiloop calculations in heavy quark effective theory , in New Computing Techniques in Physics Research, ed. by B.H. Denby , D. Perret-Gallix (World Scientific , Singapore, 1995 ), p. 217 . arXiv:hep-ph/9504400 12. A. Czarnecki , K. Melnikov , Threshold expansion for heavy light systems and flavor off diagonal current correlators . Phys. Rev. D 66 , 011502 ( 2002 ). arXiv:hep-ph/0110028 13. A.G. Grozin, Lectures on multiloop calculations. Int. J. Mod. Phys. A 19 , 473 ( 2004 ). arXiv:hep-ph/0307297 14. A.G. Grozin, Lectures in QED and QCD, Chap 10 (World Scientific, Singapore, 2007 ) 15. A.G. Grozin , A.V. Kotikov , HQET heavy-heavy vertex diagram with two velocities . arXiv:1106 .3912 [hep-ph] 16. A. Palanques-Mestre , P. Pascual , The 1/N f expansion of the γ and β functions in QED . Commun. Math. Phys. 95 , 277 ( 1984 ) 17. A.I. Davydychev , M. Yu . Kalmykov, New results for the ε - expansion of certain one-, two- and three-loop Feynman diagrams . Nucl. Phys. B 605 , 266 ( 2001 ). arXiv:hep-th/0012189 18. M. Neubert , Higher order perturbative corrections to b → c transitions at zero recoil . Phys. Lett. B 341 , 367 ( 1995 ). arXiv:hep-ph/9409453 19. D.J. Broadhurst , A.G. Grozin , Matching QCD and HQET heavylight currents at two loops and beyond . Phys. Rev. D 52 , 4082 ( 1995 ). arXiv:hep-ph/9410240 20. M. Beneke , V.M. Braun , Power corrections and renormalons in Drell-Yan production . Nucl. Phys. B 454 , 253 ( 1995 ). arXiv:hep-ph/9506452 21. M. Neubert , C.T. Sachrajda , Cancellation of renormalon ambiguities in the heavy quark effective theory . Nucl. Phys. B 438 , 235 ( 1995 ). arXiv:hep-ph/9407394 22. A.I. Davydychev , M. Yu . Kalmykov, Massive Feynman diagrams and inverse binomial sums . Nucl. Phys. B 699 , 3 ( 2004 ). arXiv:hep-th/0303162 23. A. Devoto , D.W. Duke , Table of integrals and formulae for Feynman diagram calculations . Riv. Nuovo Cim. 7N6 , 1 ( 1984 ), formulas (2.1.7)-(2.1.9) . http://functions.wolfram.com/ ZetaFunctionsandPolylogarithms/PolyLog3/17/01/ 24. S.A. Larin , J.A.M. Vermaseren , The αs3 corrections to the Bjorken sum rule for polarized electroproduction and to the GrossLlewellyn Smith sum rule . Phys. Lett. B 259 , 345 ( 1991 ) 25. S.A. Larin , The renormalization of the axial anomaly in dimensional regularization , in Quarks-92, ed. by D. Yu. Grigoriev , V.A. Matveev , V.A. Rubakov , P.G. Tinyakov (World Scientific , Singapore, 1993 ), p. 201 . arXiv:hep-ph/9302240 26. S.A. Larin , Phys. Lett. B 303 , 113 ( 1993 ) 27. T. Huber , D. Maître, HypExp, a mathematica package for expanding hypergeometric functions around integer-valued parameters . Comput. Phys. Commun . 175 , 122 ( 2006 ). arXiv:hep-ph/0507094 28. T. Huber , D. Maître, HypExp 2, expanding hypergeometric functions about half-integer parameters . Comput. Phys. Commun . 178 , 755 ( 2008 ). arXiv: 0708 .2443 [hep-ph] 29. D. Maître, HPL, a mathematica implementation of the harmonic polylogarithms . Comput. Phys. Commun . 174 , 222 ( 2006 ). arXiv:hep-ph/0507152 30. D. Maître , Extension of HPL to complex arguments . Comput. Phys. Commun . 183 , 846 ( 2012 ). arXiv:hep-ph/0703052 31. E. Remiddi, J.A.M. Vermaseren , Harmonic polylogarithms . Int. J. Mod. Phys. A 15 , 725 ( 2000 ). arXiv:hep-ph/9905237