Four-loop quark form factor with quartic fundamental colour factor

We analytically compute the four-loop QCD corrections for the colour structure $(d_F^{abcd})^2$ to the massless non-singlet quark form factor. The computation involves non-trivial non-planar integral families which have master integrals in the top sector. We compute the master integrals by introducing a second mass scale and solving differential equations with respect to the ratio of the two scales. We present details of our calculational procedure. Analytical results for the cusp and collinear anomalous dimensions, and the finite part of the form factor are presented. We also provide analytic results for all master integrals expanded up to weight eight.


Introduction
Form factors are indispensable vertex functions which enter a number of quantities in precision physics. Most prominent examples are the virtual corrections to the Drell-Yan process or inclusive Higgs boson production. Form factors are furthermore the simplest Green's function with a non-trivial infrared structure. In fact, from the pole parts of the form factors it is possible to extract universal quantities, like the cusp or collinear anomalous dimension. They enter general formulae which predict the infrared pole structure of massless on-shell multi-loop multi-leg QCD amplitudes [1,2].
In this paper we consider the quark-anti-quark-photon form factor with massless quarks which is obtained from the corresponding vertex function Γ µ q via F q (q 2 ) = − 1 4(1 − ǫ)q 2 Tr q 2 / Γ µ q q 1 / γ µ , (1.1) where we work in d = 4 − 2ǫ space-time dimensions, q = q 1 + q 2 , and q 1 (q 2 ) is the incoming quark (anti-quark) momentum. Two-loop corrections to F q have been computed for the first time more than twenty years ago [3][4][5][6] and the three-loop terms are available since about ten years [7][8][9][10][11] (for the computation of master integrals see also ref. [12]). Only two years ago first four-loop result for F q became available: in a first step the large-N c limit has been considered, where only planar Feynman diagrams contribute, and the fermionic and non-fermionic corrections have been computed in refs. [13] and [14], respectively. Fermionic corrections with three closed quark loops have been computed in ref. [15]; the complete terms proportional to n 2 f are available from [16].
Important information about QCD amplitudes is already obtained from the pole part of the form factor. Of particular interest in this respect is the cusp anomalous dimension, γ cusp [17], which can be extracted from the 1/ǫ 2 pole of F q . At three-loop order first results for γ cusp have been computed from the asymptotic behaviour of splitting functions [18] JHEP02 (2019)172 where the fractional hadron momentum tends to 1. The results have been confirmed afterwards by a dedicated calculation of the pole parts of the form factor [19]. Also at four-loop order there are two approaches to obtain γ cusp : the n 3 f terms of γ cusp has been obtained in refs. [15,20,21] and analytic results in the large-N c limit and for the (complete) n 2 f contributions have been obtained in refs. [13,14,16] and [22] from the explicit calculation of the form factor and the splitting functions in the threshold limit, respectively. The approach used in [22] could be extended to all colour structures; numerical results are presented in refs. [23,24]. Recently the abelian four-loop contribution of the linear n f term to γ cusp has been computed analytically in ref. [25]. The main focus of [25] is the cusp anomalous dimension for massive fermions in QED. The abelian n f term for massless quarks is obtained as a by-product.
We define the expansion of F q in terms of the bare strong coupling constant as The universal quantities γ cusp and γ q are conveniently extracted from the pole part of log(F q ) after renormalization of α s (see, e.g., refs. [2,8,17]). We define the corresponding loop expansions as follows In order to fix the normalization we provide the one-loop results which read γ 0 cusp = 4 and γ 0 ). In this work, we provide analytic four-loop results for γ cusp , γ q and F q for the colour structure (d abcd F ) 2 which for a SU(N c ) group is given by This colour structure appears for the first time at four-loop order and thus behaves in many respects as a leading order contribution. For example, there is no contribution from ultraviolet renormalization to this colour factor. It is furthermore responsible for the breaking of the naive Casimir scaling (i.e. the replacement of C F by C A ) used to relate the cusp anomalous dimensions of the fermion and gluon form factors (see, e.g., ref. [24] for more details).
The colour factor of eq. (1.4) arises from diagrams where four gluons connect the two external fermion lines, see figure 1(a). Note that there are also singlet diagrams with colour factor proportional to (d abcd F ) 2 , see figure 1(b). In this work we only consider non-singlet contributions.

