Matter dependence of the four-loop cusp anomalous dimension

We compute analytically the matter-dependent contributions to the quartic Casimir term of the four-loop light-like cusp anomalous dimension in QCD, with $n_f$ fermion flavours and $n_s$ scalar flavours. The result is extracted from the double pole of a form factor. The calculation is significantly simplified by our choice of basis of master integrals. We first identify, for each topology, an over-complete set of integrals of uniform transcendental weight, by requiring that their integrands can be put in a dlog form. We then perform a systematic analysis of their soft and collinear singularities, in order to build linear combinations of integrals with simpler pole structure. With this information, we build a dlog basis such that only integrals with ten or fewer propagators contribute to our result for the cusp anomalous dimension. These integrals are then computed via the method of differential equations, through the addition of an auxiliary scale. Combining our result with that of a parallel paper, we obtain the complete $n_{f}$ dependence of the four-loop cusp anomalous dimension. Finally, using known numerical results for the gluonic contributions, we obtain an improved numerical prediction for the cusp anomalous dimension in $\mathcal{N}=4$ super Yang-Mills.


INTRODUCTION
The cusp anomalous dimension is a universal quantity appearing in QCD. It governs infrared divergences of scattering amplitudes [1][2][3], and appears in the large spin limit of twist-two operators [4]. It also controls the resummation of large logarithms due to soft radiation and is therefore relevant to many observables [5][6][7][8][9]. Its colour dependence is governed by non-Abelian exponentiation, which allows for the first time for quartic Casimir terms at four loops. The latter have received a lot of attention, in particular because their presence implies a breaking of the Casimir scaling property, and since they represent the last missing four-loop ingredient in the above calculations. Further interest comes from the fact that these are the first truly non-planar terms in N = 4 super Yang-Mills (sYM). In this theory, the planar cusp anomalous dimension is known from integrability [10], and it remains an open question whether integrability extends to the non-planar sector.
The planar cusp anomalous dimension is known to four loops in QCD [11], and some of the n f dependent terms have been computed [12][13][14][15][16]. Recently, numerical results were obtained for the quartic Casimir terms in N = 4 sYM [17,18] and in QCD [19]. In this Letter, we present the first analytic result for the n f terms in QCD.
Given the complexity of such a non-planar four-loop calculation, we develop and use cutting-edge methods to achieve this goal. The latter may be of interest in their own right, as we expect they can be applied to many other situations.
We use as our starting point a form factor of composite operators inserted into two on-shell states. Thanks to the universality of the cusp anomalous dimension, we are free to choose a suitable operator, and we make a particularly simple choice, as explained below. The kinematic dependence is fixed by dimensional analysis, so that the form factor essentially depends on , the parameter of dimensional regularization in D = 4 − 2 dimensions, only.
In recent years, it has become standard to make an educated choice of Feynman integral basis [20][21][22], where the integrals are of uniform transcendental weight (UT), or so-called pure functions. A given L-loop Feynman integral with this property has the expansion I pure = −2L k c k k , where the c k are numbers of transcendental weight k. This property is particularly useful in N = 4 sYM where, conjecturally, the form factors have uniform and maximal weight. While in general the form factors are expressed as F = i r i ( )I i pure , with some rational functions r i ( ), in the latter theory the r i are just numbers, i.e. -independent. Note that this property only becomes visible when a basis of pure functions is chosen. One of the first applications of these ideas was at the level of the three-loop form factor [23]. We argue that such a basis choice will be of crucial importance also in QCD. Experience from lower loops shows that terms having at least one factor of n f have a drop of transcendental weight. In a UT basis, this means that all coefficients have the form r i ( ) = q i ( ). Making this property manifest allows us to take a calculational shortcut.
The quartic Casimir terms appear for the first time at four loops, and as a consequence of renormalizability, they come with a 1/ 2 pole (whose coefficient is the cusp anomalous dimension). Thanks to the additional factor of mentioned above, we need to know the four-loop integrals only up to (and including) the 1/ 3 pole. In order to take advantage of this fact, we classify the pure functions according to their soft and collinear divergence properties [20,24,25]. In this way, we can arrange complicated integrals (typically, integrals having many propagators) into linear combinations having only 1/ 2 or better pole structure, and hence irrelevant for the determination of the cusp anomalous dimension. In this way, only a subset of form factor integrals is needed.
Expressions for all planar four-loop form factor integrals were obtained previously [11] by an application of the differential equations method [21,26]. Here, we evaluate all required non-planar integrals using the same method. We validate these results, whenever possible, by comparing to numerical integration [27], and as crosschecks verify the expected infrared structure, as predicted by our integrand analysis. We then use these results to expand the form factor up to the 1/ 2 pole, and extract the cusp anomalous dimension analytically.

