Heavy Quark Form Factors at Three Loops in the Planar Limit

We compute the color-planar and complete light quark non-singlet contributions to the heavy quark form factors in the case of the axialvector, scalar and pseudoscalar currents at three loops in perturbative QCD. We evaluate the master integrals applying a new method based on differential equations for general bases, which is applicable for all first order factorizing systems. The analytic results are expressed in terms of harmonic polylogarithms and real-valued cyclotomic harmonic polylogarithms.

Form factors are the matrix elements of local composite operators between physical states. In perturbative Quantum Chromodynamics (QCD) these objects play a significant role in determining physical observables. In scattering cross-sections, they provide important contributions to the virtual corrections. The massive form factors are of importance for the forward-backward asymmetry of bottom or top quark pair production at electron-positron colliders and to static quantities like the anomalous magnetic moment of a heavy quark and other processes. They are also of importance to scrutinize the properties of the top quark [1,2] during the high luminosity phase of the LHC [3] and the experimental precision studies at future high energy e + e − colliders [4].
In this letter, we calculate both the color-planar and complete light quark non-singlet threeloop contributions to the massive form factors for axialvector, scalar and pseudoscalar currents. Our results for the vector current, including a detailed account of the techniques used in these calculations, will be presented elsewhere [5]. The two-loop QCD corrections to the massive vector, axialvector form factors, the anomaly contributions, and the scalar and pseudoscalar form factors were first presented in [6][7][8][9]. In [10], an independent computation led to a cross-check of the vector form factor, giving also the additional O(ε) terms in the dimensional parameter ε = (4 − D)/2. Recently, the contributions up to O(ε 2 ) for all the massive two-loop form factors were obtained in Ref. [11]. The color-planar contributions to the massive three-loop form factor for the vector current have been computed in [12,13] and the complete light quark contributions in [14]. The large β 0 limit has been considered in [15].
Our notations follow those used in Ref. [11]. We consider the decay of a virtual massive boson of momentum q into a pair of heavy quarks of mass m, momenta q 1 and q 2 and color c and d, through a vertex X cd , where X cd = Γ µ A,cd , Γ S,cd and Γ P,cd indicates the coupling to an axialvector, a scalar and a pseudoscalar boson, respectively. Here q 2 = (q 1 + q 2 ) 2 is the center of mass energy squared and the dimensionless variable s is defined by The amplitudes take the following general form whereū c (q 1 ) and v d (q 2 ) are the bi-spinors of the quark and the anti-quark, respectively. We denote the corresponding UV renormalized form factors by F I , with I = A, S, P . They are expanded in the strong coupling constant α s = g 2 s /(4π) as follows The following generic forms for the amplitudes can be found by studying the general Lorentz structure. For the axialvector current, it can be cast as where σ µν = i 2 [γ µ , γ ν ] and a Q is the Standard Model (SM) axialvector coupling constant. Likewise, for the scalar and pseudoscalar currents, one has where v = ( √ 2G F ) −1/2 is the SM vacuum expectation value of the Higgs field, with G F being the Fermi constant, s Q and p Q are the scalar and pseudoscalar couplings, respectively. Finally, to obtain the unrenormalized form factors, we multiply appropriate projectors as provided in [11] and perform the trace over the color and spinor indices. For later purpose we denote by N c the number of colors, and n l and n h are the number of light and heavy quarks, respectively. We will set n h = 1 in the following.
Since we use dimensional regularization [16], one important point is to define a proper description for the treatment of γ 5 . Both the color-planar and complete n l contribution belong to the so-called non-singlet case, where the axialvector or pseudoscalar vertex is connected to open heavy fermion lines. Hence, both γ 5 -matrices appear in the same chain of Dirac matrices, which allows us to use an anti-commuting γ 5 in D space-time dimensions, with γ 2 5 = 1. This is implied by the well-known Ward identity, which in terms of the form factors, takes the form Here ns denotes the non-singlet contributions. For convenience, we introduce the kinematic variable [17] x which we use in the following. In particular, we focus on the Euclidean region, q 2 < 0, corresponding to x ∈ [0, 1[.  The Feynman diagrams for the different form factors are generated using QGRAF [18], the color algebra is performed using Color [19], the output of which is then processed using Q2e/Exp [20,21] and FORM [22,23] in order to express the diagrams in terms of a linear combination of a large set of scalar integrals. These integrals are then reduced using integration by parts identities (IBPs) [24,25] with the help of the program Crusher [26] to obtain 109 master integrals (MIs), out of which 96 appear in the color-planar case. In the color-planar limit, the families of integrals can be represented by eight topologies, shown in Figure 1, whereas for the complete light quark contributions, three more topologies, cf. Figure 2, are required 1 .
Finally, the master integrals have to be computed. For this we use the method of differential equations, see also [27][28][29][30]. The corresponding differential equations are obtained from the IBP relations. Here a central question is whether the corresponding linear system of differential equations is first order factorizable or not. Using the package Oresys [31], based on Zürcher's algorithm [32,33], we have proved that the present system is indeed first order factorizable in x-space. Without any need to choose a special basis, one is therefore in the position to solve the system in terms of iterated integrals of whatsoever alphabet, cf. Ref. [5] for details. The differential equations are solved order by order in ε successively, starting at the leading pole terms ∝ 1/ε 3 . The successive solutions in ε contribute to the inhomogeneities in the next order. We compute the master integrals block-by-block, where for an m × m system a single inhomogeneous ordinary differential equation of order m or less is obtained, which we solved using the variation of constant. The other m − 1 solutions result from the former solution immediately. The boundary conditions can be determined by a separate calculation at x = 1. The calculation is performed by intense use of HarmonicSums [34][35][36][37][38][39][40], which uses the package Sigma [41,42]. We finally have checked all master integrals numerically using FIESTA [43][44][45].
In the present case, the emerging harmonic polylogarithms stem from the inhomogeneities, adding further letters which result from the rational coefficients in the differential equations. They are obtained by partial fractioning as the k-th powers of letters, k ∈ N, which have to be transformed to the letters by partial integration in case. This method has some relation to the method of hyperlogarithms [46,47]. One obtains up to weight w=6 real-valued iterated integrals over the alphabet i.e. the usual harmonic polylogarithms (HPLs) [48] and their cyclotomic extension [34], including the respective constants in the limit x → 1, i.e. the multiple zeta values (MZVs) [49] and the cyclotomic constants [34,50]. In case of the iterated integrals we apply the linear representation. For a numerical implementation the use of the shuffle algebra [51] implemented in HarmonicSums reduces the number of functions accordingly. In the MZV and cyclotomic case there are proven reduction relations to weight w = 12 [49] and w = 5 [50], respectively, which we have used. The 64 cyclotomic constants which appear up to w = 5 reduce to 18. At w = 6 124 cyclotomic constants remain at the moment. Note that there are more conjectured relations, cf. [52], based on PSLQ [53]. If these conjectured relations are used, only multiple zeta values remain as constants in all form factors using our real representation for the cyclotomic harmonic polylogarithms. The analytic result for the different form factors in terms of HPLs and cyclotomic HPLs [34,48] can be analytically continued outside x ∈ [0, 1[ by using the mappings x → −x, x → (1 − x)/(1 + x) on the expense of extending the cyclotomy class in cases needed. The UV renormalization of the form factors has been performed in a mixed scheme. We renormalize the heavy quark mass and wave function in the on-shell (OS) renormalization scheme, while the strong coupling constant is renormalized in the MS scheme, which is given by setting the universal factor S ε = exp(−ε(γ E − ln(4π)) for each loop order to one at the end of the calculation. The required renormalization constants are well known and are denoted by Z m,OS [54][55][56][57][58], Z 2,OS [54][55][56]59] and Z as [60,61] for the heavy quark mass, wave function and strong coupling constant, respectively. For all the cases, the renormalization of the heavy-quark wave function and the strong coupling constant are multiplicative, while the renormalization of massive fermion lines has been taken care of by properly considering the counter terms. For the scalar and pseudoscalar currents, the presence of the heavy quark mass in the Yukawa coupling employs another overall mass renormalization constant.
The infrared (IR) singularities of the massive form factors can be factorized [62] as a multiplicative renormalization factor. Its structure is constrained by the renormalization group equation (RGE), as follows, where F fin I is finite as ε → 0. The RGE for Z(µ) reads where Γ is the corresponding cusp anomalous dimension, which is by now available up to threeloop order [63,64]. Notice that Z does not carry any information regarding the vertex. Both Z and Γ can be expanded in a perturbative series in α s as follows and one finds the following solution for Eq. (11) Eq. (13) correctly predicts the IR singularities for all massive form factors at the three-loop level.
We finally obtain the color-planar and the complete light quark non-singlet (n l ) contributions for the three-loop massive form factors for the same currents as before. Since the expressions are very long, we provide them as supplemental material along with this publication only. In Figures 3-4 we illustrate the behaviour of the O(ε 0 ) parts of the different form factors as a function of x ∈ [0, 1]. We also show their small-and large-x expansions. The latter representations are obtained using HarmonicSums. The different limits are characterized as follows : Low energy region (x → 1): In the space-like case (q 2 < 0) we expanded the form factors, redefining x = e iφ , φ = 0. High energy region (x → 0): Here we expand the form factors up to O(x 4 ). The chirality flipping form factors F V,2 and F A,2 vanish and the effect of γ 5 gets nullified in this limit implying F V,1 = F A,1 and F S = F P .
Threshold region (x → −1): Here expansions of the form factors in β = 1 − 4m 2 q 2 describe the dominant terms. For the numerical evaluation of the HPLs and the cyclotomic HPLs we use the GiNaC-package [65,66]. We performed a series of further checks. The Ward identity (7) has been checked by an explicit calculation. By maintaining the gauge parameter ξ to first order, a partial check on gauge invariance has been obtained. After α s -decoupling the UV renormalized results satisfy the universal IR structure, confirming again the correctness of all pole terms. Finally, we compared our results with those of Ref. [67], which has been obtained using partly different methods, and agree by adjusting the respective conventions.