Massive three-loop form factor in the planar limit

We compute the three-loop QCD corrections to the massive quark-anti-quark-photon form factors F1 and F2 in the large-Nc limit. The analytic results are expressed in terms of Goncharov polylogarithms. This allows for a straightforward numerical evaluation. We also derive series expansions, including power suppressed terms, for three kinematic regions corresponding to small and large invariant masses of the photon momentum, and small velocities of the heavy quarks.


Introduction
Massive form factors are important building blocks for various physical quantities involving heavy quarks. Among them are static quantities like anomalous magnetic moments but also production cross sections and decay rates. Furthermore, form factors are the prime examples for studying the infrared behaviour of QCD amplitudes.
We consider QCD corrections to the quark-photon vertex. The latter can be parametrized as follows, V µ (q 1 , q 2 ) =ū(q 1 )Γ µ (q 1 , q 2 )v(q 2 ) , (1.1) where the colour indices of the quarks are suppressed andū(q 1 ) and v(q 2 ) are the spinors of the quark and anti-quark, respectively. The momentum q 1 is incoming and q 2 is outgoing with q 2 1 = q 2 2 = m 2 . The vertex function Γ µ (q 1 , q 2 ) can be decomposed into two scalar form factors which are usually introduced as where q = q 1 −q 2 is the outgoing momentum of the photon and σ µν = i[γ µ , γ ν ]/2. Q q is the charge of the considered quark. Sample Feynman diagrams can be found in figure 1. Two-loop QCD corrections to the electric and magnetic form factors for the heavy quark vector current have been computed for the first time in ref. [1] (axial vector and anomaly contributions have been considered in [2,3]) where analytic results have been obtained. An independent cross check of the twoloop results for F 1 and F 2 has been performed in [4] where also O(ǫ 2 ) and O(ǫ) terms have been added to the one-and two-loop results, respectively. The results have been used to obtain predictions for the three-loop form factor F 1 in the high energy limit, by exploiting evolution equations and the exponentiation of infrared divergences (see also ref. [5] for earlier considerations).
In this paper we compute the three-loop form factor in the planar limit, keeping the exact mass dependence. After expanding our exact result for small quark masses we can compare to the high-energy results of [4] mentioned above, and complete them by determining the unknown constants in the 1/ǫ and ǫ 0 part. We furthermore provide power-suppressed terms.
Massive form factors have infrared divergences that are well understood. After the ultraviolet renormalization, all poles in dimensional regularization are given in terms of the cusp anomalous dimension [6,7], and the beta function. The three-loop cusp anomalous dimension was computed in refs. [8,9]. By verifying the infrared pole structure at the three-loop order, we provide a first independent check of the result of refs. [8,9] (in the planar limit).
In the static limit, the infrared divergences disappear, and F 1 and F 2 are finite. In fact, F 1 vanishes and F 2 determines the anomalous magnetic moment of a heavy quark which has been considered at two-loop order in ref. [10]. A dedicated calculation at three loops has been performed in ref. [11] which serves as a welcome check for our exact result expanded for q 2 → 0.

JHEP01(2017)074
The remainder of the paper is organized as follows. In section 2 we provide technical details on the calculation of the amplitudes. In particular we briefly describe the renormalization procedure. The infrared structure of the form factors is presented in section 3. Our results for F 1 and F 2 are discussed in section 4 including the three-loop results for the static limit, the high-energy limit and for small quark velocities. We conclude in section 5.

