The full four-loop cusp anomalous dimension in $\mathcal{N}=4$ super Yang-Mills and QCD

We present the complete formula for the cusp anomalous dimension at four loops in QCD and in maximally supersymmetric Yang-Mills. In the latter theory it is given by \begin{equation} {\Gamma}^{\rm}_{\rm cusp}\Big|_{\alpha_s^4} = -\left( \frac{\alpha_s N}{\pi}\right)^4 \left[ \frac{73 \pi^6}{20160} + \frac{ \zeta_{3}^2}{8} + \frac{1}{N^2} \left( \frac{31\pi^6}{5040} + \frac{9 \zeta_3^2}{4} \right) \right] \,.\nonumber \end{equation} Our approach is based on computing the correlation function of a rectangular light-like Wilson loop with a Lagrangian insertion, normalized by the expectation value of the Wilson loop. In maximally supersymmetric Yang-Mills theory, this ratio is a finite function of a cross-ratio and the coupling constant. We compute it to three loops, including the full colour dependence. Integrating over the position of the Lagrangian insertion gives the four-loop Wilson loop. We extract its leading divergence, which determines the four-loop cusp anomalous dimension. Finally, we employ a supersymmetric decomposition to derive the last missing ingredient in the corresponding QCD result.

Abstract: We present the complete formula for the cusp anomalous dimension at four loops in QCD and in maximally supersymmetric Yang-Mills. In the latter theory it is given by

Introduction
The cusp anomalous dimension is an important quantity in four-dimensional Yang-Mills theories ranging from QCD to maximally supersymmetric N = 4 Yang-Mills theory. It controls the leading ultraviolet divergences of Wilson loops evaluated along a closed contour in Minkowski space-time containing cusps formed by two light-like tangent vectors [1,2]. Two examples of such contours relevant for our discussion are a wedge formed by two semi-infinite lines and a null polygon with edges along different light-cone directions. The former example plays an important role in the study of infrared asymptotics of onshell scattering amplitudes and form factors [3][4][5]. It also naturally appears in the analysis of DGLAP splitting functions in the semi-inclusive limit [6,7], and in the resummation of large Sudakov corrections due to soft and collinear emissions [8,9]. In all these cases, the contribution of soft particles is in a one-to-one correspondence with ultraviolet divergences of semi-infinite cusped Wilson loops [3], allowing us to find the asymptotic behaviour of the corresponding physical quantities in terms of the cusp anomalous dimension [10].
The study of infrared and collinear singularities in Yang-Mills theory is of considerable interest (see for example refs. [11][12][13][14]) and the cusp anomalous dimension plays a key role in it. In particular, the four-loop cusp anomalous dimension is the only missing contribution to understand the double pole in the dimensional regulator of four-loop scattering amplitudes. Currently, form factors relevant for the precise determination of hadron collider observables such as the Drell-Yan or Higgs boson production cross section [15,16] are known to third loop order [17] and the advent of the next order is apparent [18][19][20][21]. Knowing the four-loop cusp anomalous dimension, and consequently more of the singularity structure of these quantities serves as stringent check of such a computation. In combination with the understanding of how scattering amplitudes factorise as particles approach kinematic limits the understanding of infrared singularities allows to resum parts of the perturbative expansion to all orders. The four-loop cusp anomalous dimension is a necessary ingredient for resummation at next-to-next-to-next-to-leading logarithmic accuracy. Currently, such computations exist with numerical approximations for the four-loop cusp anomalous dimension for a variety of observables. This applications range from the extraction of the strong coupling constant from e + e − event shapes [22,23] to precision observables at hadron colliders like the transverse momentum distribution of electro-weak gauge bosons [24][25][26][27].
Furthermore, null polygon Wilson loops and the cusp anomalous dimension are of special interest in N = 4 super Yang-Mills theory (sYM). This theory has a number of remarkable properties, e.g. the celebrated AdS/CFT correspondence, and is believed to be integrable, at least in the planar limit [28]. Integrability has been successfully exploited to predict the cusp anomalous dimension in planar N = 4 sYM theory for any value of the 't Hooft coupling λ = g 2 N [29]. Extending the integrability approach to predicting non-planar corrections to the cusp anomalous dimension is currently under active investigation.
Another remarkable property of N = 4 sYM theory is that light-like Wilson loops describe the asymptotic behaviour of off-shell correlation functions in the light-like limit, when the operators approach the position of vertices of null polygon, and they are dual in the planar limit to the so-called MHV on-shell scattering amplitudes. We shall use this relationship below.
At present, the cusp anomalous dimension is known in full in QCD at three loops [2,30]. At four loops, it has been computed analytically up to one special term that we specify presently. The cusp anomalous dimension is known to have the so-called Casimir scaling up to three loops [6]. Namely, the dependence of Γ cusp on the representation of the SU (N ) gauge group only enters through a quadratic Casimir. Starting from four loops this property is violated due to the appearance of a new SU (N ) color factor -the quartic Casimir built out of two completely symmetric invariant d−tensors with all indices contracted [31,32].
As mentioned above, in N = 4 sYM theory the cusp anomalous dimension is known in the planar limit to any loop order. The non-planar correction first appears at four loops and it arises precisely from the quartic Casimir mentioned above. At large N , it is given by the sum of two terms with one of them suppressed by the factor of 1/N 2 . It is this term that generates a nonplanar correction to the cusp anomalous dimension in N = 4 sYM theory at four loops. 1 The new color factor arises from pure gluonic diagrams. Their contribution to Γ cusp is not sensitive to the matter content of the theory and, therefore, is the same in N = 4 sYM theory and in QCD. Thus, computing the four-loop correction to the cusp anomalous dimension proportional to this color factor in the former theory, and adjusting the color factors, we can predict the analogous contribution in QCD. This is the only missing term in the existing four-loop expression for Γ cusp in QCD mentioned above.
In this paper, we compute for the first time analytically the non-planar correction to the four-loop cusp anomalous dimension, both in N = 4 sYM and in QCD. We do so using an approach based on a relationship of Wilson loops with correlation functions of local operators. It allows us to avoid complicated Feynman graph calculations. We validate our results by comparing them against previously known analytic (in the planar case) and numerical (in the non-planar case) values, finding perfect agreement. Our result for the four-loop cusp anomalous dimension in N = 4 sYM in the adjoint representation of SU (N ) is where α s = g 2 /(4π) is the fine structure constant. The corresponding QCD result can be found in eqs. (6.2) and (6.3) below. The outline of the paper is as follows. In section 2, we begin by reviewing correlation functions in N = 4 super Yang-Mills and their relation to light-like polygonal Wilson loops. We focus in particular on the case of the latter with a Lagrangian insertion, and review their relationship to the cusp anomalous dimension. Finally, we discuss the integrand, F , for the three-loop correlation function of the null rectangular Wilson loop with a Lagrangian f for a scalar or fermion in the loop, respectively; diagram (b) contributes to Γ (4) ; diagram (c) contributes to Γ (4) and Γ (4) g .
insertion. Section 3 is dedicated to the analytic calculation of the relevant three-loop Feynman integrals. After reviewing the definitions in 3.1, we explain in section 3.2 our method for choosing a basis of integrals possessing simple transcendental weight properties. We make use of a convenient choice of loop variables to simplify this analysis. In section 3.3 we apply the differential equations method to compute all basis integrals. We present the integrated results for F to three loops in section 4. In section 5, we derive general formulae for performing the integration over the Lagrangian insertion, to extract the cusp anomalous dimension at the next loop order. Section 6 contains the main results of this paper -the expressions for the four-loop cusp anomalous dimension in N = 4 super Yang-Mills and in QCD, for an arbitrary representation of the Wilson line. Finally, we conclude and give an outlook in section 7.