SETUP AND DEFINITIONS
We work in massless QCD with gauge group SU (N c ) and n f fermion flavors. For convenience, we couple the theory canonically, i.e. through covariant derivatives, to n s complex scalar fields, with canonical kinetic term φ φ . This allows us to consider a composite operator O = φφ, inserted into on-shell scalar states, i.e. with p 2 1 = p 2 2 = 0, Here the scalar fields are considered to be in the representation R of SU (N c ), which we take either to be the fundamental (F), or adjoint (A). In the following, we will set the only kinematic scale 2p 1 · p 2 = −1, and the dimensional regularization scale µ 2 = 1, without loss of generality. The fact that O has spin zero means that in momentum space, no additional momentum operator is inserted into the diagram at the cusp. As a consequence, the corresponding Feynman diagrams contain one numerator factor less compared to what one would have obtained for a (spin-one) fermion current O µ =Ψγ µ Φ, for example. This turns out to be an important technical simplification.
The cusp anomalous dimension is universal, that means it does not depend on the types of external particles, in this case scalars. There will be a difference at the subleading order in , however, for two reasons: the φ 2 operator has a UV anomalous dimension, and the collinear anomalous dimension depends in general on the particle type.
We are interested in the four-loop contribution of F to the quartic Casimir structure [28] where d abcd R , R denotes the SU (N c ) representation of the external scalars and X the representation of the internal matter fields (n f fermions and n s scalars). The quartic Casimir n f and n s contributions originate from a small set of four-loop Feynman diagrams with an internal fermion box, as shown in Fig. 1, and internal scalar box, triangle, and bubble subdiagrams, respectively. There is also a corresponding gluonic quartic Casimir term, which however is beyond the scope of the present paper.
The general structure of infrared divergences of the form factor, together with the fact that the quartic terms appear for the first time at this loop order, implies that where α s is the strong coupling. Our goal is to determine K 4 | n f , ns quartic . We perform the calculation in a general covariant gauge with parameter ξ, and verify that the linear terms in ξ disappear from the result.

INTEGRAL REDUCTION
The form factor F | n f , ns quartic is expressed in terms of scalar Feynman integrals. A first step in its calculation consists in exploiting integration-by-parts identities (IBP) [29] in order to reduce the expression to a minimal number of so-called master integrals (MI). This step is an important bottleneck in our calculation. We find that even on state of the art compute servers, publicly available IBP reduction programmes run into difficulties.
For this reason, we use novel techniques pioneered in [30,31], based on reconstructing analytic results from numerical evaluations over finite fields. More in detail, we generate the system of IBP equations and mappings for each topology, with the help of LiteRed [32]. We then solve it modulo prime integers p, for several numerical values of and p using a custom generic linear solver for sparse systems over finite fields. In this way we avoid introducing complex expressions at intermediate stages of the calculation. The full analytic result is then reconstructed using the techniques illustrated in [31]. We also identify equivalences between integrals appearing in different integral families, such that the result is expressed in terms of a minimal number of master integrals for the combined set of integral families.
For the reduction of the diagrams contributing to the form factors, rather than producing large and complex IBP tables, we directly apply the numerical reduction to the diagram and only reconstruct the reduced final expression in terms of a minimal set of master integrals. We make a special basis choice, where the coefficients are expected to be relatively simple. We then add this definition of the basis integrals to the system of equations  and select them as preferred master integrals during its solution, thereby obtaining directly the simpler integral coefficients. The next two sections are dedicated to explaining how an educated basis choice is made. Table I gives an overview of the number of MI for each of the integral families. The second to fourth column state the number of MI per family, grouped according to the number of propagators. Σ is the total number of MI of the corresponding family, while Σ * is the number of MI excluding all integrals that can be related to integrals of an integral family previously considered (i.e. that appears above in the same table). We ordered the families such that the first two families are the planar topologies and then we have the non-planar topologies. In total, adding all entries for Σ, we have 272 MI. After considering all relations between MI this number is reduced to 100 (the sum of all entries in the Σ * column).