Setup and calculation
The form factors F 1 and F 2 appearing in eq. (1.2) are conveniently computed with the help of projectors which are applied to Γ µ (q 1 , q 2 ). Using the kinematics defined in eq. (1.2) we have (i = 1, 2) 1 and s = q 2 . It is convenient to introduce the dimensionless variable Then the low-energy, high-energy and threshold limits correspond to x → 1, x → 0 and x → −1, respectively. Note that for x > 0 we have s < 0 and thus the form factors do not have imaginary parts. The same is true for x ∈ C with |x| = 1. For 0 < s < 4m 2 we have that x is on the upper half of the unit circle. It is convenient to write the perturbative expansion of F i (i = 1, 2) in the form with F (0) 1 = 1 and F (0) 2 = 0. In the large-N c limit we furthermore have that F where n l counts the number of closed massless quark loops and N c is the number of colours. Note that we do not consider contributions with massive closed fermion loops. In eq. (2.4) we suppress the scale dependence of α s and F (n) i . The calculations performed in this paper use the groundwork performed in ref. [12] where all scalar integral families up to three loops, which are needed for the massive form factors F 1 and F 2 in the large-N c limit, have been classified and the corresponding master integrals have been computed analytically in terms of Goncharov polylogarithms [13]. 2 We JHEP01(2017)074 use in particular the information from figure 1 of ref. [12] where eight three-loop families are defined. This information is used to generate with the help of the programs qgraf2.0 [14] and q2e/exp [15][16][17] amplitudes for F 1 and F 2 which are expressed in terms of linear combinations of integrals from the eight three-loop families. We also use formulae for reduction of Goncharov polylogarithm values at sixth roots of unity derived in [18].
For the reduction to master integrals we use the program FIRE5.2 [19][20][21][22] in combination with LiteRed [23,24]. Once the reduction for each family is complete we use the program tsort, which is part of the latest FIRE version [21] (implemented in the command FindRules) and based on ideas presented in ref. [20], to obtain relations between primary master integrals, and to arrive at a minimal set. This leads to 89 master integrals needed for the large-N c limit of F 1 and F 2 .
In our calculation we allow for a general QCD gauge parameter ξ but set ξ 2 terms (ξ = 0 corresponds to Feynman gauge) to zero before performing the reduction to master integrals. The bare form factors still contain linear ξ terms which only drop out after renormalization. This serves as a welcome check for our calculation.
The ultraviolet renormalized form factors are obtained by renormalizing the strong counpling constant α s in the MS scheme and the heavy quark mass on-shell. Both counterterms are needed to two-loop accuracy and are well-known in the literature. Note, however, that for the on-shell mass counterterm higher order ǫ terms are needed. The latter can be found in ref. [25].
In this context we would like to mention that in ref. [1] an non-standard version of the MS scheme has been employed as can bee seen from eq. (24) of that paper. The quantity C(ǫ), which enters the definition of the renormalization constant, induces π 2 terms which enter the ǫ 0 part of the two-loop form factor. See also the discussion in ref. [5] on this subject.
A further ingredient to the renormalization procedure is the on-shell wave function renormalization constant for the external heavy quarks which is needed to three-loop order and can be found in refs. [25,26].

Infrared divergences of massive form factors
Form factors of massive particles have infrared divergences originating from exchanges of soft particles. The latter can be described in the eikonal approximation. In this way, the infrared divergences of the form factors can be mapped to ultraviolet divergences of Wilson lines [27]. The relevant Wilson line has the geometry of a cusp formed by the particle momenta. It obeys a renormalization group equation that is governed by the cusp anomalous dimension [6,7,28].
Applying this correspondence to the original form factors, one has where Z is an infrared renormalization factor (in minimal subtraction), F is the ultravioletrenormalized form factor, and F f is finite both in the ultraviolet and infrared. In other words, all infrared poles of F are reproduced by Z.

JHEP01(2017)074
Z satisfies the following renormalization group equation where α s is the renormalized strong coupling and β D is the D-dimensional β function, Here C F = (N 2 c −1)/(2N c ) and C A = N c are the quadratic Casimir operators of the SU(N c ) gauge group in the fundamental and adjoint representation, respectively, n l is the number of massless quark flavors, and T = 1/2.
The perturbative expansions of Γ cusp and Z have the form Solving eq. (3.2) to three loops, one finds The cusp anomalous dimension in QCD was computed to three loops in [6][7][8][9].
In this way, one can see explicitly the poles generated by the right-hand side of eq. (3.1). We have verified that this equation correctly predicts all infrared poles in F 1 and F 2 to three loops.

Structure of results for form factors
Before presenting explicit results we briefly discuss the general structure of our analytic expressions.
All relevant master integrals were computed analytically in ref. [12]. From this it is clear that the form factors are given in terms of iterated integrals, with certain rational prefactors. The required set of integration kernels are We sometimes refer to the arguments of the logarithms x, Up to two-loop order and for the three-loop fermionic contributions (i.e. the n 1 l and the n 2 l terms) we observe only master integrals with letters x, 1 − x and 1 + x. This means that all of them can be expressed in terms of usual harmonic polylogarithms [29,30].
On the other hand, the non-fermionic three-loop part has the additional letter 1 − x + x 2 . Introducing the complex roots of this polynomial, In this way, all results can be written in terms of Goncharov polylogarithms. See ref. [12] for more details.
At three-loop order we observe that r 1 = e iπ/3 plays a special role for the form factors F 1 and F 2 since the coefficients of the Goncharov polylogarithms develop poles up to sixth order in x − r 1 . We could show that these poles are artificial by expanding the Goncharov polylogarithms around x = e iπ/3 . The analytic expressions for the finite result for F 1 and F 2 for x = r 1 are quite lengthy and can be found (for µ 2 = m 2 ) in the ancillary file.