Calculation
There are 18 Feynman diagrams with a closed fermion loop which is connected to the external fermion line via four gluons. A representative diagram is shown in figure 1  other diagrams are obtained by the various possibilities to connect the four gluons to the external fermion lines. We can map the 18 (six planar and twelve non-planar) diagrams to six integral families, two planar and four non-planar ones. They are illustrated in figure 2 where thin solid lines represent massless propagators. 1 The thick external line carries the virtuality q 2 . The planar families have been studied in refs. [13,14] where in particular all master integrals have been computed. Analytic results for the non-planar families in figure 2 are not yet available in the literature. In the following we concentrate our discussion on them.
With the help of a suitably chosen projector to obtain F q (introduced in eq. (1.1)) we can express the amplitude as a linear combination of scalar functions, which correspond to the family definitions of figure 2 propagators and six for irreducible numerators. We use FIRE [26][27][28] in combination with LiteRed [29,30] for the reduction to master integrals. In table 1 we present some information about the individual (non-planar) families. Altogether we have to compute about 50 000 integrals which can be reduced to almost 200 master integrals. We refrain from minimizing the master integrals among the various families since our approach (see below) is applied to a whole family and provides simultaneous results for all master integrals. We nevertheless establish relations between master integrals of different families and use them as cross checks for our results. For example, 36 of the 41 master integrals from df2-5 can be mapped to master integrals of df2-2. Note that we have performed the calculation in Feynman gauge. For the computation of the master integrals we use the idea suggested in [31] and used in our previous works for the planar [13,14] and n 2 f calculation [16]: we introduce a second mass scale q 2 2 = xq 2 as the virtuality of one of the external quarks. This increases, of course, the complexity of the problem. We encounter a more difficult reduction problem and there are significant more master integrals present in the individual families (compare "# 1-scale MIs" and "# 2-scale MIs" in table 1). However, the introduction of the second mass scale has the advantage that we can use the powerful method of differential equations. In fact, the basic idea is to choose x = 1 in order to fix the boundary conditions, since in this limit one has to compute massless two-point functions which are well studied in the literature [32,33]. The differential equations are then used to transport the information to the point x = 0. The method has been described in some details in ref. [16] where for the first time non-planar four-loop families have been considered. For the integral families considered in this paper the method had to be further refined. Note that in ref. [16] no non-planar master integrals had to be computed in the top sector where the indices of all twelve propagators are positive.
For each family we can introduce a system of differential equations of the form 2 where j(x) is a vector of (two-scale) master integrals in the primary basis chosen by FIRE and m(x) is a square matrix. We use the idea suggested in refs. [34,35] to turn to a so-called ǫ or canonical basis where the right-hand side of the differential equations is proportional to JHEP02(2019)172 ǫ and singularities with respect to the variables of the differential equations are Fuchsian, i.e., of the form 1/(x − a). To arrive at a canonical basis, we use the algorithm of ref. [36] 3 and its private implementation. We apply this procedure to each family separately and arrive at an ǫ form given by where J are the master integrals in the canonical basis, which are connected to the ones in the primary basis via j( with constant matrices M a . In our case the sum only includes two terms, a = 0 and a = 1, which correspond to the physical point and the point where we want to fix the boundary conditions, respectively. Next, we introduce, as in [16], the path-ordered exponent and define the quantities (with a slight abuse of notation) 4 U (x, 0) = lim which have the properties Note that U (x, 0) and U (x, 1) can be obtained in a straightforward way as an expansion in ǫ in terms of Harmonic polylogarithms (HPLs) [42] with arguments (1 − x) and x, respectively. Furthermore, both U (x, 0) and U (x, 1) solve the system (2.2) and are thus related by a matrix U 01 which only depends on ǫ but not on x: We will call the matrix U 01 ≡ U 01 (ǫ) the associator. It can be constructed by multiplying eq. (2.7) by x −ǫA 0 from the left and taking the limit x → 0 which leads to In practice, the right-hand side of eq. (2.8) is evaluated by extracting all log(x) terms contained in U (x, 1) with the help of shuffle relations to eliminate the leading letter "1"

JHEP02(2019)172
from the HPLs. 5 They have to cancel against the log(x) terms from x −ǫA 0 such that the limit x → 0 can be taken. Let us in a next step discuss the boundary conditions which we compute for x = 1. Note that in this limit our integrals are analytical and thus we do not have contributions of the form (1 − x) −kǫ with k = 0. In the canonical basis we can thus write where C 1 is a vector with ǫ-dependent components. Similarly we have (2.10) Note that in this limit the integrals in J have a logarithmic dependence on x. We are only interested in the so-called hard part which means that from the various contributions of the form x −kǫ we only take those with k = 0.
Next we want to relate the constants C 0 and C 1 to coefficients of integrals from the primary basis evaluated near x = 0 and x = 1, respectively. These relations have the form where L 0,1 are matrices depending on ǫ, and c 0,1 are the column vectors of the specific coefficients in the asymptotics x → 0 and x → 1, respectively. Note that the vector c 1 is obtained from the boundary conditions, and the aim of our calculation is the hard part of c 0 . In the following we present details about how we determine which set of coefficients c 0 suffices and calculate the matrix L 0 . L 1 and c 1 are calculated in analogy. We start with the generalized series expansion of T (x, ǫ) U (x, 0) which can be cast in the form where α = n 1 + ǫn 2 with integer n 1 and n 2 , and u (α, k) are matrices which depend on ǫ.
The key point is that, using the approach of ref. [44], we can calculate plenty of terms in the above expression, keeping the exact ǫ dependence. After applying eq. (2.12) to C 0 we have where c (α, k) = u (α, k) C 0 . (2.14) Each c (α, k) is a column vector of the form (c 1 (α, k) , . . . , c N (α, k)) ⊺ , where N is the number of two-scale master integrals of the considered family.

JHEP02(2019)172
In a next step we select from the coefficients c i (α, k) (for various i, α, and k) the minimal set, which is sufficient to determine all constants in C 0 = (C 01 , . . . C 0N ) ⊺ . Let this set be Here sufficient refers to the rank of the matrix which has to be greater or equal to the number of master integrals N , and minimal means that M = N . In other words, R is a square matrix, which is invertible and we have Of course, this procedure does not lead to unique quantities c 0 and L 0 , which, however, is not a problem since the arbitrariness cancels after performing the matching to the one-scale master integrals. As a rule of thumb we first try to pick coefficients only among the leading coefficients of the asymptotic expansion of the integrals j(x) and then extend the search to subleading terms in x, if necessary. Using eqs. (2.7), (2.9), (2.10) and, (2.11) we finally arrive at which is used to obtain the coefficients at x = 0 from the ones at x = 1. Note that L 0 and L 1 are exact in ǫ but U 01 is usually known as an expansion for ǫ → 0. The number of components of c 0 is the number of the two-scale master integrals. For example, for df2-2, it is 337. Our goal is the determination of the coefficients in the naive part of the expansion, i.e. the part of the expansion with non-negative integer powers of x. For df2-2, c 0 contains 116 coefficients corresponding to the naive limit. One can expect that this number is equal to the number of one-scale master integrals, which is, however, not the case. The reason is the additional symmetry of the one-scale integrals, related to the permutation of two massless legs. This symmetry reduces the number of one-scale master integrals to 71. Therefore, there are 116−71 = 45 redundant relations which we use as a check once we have satisfied 71 relations using explicit results for the one-scale master integrals. In practice, most of the one-scale master integrals have the same indices as JHEP02(2019)172 the corresponding two-scale master integrals so that the results for these one-scale master integrals are obtained directly from the naive part of the two-scale master integrals. For the remaining one-scale master integrals (where an index equal to two is chosen in another place), results are obtained after solving simple linear systems of equations.
Let us stress that the basic ideas of the described procedure have already been discussed in ref. [16], However, the approach presented here is more algorithmic and has now reached a state where it can be applied to highly non-trivial non-planar integral families, as it is demonstrated in this paper.
Note that in our case, we had to expand U 01 up to ǫ 9 (weight 9) for df2-2 and df2-3 since the property of uniform transcendentality is destroyed when mapping the two-scale master integrals to one-scale master integrals in the limit x → 0. In the final result for the form factor all weight-nine constants drop out. This happens separately for each family. In principle it is possible to adapt the basis of the one-scale master integrals such that only an expansion of U 01 up to ǫ 8 is necessary. However, our approach is powerful enough such that an expansion up to ǫ 9 did not pose any serious technical problems. For df2-5 and df2-6 an expansion up to weight eight is sufficient.
The reduction of one-scale as well as of two-scale integrals, needed for the derivation of differential equations for the (two-scale) master integrals, took several months for each of the four non-planar families. Using the standard version of FIRE we have failed to reduce the two-scale integrals of family df2-2 in the top sector. However, following the ideas of ref. [15], based on modular arithmetics, we managed to improve the performance of FIRE [45]. The new version can be used in a massive parallel mode on supercomputers which allows us to obtain the missing reductions.
In ref. [46] many (planar and non-planar) four-loop vertex integrals have been computed numerically. Among them are uniformly transcendental integrals in the top sectors of df2-2 and df2-3. Reducing these integrals to our primary bases and using our analytic results we can confirm the results (A.4)-(A.7) of ref. [46].
Let us finally mention that we have performed numerical cross checks of all master integrals of families df2-2, df2-3, df2-5 and df2-6 with up to ten positive indices expanded up to order ǫ 0 using FIESTA [47].
Analytic results for all master integrals can be downloaded in electronic form from [48]. For illustration we show for families df2-2 and df2-3 the master integrals with twelve lines in the appendix. Families df2-5 and df2-6 have no twelve-line master integrals.

Results
After inserting the analytic results for the master integrals into the amplitude for the form factor we observe that all poles higher than 1/ǫ 2 cancel. This is expected since the coefficients of the 1/ǫ 8 , . . . , 1/ǫ 3 poles are determined by lower-loop contributions. Since the colour structure (d abcd F ) 2 appears for the first time at four-loop order it can at most have 1/ǫ 2 poles. For the same reason there are no renormalization contributions to the (d abcd F ) 2 contribution.

JHEP02(2019)172
Our result for F (4) q (see eq. (1.2)) reads where N F = N c = 3 and ζ n is Riemann's zeta function evaluated at n. Note that the coefficient of the 1/ǫ 2 pole contains constants of at most weight 5 although in principle also weight 6 terms could appear. A similar weight-drop is also observed for the 1/ǫ term and the constant contribution. The cusp and collinear anomalous dimension can be extracted from the 1/ǫ 2 and 1/ǫ poles, respectively. For convenience of the reader we present the corresponding results separately. They are given by In refs. [23,24] the quark and gluon splitting functions at four-loop order have been considered. As a by-product numerical results for cusp anomalous dimensions have been obtained, in particular for C F γ 3 cusp | (d abcd For the case of a massive form factor, the corresponding cusp anomalous dimension Γ cusp is a function of the scalar product of the velocities of the two heavy quarks, v µ 1 and v µ 2 . Usually one introduces the variable φ = iϕ via cosh ϕ = v 1 · v 2 . It is interesting to note that eq. (3.2) corresponds to Γ cusp in the limit ϕ → ∞. Note that the limit ϕ → 0 has been considered in ref. [49].

Conclusions
We perform the next step towards the computation of massless four-loop form factors and compute the contribution of the quartic colour structure (d abcd F ) 2 to the photon-quark form

JHEP02(2019)172
factor. We have to consider two planar and four non-planar integral families which are shown in figure 2. We want to stress that this is the first time that master integrals with twelve propagators corresponding to non-planar graphs have to be considered. Our main results are shown in eqs. (3.1), (3.2) and (3.3). Furthermore, we provide analytic results for all master integrals in a supplementary file to this paper. We have used this calculation to further refine our method, which is used to obtain analytic results for the master integrals. The new element is the construction of the socalled associator which directly relates the coefficients in the boundary condition to the coefficients of the integrals in the physical limit. We are confident that the remaining contributions can be computed along the same lines. However, one has to keep in mind that much more families have to be considered and that the reductions to master integrals (both with one and two mass scales) require a significant amount of CPU time.

A Explicit results for twelve-line non-planar master integrals
In this appendix we present explicit results for the most complicated master integrals of the families df2-2 and df2-3 with twelve lines. We provide the ǫ expansion up to the constant term. Our results read The subscripts denote the exponents of the propagators, where the order is defined in figure 2. The six indices for the numerators are not shown; they are zero. Furthermore, we have s 8a = ζ 8 + ζ 5,3 ≈ 1.0417850291827918834 .
Note that all constants of weight 8 cancel in the combination of the master integrals which leads to the (d abcd F ) 2 part of the photon quark form factor, see eq. (3.1).