MASTER INTEGRALS AND PURE FUNCTIONS
As mentioned in the introduction, it is advantageous to select a basis of pure integrals. Conjecturally, the latter can be identified by checking that their four-dimensional loop integrands can be put into a so-called dlog form [33].
(See also [34] for recent developments on the topic of identifying UT integrals.) We systematically find such integrals using the algorithm [22]. The last column of Table I gives the number of dlog integrals we found in the different families, based on an ansatz with heuristic power counting constraints. We find that it is possible to choose a subset of these dlog integrals that can be used as a complete basis of MI.
Using this basis (denoted by f ) to express the result, we find where the C(n f , n s ) denotes the overall normalization for the case of internal fermions or scalars and the q i are IBP coefficients of O( 0 ). Remarkably, all integral coefficients are proportional to ! This confirms our expectation that F | n f ,ns quartic has a transcendental weight drop. Note that it is essential to use a basis of pure functions to observe this property prior to computing the integrals.

CHOOSING AN INTEGRAL BASIS WITH GOOD INFRARED PROPERTIES
The pure functions we found in the previous section may have several (in general, nested) regions of soft and collinear divergence, due to the on-shell light-like kinematics. At four loops, these regions lead to poles of up to −8 . On the other hand, we expect the quartic Casimir contribution to be given by a −2 pole only, see eq. (3).
As this term appears for the first time in perturbation theory, its singular behaviour must emerge in the sum of the four-loop diagrams shown in Fig. 1, without the need of exponentiating lower-loop diagrams, as would be the case for other color structures. Indeed, the physical origin of the two poles is an overall (i.e., involving all loop momenta) soft and collinear divergence, respectively. This motivates the question whether the infrared structure of the four-loop integrand can be made manifest. Similar questions were pursued in [20,24,25,35] for scattering amplitudes and correlation functions (see also [36]). The physical idea is to look for linear combinations of loop integrands for which certain soft and collinear regions (and hence the associated divergences) are suppressed. In the above references, integrands without one-loop soft/collinear subdivergences were constructed. Here, we perform a dedicated analysis, analysing all Lloop soft/collinear regions of loop integration algorithmically. This information then allows us to construct loop integrals for which we can give an upper bound on the degree of divergence.
In the following we briefly sketch the basic ideas behind our IR analysis algorithm. In order to test the region, where the loop momentum k i becomes collinear to p 1 (γ 1,i → 0) and/or soft (β i → 0), we parameterize with p 1 ·k ⊥i = p 2 ·k ⊥i = 0, and analogously for k||p 2 . We also consider consecutive p 1 -and p 2 -collinear limits of k i (γ 1,i → 0 and γ 2,i → 0, respectively) using the parametrization Choosing region momenta for planar integrals, we can now take soft and collinear limits of each loop momentum separately in arbitrary order. We do this by Laurent expanding the integrand in the soft and collinear parameters β i and γ 1/2i , respectively. If we find a single pole of the form dβ i /β i or dγ 1/2i /γ 1/2i we conclude that the corresponding limit (potentially) contributes a 1/ pole to the integral. We then proceed with the residue of this pole and test the next limit, and so on. Note that the dlog property gaurantees that we never encounter more that single poles in this procedure.
We also test simultaneous soft or collinear limits involving several loop momenta. For example, consider the case where the loop momenta k 1 , k 2 and k 3 become soft simultaneously. We then set β i = β 123 b i for i = 1, 2, 3 and send β 123 → 0 in order to look for a dβ 123 /β 123 pole. An analogous procedure is performed for joint collinear limits. Our code systematically checks all consecutive (joint) soft/collinear limits of the k i and keeps track of the origin of singularities found in this way.
As an example, consider the Feynman integral shown in Fig. 2. According to our algorithm, the maximal singular behaviour comes from first taking the joint p 1collinear limit of loop momenta {1, 2, 4} (γ 1,124 → 0), then γ 1,3 → 0 followed by the joint limit γ 2,234 → 0 and finally γ 2,1 → 0. (Obviously, the analogous behaviour is found in the p 1 ↔ p 2 symmetric version of the limits.) Hence, we expect the maximal degree of divergence to be 1/ 4 . This is indeed confirmed by the available analytic result [11]: π 4 /(5184 4 ) + O( −3 ). The above analysis was facilitated by the existence of a canonical momentum routing for planar integrands. For nonplanar integrals we therefore adopt a pragmatic approach and run the algorithm for all momentum routings where four of the twelve propagators of the nonplanar topologies coincide with 1/k 2 i . In general, our method gives an upper bound on the degree of divergence of an integral only, as there may be cancellations, or some regions may give zero contributions, e.g. due to scale-less integrals. Note that we make the physical assumption that only soft and collinear regions are relevant to this analysis, so that the potential presence of other scaling regions could alter the conclusions. For all integrals that we explicitly computed, we verified that our bound was satisfied, thereby validating the procedure.
As was explained above, we expect the sum of all Feynman diagrams to have one overall soft and one overall collinear divergence only. On general grounds one expects this property to appear in the sum of all diagrams, so that it is a priori not clear to what degree one can make this property manifest at the level of individual basis elements, within the different integral families. Here, we will profit from the insights on the infrared properties of the basis integrals in another way.
We use the additional information on the infrared behaviour of the individual integrals to assemble a dlog basis where integrals having more than 10 propagators have at most 1/ 2 poles. As all integrals appearing in F with a coefficient that is O( ), cf. eq. (4), this implies that we do not need integrals more than 10 propagators in order to extract K 4 | n f , ns quartic . We compute these integrals in the next section.