Cusp anomalous dimension from a correlation function
The connection between infrared asymptotics of on-shell scattering amplitudes and form factors and ultraviolet divergences of semi-infinite cusped Wilson loops has been previously used to obtain the cusp anomalous dimension. In this paper, we follow another approach to computing the cusp anomalous dimension that relies on the relation between off-shell correlation functions of local operators and light-like Wilson loops.

The quartic Casimir terms in different gauge theories
As was mentioned above, the cusp anomalous dimension governs the ultraviolet divergences of light-like Wilson loops. The simplest example of such an object is a null rectangular Wilson loop where the contour C is a rectangle with vertices located at four points x i that are light-like separated, x 2 i,i+1 = 0 (with i = 1, . . . , 4 and i + 4 ≡ i). We took the representation to be the adjoint, R = A, the reason for such a choice will be clear in a moment. We normalized the Wilson loop by N A = N 2 − 1, so that its perturbative expansion starts with 1. The Wilson loop (2.1) has ultraviolet (UV) divergences due to the presence of the four cusps. In dimensional regularization, with D = 4 − 2ǫ, its leading divergence is [34,35] cusp,A are expansion coefficients of the cusp anomalous dimension in the adjoint representation 3) The relation (2.2) holds in a conformal Yang-Mills theory, otherwise it is valid up to terms proportional to the beta function.
In a generic Yang-Mills theory with SU (N ) gauge group, containing n f fermions and n s scalars, the cusp anomalous dimension takes the following general form at four loops where the subscript R refers to the SU (N ) representation in which the Wilson lines are defined, and N R = tr R 1 is the dimension of the representation. In the case of QCD and N = 4 sYM the relevant representations are fundamental (R = F ) and adjoint (R = A). The four terms inside the brackets in (2.4) contain four different color factors depending on R. Sample diagrams contributing to these terms are shown in Figure 1. The first term in (2.4) is proportional to the quadratic Casimir C R = T a R T a R . The corresponding expansion coefficient Γ (4) is independent of R. It only depends on the quadratic Casimir in the adjoint representation, C A = N , as well as on the number of fermions (n f ) and scalars (n s ) and on the quadratic Casimirs in the representations in which these particles are defined. This proportionality property is usually referred to as the Casimir scaling. The coefficient function Γ (4) has been computed in refs. [18,19].
The three remaining terms in (2.4) involve quartic Casimirs built out of completely symmetric tensors in different representations. The above sum runs over all permutations of color indices. Such color factors appear first at four loops and they violate the Casimir scaling. The last two terms in (2.4) come from Feynman diagrams containing fermion and scalar loops coupled to four gluons. The corresponding expansion coefficients Γ (4) f and Γ (4) s have been 2 Here the additional factor of 2 is inserted to take into account that the Wilson loop is defined in the adjoint representation (and not in the fundamental representation as in the mentioned papers). computed in refs. [36,37]. The main result of this article is the analytic calculation of the last missing coefficient Γ (4) g . Numerical results for this quantity were obtained in refs. [37][38][39][40]. Since the matter dependence is known, the gluonic coefficient can be computed in maximally supersymmetric Yang-Mills theory. 3 This theory has a number of remarkable properties that simplify the calculation significantly.