Analytical results
We refrain from providing the results for the full three-loop form factors since the analytic expressions are too lengthy. All results which are discussed in this section can be downloaded from https://www.ttp.kit.edu/preprints/2016/ttp16-053/.
It is instructive to consider the form factors F 1 and F 2 in various kinematical regions which have already been mentioned in the Introduction. They are discussed in the remaining part of this section. In section 4.3 they are numerically compared to the exact result.
We start with the limit s ≪ m 2 which we obtain by expanding the Goncharov polylogarithms in the master integrals for x → 1. The expansion has to be carried out carefully since there are higher order poles in 1/(1 − x) in the prefactor. In fact, we expand all master integrals up to order (1 − x) 9 and obtain F 1 and F 2 up to order (1 − x) 4 . For the presentation in the paper we write x = e iφ and we restrict ourselves to expansion terms up JHEP01(2017)074 to order φ 2 , which for µ 2 = m 2 are given by Note that F 2 (x = 1) is finite and agrees with eqs. (54) and (55) of ref. [11] after adapting the large-N c limit. F 1 (x = 1) = 0 is a welcome check of our calculation, see also [31]. We expand all master integrals for x → 0 up to order x 6 which is sufficient to obtain F 1 and F 2 up to order x 4 . For illustration we show the one-, two and three-loop results including the first power-suppressed corrections of order x 1 . It is convenient to write the n-loop component of F i in the high-energy limit as follows Our results for F 1 read (for  (3)   with l x = log(x). It is interesting to note that the leading term f 1,lar the highest power is five, in agreement with the investigations of ref. [32] in the framework of QED. In that reference the leading behaviour of F 1 in the limit x → 0 has been considered and the coefficient of l 6 x has been determined. Translated to QCD, this correspond to sub-leading colour factors C 3 F . For F 2 we get the following expansion coefficients Note that the coefficients f (n,k) i,lar contain logarithmic terms in x which leads to a divergent behaviour of F 1,lar have been predicted using evolution equations. However, the ǫ 0 term, the sub-leading terms of order x n with n ≥ 1, and the results for f (3,n) 2,lar are new (see also subsection 4.4). Finally, we want to remark that higher order ǫ terms for the oneand two-loop coefficients can be found in the ancillary file.

Threshold
To obtain the threshold limit we expand the master integrals up to order (1 + x) 6 . After inserting the expanded results into the expressions for the form factors it is convenient to use is the velocity of the produced heavy quarks. Note that the ultravioletly renormalized form factors develop poles up to order 1/β n where n = 1, 2, 3 is the number of loops. On the other hand, the bare form factors have poles up to 1/β 2n (cf. ref. [4] where bare two-loop results are presented). Since the resulting expressions are quite large we refrain from displaying them in the paper but refer to the ancillary file which comes together with this paper. It is, however, instructive to look into the cross section σ(e + e − → QQ), where Q is a heavy quark. Close to threshold it is determined by the virtual correction, i.e., the form factors F 1 and F 2 , since the contributions from real radiation are suppressed by a relative factor β 3 . In fact, we can write where σ 0 = 4πα 2 Q 2 Q /(3s). Our calculation of F 1 and F 2 determines the first three terms for each ∆ (n) in the expansion for β → 0. Note that individually F 1 and F 2 still contain poles in ǫ, however, the combination given in eq.
where the ellipses refer to higher order terms in β. The one-and two-loop expressions agree with the large-N c limit of refs. [33][34][35] and the three-loop terms agree with refs. [36][37][38]. 3 At n-loop order the leading term of ∆ (n) behaves as (α s /β) n which is determined by the Sommerfeld factor 4 S = z/(1 − e −z ) with z = C F α s π/β. It is interesting to note that the series expansion of S has no term of order α 3 s and thus ∆ (3) starts at order 1/β 2 which is confirmed by our explicit calculation.
In the context of effective theories an important quantity derived from the form factor F 1 is the matching coefficient between QCD and non-relativistic QCD of the vector current. It is obtained by considering the on-shell photon-quark vertex for q 2 = 4m 2 . Due to the singularities in 1/β (see above) it is not possible to obtain the matching coefficient from the general result for F 1 . Rather a dedicated calculation is necessary which has been performed in [40] to three-loop order using semi-analytical methods. The planar master integrals of [40] have been computed in [12] as by-product of the calculation of all master integrals used in this calculation.

Numerical results
This subsection is devoted to the numerical evaluation of the form factors which we perform with the help of ginac [41,42]. In figures 2 and 3 F 1 and F 2 are shown as a function of x where the leading term of eq. (4.3) is subtracted to obtain a regular behaviour for x = 0 (which corresponds to s = ∞). From left to right the one-, two-and three-loop results are shown and the upper plots correspond to the real and the lower ones to the imaginary parts. Note that the latter are zero for x > 0. One observes that the expansions for s ≫ m 2 (which include terms up to order x 4 ) provide a good approximation to the exact result in the interval −0.3 x 0.3 which corresponds to 0.18 m 2 /s −0.61. On the other hand, the approximations obtained for s ≪ m 2 (which include terms up to order (1 − x) 4 ) agree with the exact result for x 0.4. Figure 4 shows the dependence of F 1 (top plots) and F 2 (bottom plots) as a function of φ where x = e iφ . In this region the form factors are real. One observes good agreement of the expanded and exact result up to φ ≈ 0.5 which corresponds to s/m 2 0.25.

Checks
Our result has passed several cross checks and consistency relations which we describe in this subsection.
We have successfully compared our bare and UV-renormalized one-and two-loop results (expanded up to O(ǫ 0 )) to the expressions provided in refs. [1,4] after taking the large-N c limit. Note that in [1] a different renormalization scheme has been used which leads to a difference in the finite contribution proportional to π 2 . This is due to the factor

JHEP01(2017)074
Γ(1 + ǫ) which is present in the counterterm for the strong coupling constant in eq. (24) of ref. [1] (see also discussion above). For the UV-renormalized two-loop form factor F 2 we agree with ref. [4] including O(ǫ 1 ) terms. For F 1 we disagree with the order ǫ term of ref. [4] in a term which is independent of Goncharov polylogarithms. The difference of our result and the one of [4] reads (4.10) In our expression there is no 1/(1 + x) 6 term at all. Such a term leads to a different lowenergy and threshold behaviour. In particular, the O(ǫ 1 ) term of the renormalized two-loop form factor would have a stronger divergence than the expected 1/β 2 behaviour, cf. the ancillary file to this paper. Furthermore, a term as in eq. (4.10) influences via renormalization the ǫ 0 terms of the three-loop F 1 which would lead to different low-energy and threshold expansions than the ones discussed in section 4.2. In particular, F 1 (x = 1) would be different from zero and the agreement of ∆ (3) in eq. (4.9) with the literature would be destroyed. As a further cross check we also compared with predictions of three-loop corrections to F 1 in the high-energy limit which have been obtained in ref. [4] on the basis of evolution equations. We find agreement including the log(x)/ǫ terms. The remaining 1/ǫ and the ǫ 0 terms cannot be predicted using the method of ref. [4]. However, these terms are contained in our result.
From the 1/ǫ pole of our result we can extract with the help of eq. (3.6) the cusp anomalous dimension Γ cusp up to three-loop order in the large-N c limit. Up to two-loop order we find agreement with refs. [7,43] and at three loops we can reproduce the results of [8,9]. This is the first independent check of (part of) the results obtained in [8,9] using a completely different method.
We have checked that the renormalized form factors have the correct static limit. In particular, F 1 (0) vanishes and F 2 (0) agrees with the explicit calculation of the three-loop corrections to the anomalous magnetic moment of a heavy quark which was performed in ref. [11].
For x ∈ (0, 1] we have that s ≤ 0. Thus the results for the form factors have to be real. Since the individual Goncharov polylogarithms are complex-valued this is a useful cross check. Similarly, if x = e iφ with either φ ∈ [0, π] or φ ∈ [−π, 0] (i.e. x is on the upper or lower semi-circle) we have that s is below threshold with 0 ≤ s/m 2 ≤ 4. Again, the form factors must be real-valued.

Conclusions and outlook
In this paper we evaluated for the first time massive three-loop form factors, in the planar limit. As a byproduct, we confirmed the recent result for the three-loop cusp anomalous dimension in the large-N c limit, which describes the infrared divergences of the form factors. We expressed the results analytically in terms of Goncharov polylogarithms. The latter allow for a straightforward numerical evaluation.

JHEP01(2017)074
We investigated analytically the low-energy, threshold, and high-energy limits, and derived expressions containing logarithmically enhanced as well as power suppressed terms. It would be interesting if some of these expansions could be obtained from effective field theory methods. See for example refs. [44,45] for work on power-suppressed terms.
Our results can be used to predict infrared divergent terms at higher loop orders, via renormalization group equations, along the lines of refs. [4,5].