Quark and gluon form factors in four loop QCD: the $N_f^2$ and $N_{q\gamma} N_f$ contributions

We calculate the four-loop massless QCD corrections with two closed quark lines to quark and gluon form factors. The results for the gluon form factor and the singlet part of the quark form factor are given for the first time. From our analytic expressions for the form factors, we determine the corresponding cusp anomalous dimensions. The relevant Feynman integrals are obtained with novel integral reduction techniques and direct integration methods.


INTRODUCTION
In this Letter, we continue our ongoing study of the four-loop corrections in massless Quantum Chromodynamics (QCD) to the basic quark and gluon form factors for the photon-quark-antiquark vertex and the effective Higgs-gluon-gluon vertex. Despite a significant amount of recent attention, only partial results are available in the literature so far. The N 3 f contributions to the quark and gluon form factors with three closed fermion lines were calculated, respectively, in [1] and [2]. For the quark form factor, the contributions of order N 2 f were computed in [3] and the leading color limit was obtained in [1,4]. Very recently, the quartic Casimir color structure of the N f contributions to the quark form factor was calculated in [5] and analytic results for the corresponding contributions to the quark cusp anomalous dimension were also derived in a complementary approach in [6].
Quartic Casimir structures in the cusp anomalous dimensions are of particular interest since they violate the Casimir scaling principle proposed in references [7] and [8]. The existence of such Casimir scaling violations was demonstrated by explicit numerical calculation in N = 4 super Yang-Mills theory [6,9,10] and in QCD [11,12]. We wish to remark that, using approximation techniques [13][14][15], references [11] and [12] actually delivered complete numerical results for both the quark and gluon cusp anomalous dimensions and suggested a generalization of the Casimir scaling principle.
In this work, we substantially improve upon the current state-of-the art by providing complete analytic results for all contributions to the quark and gluon form factors with two closed fermion lines. While the N 2 f contributions to the quark form factor were already known, the singlet contributions to the quark form factor of order N qγ N f , where the photon couples to one of the two closed fermion lines, are new. We furthermore calculate the complete set of N 2 f contributions to the gluon form factor for the first time, including a determination of the N 2 f terms of the gluon cusp anomalous dimension. Due to the complexity of the relevant Feynman diagrams (see Fig. 1), it was necessary to build upon specialized computational technology studied by us in earlier works [16][17][18].
One of our main calculational challenges was the re-duction of four-loop, three-point Feynman integrals, for which we employ finite field-based techniques [16,19] implemented in the private program Finred. To evaluate the master integrals, we also make use of integrals with a significant number of dots (higher powers of the propagators). For their reduction, a newly developed algorithm based on the Lee-Pomeransky parametric representation [20] was used, building upon ideas discussed in [21,22]. The finite integral method of reference [18] was applied to calculate all of the 163 four-loop master integrals which remain after integration by parts reduction [16,[23][24][25]. Nineteen of those which we computed for the first time are of the more challenging non-planar type, free of massless one-loop bubble insertions. These include in particular master integrals for the first three top-level topologies shown in Fig. 1.
To carry out the calculation, we pass to a finite integral basis and analytically integrate the ǫ expansions of the resulting integrals starting from their Feynman parametric representations. For this purpose, we employ the program HyperInt [26]. This direct integration in terms of multiple polylogarithms requires the integrals to be linearly reducible [27,28], a criterion which was satisfied by most of our finite integrals. Although the first two integral topologies of Fig. 1 were not linearly reducible initially, we were able to make simple changes of variables which rendered them so. As will be discussed below, the integration is sometimes greatly simplified by making a judicious choice of finite integral basis.

SETUP AND INTEGRAL REDUCTION
We consider the perturbative amplitudes for the decays of photons and Higgs bosons into massless partons, γ * (q) → q(p 1 )q(p 2 ) and h(q) → g(p 1 )g(p 2 ), respectively, with p 2 1 = p 2 2 = 0 and q 2 = (p 1 +p 2 ) 2 . Interfering the bare amplitude with the tree amplitude, summing over polarizations and color, and normalizing to the corresponding tree-level expressions, we obtain the form factors where we take the L-loop massless QCD corrections into account. We work in conventional dimensional regularization with the bare strong coupling constant α bare s , the 't Hooft scale µ ǫ , Euler's constant γ E , and the parameter of dimensional regularization ǫ = (4 − d)/2. The amplitudes are calculated using a general R ξ gauge for the internal gluons, with up to one power of 1 − ξ, and arbitrary reference vectors for the external gluons. We denote the number of light quark flavors by N f and the chargeweighted sum of the N f quark flavors normalized to the charge of the external quark q by N qγ ≡ q ′ e q ′ /e q .
We generate the four-loop diagrams with the program QGraf [29] and consider the gauge-invariant subset with two closed fermion lines. We match the planar and nonplanar diagrams to nine complete sets of eighteen denominators (integral families) with Reduze 2 [30][31][32], where one such set of denominators may cover several twelveline top-level topologies and equivalent subtopologies are identified. In total, we encounter forty-six twelve-line top-level topologies, out of which twenty-two are nonplanar, see Fig. 1. We evaluate the color algebra for a general compact simple Lie group with Color.h [33] and evaluate the Lorentz and Dirac algebra with Form 4 [34].
We employ the in-house program Finred to reduce the resulting Feynman integrals to master integrals. For the reduction of the amplitude, we employ conventional momentum space integration-by-parts, Lorentz, and sector symmetry identities. For the basis change to our finite integrals, we employ first-and second-order annihilators [21,22] in the Lee-Pomeransky representation. Instead of resorting to an external computer algebra system, we calculate the required syzygies using linear algebra methods [35,36] with Finred as a linear solver.
Finite field sampling allows us to easily discard redundant equations, which reduces the number of equations and speeds up the reduction process considerably. Using different finite fields and different samples for d allows us to solve linear systems with 64 bit integers as coefficients and to reconstruct the rational functions from these finite field solutions. Robust vetos of fake reconstructions allow us to work without a priori assumptions concerning the required number of samples. We typically need O(10 1 ) finite fields and O(10 2 ) values for d for a successful reconstruction. We verify the reconstructed solution using five independent samples with unrelated values for d and the modulus. The reduction is run in a distributed manner and our final integral tables amount to several terabytes of compressed data.
Our unreduced amplitudes contain a total of 21286021 integrals with up to twelve propagators and six inverse propagators (irreducible numerators), for which we reconstructed a total of O(10 9 ) reduction identities in 1863 non-equivalent non-zero topologies, see Tab. I for morē  statistical data. After insertion of the reduction identities in the amplitudes, we find 163 irreducible master integrals and observe a non-trivial cancellation of all gauge parameter-dependent terms, something which constitutes an important sanity check on our calculation.

MASTER INTEGRALS
To evaluate the master integrals using direct analytical integration, we proceed as follows. We apply the algorithm of reference [17] to generate a reasonably long (over-complete) list of finite integrals in higher dimensions for a given sector. For each of these integrals, we evaluate the leading term by direct integration of the Feynman parametric representation using HyperInt. To avoid the evaluation of many terms in the ǫ expansion, it is often desirable to select a basis of finite integrals with high maximal transcendental weights at leading order in ǫ. Using dimensional recurrence [37,38] and dedicated reduction identities, we express the conventional master integrals in 4 − 2ǫ dimensions in terms of the new finite master integrals and their subsectors, which should be thought of as being known in a bottom-up approach.
One advantage of the finite master integrals constructed in this way is that they enter the ǫ expansion of amplitudes at higher orders than their conventional counterparts. Let us take, for example, the first integral topology of Fig. 1. In the conventional basis, we find which contributes to all terms in the ǫ expansions of the bare form factors through to O ǫ 0 . Note that, in Eq.
(2) and throughout this work, we adopt the conventions and definitions of reference [18] for our Feynman integrals and multiple zeta values. Now, if one replaces the above master integral with a suitable finite representative, one can make explicit that no evaluation of the new, finite master is required to extract the finite parts of the form factor contributions discussed in this work. This is simply because the finite parts of the contributions of order N 2 f and N qγ N f turn out to have a maximal weight of six (see Eqs. (6), (7), and (8) below), but the leading term in the ǫ expansion of the finite master integral for this sector may be consistently chosen to have a maximal weight of seven. This feature is particularly compelling for the example considered, because its finite integrals were only accessible to HyperInt after we made a change of variable of the form i /x ℓ for one of the Feynman parameters, x i . Let us elaborate further on what we require from our finite basis integrals. In most cases, it is very useful to pick finite integrals which faithfully map the weights of the underlying multiple zeta values when passing from the preferred finite basis back to the conventional basis. That is to say, we try whenever possible to ensure that a maximal weight of w at leading order in ǫ on the finite integral side implies complete weight w information on the conventional integral side. This is not automatic and is closely related to the problem of finding an integral basis for a given sector which is simultaneously finite and uniform weight [39]. As a dramatic example, consider the conventional master integrals  With the master integrals above, the reduced integrand for the N 2 f four-loop gluon form factor would actually require the unknown O (ǫ) terms of Eqs. (3) and (4) to extract final results through to O ǫ 0 . On the other hand, a suitable finite integral basis including offers the possibility to completely avoid evaluating anything but the weight six, leading-order term of Eq. (5).
For the majority of the 151 irreducible topologies relevant to the N 2 f gluon form factor, we were able to find finite basis integrals which allowed us to avoid computing spurious orders of their ǫ expansions.

RESULTS AND CROSS-CHECKS
With the techniques described above, we find for the non-singlet quark form factor F q 4 (ǫ) and for the gluon form factor The color factors above are defined in Eqs. (185) -(188) of reference [33] for a SU (N c ) gauge group (T F = 1/2). Our cross-checks were as follows. First of all, we observed that the gauge-dependent parameters of our calculation dropped out of our final results. Most of our master integrals were computed twice using different choices for the finite integrals and it was gratifying to see that we could often produce uniform weight integrals by mapping our finite integrals back to conventional ones (see e.g. Eq. (3)). Due to the simplicity of finite Feynman integrals [40], we found that it was possible to check all our non-planar master integrals to a relative precision of 10 −4 numerically using FIESTA 4 [41]. We were also able to successfully check some of our integrals analytically against the results of [3,5]. Furthermore, we agreed with the higher-order pole predictions of [42] and the known result for the N 2 f part of the bare quark form factor [3]. Finally, we extracted the cusp anomalous dimensions from the ǫ −2 poles: As predicted by the Wilson loop picture, the color structure C 2 F N 2 f in the gluon form factor enters the cusp with zero coefficient, as does the quartic Casimir term. In fact, we have confirmed by direct calculation that Casimir scaling holds for the N 2 f contributions. Eq. (9) is in agreement with the analytic result of reference [43] and Eq. (11) represents the first direct extraction of Γ g 4 | N 2 f .

SUMMARY
In this Letter, we calculated previously-unknown contributions to four-loop quark and gluon form factors with two closed fermion loops. We also derived the corresponding contributions to the cusp anomalous dimensions. The master integrals were computed by direct integration of finite integrals in the Feynman parametric representation. For two integral topologies, we had to perform simple changes of variables to render the representation linearly reducible and thus accessible to an integration with HyperInt.