The cusp anomalous dimension from a finite ratio of Wilson loops
In principle, one could attempt to compute W A (x 1 , x 2 , x 3 , x 4 ) in N = 4 sYM directly in perturbation theory in order to extract Γ cusp,A . The main complications would be a proliferation of Feynman diagrams and the evaluation of the corresponding complicated four-loop Feynman integrals. Instead of attempting this direct calculation, we use a method that avoids to a large degree the need for Feynman diagrams, and that only requires the evaluation of a finite three-loop quantity.
The key insight is that we do not need to evaluate log W A (x 1 , x 2 , x 3 , x 4 ) fully, as the cusp anomalous dimension appears in (2.2) as the leading double pole in dimensional regularization. We can imagine log W A (x 1 , x 2 , x 3 , x 4 ) as being represented by multiple integrals. Physically, it is clear that the cusp divergences arise from the integration over particles that propagate at short distances along the light-like edges adjacent to the cusps. The idea is that, upon rescaling the distances, we can isolate a divergent integral over the overall scale and express the cusp anomalous dimension as the finite result of the remaining integration. This is similar in spirit to [42]. Here we follow the approach of [43][44][45].
Let us consider a slightly more general object, namely the correlation function of the Wilson loop with the insertion of the Lagrangian normalized by W A [43,46,47] (see footnote 2) It will serve us as the finite integrand mentioned above. It is possible to show following [35] that UV divergences cancel in the ratio of correlation functions on the left-hand side of (2.6), so that it remains finite for ǫ → 0. As a consequence, the first term on the right-hand side of (2.6) can be defined at D = 4. Conformal symmetry fixes it up to an arbitrary function F (x) of cross-ratio [43] x = It depends on four points x 1 , x 2 , x 3 , x 4 defining vertices of a null rectangle and the Lagrangian insertion point x 5 . The invariance of the Wilson loop under cyclic permutation of points x 1 , . . . , x 4 leads to the relation F (x) = F (1/x). In addition, F (x) depends on the rank of the gauge group N , and on the Yang-Mills coupling g 2 . As customary in the QCD literature 4 , we expand F (x) in powers of the fine structure constant α s /π = g 2 /(4π 2 ), (2.8) Note that the expansion starts at order α s . The function F (0) defines the correlator (2.6) at Born level and the functions F (L) describe corrections at L loops.
Let us now explain how we can combine eqs. (2.6) and (2.2) to obtain the cusp anomalous dimension from F . On the one hand, we observe that once we integrate the l.h.s. of eq. (2.6), this yields a derivative of log W A (x 1 , x 2 , x 3 , x 4 ) w.r.t. the coupling. On the other hand, comparing to the r.h.s. of eq. (2.6), we see that the leading divergence of that quantity is controlled by the cusp anomalous dimension: As a consequence, once we know the finite ratio in eq. (2.6) at (L − 1) loops, integrating over x 5 and extracting the leading divergence allows us to compute the cusp anomalous dimension at L loops. The detailed derivation was performed in ref. [45] using two different regularizations, with the same result. The effect of the integration over x 5 is given by a functional I relating the function F (x) to the cusp anomalous dimension. We have where I acts on individual terms as in In order to apply the last relation we write F (x) in the form of a small x expansion.
The reason why small x asymptotics is related to the cusp anomalous dimension can be understood from eq. (2.9) -the cusp singularities arise from x 5 approaching the cusp points x i . Using the definition (2.7), it is easy to see that in this limit x → 0 or x → ∞. Due to the symmetry x → 1/x we can map both regions to small x. Eq. (2.10) allows us to determine Γ cusp at L loops from F (x) at (L − 1) loops, cf. eq. (2.8). This reduces the complexity of the calculation considerably.