COMPUTATION OF MASTER INTEGRALS
In recent years, the method of computing Feynman integrals via differential equations in canonical form [21] has become very popular. The scale dependence of the FIG. 3: Typical form factor integral, which we compute for p 2 1 = 0. We then extract its value at p 2 1 = 0.
form factor integrals is trivial, so that a direct application of the method is not possible. However, one can use the trick employed in [26] (and in subsequent work [12,15]) of introducing another scale, e.g. p 2 1 = 0. One then sets up the canonical differential equations for the form factors integrals with two off-shell legs. Note that as p 2 → 0, the integrals degenerate to propagator-type integrals that are all known. Therefore one can use the differential equations to determine the desired on-shell form factor integrals by relating them to the known propagator-type integrals [37].
We follow this approach and obtain all required integrals analytically, up to transcendental weight 8, which is sufficient for a four-loop computation. The necessary integral reductions, this time for form factor integrals with two off-shell legs, are again performed using the algorithm [31]. It is interesting to note that, as a consequence of the simple form of the differential equations [26], it follows that the equations relating the results of the propagator-type integrals to our form factor integrals involve involve only harmonic polylogarithms [38,39] with indices 0 and 1, evaluated at 1. As a result, only multiple zeta values appear in these equations, to any order in . The results are provided in ancillary files.
As an example, consider the integral shown in Fig. 3. We find that in the on-shell case p 2 1 = p 2 2 = 0, it is given by This agrees with our infrared analysis.

MAIN RESULTS
With the result of the previous sections, we can expand our result for the form factor to order −2 , and obtain K 4 | n f ,ns quartic . For the n f results quoted below we assume the fermions to live in the fundamental representation (X = F ), for the n s and N = 4 sYM results we put all fields in the adjoint representation (R = X = A). We find and The n f term is in perfect agreement with the numerical result of [16,19].
We can now use our novel analytic result in eq. (8), together with the analytic expression for the planar n f term in K 4 from [12], and a conjectured new result for the n f T F C R C F C A term from a parallel paper [40], to obtain the complete analytic (linear) n f term of the light-like QCD cusp anomalous dimension: Next, we may use the numerical result of [19] for the purely gluonic quartic Casimir term, together with the analytic matter contributions computed here, to obtain the result in N = 4 sYM, This agrees perfectly with the result of [17], and improves the numerical precision by two decimal places.

DISCUSSION AND OUTLOOK
In this Letter, we computed the matter-dependent contributions to the quartic Casimir term of the four-loop light-like cusp anomalous dimension in QCD. Combining this with other results, we obtained the full four-loop n f dependence of the cusp anomalous dimension in QCD. We also obtained a more precise numerical result for the cusp anomalous dimension in N = 4 sYM.
Our calculation was considerably simplified by using a basis of master integrals of uniform transcendental weight with improved soft and collinear properties. In this way, we did not require any master integrals with 11 or more propagators.
Let us comment on the extension of this idea to more general types of integrals. We note that while for the calculation of form factor integrals, only soft and collinear scaling regions are expected to be relevant, in general, one may have to take into account other regions as well (such as Glauber regions). In principle, one can verify the predictions of this analysis by sector decomposition [41]. This check can be performed purely analytically and it does not require to perform the numerical integration step of the sector decomposition method implemented in [27,42]. Finally, while we mostly benefitted from the analytic improvements of our integral basis, we expect that the integrals constructed in this way may also be more stable numerically [17,18].