The Wilson loop integrand from a correlation function
Our goal is therefore to obtain the function F (x) at three loops. In the conventional Feynman diagram approach one would compute W A L(x 5 ) and W A in dimensional regularization and then find their ratio (2.6). Here we follow a different strategy that is based on a relationship between correlation functions and light-like Wilson loops. It will allow us to obtain the integrand for F (x), bypassing Feynman diagram computations.
Our starting point is the four-point correlation function of half-BPS operators in N = 4 sYM theory is bilinear in scalar fields Φ I (with I = 1, . . . , 6) defined in the adjoint representation of the SU (N ). The operator's scaling dimension ∆ = 2 is protected from quantum corrections by supersymmetry.
At weak coupling, the perturbative expansion of G 4 in powers of the coupling α s can be obtained using the method of Lagrangian insertions. This method relies on the observation that a derivative of an n-point correlation function with respect to the coupling is given by an (n + 1)-point correlation function involving the insertion of Lagrangian Differentiating repeatedly with respect to the coupling, we can express the L−loop correction to G 4 in terms of (4 + L)-point correlation functions involving L insertions of N = 4 sYM Lagrangians Here the superscript '(0)' was introduced to indicate that the correlation functions are computed in the Born approximation. The relation (2.14) allows us to interpret G 4,L as Feynman integrands at order g 2L . They are rational functions of the distances . The form of the functions is heavily constrained by the symmetries of N = 4 sYM theory. Namely, the fact that the half-BPS operators O and the Lagrangian L belong to the same supermultiplet leads to a hidden symmetry of (2.14) under permutation of points x 1 , . . . , x 4 and x 5 , . . . , x 4+L [48]. At four loops, this symmetry, combined with conformal symmetry and the requirement for (2.14) to have correct behavior in various OPE limits, fixes the leading, planar part of G (0) 4,L uniquely whereas the non-planar part of the integrand can be determined up to four arbitrary coefficients [49].
These coefficients were determined in the recent paper [50] by matching the non-planar part of G (0) 4,L=4 computed in ref. [51] using the reformulation of N = 4 sYM in twistor space to the analogous expression for the same correlation function obtained in refs. [48,49]. Thus, the construction of the four-loop integrand of the correlation function (2.12) is now completed.
To make a connection to the Wilson loops considered in the previous subsection, we examine the limit x 2 12 , x 2 23 , x 2 34 , x 2 41 → 0 when four operators in (2.12) become light-like separated in a sequential manner. As was shown in ref. [52], the correlation function simplifies in this limit and its leading asymptotic behavior is given by the product of Bornlevel contribution G (0) 4 and light-like rectangular Wilson loop W A , lim This relation has a transparent physical meaning. In the first quantized picture, G 4 describes the propagation of a scalar particle along a closed contour that goes through the point x i . In the light-like limit this particle has infinite energy. As a consequence, it propagates along a classical trajectory that coincides with C and its interaction with an induced radiation gives rise to an eikonal phase. The latter is given by a Wilson loop in the same representation in which the scalars are defined. This explains the choice of the representation in (2.1). We note that eq. (2.15) is somewhat formal, as the ratio of four-dimensional correlation functions on the left-hand side diverges in the light-like limit x 2 i,i+1 → 0 and requires a regularization. However, for the purposes of the present paper, we can use the relationship (2.15) at the level of the loop integrands. In other words, we can profit from the known integrands for the correlation function (2.14) in order to obtain the four-dimensional integrand for the function F (x). In the following, for simplicity of notation, we write the relations between G 4 and F (x) in an integral form, but keeping in mind that they hold at the level of the four-dimensional integrands.
Combining together (2.6) and (2.15) we obtain Replacing G 4 with its expression (2.14) we can now construct the function F (x) at three loops. Indeed, as was shown in [52], the ratio of correlation functions take the following form in the light-like limit where the integrand I L is obtained from G where I L are expressed in terms of the functions I L according to Here P (L) are homogeneous polynomials in x 2 ij of degree (L − 1)(L + 4)/2 whose explicit expressions can be found in [48,49].
Up to three loops, i.e. for L ≤ 3, the polynomials P (L) do not depend on N . At four loops, P (4) receives a non-planar correction (2.21) As was mentioned above, the planar part can be fixed uniquely whereas the general expression for the non-planar part involves a few arbitrary coefficients where the explicit expressions for Q i and R j can be found in [49] (see Eqs. Combining together relations (2.16) -(2.21), we obtain the following expression for the non-planar correction to the function (2.8) (2.23) The planar correction to F (x) admits a similar representation with P non−planar replaced by a lengthy expression involving a linear combination of P (4) planar and the product of lower loop polynomials P (ℓ 1 ) P (ℓ 2 ) . . . subject to ℓ 1 + ℓ 2 + · · · = 4. To save space, we do not present the corresponding expression.
We can use the relation (2.23) together with (2.10) to compute the non-planar correction to the cusp anomalous dimension at four loops. Notice that applying (2.10), we have to treat x µ 5 as a D−dimensional vector. In this case, the Gram determinants do not vanish, R j = 0, and may contribute to (2.22). As was shown in [48,49], the consistency of (2.17) with the OPE expansion of the correlation function (2.12) for x 2 i,i+1 → 0 implies that the polynomial P (4) nonplanar has to vanish when one of the integration points x i (with i = 5, . . . , 8) approaches the edges of light-like rectangle with vertices at the points x 1 , x 2 , x 3 , x 4 . In D = 4 dimension, this leads to the general expression (2.22). Repeating the same analysis in D = 4 we find that the above mentioned condition is satisfied provided that We use this relation and replace the coefficients c i with their values to get from (2.22) with d 1 and d 3 arbitrary. The terms proportional to d 1 and d 3 satisfy separately the OPE condition mentioned above. As we discussed, this implies that the regions of loop (a) x 1 x 2 x 3 x 4 x 2 x 6 x 7 x 8  integration that potentially produce divergences are suppressed. For this reason, we expect the corresponding terms to have a smooth and vanishing four-dimensional limit. To verify this, we show in the next section that the d 1 term on the right-hand side of (2.25) produces a vanishing contribution to (2.23) and, therefore, can be safely omitted. For the term proportional to d 3 , we assume the same based on the above argument, and thus set d 3 = 0 in our calculation. In summary, we obtained the loop integrand for F (x) at three loops from the corresponding correlation function G 4 at four loops.
Before proceeding to computing the Feynman integrals contributing to F (x), we make two further simplifications. First, we make use of the conformal invariance of F (x) to send the point x 5 to infinity. This leads to the expression x = x 2 13 /x 2 24 for the cross-ratio (2.7), and likewise removes the x 5 −dependence from eq. (2.23). Second, the resulting expressions for the integrals in (2.23) depend on four points x 1 , x 2 , x 3 , x 4 that are null separated x 2 i,i+1 = 0. Introducing the dual variables p i = x i+1 − x i , we notice that these integrals resemble momentum integrals contributing to on-shell four-particle scattering amplitudes. This will be discussed in the next section.

Three-loop integrals from differential equations
In this section we compute the three-loop Feynman integrals contributing to F (x).

Definition of the integral family
The Feynman integral in (2.23) involves the conformal polynomial P (4) non−planar defined in (2.25). It can be expanded into a sum of terms given the product of distances x 2 ij . In spite of the fact that the integral (2.23) is finite, each individual term may produce divergences.
For this reason we set up the calculation in D = 4 − 2ǫ dimensions, and take the limit ǫ → 0 at the end.
After sending x 5 to infinity, as explained at the end of the previous section, we find that the Feynman integrals contributing to (2.23) and to the analogous expression for F (3.1) We will see presently that this integral resembles momentum integrals computed previously for three-loop two-to-two scattering processes [53]. In order to do so, we interpret the x's as dual variables, according to Here, we consider four p i corresponding to a four-particle scattering process. We take p i as ingoing, ordered according to p 1 , p 2 , p 3 , p 4 6 . Furthermore, the sum over momenta is conserved and they are light-like.
We define scalar products of the momenta by Using the definition (3.1), one may represent all planar three-loop two-to-two scattering integrals. For example, the triple ladder integral is given by G 1,1,1,1,1,1,1,1,1,1,0,0,0,0,0 , see Figure 2(a). However, the use of dual coordinates has further advantages. It allows us also to write down products of lower-loop integrals. For example, the product of a one-loop box integrals with a two-loop ladder integral is given by G 1,1,1,1,1,1,1,1,0,1,0,0,1,1,0 . What integrals occur in F ? At leading color and up to three loops, we have precisely planar three-loop integrals, and products of planar lower-loop integrals. Up to the unusual fact that products are considered, these integrals have a clear interpretation in terms of planar scattering processes. The non-planar corrections to F do not have such a simple interpretation, but they can still be expressed in our notation. We found the graphical representation shown in Fig. 2(b) useful.
Some readers might wonder wether the integrals under consideration are related to non-planar three-loop integrals considered in refs. [54,55]. This is not the case. Moreover, there is one important difference, which is somewhat surprising. Although the integrals defined in eq. (3.1) cannot in general be drawn as planar two-to-two scattering diagrams, they do share important analytic properties with planar integrals. It turns out that their Feynman parametrization is positive definite if s, t < 0, such that in this region the integrals are real-valued. Moreover, it can be seen from Fig. 2(b) that the integrals only have cuts in the s-and t-channel, just like planar integrals. For these reasons we expect them to be free of u-channel cuts [56]. This will be important when fixing boundary conditions for differential equations in section 3.3.
In practice, in order to analyze the integrals containing polynomials Q 3 , R 1 and R 2 (see Eqs. (2.23) and (2.25)), it turns out that we need integrals with up to thirteen propagators, i.e. thirteen non-negative indices a i in eq. (3.1). Note that this is three more, and hence considerably more complicated, compared to the integrals that were computed in [53].
We performed the reduction of all integrals needed using the program FIRE5 [57]. As a result, we obtained 257 master integrals, meaning that any integral that we need can be expressed in terms of the latter. 7 In the next section, we explain how we chose a particularly convenient integral basis.

Integrand analysis via improved Baikov variables
In the previous section we found that the integral family of eq. (3.1) is described by 257 master integrals. Before proceeding with their calculation, we will choose a convenient integral basis for the latter. In order to do so, we rely on insights into the structure of Feynman integrals from [58,59]. The main idea is that the Feynman loop integrand contains the key information on the structure of the functions appearing after carrying out the space-time integrations. In particular, the singularity structure of the integrand, and its residues, allow one to predict key features of the outcome of the integration, such as the transcendental weight properties of the answer.
In practice, following [58,59], we consider loop integrands of the type appearing in eq. (3.1), and search for integrands that are free of double poles, and whose leading singularities are normalized to be independent of the external kinematics. In order to do so, we make use of the algorithm [60] (and an independent implementation of it) that is based on a particular parametrization of the integration variables, combined with partial fractioning. Although this algorithm has been applied in many similar situations, the present case is challenging due to the large number of propagator factors. For this reason we made certain improvements that we discuss below.
The first improvement is based on the observation that certain loop integrands can be written as dlog differential forms, in such a way that the integrand is expressed as a single term [61,62]. The key example we need is that of the one-loop on-shell box integrand, Here the new variables α i are closely related to the original propagators, where k * ± corresponds to the two solutions of the maximal cut condition. We note that this is closely related in spirit, but different to the Baikov representation. There, one chooses as new integration variables a set of propagator factors (and possibly irreducible scalar products), in D dimensions. Here, we consider the four-dimensional part of the loop integrand, and our new variables are normalized by (k − k * ± ) 2 . As a result, the change of variables is rational in our case.
Let us now consider the integrand in eq. (3.1) for a i = 1, and with an arbitrary numerator P (x 2 ij ). The denominator is naturally written in the dual x coordinates, and contains three factors of the type (3.5). Therefore we find it natural to perform this change of variable for each loop integration variable. We do so with α i , β i and γ i corresponding to x 6 , x 7 and x 8 , respectively. The only new calculation w.r.t. (3.5) we need to do is for factors of the type x 2 67 . The latter is a rational function of α i , β i , s, t. As a result, we obtain for the integrand where P stands for a polynomial in the x 2 ij variables. We call this representation improved Baikov representation.
Our goal is to find as many numerators P as possible for which the integrand is free of double poles and has constant maximal residues. In principle we could start with a general ansatz for P , subject to certain power counting constraints, and compute residues. As already mentioned, this turns out to be prohibitive due to the size of the integrand. For this reason,we use a divide-and-conquer approach. First of all, most integrals we are interested in contain only a subset of the 15 propagators, so that we can perform a dedicated analysis for each of them. Each propagator structure defines what we call an integral sector. Second, instead of running the algorithm directly for each of these integral sectors, we analyze the integrand on cuts. In this way, each run of the algorithm is faster, and we can easily combine the constraints on P arising from each cut. One of the main advantages of the improved Baikov representation (3.7) is that it makes taking cuts simple. In practice, we find that only a subset of cuts is needed to obtain candidates for dlog integrands. Finally, once a candidate is found, we can test it using our algorithm (this is considerably less complicated than running the algorithm for the whole ansatz).
In this way, we arrive at a large set of dlog integrands. We complement them with information on uniform weight integrals from ref. [53]. We then perform the integral reduction of the corresponding (D-dimensional) integrals, which allows us to relate the candidate integrals to a basis of master integrals. We then select a linearly independent set of candidate integrals, and choose this as our new integral basis. In the next section we derive and solve the differential equations satisfied by these integrals.

Analytic computation of the master integrals
Here we discuss the computation of the three-loop integrals via differential equations. We denote by f the basis of integrals that was found with the methods discussed in the last section. We find that the differential equations (for reviews, see [63,64]) for this basis takes the form where d = ds∂ s + dt∂ t and A i are some matrices. These equations are in the canonical form proposed in [59], as expected for uniform weight integrals. Indeed, it is easy to read off what properties the perturbative solution in ǫ has, as we discuss presently. As one can always choose one overall scale (by dimensional analysis), we can set t = −1 without loss of generality. In this way it is clear that the solution to (3.8) can be written, to any order in ǫ, in terms of iterated integrals of the alphabet where x = s/t. This is slightly more complicated compared to the differential equations for the three-loop integrals derived in [53][54][55], which are described by the alphabet {x, 1 + x}. We note that it is possible to rewrite the alphabet (3.9) in terms of linear letters, with singularities corresponding to fourth roots of unity, 1, i, −1, −i, simply by changing variables to x = −z 2 . In this way, one can rewrite the solution in terms of Goncharov polylogarithms. As we will see below, while the additional letter (1 Feynman integrals, remarkably, we find that it drops out of the results for F . Let us now solve eq. (3.8) for fixed t = −1 and x = s/t. As discussed in section 3.1, f is real-valued for x > 0, and we expect a branch cut along the negative real x-axis, so that one needs to add a small imaginary part to x when analytically continuing to that region. This is relevant when implementing the finiteness of the integrals at u = 0, as discussed above, which corresponds to x = −1.
Our strategy for solving the differential equation is a follows: we write down the general solution in terms of iterated integrals, with base point x = 0 + . Here the superscript '+' indicates that we approach 0 keeping x > 0. This specification is necessary since the Feynman integrals under consideration in general diverge as x → 0, with the rate of divergence being controlled by the matrix A 1 . We parametrize the solution in terms of real-valued integration constants at x = 0 + .
Then, we analytically continue to negative x, and impose the finiteness condition at x = −1. This constraint, together with the real-valuedness of the integration constants, turns out to fix all integration constants, except for a trivial overall normalization. 8 This useful feature that the integration constants can be obtained effortlessly from physical consistency conditions, has already been observed and stressed in a number of calculations, see e.g. [53,56,65]. Finally, we fix the overall normalization by explicitly evaluating one trivial integral in terms of Gamma functions.
In this way, we can obtain the analytic solution for all master integrals to any desired order in ǫ. We are interested in computing F , a finite quantity. On the other hand, typical on-shell three-loop integrals contain nested soft and collinear divergences that lead to up to 1/ǫ 6 poles. As a consequence, we need to expand the solution of the integrals up to finite order, which corresponds to up to transcendental weight 6.
With this we can compute all integrals contributing to F (x) to three loops. In particular, an important check is the finiteness of the answer. Note that this property is highly non-trivial, since individual three-loop integrals have poles up to 1/ǫ 6 . We verify that this is the case for the planar contributions to F .
Proceeding to the non-planar contribution, we are particularly interested in the nonplanar integral Q 3 and the Gram determinants R 1 and R 2 , as discussed in the previous section. We found there that Q 3 by itself does not vanish in the D-dimensional OPE limit, but the linear combination of eq. (2.25) does. Indeed, we find that Q 3 has a 1/ǫ pole in dimensional regularization, while the combination appearing in eq. (2.25) is finite. Moreover, we checked explicitly that the finite result is independent of the free parameter d 1 . We take this as supporting evidence that our procedure is consistent and produces an unambiguous finite answer.

Wilson loop with Lagrangian insertion at three loops
Here we present our novel three-loop results for F (x) defined in eq. (2.8). The tree-level, one-loop, and two-loop contributions to F have been computed in [43][44][45]. The first nonplanar corrections can appear in eq. (2.8) at three loops, i.e. at order α 4 s . They are the main objectives of our work.
As a warm-up we begin by reproducing the lower-loop results. This serves as a welcome validation of our procedure, and of the integrals computed in the previous section. We give the formulas for completeness. We express them in terms of harmonic polylogarithms (HPL) [66], since the latter will be appropriate when giving the three-loop results below. We have correction (unlike the planar one) depends on more than one rational structure. This corresponds to the fact that the non-planar integrand has two leading singularities, namely 1 and 1/(1 + x), while the leading singularities of the planar integrand are x−independent. This is reminiscent of the situation for non-planar three-loop two-to-two scattering amplitudes, where the same leading singularities occur [55].

Cusp anomalous dimension from integration over Lagrangian insertion
The knowledge of F (x) at (L − 1) loops allows one to obtain the cusp anomalous dimension at L loops. The detailed relation was found in [45], and we briefly reviewed it in section 2, see eqs. (2.10) and (2.11). The net result is that we need to perform a certain integral transformation of F (x). The action of the latter on a power of x was given in eq. (2.11), and we repeat it here for convenience, In this section, we show how to apply the functional I to the terms appearing in F (x).
In particular, we are interested in the scenario where the function F (x) is a linear combination of harmonic poly-logarithms of the form of our results F (n) . First, we observe that the r.h.s. of eq. (5.1) vanishes if p is a non-zero, positive integer. Consequently, only the first, O(x 0 ) term of a convergent Taylor series expansion of a function would contribute to this integral. However, the functions we are interested in contain also logarithms log a (x) for positive integers a. We find that The above expression is valid for positive integer p. Next, we want to apply this result to the type of functions we are interested in. Note, that a HPL that admits a convergent power series expansion for x ∈ [0, 1] can be written as Our integration rule of eq. (5.2) may now be easily applied to each term of such a series expansion. The integral over HPLs that are regular at x = 0 consequently vanishes. If they are multiplied by powers of log(x) we obtain a closed form expression by manipulation of the indices appearing in the definition of a HPL (eq. (5.3)). With this we find the following two integration rules. To derive the result of the second integration rule above one can simply make use of the series representation of 1/(1 + x) around x = 0 and re-arrange summation indices. Notice that the right-hand side of the above integration rules is given by linear combinations of HPLs with argument −1, which can be evaluated in terms of multiple-zeta values using for example the package HPL [67]. Since the functions F are of the form considered above we are now equipped to apply eq. (2.10).
6 The four-loop cusp anomalous dimension in N = 4 sYM and QCD Applying the method of the previous section we evaluate I[F (L) ] for L = 0, 1, 2, 3 and use eq. (2.10) in order to determine the cusp anomalous dimension in N = 4 super Yang-Mills. We find This is the full result to four loops for the cusp anomalous dimension in N = 4 super Yang-Mills, for a Wilson loop in the adjoint representation. 9 The expression is valid in the supersymmetric DR scheme. The quartic Casimir term, and equivalently the first subleading color term in eq. (6.1), is independent of the scheme choice at four loops. We note that the non-planar four-loop contribution is negative. We use (6.1) to derive the quartic Casimir part for pure Yang-Mills theory, following [37]. Together with the known n f -dependent terms, this completes the full four-loop cusp anomalous dimension in QCD. We present the full result here, in the MS scheme. First, we recall the expressions up to three loops [2,30]: At four loops, depending on the color factor, various terms were known either analytically or numerically [18-20, 36, 37, 39, 40, 68-72]. 10 We present for the first time the fully analytic result, which is 3 .

(6.3)
For convenience of the reader we also print the above formula to six significant digits.
This result is given for an arbitrary representation R of the Wilson loop. In the QCD context there two choices for the representation, namely R = A (gluons) and R = F (quarks). For the normalization of the SU (N ) generators, we follow the conventions of [73,74], which are We remark that if one retains only the leading transcendental pieces of the QCD result (6.3), i.e. the transcendental weight 2(L − 1) pieces at L loops, and switches to the adjoint representation R → A, then one recovers (6.1), as expected [75,76].

Conclusion and outlook
In this paper we computed analytically the non-planar part of the cusp anomalous dimension in N = 4 super Yang-Mills at four loop. Our result agrees with the previous numerical results of refs. [37,38]. N = 4 super Yang-Mills is expected to be integrable, and there is considerable interest in finding ways to solve the theory. The planar part of the cusp anomalous dimension is successfully described by integrability [29], with important input from perturbative calculations [77]. We expect that our result, which constitutes the leading non-planar contribution to this quantity, will be an important reference value for future integrability studies. As was mentioned in the Introduction, the cusp anomalous dimension controls behavior of the DGLAP splitting functions close to the end-point, or equivalently the large spin asymptotics of twist-two anomalous dimensions, γ(S) = 2Γ cusp log S + O(S 0 ) [6]. Using the obtained result, we can predict nonplanar correction to these anomalous dimensions. At present, nonplanar corrections to γ(S) are known for lowest value of spins S = 0, 2, 4, 6, 8 both in N = 4 sYM and in QCD [50,[78][79][80]. These expressions exhibit an interesting structure and several conjectures have been formulated about the possible form of nonplanar corrections to γ(S) for arbitrary spin S. It would be very interesting to find such a formula. Our result provides the large spin asymptotics of γ(S) and it can be used to constrain an ansatz for this function.
Having the result in N = 4 super Yang-Mills, we derived the value for the purely gluonic quartic Casimir contribution to the cusp anomalous dimension in QCD. The analytic value we computed agrees with the numerical result of ref. [40]. Our result provides the last missing ingredient for the full four-loop result in QCD. By assembling the known terms from the literature, together with our new result for the gluonic quartic Casimir term, we presented for the first time the complete result for the four-loop cusp anomalous dimension in QCD, for an arbitrary color representation of the Wilson lines.
For convenience of the reader we provide ancillary files with our paper that contain the results for F (x) given in eqs. (4.1) -(4.5), as well as the result for the cusp anomalous dimension in N = 4 sYM, eq. (6.1), and QCD, eqs. (6.2) and (6.3).