The three-loop cusp anomalous dimension in QCD and its supersymmetric extensions

We present the details of the analytic calculation of the three-loop angle-dependent cusp anomalous dimension in QCD and its supersymmetric extensions, including the maximally supersymmetric $\mathcal{N}=4$ super Yang-Mills theory. The three-loop result in the latter theory is new and confirms a conjecture made in our previous paper. We study various physical limits of the cusp anomalous dimension and discuss its relation to the quark-antiquark potential including the effects of broken conformal symmetry in QCD. We find that the cusp anomalous dimension viewed as a function of the cusp angle and the new effective coupling given by light-like cusp anomalous dimension reveals a remarkable universality property -- it takes the same form in QCD and its supersymmetric extensions, to three loops at least. We exploit this universality property and make use of the known result for the three-loop quark-antiquark potential to predict the special class of nonplanar corrections to the cusp anomalous dimensions at four loops. Finally, we also discuss in detail the computation of all necessary Wilson line integrals up to three loops using the method of leading singularities and differential equations.

1 Introduction and summary The predictive power of QCD as a theory of strong interaction relies on the possibility to predict the scale dependence of various observables in terms of anomalous dimensions as a function of the strong coupling constant and various kinematical invariants. Well-known examples include the anomalous dimensions of twist-two operators, which govern the scale violation of structure functions of deep inelastic scattering. In this paper we study another important anomalous dimension that appears in many physical quantities involving heavy quarks, the so-called cusp anomalous dimension [1][2][3][4][5][6][7].
The simplest physical process that leads to the appearance of this anomalous dimension is the scattering of a heavy quark off an external potential (see e.g. [8][9][10]). In the infinite mass limit, m Q → ∞, the quark behaves as a classical charged particle -it moves with velocity v µ 1 that changes to v µ 2 after scattering off the external source with the momentum transferred q µ = m Q (v 1 − v 2 ) µ . Due to the instantaneous acceleration, the heavy quark starts emitting gluons with arbitrary momenta. The gluons with small momenta generate infrared divergences (IR), whereas the gluons with large momenta introduce a dependence on the ultraviolet (UV) cut-off. As was shown in [11,12], the dependence of the scattering amplitude on both IR and UV cut-offs is controlled by the cusp anomalous dimension Γ cusp (φ, α s ), which depends on the Minkowskian recoil angle of the heavy quark, (v 1 v 2 ) = cosh φ M , where φ = iφ M .
The heavy quark scattering and its cross channel, the heavy quark production, enter as partonic subprocesses in various important physical applications, e.g. heavy meson form factors in QCD [13] and top quark pair production [14]. In these processes IR and UV cut-offs are replaced by relevant physical scales leading to the appearance of large perturbative corrections enhanced by powers of logarithms of the ratios of these scales. Such logarithmic corrections can be resummed to all orders in the QCD coupling constant. The cusp anomalous dimension is an important ingredient of the resulting resummation formulas which have numerous phenomenological applications (see e.g. [8,[14][15][16]).
The cusp anomalous dimension has a simple interpretation in terms of Wilson loops [12]. The heavy quark couples to gluons through an eikonal current and, as a consequence, the heavy quark scattering amplitude reduces to an eikonal phase. This phase is given by an expectation value of a path-ordered exponential of the gauge field integrated along the classical trajectory of heavy quark. The latter consists of two semi-infinite rays separated by a relative angle φ, i.e. a Wilson loop with a cusp. Due to the presence of the cusp on the integration contour, the vacuum expectation value of the Wilson loop develops specific ultraviolet divergences. The cusp anomalous dimension Γ cusp (φ, α s ) governs its dependence on the ultraviolet cut-off.
The two-loop result for the cusp anomalous dimension has been known for more than 25 years [7] (see also [17]). In this paper, we describe details of the three-loop calculation of this fundamental quantity in QCD. The result was previously reported in [18,19]. The obtained expression for the cusp anomalous dimension has an interesting dependence on the cusp angle. The following three limits are of physical importance: • In the small angle limit φ → 0 the cusped Wilson loop reduces to a straight Wilson line. In this limit the cusp divergences disappear and the cusp anomalous dimension vanishes as −B(α s )φ 2 , with B(α s ) being a positively definite function of the coupling constant.
• In the large (Minkowskian) angle limit, φ = iφ M , with φ M → ∞, the cusp anomalous dimension scales linearly with the angle, K(α s )φ M , with K(α s ) being the lightlike cusp anomalous dimension, which also governs the large-spin asymptotics of the anomalous dimension of twist-two operators [20].
• In the limit of a backtracking Wilson line, φ → π, the three-loop cusp anomalous dimension develops a pole V (α s )/(π − φ). In a gauge theory with exact conformal symmetry, V (α s ) coincides with the analogous function defining quantum corrections to the static quark-antiquark potential. We demonstrate that this relation holds in QCD to up to conformal symmetry breaking corrections proportional to the beta function.
In addition to QCD, we also compute Γ cusp (φ, α s ) in supersymmetric extensions of QCD. There are several reasons for doing this. The cusp anomalous dimension depends on the particle content of the theory. We show that, surprisingly enough, the latter dependence can be eliminated by expressing Γ cusp (φ, α s ) in terms of an effective coupling constant a ∼ K(α s ) closely related to light-like cusp anomalous dimension mentioned above. We find that the cusp anomalous dimension viewed as a function of the cusp angle φ and the new coupling a reveals a remarkable universality property -it takes the same form in QCD and its supersymmetric extensions, to three loops at least. Among various supersymmetric gauge theories, the maximally supersymmetric N = 4 Yang-Mills theory plays a special role. This theory is believed to be integrable in the planar limit (see e.g. [21]), which opens up the possibility of determining the above-mentioned universal function for an arbitrary coupling constant (in the planar limit at least).
The coefficients of the expansion of the cusp anomalous dimension in powers of the coupling constant depend on various color factors of the SU (N ) gauge group. Up to three loops, the latter are given by quadratic Casimirs of the SU (N ) gauge group, whereas starting from four loops new color factors proportional to higher Casimirs can appear [22,23]. A distinguished feature of such color factors is that they generate non-planar corrections to Γ cusp (φ, α s ). We exploit the above mentioned universality property of the cusp anomalous dimension and make use of the known result for the three-loop quark-antiquark potential to predict (conjecturally) this special class of nonplanar corrections to Γ cusp (φ, α s ) at four loops.
In order to compute the cusp anomalous dimension, we need to separate the IR and UV divergences of the cusped Wilson loop mentioned above. We use an infrared suppression factor to remove the IR divergences coming from the integration region at large distances, and employ dimensional regularization (dimensional reduction in the supersymmetric case) to regulate the UV divergences. The cusp anomalous dimension is obtained from the latter in the usual way via a renormalization group equation.
We carry out the calculation in momentum space, where the Wilson lines are replaced by eikonal propagators. As a technical trick, we use eikonal identities to relate all nonplanar integrals appearing in our calculation to (sums and products of) planar integrals. We classify all planar three-loop vertex diagrams of this type, and relate them to master integrals using integral reduction programs. All Feynman diagrams are generated in an automatic way, in an arbitrary covariant gauge, and expressed in terms of the master integrals.
We compute the master integrals applying the differential equations method [24][25][26][27][28] and using the new ideas of [29], recently reviewed in [30]. It was proposed in that paper that the differential equations can be cast into a canonical form that makes properties of the answer manifest, and that can be easily solved. The canonical form is achieved by writing the differential in a certain basis that can be found systematically using the criteria described in [29]. In particular, integrals having constant leading singularities [31] play an important role. The leading singularities [32] of a Feynman integral essentially correspond to residues at certain poles of the integrand, and hence are easily computed. 1 We give a pedagogical introduction to this method, presenting the two-loop computation in full detail, and giving three-loop examples.
We present the analytic result for the three-loop cusp anomalous dimension in terms of harmonic polylogarithms [34]. The latter can be readily evaluated numerically [35,36], analytically continued, or expanded [35][36][37][38][39][40][41] around the above-mentioned interesting physical limits. In this way, we reproduce the known result for Γ cusp (φ, α s ) in the lightlike limit [42][43][44][45][46] and provide new insights into on the relation to the quark-antiquark potential [47][48][49][50] in the backtracking Wilson line limit. We carry out a number of checks of our results. At the level of the Feynman integrals, we reproduce correctly previously known results, including (the gauge-dependent) heavy quark wave function renormalization [51,52]. Another important check is the gauge independence of the final result.
The paper is organized as follows. In section 2 we discuss the main properties of the cusp anomalous dimension, using the one-loop case as an example. In section 3 we discuss our three-loop Feynman diagram calculation, while section 4 is devoted to the calculation of the Feynman integrals. It contains a detailed discussion of the two-loop case. In section 5 we summarize our main results, and in section 6 we discuss their properties. Section 7 contains concluding remarks. There are four appendices. Appendix A summarizes our conventions, Appendix B discusses the heavy quark effective theory (HQET) approach to computing the cusped Wilson loop, Appendices C and D contain a calculation of certain infinite classes of large n f terms of the cusp anomalous dimension and quark-antiquark potential.

Cusped Wilson loop
In a generic four-dimensional Yang-Mills theory a cusped Wilson loop is defined as where the gauge field A µ (x) = A a µ (x)T a is integrated along a closed integration contour C. The latter is smooth everywhere, except for a single point where is has a cusp. Here T a are the generators of the SU (N ) gauge group in an arbitrary representation R and the normalization factor N R = tr R 1 is introduced to ensure that W = 1 + O(g 2 ). The cusped Wilson loop depends on the choice of the representation R, which we take to be an arbitrary irreducible representation of SU (N ). Later in the paper we shall discuss the dependence of the cusp anomalous dimension on R.
For our purposes we shall choose the integration contour C in (2.1) to consist of two semi-infinite rays running along two directions v µ 1 and v µ 2 (with v 2 1 = v 2 2 = 1), with the Euclidean cusp angle φ (see figure 1) In Minkowski space-time the analogous angle is defined as cosh φ M = v 1 · v 2 and it is related to (2.2) as φ = iφ M . The reason for such a choice of the integration contour is twofold. Firstly, the corresponding cusped Wilson loop W has a clear physical meaning in the context of heavy quark effective theory (after analytical continuation from Euclidean to Minkowski space). Namely, it describes the amplitude for a heavy quark with velocity v 1 to undergo the transition into the final state with velocity v 2 . We shall make use of this interpretation later in this section. Secondly, the above choice of the contour facilitates significantly the evaluation of perturbative corrections to W . In particular, it allows us to apply a powerful technique for computing higher loop Feynman integrals, as will be discussed in section 4.

Case of study
Using the definition (2.1), we can expand W in powers of the coupling constant where D µν (x) is the free propagator of the gauge field and C R = T a T a is the quadratic Casimir of the SU (N ) in the representation R. As follows from this relation, the lowest order correction to W does not depend on a particle content of Yang-Mills theory. The latter dependence arises at order O(g 4 ). Going beyond the leading order, we have to specify the underlying gauge theory. In what follows, we shall consider two special cases: (i) gauge field coupled to n f species of fermions all in the fundamental representation of the SU (N ); (ii) gauge field coupled to interacting n s scalars and n f fermions all in the adjoint representation of the SU (N ).
The corresponding Lagrangians are specified in Appendix A. The interaction terms are chosen in such a way that, fine tuning the number of fermions and scalars, we can use these two cases to describe QCD and supersymmetric Yang-Mills theories, respectively. In particular, the maximally supersymmetric N = 4 SYM theory corresponds to n s = 6 and n f = 4. For arbitrary number of fermions and scalars the above mentioned theories are neither conformal, nor supersymmetric. The N = 4 SYM theory plays a special role since both symmetries are exact for any value of the coupling constant. In addition, we can define in this theory a supersymmetric extension of the cusped Wilson loop [53,54] where we introduced the parameterisation of the integration contour x µ = x µ (t) witḣ x µ = ∂ t x µ (t). As compared with (2.1), it has an additional coupling to six scalars φ I that depends on a unit vector n I = n I (t) in the internal S 5 space. In analogy with the previous case, we take this vector to be constant along two semi-infinite rays, n I 1 and n I 2 , except the cusp point where it forms an additional internal cusp angle cos θ = I n I 1 n I 2 . The vacuum expectation value of this Wilson loop operator has been studied in many papers. Perturbative results are available up to four loops in the planar case (and in part in the non-planar case) [48,49,55,56]; results at strong coupling are available via the AdS/CFT correspondence [48,55,57]; the small angle asymptotics is known exactly [58]. Finally, the system is governed by integrability [59,60], for further work in this direction see e.g. [61].

Cusp anomalous dimension
To explain the general framework of our analysis, let us revisit the one-loop calculation of the cusped Wilson loop (2.3). Anticipating the appearance of divergences in W , we introduce the dimensional regularization with D = 4 − 2ǫ. Then, we use gauge invariance of W to choose the Feynman gauge, D µν (x) = g µν D(x), with To the lowest order in the coupling, we find from (2.3) that W is given by a gluon propagator integrated over the position of its end points on the integration contour. Parameterizing points on two semi-infinite rays as −v µ 1 s and v µ 2 t with 0 ≤ s, t < ∞, we arrive at the following integral , (2.6) where in the second relation we switched to the momentum representation. Then, where the second term inside the square bracket comes from the contribution of the gluon propagator attached to the same semi-infinite ray. It is easy to see from the second relation in (2.6) that I(φ) develops poles in ǫ both in the infrared, k µ → 0, and in the ultraviolet, k µ → ∞. In the configuration space, for ρ = s + t, the same poles arise from integration over ρ → ∞ and ρ → 0, respectively. Moreover, since the integral in (2.6) does not involve any scale, it vanishes in the dimensional regularisation, I(φ) = 0, thus indicating that infrared (IR) and ultraviolet (UV) poles of W are in a oneto-one correspondence with each other [11].
In order to compute the cusp anomalous dimension, we have to disentangle UV and IR divergences of W . This can be done by introducing inside the first integral in (2.6) the additional factor exp(−iδ(s + t)) with Im δ < 0. It suppresses the contribution of large (s + t) and introduces the dependence on the IR cut-off δ. In this way, we obtain from (2.6) where we introduced the subscript δ to indicate the dependence on this scale. 2 Changing the integration variables in the first relation to y = s/(s + t) and ρ = s + t, we obtain where the cusp angle φ was defined in (2.2). As expected, I δ (φ) develops a UV pole 1/ǫ. Together with (2.7) this leads to the well-known result for one-loop cusp UV divergence [1] log W = − 1 2ǫ The coefficient in front of 1/ǫ is gauge invariant, it does not depend on the IR regulator δ and defines the one-loop correction to the cusp anomalous dimension. The properties of cusp singularities of Wilson loops are well understood to all loops [1][2][3][4][5][6]. The cusped Wilson loop can be made finite by subtracting UV poles and expressing the resulting quantity log W − log Z in terms of renormalized coupling constant. In the MS scheme, the renormalization Z−factor has the following form The QCD beta-function is well known [62] (we need it to two loops), while for the case of theory (ii) renormalization was discussed in [63][64][65]. The expansion coefficients Γ (i) carry the dependence on the cusp angle φ and define the cusp anomalous dimension Matching (2.10) into (2.11) we obtain the one-loop cusp anomalous dimension 14) where the notation was introduced for x = e iφ and ξ = (1 + x 2 )/(1 − x 2 ).

Regularization
Comparing (2.8) with (2.6) we observe that the net effect of the IR cut-off is to shift the position of poles of the eikonal propagators. This transformation has a simple interpretation in the context of heavy quark effective theory (see Appendix B for more details). As was mentioned earlier, the cusped Wilson loop (2.1) can be interpreted as an amplitude for a heavy quark with the velocity v 1 to undergo the transition into the state with the velocity v 2 (see figure 2). Indeed, the heavy quark propagates along the straight line in the direction of its velocity and the effects of its interaction with gauge fields is described by a Wilson line evaluated along the classical trajectory of a heavy quark. The heavy quark transition amplitude suffers from both IR and UV divergences. The former can be regularized in the momentum representation by slightly shifting the heavy quark from its mass shell with the IR cut-off δ having the meaning of the residual energy of heavy quark. As was already mentioned in the previous subsection, the corresponding Feynman integrals are homogenous functions of loop momenta k. This fact allows us to assign to δ an arbitrary real value. It proves convenient to choose δ = 1/2. Applying the regularization (2.15) with δ = 1/2, we can make use of a very efficient diagram technique for computing W in the momentum representation beyond the leading order (see figure 3 for the corresponding Feynman rules). 3 v . . . To regularize UV divergences we employ dimensional regularization. The cusp divergences come both from the one-particle irreducible vertex corrections V (φ) and from self-energy corrections Σ to the heavy quark propagators (see figure 2). In virtue of Ward identities, the latter contribution is related to the vertex correction at zero recoil angle  V (0) leading to [7] log This relation allows us to compute log Z from the subset of Feynman diagrams corresponding to vertex corrections V (φ), i.e. with non-trivial angular dependence.

Nonabelian exponentiation
The calculation of the cusp anomalous dimension can be significantly simplified by making use of the nonabelian exponentiation property of the Wilson loop [22,23,66]. It allows us to express a logarithm of the Wilson loop, log W , in terms of a special class of 'maximally nonabelian' diagrams, the so-called webs.
In the special case of gauge theories in which all fields are defined in the adjoint representation of SU (N ), this leads to the following general expression where C A = N is the quadratic Casimir operator of SU (N ) in the adjoint representation, f abc f and = C A δ cd , and V n (φ) stands for the sum of certain Feynman integrals defining n−loop corrections to the (one-particle irreducible) vertex function (see figure 4). Notice that the expression on the right-hand side of (2.17) only depends on the quadratic Casimirs. In addition, it is proportional to C R that depends on the representation in which the Wilson loop (2.1) is defined, the so-called Casimir scaling. It is expected that both properties are violated at four loops since the color factors start to depend on higher Casimirs of SU (N ). The power of the nonabelian exponentiation (2.17) is that it allows us to discard the diagrams whose color factor does not contain terms of the maximally nonabelian form. Moreover, we can use (2.17) to express their contribution in terms of Feynman integrals V n that appear on the right-hand side of (2.17). To illustrate this point consider the Feynman diagrams shown in figure 4. The one-loop diagram shown in figure 4(a) has the color factor C R and the corresponding Feynman integral defines V 1 (φ). The two-loop diagrams shown in figures 4(b) and (c) have the color factors C 2 R and C R (C R − C A /2), respectively.
Since the second color factor contains the maximally nonabelian term C R C A , only diagram shown in figure 4(c) contributes to (2.17) at two loops. At the same time, the nonabelian exponentiation implies that the two-loop contribution to W proportional to C 2 R should be related to one-half of the square of one-loop contribution. This leads to the following relation between the Feynman integrals corresponding to diagrams shown in figure 4: 4 Indeed, in configuration space the diagrams shown in figure 4(b) and (c) only differ in the ordering of gluons attached to two semi-infinite rays and the relation (2.18) follows from the identity θ(t 1 − t 2 ) + θ(t 2 − t 1 ) = 1.
Notice that the diagram shown in figure 4(c) is nonplanar. We can then use (2.18) to turn the logic around and express the contribution of this diagram to (2.17) in terms of planar integrals only. The same is true at higher loops. Namely, up to three loops, the vertex function V n (φ) on the right-hand side of (2.17) can be expressed in terms of planar Feynman integrals only. To see this we observe that the sum in (2.17) only depends on C A = N and does not contain nonplanar corrections. Therefore, computing log W in the planar limit we can unambiguously determine V n (φ) up to three loops. Starting from four loops, log W depends on higher SU (N ) Casimirs that generate subleading (nonplanar) corrections suppressed by powers of 1/N 2 (see section 6.1 below). They are accompanied by Feynman integrals that are not necessarily planar.
The fact that only planar Feynman integrals are needed up to three loops is a technical simplification that will be helpful (but not essential) in the calculation described in section 4. The main advantage of planar integrals is that we can define canonical region (or dual) coordinates that make it easy to deal with the loop integrand, without having the ambiguity of redefinitions of the loop momenta. This also makes the classification of all required integrals rather straightforward.
An immediate consequence of nonabelian exponentiation (2.17) is that the cusp anomalous dimension (2.13) has a similar dependence on the SU (N ) Casimirs, with γ, γ A and γ AA depending on the cusp angle φ. We recall that this relation only holds in gauge theories with all fields defined in the adjoint representation. If some of the fields are defined in the fundamental representation, as it happens in QCD, the color factors of maximally nonabelian diagrams have more complicated form and depend on quadratic  anomalous dimension in QCD with n f fermions in the fundamental representation has the following form where T F defines the normalization of the SU (N ) generators in the fundamental representation, tr F (T a T b ) = T F δ ab , and the coefficient functions are different, in general, from those in (2.19). As compared with (2.19), the cusp anomalous dimension in QCD contains the additional terms proportional to powers of T F n f . They come from diagrams involving fermion loops (see figure 5).

Dependence on the cusp angle
In order to discuss the dependence of Γ cusp (φ, α s ) on the cusp angle, it proves convenient to introduce auxiliary (complex) variables In Euclidean space, for 0 ≤ φ ≤ π, we have |x| = 1. In Minkowski space, for φ = iφ M with φ M real, the variable x = e −φ M can take arbitrary nonnegative values. Moreover, due to the symmetry of the definition (2.21) under x → 1/x we can assume 0 < x < 1 without loss of generality.
We can use the one-loop result (2.14) to illustrate interesting asymptotic behaviour of the cusp anomalous dimension in three different limits. For φ → 0, or x → 1, the integration contour in figure 1 transforms into a straight line leading to the vanishing of the cusp anomaly with B = C R α s /(3π) + O(α 2 s ) the so-called bremsstrahlung function. For φ → π, or x → −1, the integration contour degenerates into two antiparallel lines and the cusp anomalous dimension develops a pole with V (α s ) = C R α s + O(α 2 s ) being closely related to the heavy quark-antiquark potential (we shall clarify this relation in section 6.5). In Minkowski space, for large cusp angle, φ M → ∞, or x → 0, the cusp anomaly scales logarithmically with K(α s ) = C R α s /π + O(α 2 s ) being the light-like cusp anomalous dimension. Finally, the Wilson line integrals are naively invariant under the crossing transformation v 2 → −v 2 , or equivalently x → −x (see e.g. one-loop integral (2.6)). This invariance is broken by the Feynman '+i0' prescription and, therefore, we expect it to be valid only up to terms picked up from crossing the branch cut on the negative real axis, where 0 < x < 1 and Disc Γ(−x) := Γ(−x + i0) − Γ(−x − i0) denotes the contribution originating from crossing the branch cut. For example, at one loop we have from (2.14),

Setup of the three-loop calculation
As explained in section 2.3, the cusp anomalous dimension can be calculated within the framework of heavy quark effective theory (HQET). More precisely, we have to compute three-loop corrections to the vertex function V (φ) (see figure 2) and, then, apply (2.16) and (2.11). The corresponding HQET diagrams contributing to V (φ) contain two external heavy quarks with velocities v 1 and v 2 . Note that, by definition, the heavy quarks do not propagate within loops and therefore, we only have to consider massless particles (gluons, fermions and scalars) inside the diagrams. The interaction between massless particles is described by the Lagrangians specified in Appendix A. We use QGRAF [67] to generate all (one-heavy-quark-irreducible) vertex diagrams in HQET. In total there are 315 three-loop diagrams involving gluons and fermions plus 100 additional diagrams involving scalars. As explained in Section 2.4, due to nonabelian exponentiation we only need to calculate the planar diagrams. We find that in the planar limit there are only 120 diagrams, plus 32 diagrams with scalars.
Computing the contribution of three-loop planar diagrams to the vertex function, we performed the numerator algebra using Form [68], TForm [69] or Reduce [70]. In this way, we obtained scalar HQET integrals which can be mapped onto the 8 generic topologies 5 shown in Fig. 6 either manually or using q2e and exp [71,72]. Applying integration-by-parts identities [73], the three-loop integrals are then reduced to a set of 71 master integrals with the help of Crusher [74], Fire [75][76][77] or LiteRed [78,79]. The evaluation of the master integrals, which plays a central role in the calculation, is described in detail in Section 4. Matching the divergent part of the vertex function to the expected form (2.16) and (2.11), we obtain the three-loop cusp anomalous dimension given in Section 5. Its properties, which also serve as indispensable checks of our calculation, are discussed in Section 6.
In addition, in the case of QCD, we compute corrections to the cusp anomalous dimension enhanced by the number of fermions, The details of the calculation can be found in Appendices C and D.

Three-loop calculation of HQET integrals
In this section, we describe our choice of the basis of Feynman integrals that contribute to the cusp anomalous dimension at three loops and present their calculation. An unusual feature of these integrals as compared with the conventional Feynman integrals is that they involve eikonal or heavy quark propagators (see figure 3). In what follows we refer to them as HQET integrals.
In section 4.1, we start by introducing the generalized polylogarithm functions required in our calculation. In section 4.2 we discuss their weight properties and relation to Feynman integrals. To explain the procedure, we first explain our method for computing the master integrals using differential equations. The two-loop case is reviewed as a pedagogical example in section 4.3, Next, in sections 4.4 and 4.5, we explain in detail our choice of integral basis. We give there two complementary points of view, the first being based on analyzing the Wilson line integrals in position space, and the second analyzing generalized cut properties of the same integrals in the momentum-space. Finally, in section 4.6, we perform the three-loop calculation of the master integrals.

Iterated integrals
We will see below that the HQET integrals required in our calculation can be expressed in terms of certain iterated integrals studied in the mathematical literature [80,81]. More precisely, a particular subclass of such integrals, known in the physics literature as harmonic polylogarithms (HPL) [34,82], is sufficient to express all results.
The harmonic polylogarithms H a 1 ,a 2 ,...,an (x) depend on the set of indices a 1 , . . . , a n taking values {−1, 0, +1}. They are defined iteratively with respect to their weight n. The iteration starts with the weight-one functions For all indices being different from zero, a i = 0 for i = 1, . . . , n, higher weight functions are defined as with the integration kernels In the case of all indices being zero, we have The weight of H a 1 ,a 2 ,...,an (x) refers to the number of integrations with logarithmic kernels dx/x, dx/(x + 1), dx/(x − 1) and equals the length of the index vector a = (a 1 , a 2 , . . . , a n ).
In the physics literature, it is sometimes colloquially referred to as "transcendentality". Iterated integrals satisfy a shuffle algebra, which expresses the product of a weight n and a weight m function as a sum over weight k = n + m functions, where the list c of length n + m arises from "shuffling" the lists a = (a 1 , . . . , a n ) and b = (b 1 , . . . , b m ), like a deck of cards. Special values of harmonic polylogarithms at x = 1 and x = −1 are related to nested sums, called Euler sums. The latter satisfy additional relations, see e.g. [34,83], that allow us to reduce them to a minimal number of constants. It turns out that in our calculation, only zeta values ζ n = k≥1 1/k n are needed.

Pure functions of uniform weight
In our calculation, functions of uniform weight play an important role. The latter are defined as linear combinations of iterated integrals of the same weight with coefficients rational in x. For example, is a function of uniform weight 3. We can go further and define so-called pure functions, which are linear combinations of uniform weight functions with rational coefficients, e.g.
Pure functions have nice properties that make them easy to compute. Most importantly, their differential has a simple form. If f (k) (x) is a pure function of weight k, then where c i ∈ Q, the α i (x) are at most algebraic functions, and g For k = 1, the expression on the right-hand side of (4.8) contains only one term, since there is only one (independent) weight zero function f (0) (x) = 1. As a consequence, the relation (4.8) allows for a simple recursive way of defining a weight k function, through differential equations. This is precisely the route that we take in section 4.3 below. Let us imagine we have a set of Feynman integrals that can be evaluated in terms of iterated integrals. Taking certain linear combinations of these integrals, we may try to express them in terms of pure functions. But is there a way to tell in advance which Feynman integrals will evaluate to pure functions?
A proposal in this direction was made in [31], based on ideas related to generalized unitarity [32,84], and relying on a large body of evidence from computations in N = 4 super Yang-Mills. To understand this, let us imagine a Feynman integral depending on many kinematic variables. Iterated integrals are multivalued functions in these variables. The idea is that taking generalized unitarity cuts of Feynman integrals should in some way correspond to taking discontinuities of these functions (the precise correspondence between the two objects is an open problem). Then, taking different discontinuities should project onto different terms in the expression of an integral in terms of iterated integrals. The conjecture of [31] is that if all leading singularities (corresponding to a series of cuts that completely localize a Feynman integral) are rational numbers, the answer is a pure function.
This conjecture was verified in a number of non-trivial examples, see e.g. [85,86]. Moreover, it turned out that choosing such integrals as a basis rendered the physical answer much simpler already at the integrand level, before carrying out the integrations. (This is in part due to the close relationship between certain unitarity cuts and infrared divergences, whose appearance is clearer in the new basis choice.) The understanding of the relationship between Feynman integrals and uniform weight functions was put on a firmer footing in [29], by providing a way of proving the conjecture with the help of differential equations. It is known that Feynman integrals viewed as functions of kinematical invariants satisfy a system of first-order partial differential equations (see e.g. [28,30] for reviews). Denoting the set of such integrals by f (x) we have in complete generality whereÃ(x, ǫ) is a nontrivial matrix depending on the dimensional regularization parameter ǫ and the differential on both sides is taken with respect to x. In [29] it was suggested that by changing the basis f → T (x, ǫ)f to uniform weight functions, the differential equation (4.10) should take a simple canonical form, with the matrixÃ(x) being ǫ independent and given by a linear combination of logarithms with rational coefficients. We can write a formal solution to (4.11) as a path-ordered exponential with some contour C connecting the base point at x = 1 with the function argument, and g(ǫ) = f (1, ǫ) being the boundary values.
In fact, in the canonical form (4.11), each term in the ǫ−expansion of f = ǫ k f (k) satisfies an equation of the form (4.8), which allows one to prove that f (k) has uniform weight k. 6 Moreover, assigning weight (−1) to ǫ, each f has uniform weight, in the sense of the ǫ expansion.
How does one find an appropriate basis f ? Given the conjecture of [31], integrals having constant leading singularities in the sense of that paper are a natural choice. The differential equations then allow one to prove the uniform weight properties of those functions. We remark that generalized unitarity cuts are also very natural in the context of differential equations. The reason is that the cut integrals satisfy the same differential equations as the original integrals, but with different boundary conditions. The cuts allow one to focus on a subset of integrals that share a common propagator structure. This can be used e.g. for making consistency checks before the whole system of differential equations is considered.
Another strategy put forward in [29] is to try to deduce the weight properties from an integral representation, e.g. in Feynman parametrization. This works particularly well in cases with few propagators, and for Wilson line integrals. We will present various examples below. 7 In subsections 4.4 and 4.5 we will see various examples of analyzing weight properties of Feynman integrals before integration, either based on parametric representations, or based on generalized unitary cuts.

Two-loop master integrals and differential equations
Let us present as an example the full set of differential equations (4.11) for the two-loop case. The analysis at three loops will be almost identical, except that the system will be much larger. Going through the simple two-loop example therefore allows us to be more explicit.
In order to discuss all planar two-loop integrals (recall that non-planar integrals can be obtained via eikonal identities), we introduce the following notation where a 1 , . . . , a 7 are arbitrary integer indices and (4.14) Examining (4.13) for various values of indices, we identify two families of integrals that match topology of planar Feynman diagrams contributing to the cusped Wilson loop at two loops. They are shown in figure 7(a), (b) and are given by G a 1 ,a 2 ,a 3 ,a 4 ,a 5 ,a 6 ,0 , G 0,a 2 ,a 3 ,a 4 ,a 5 ,a 6 ,a 7 , (4. 15) respectively. All other planar two-loop integrals can be obtained by pinching lines (setting some of the a i to zero), or adding numerators (setting some a i to negative values).
Integral reduction using integration-by-parts (IBP) identities [73] shows that there are nine master integrals that can come from the two families of integrals shown in (4.15). Notice that integral reduction programs automatically choose a particular integral basis f according to certain criteria. Such a basis typically does not have the uniform weight properties discussed above, and hence leads to a complicated form of differential equations (4.10). In order to bring the differential equations to a simple canonical form (4.11) we make the following choice of master integrals, where χ = (1 − x 2 )/x. Here dot denotes a propagator squared in momentum space. A distinguished feature of this basis is that all functions f 1 , . . . , f 9 have a uniform weight. This property is by no means obvious and can be established using the methods discussed in subsections 4.4 and 4.5.
All integrals depend on dimensionless kinematical variable x with v 2 1 = v 2 2 = 1, and are normalized in such a way that f 1 , . . . , f 9 are expected to be pure functions of weight zero. This can be verified by computing their differential with respect to the kinematic variable x. Using the definition (4.25), we can implement this by differentiating the defining Feynman integrals with respect to v 1 , In this way, we find that the set of nine basis integrals f = (f 1 , . . . , f 9 ) satisfies the differential equation (4.11) with b 2 = diag(4, 2, 4, 0, 0, 2, 2, 0, 2) and In order to solve the differential equation (4.27) we also need boundary conditions. The latter can be easily fixed for x = 1, or equivalently v µ 1 = v µ 2 , where no singularities are expected from the Feynman integrals. Since χ = 0 in this limit, the only non-vanishing integrals in (4.16) -(4.24) are f 4 , f 5 , and f 8 . For x = 1 they are reduced to integrals with bubble insertions and can be easily evaluated (see relations (4.43) below). In this way, we find with all other integrals vanishing at x = 1. Making use of it is easy to see that the above expressions give rise to uniform weight ǫ expansions. 8 Returning to the differential equation (4.27), we can write its solution as Matching the coefficients in front of powers of ǫ on the both sides of (4.27), we find that f (k) (x) are given by a Q-linear combination of harmonic polylogarithms of weight k defined in section 4.1. It is then straightforward to expand f (x) to any desired order in ǫ. For instance, we have (4.32) In agreement with our expectations, the coefficients in front of powers of ǫ are pure functions.
We should mention that the basis choice of f (x) is not unique. As we show in the next subsection, we can introduce two other integrals that are also pure functions of weight zero. They are related however, via IBP identities, to the nine basis functions f (x). In order for g 1 and g 2 to be pure functions, they should be given by a Q-linear combination (independent of x and ǫ) of the basis integrals. Indeed, we find that This also means that, replacing e.g. f 7 by g 2 would have lead to an equally nice set of differential equations.
Of course, not all integrals have such nice properties. As an example, consider the following integral that can appear in the Feynman diagram calculation This integral is obviously not a pure function. Choosing it as a basis integral would lead to an unnecessarily complicated dependence of the differential equations on ǫ and x .

Wilson line diagrams in position space and uniform weight integrals
In this and in the following subsection we explain the method that we use to identify integrals that can be evaluated in terms of pure functions. As a warm up example we revisit the calculation of Wilson line integral (see figure 8(a)) that contributes to the one-loop cusp anomalous dimension. It is given in position space by an integral of a scalar propagator connecting two points −sv µ 1 and tv µ 2 , with s and t being the line integration parameters, . (4.37) Using this identity, the integrand can be written in the so-called "d-log"' form [56] (4.38) In this form, it is manifest that the integral, multiplied by (1 − x 2 )/x, has a differential of the form (4.8), with n = 2, and is hence a pure function of weight two. Likewise, Wilson line integrals with more propagators stretched between two (or more) semi-infinite rays are seen to be pure functions of higher weight. (An algorithm for computing all such contributions was given in ref. [56].) The above analysis is rather formal since the Wilson integral is divergent and requires regularization both in UV and IR. As we will see in a moment, regularization does not affect the uniform weight properties of the integral. For example, at one loop the regularized Wilson line integral (2.8) is given, up to overall factor (4.39) Here we changed variables according to s = ρz, t = ρz and introduced notation for The ρ integral in (4.39) gives UV divergence 1/ǫ. More generally, introducing ρ as an overall scale in a given Wilson line integral, we can always separate the ρ integration from the rest of the calculation. Moreover, the ρ integral can always be evaluated in terms of gamma function, typically Γ(Lǫ) at L loops, which does not change the weight properties of the answer, except for an overall offset. Therefore, in the examples below, we will not discuss further the ρ integration. Let us examine the integral I 1 (x, ǫ). The integrand in (4.40) can be Taylor expanded in ǫ. At ǫ = 0, the integral I 1 (x, 0) is obviously evaluates to a logarithm, i.e. a pure weightone function. It is easy to see that expanding (4.40) at higher orders in ǫ will increase the weight of the resulting function accordingly. With the convention that ǫ has weight (−1), we can therefore see that I 1 (x, ǫ) has uniform weight one. We could proceed along the lines of section 4.2 (see also refs. [56,88]) and evaluate the integral, at a given order in ǫ. Instead, in this paper, we evaluate all such integrals using differential equations, with ǫ as a parameter.
We can show in a similar manner that the integrals with L propagators attached to two (or more) semi-infinite rays are expressed in terms of pure functions of weight 2L. Indeed, we can apply the identity (4.38) to each propagator to deduce that the integral is given by (1/χ) L (with χ = (1 − x 2 )/x) times a pure function of weight 2L. As in the one-loop case, introducing regularization does not affect this result. 10 For example, at two loops, the ladder integral that enters into the definition (4.16) of the basis function f 1 , is given by the product of χ −2 and a pure function of weight 4. Multiplying the ladder integral by ǫ 4 χ 2 we therefore obtain a pure function of weight 0. This explains the origin of the normalization factors in the definition of f 1 .
The above analysis can be generalized to integrals where propagators are raised to some power. For example, consider the integral of figure 8(b) where a dot denotes a propagator squared (in momentum space). Parametrizing the end-point of propagators according to s 1 = ρz, t 1 = ρzy and t 2 = ρz (withz = 1 − z and similar forȳ), this leads to (up to an inessential overall factor and terms suppressed by powers of ǫ), Here the UV pole 1/ǫ comes from ρ−integration and the additional factor of z(1 − z) on the left-hand side comes from the Jacobian of the change of variables and the doubled eikonal propagator. This shows that the integral is given by a function of weight three. We can convert it into a pure function of weight zero by multiplying the integral by the normalization factor ǫ 3 χ 2 . The resulting function coincides with g 1 (x) defined in (4.33). Another example is the integral shown in figure 8(c). The Fourier transform of the doubled propagator gives ∼ (−x 2 ) −ǫ /ǫ, so that this factor is irrelevant at the level of the integrand and can be replaced with 1/ǫ. Parametrizing the line integrals as in the previous example, we obtain We conclude that this integral multiplied by ǫ 3 χ yields a pure function of weight zero. It coincides with the basis function f 6 defined in (4.21). These examples might mislead the reader in thinking that the uniform weight property is rather trivial. However, this is not the case. For instance, just moving the dot in the above examples to another propagator destroys this property. In our analysis, it would result in the impossibility of rewriting the integrand in a "d-log" form.
For integrals with fewer propagators, bubble subintegrals can appear. Whenever this happens, the latter can be integrated out, leaving one with an integral that effectively has one loop less, up to some gamma functions coming from the integration. This means that many integrals can be chosen based on the knowledge of pure functions at the lower loop order. The relevant formulas are obtained by elementary integrations in Feynman parameter space. For a bubble on an eikonal line (see (4.17)) and for a scalar bubble (see (4.18)) we have, respectively, The momentum dependence of these integrals that is important for the present analysis can be simply obtained by power counting. We can use the relations (4.43) to express the two-loop integrals entering the definition of basis functions (4.17), (4.18), (4.20) and (4.24) in terms of one-loop integrals. Moreover, the bubble-type integrals entering (4.19) and (4.23) can be entirely evaluated in terms of Γ−functions. In this way, we verify that f 2 , f 3 , f 4 , f 5 , f 8 and f 9 are indeed pure functions of weight zero. 11 Let us now discuss the two-loop master integrals with an internal interaction vertex, cf. (4.22) and (4.34). It is convenient to analyze them in position space as well. For simplicity, we will carry out the analysis in four dimensions. Let us begin with the integral in (4.34) and denote by x 1 , x 2 , x 3 the points the three-point vertex is attached to, with x 2 , x 3 lying on the same Wilson line segment. (For the integral in (4.22) we can set x 2 = 0.) These points can be parametrized by with s 1 > 0 and t 2 > t 1 > 0. Consider carrying out the integration over the internal vertex. The integral involves three scalar propagators attached to this vertex and it gives rise in four dimensions to [90] 1 iπ 2ˆd where is a known pure function of weight two. Its explicit expression is not relevant for the present analysis. The latter focuses on the question whether the integrand can be put in "d-log" form.
Just as in the case of propagator exchanges, there are simplifications due to the fact that the Wilson lines lie in a plane, which leads to simplifications. After some algebra, we find for the integrand (for x < 1) where in the last relation we changed variables as t 1 = ρ(1 − z)y, t 2 = ρ(1 − z), s 1 = ρz. The ρ integration just gives an overall normalization, while the remaining integrand can be put into a "d-log" form. Remembering the weight-two functionΦ (1) , we expect that the integral, normalized by (1 − x 2 )/x, gives a pure weight four function. Then, we multiply it by ǫ 4 to obtain a pure function (4.34) of weight zero. We can use the calculation above in order to also analyze the integral (4.22) where x 2 = 0. This is simply achieved by setting t 1 = 0 and dropping the t 1 integration in (4.47). In this case, after changing variables according to t 2 = ρ(1 − z), s 1 = ρz we obtain Notice that the ρ−integral is divergent at ρ = 0. If we introduced the dimensional regularization from the beginning, the integrand (4.48) would be modified by the factor ρ 2ǫ leading to a 1/ǫ pole upon integration over ρ. Therefore, as in the previous case, we expect the integral in (4.22) to be a uniform function of weight four and, as a consequence, the basis function f 7 to be a pure function of weight zero. It is clear that the method discussed in this subsection does not rely on a particular loop order and it proves to be very useful in selecting uniform weight integrals at the three-loop order.
The attentive reader may have noticed that the above analysis relied mainly on the properties of the denominator factors and not those of the functionΦ (1) defined in (4.46). In fact, ignoring this function corresponds to taking a generalized cut, making contact with the conjecture of [31]. Similarly, and perhaps more easily, we could have taken the maximal cut of this integral in momentum space, with the conclusion that it has a unique normalization factor x/(1 − x 2 ). As we will demonstrate in the next subsection, the approach based on generalized cuts is especially useful for more complicated integrals with many propagator factors that can be cut.

HQET integrals in momentum space and maximal cuts
In this subsection, we perform an analysis of maximal cuts of HQET integrals in momentum space. The objective is to determine whether a given integral has a unique overall normalization factor, consistent with being a pure function. We will start by reviewing some of the integrals of the previous subsection, and then turn to an example occurring in three-loop computation.
Let us start by verifying the normalization factor of the one-and two-loop ladder integrals. We work in four dimensions but keep IR regularization with δ = 1/2. The maximal cut of the one-loop integral (2.8) is given by (here and in the remainder of this subsection we will neglect inessential x−independent normalization factors) There are various ways of evaluating this integral. To solve the massless on-shell condition for the loop momentum k 2 = 0 we make use of spinor variables (see, e.g., [91]) 12 50) or simply k = ρ|λ [λ|. Together with 2k · v i = λ|v i |λ], this leads to This is indeed the correct normalization factor, cf. eq. (4.39). Similarly, a short calculation shows that the maximal cut of the double ladder integral is given by (4.52) 12 Another way could be to use Sudakov decomposition k µ = αv µ 1 + βv µ 2 + k µ ⊥ and carry out integration over α, β and k ⊥ . which is consistent with eq. (4.16).
Let us now consider a more complicated example of three-loop integral containing 9 propagators shown in figure 9(a). Its maximal cut is best understood by first evaluating Figure 9. Three-loop integral and one-loop subintegral, whose maximal cuts are considered in the main text.
the maximal cut of the one-loop subintegral shown in figure 9(b) (with all external legs cut, i.e. k 2 . The latter is given bŷ where four delta-functions localize the k 1 −integral. Applying (4.53), we effectively reduce the integral of figure 9(a) to a two-loop integral. It contains however two additional propagators coming from the right-hand side of (4.53) and does not produce a function of uniform weight. We can improve the situation by inserting into the integral of figure 9(a) a numerator factor depending on loop momenta. The latter can be chosen, e.g., to cancel part of the factors coming from (4.53). In this way we can obtain two-loop integrals that are expected to be of uniform weight based on the analysis of the previous subsection. Explicitly, inserting the numerator factors (−k 2 2 ) and (−2k 3 · v 1 + 1), we evaluate the maximal cut as (4.55) With this choice of numerator factors, the three-loop integrals have a unique normalization factor and, therefore, they are good candidates for pure functions. This result is not too surprising, given that very similar results were obtained in [89] for massless two-to-two amplitudes.
The techniques described in this and the previous subsection allow us to easily and quickly assemble a list of candidate integrals that give rise to pure functions. In the case of the generalized unitarity cut or leading singularity analysis, this is expected based on the conjecture of [31]. The differential equation method allows us to prove the uniform weight property in the cases where it was only conjectured.

Three-loop master integrals and differential equations
In this subsection we extend the calculation of HQET master integrals to the three-loop level. As discussed in section 2.4, thanks to eikonal identities we need to calculate only planar HQET integrals. To this end, we define all planar integral families at three loops, describe the choice of master integrals and their computation via differential equations. Due to the size of the matrices involved, unlike the two-loop case, we select not to present the latter in this paper, but provide them and other results in the form of ancillary text files.

Definition of master integrals
All planar three-loop HQET integrals can be viewed as some special cases of the integral families shown in figure 10. Thanks to planarity it is possible to describe all of them using a global parametrization of the loop momenta k 1 , k 2 , k 3 . In order to do so, we define the following factors, (4.56) We then introduce the following notation for the HQET integrals, The integral families shown in figure 10 correspond to the following expressions in the notation of (4.57) (more generally, the a−indices can of course be different from 1), Numerator factors can be accommodated by negative values of the indices a i . It is worth pointing out that the labeling in G is not unique, in the sense that the same integrals can be represented by different index vectors. This is due to invariance under relabeling of loop momenta, due to symmetry of some graphs and due to a v 1 ↔ v 2 symmetry of the integrated results.
We remark that with the above setup we can also discuss factorized integrals. In particular, one-loop integrals multiplying generic two-loop integrals can be treated as a subset of the three-loop integrals. This is a useful check, and also allows for a convenient calculation of, say, higher orders of their ǫ expansion, within the same setup.
Solving the IBP relations for integrals shown in figure 10, we find that there are 71 master integrals in total. We choose the master integrals according to the uniform weight criteria explained in detail in subsections 4.4 and 4.5, following [29]. We denote the basis integrals by f = (f 1 , . . . , f 71 ), hoping that using the same letter f that we previously used to denote two-loop basis integrals with will not lead to confusion. As in the two-loop case, all basis three-loop integrals f are pure functions of x of weight zero. All except a handful of integrals could be chosen to be given by a single master integral (4.57), with certain powers of propagators, and normalized appropriately. Only in a few cases it turned out to be necessary to consider linear combinations of integrals (4.57). 13

Integral subsector at three loops
As an example of the basis integrals at three loops, let us to return to three-loop integrals discussed in subsection 4.5 (cf. eqs. (4.54) and (4.55)). In the notation of eq. (4.57) they read f 70 = ǫ 6 χ 2 G 1,1,0,0,1,1,1,1,1,1,−1,1 , We can use them to illustrate the relationship between generalized cuts and projections onto sectors of the differential equations. Let us consider the maximal cut of these integrals, i.e. replace all scalar and eikonal propagators with their cut version, and denote the resulting integrals byf 70 andf 71 . The latter satisfy a closed system of differential equations. This system of two equations is a subset of the full system of 71 equations. This follows from the fact that cut integrals satisfy the same IBP relations as standard ones [92]. Another way of saying this is that cut integrals satisfy the same differential equations as the standard integrals, but with different boundary conditions (in particular, the remaining basis integrals vanish upon taking the above mentioned cut,f i = 0 for i < 70). This means that the subsystem of basis integrals (4.58) and (4.59) is relevant for the full calculation. In particular, it can serve as a check of whether the choice (4.58) and (4.59) is consistent with the canonical form (4.11) of the differential equations.
Indeed, we find that the integrals (4.58) and (4.59) satisfy the system of differential equations, which is consistent with (4.11). Of course, removing the cut, it could be that terms violating the form (4.11) are present in off-diagonal terms. If this is the case, one can attempt to remove them using the methods discussed in ref. [93] and more recently in [30,94]. It turns out that this is not needed for the two integrals under discussion. In our calculation, we resorted to such "brute-force" methods only in the case of a handful of integrals.

Full system of differential equations
With the basis f given in ancillary files included in the arxiv submission of this article, the differential equations take the form with constant 71×71 matrices a 3 , b 3 and c 3 given in the ancillary file HQET 3loop mAtilde.m. We see that eq. (4.61) has four regular singular points, 0, 1, −1, ∞. Due to the x ↔ 1/x symmetry of the definition (4.25), only the first three are independent. They correspond, in turn, to the light-like limit (infinite Minkowskian angle), zero angle limit and backtracking limit.
As before, we can solve eq. (4.61) in a Laurant expansion (4.31). Then we can express f (x) order by order in ǫ in terms of harmonic polylogarithms. We use the value of basis integrals at x = 1 as boundary condition [52,95] (see [96] for a summary).
Most of the basis integrals can be evaluated trivially for x = 1 in terms of Gamma functions. Boundary conditions for Feynman integrals can often be obtained without additional work, by imposing physical properties. In reference [97] this was used e.g. in a bootstrap approach to compute single-scale integrals from differential equations. In the present case, we can use finiteness of the limit x → 1 as our main condition. It turns out that only one non-trivial integral is needed at x = 1. It is known up to weight five [98] (but this is the order we are interested in), G 1,1,1,0,0,0,1,1,1,0,1,1 (x = 1) = 12ζ 2 ζ 3 − 5ζ 5 + O(ǫ) , (4.62) which is exactly the order we need for our calculation. It is likely that also this integral could be obtained by inspecting the differential equations more closely, or applying bootstrap ideas as in [97].
There are a number of analytic checks. Out of 71 master integrals, 7 are straightline ones (studied in [52,95]), 8 can be chosen as products of lower-loop integrals, and 10 correspond to the one-loop triangle integral with ǫ-dependent powers of denominators (studied in [102]). One non-trivial integral at x = 1 was obtained in our approach from the finiteness of the x → 1 limit for all integrals entering the differential equations. Previously, it was computed in terms of a hypergeometric function in ref. [103] (another hypergeometric representation was derived in [104]). We can expand it in ǫ using the Mathematica package HypExp [105]. The result, written up to weight five, is We found perfect agreement with our result.
Using the obtained results, we reproduce results of the three-loop computation performed in [49]. For example, the following three-loop integral was computed there (taking into account the conversion between the different regulators), (4.69) In terms of the three-loop basis integrals defined above it reads We can use the above results to compute the three-loop cusp anomalous dimension for the supersymmetric Wilson loop in N = 4 super Yang-Mills theory (4.72) 16 Strictly speaking, the calculation in ref. [49] was performed for θ = 0 in which case ξ0 = (1 − x)/(1 + x).
In perfect agreement with findings of ref. [49], we find with y = 1 − x 2 and ξ 0 given by (4.67). This is a highly nontrivial test of the calculation of the master integrals. Within the differential equations method, the calculation of a given integral requires the knowledge of all integrals appearing in sub-topologies (obtained by removing propagator factors). Since the integrals needed here have a maximal number of propagator factors, this calculation is also a consistency check of many other integrals appearing for example in QCD.

Results
The basis integrals defined in the previous section allow us to compute the cusped Wilson loop (2.16). For example, we can express the three-loop correction to W as where C i are coefficient functions rational in x and depending on ǫ. Their explicit form is not particularly enlightening. We can write similar formulas for W (1) and W (2) . 17 Then, we extract divergent part of log W and match it into expected form (2.11) of log Z.
In this way, we verify gauge independence of Z−factor, reproduce well-known result for β−function and extract the three-loop cusp anomalous dimension.

Coefficient functions
To express our results for three-loop cusp anomalous dimension we introduce the following functions H 0,1 (y) + π 2 6 H 1,1 (y) + 2H 1,1,0,1 (y)  The three-loop cusp anomalous dimension involves particular linear combinations of these functionsÃ As follows from the definition, these functions vanish for zero cusp angle, or equivalently x = 1.

Three-loop cusp anomalous dimension
In QCD with n f fermion flavours, we obtained the following result for the three-loop cusp anomalous dimension (2.13) in the MS scheme In the gauge theory with n f fermions and n s scalars in the adjoint representation of the SU (N ), the three-loop cusp anomalous dimension is given by In section 6.2, we will also give the result for Γ N =4 in the dimensional reduction scheme. A close examination of (5.6) and (5.4) shows that the coefficients γ f f , . . . , γ s describing n f and n s -dependent contribution at three loops, involve the same functionsÃ 1 ,Ã 2 andÃ 3 that already appeared at two loops. This suggests that these coefficients are not independent. Indeed, we show in the next section that the cusp anomalous dimension has an interesting hidden structure that allows us to predict all n f and n s -dependent terms at three loops at least.
Notice that all functions in (5.2) except B 5 (x) depend on y = 1 − x 2 and, therefore, they are formally invariant under x → −x. However, due to the presence of the cut that runs along negative x, these functions acquire an additional contribution under x → −x proportional to their discontinuity across the cut (see (2.25)). For the function B 5 (x) the situation is slightly different. The linear combination of harmonic polylogarithms inside the brackets in B 5 (x) formally changes the sign under x → −x. It is compensated however by the odd prefactor x/(1 − x 2 ), so that B 5 (x) has the same parity properties as the other coefficient functions. In this way, we verify that our results for three-loop cusp anomalous dimension (5.5) and (5.6) satisfy the relation (2.25).
6 Properties of the cusp anomalous dimension We observe from (5.5) that the dependence of the three-loop cusp anomalous dimension on the representation R enters through an overall factor given by the quadratic Casimir of this representation C R = T a T a , the so-called Casimir scaling with γ(φ, α s ) being independent of R. As was already mentioned in section 2.4, we expect this scaling to be broken at four loops due to appearance of higher Casimirs. To understand this property, let us examine possible color factors that can appear in the perturbative expansion of Wilson loop (2.1) up to four loops. To simplify the analysis we first examine supersymmetric Yang-Mills theory. The color factor in this case consists of terms having the form tr R [T a 1 · · · T an ]C a 1 ...an with each T a i corresponding to a gluon attached to the integration contour. The tensor C a 1 ...an is a product of δ a i a j and if a i a j a k factors. 18 So, there always exists a if a i a j a k factor directly contracted to the trace of Figure 11. Graphical representation of the color factors C 1 , . . . , C 4 . Double line represents tr R [T a1 · · · T an ], solid line denotes δ aiaj .
T a i . Substituting if abc T c = [T a , T b ], we transform such terms into the sum of two terms having less if a i a j a k factors (but more T a i factors inside the trace). Applying this procedure recursively, we finally reduce any color factor to a linear combination of terms of the same form where all C−tensors are products of δ a i a j only. In this way, we obtain the basic color factors shown in figure 11.
The remaining color factors can be reduced to products and sums of the basic ones. Going through the calculation we find where N R = tr R 1 is the dimension of the representation and f abc f abd = C A δ cd with C A = N being the quadratic Casimir of the adjoint representation of SU (N ). An important difference of C 4 compared to (6.2) is that it cannot be expressed in terms of quadratic Casimirs only. More precisely, its takes the form where the ellipsis denotes terms involving quadratic Casimirs C R and C A . Here d abcd R and d abcd A are fully symmetric tensors with the generators T a defined in two different representations.
The color factors C n appear in the expression for the cusp anomalous dimension starting from n loops. The very fact that (6.3) is not proportional to the quadratic Casimir for the generic SU (N ) representation R implies that the Casimir scaling (6.1) should be violated at four loop unless some miraculous cancellation happens leading to the vanishing (angle dependent) coefficient function accompanying C 4 .
Notice that the color factors (6.2) contain higher power of C R . As we explained in section 2.4, in virtue of nonabelian exponentiation, the cusp anomalous dimension should involve maximally nonabelian factors only. Up to three loops they take the form C R , A . This means that the cusp anomalous dimension depends on particular combinations of the color factors, C 1 , C 2 − C 2 1 and C 3 + 2C 3 1 − 3C 1 C 2 . At four loops, the maximally nonabelian color factors are of two kinds, C R C 3 A and d abcd R d abcd A /N R . The latter color factor leads to a violation of the Casimir scaling (6.1) at four loops and induces a nonplanar correction to the cusp anomalous dimension.
Let us now consider the color factors in QCD. An important difference with the previous case is that the fermions are defined in the fundamental representation. This leads to the appearance of additional color factors proportional to the number of fermion flavours n f . Each fermion loop produces a factor of n f and the maximal power of n f scales with the loop order. In particular, the color factors linear in n f have the form n f tr R [T a 1 · · · T an ] tr F [T b 1 . . . T bm ]C a 1 ...an;b 1 ...bm , with the C tensor being given by a product of Kronecker symbols. As in the previous case, up to three loops n f -dependent color factors can be expressed in terms of quadratic Casimirs C R , C A and C F , where T a T a = C F is the quadratic Casimir of the fundamental representation of SU (N ). Most importantly, the additional n f dependence does not affect the Casimir scaling (6.1) at three loops but it modifies the form of the function γ(φ). At four loops, we encounter the color factor n f d abcd R d abcd F /N R analogous to (6.3), with the completely symmetric d F tensor given by (6.4) in the fundamental representation. As before, it is not proportional to C R and, therefore, leads to violation of the Casimir scaling.
To summarize, the general expression for the four-loop contribution to the cusp anomalous dimension violating the Casimir scaling is where f A (φ) and f F (φ) are some functions of the cusp angle depending on the choice of the gauge theory. Here the second term inside the brackets is present only if fermions are defined in the fundamental representation, e.g. f F (φ) = 0 in N = 4 SYM. In the special case of R being the fundamental representation of the SU (N ), we have with C F = (N 2 − 1)/(2N ). Since these color factors involve various powers of N , the expression on the right-hand side of (6.5) generates nonplanar corrections to the cusp anomalous dimension.

Renormalization scheme change
We recall that the three-loop calculation of the cusp anomalous dimension has been performed using dimensional regularization (DREG). However supersymmetry is broken in DREG since for D = 4 − 2ǫ the number of bosonic and fermionic degrees of freedom do not match for ǫ = 0. To restore the supersymmetry, we can employ dimension reduction (DRED) [106]. In this scheme the gauge fields have four components in D dimensions and the difference with DREG comes from the contribution of additional (4 − D) components of the gauge field, the so-called ǫ-scalars. Since the number of scalars n s is a free parameter in our calculation, we can easily accommodate the contribution of ǫ-scalars by replacing n s → n s + 2ǫ. Additional complications arise due to necessity to introduce evanescent coupling constants describing the self-interaction of ǫ-scalars and their coupling with fermions. In a generic gauge theory, the renormalization group evolution of the evanescent couplings differs from that of the gauge coupling and, therefore, they have to be treated differently. However, in a supersymmetric theory the beta-functions of these two sets of coupling necessarily coincide allowing us to identify them at any scale. In this case, to compute the cusp anomalous dimension in the DR scheme it suffices to replace n s → n s + 2ǫ in expression (2.11) for the Z factor in the MS scheme, identify the residue at the pole 1/ǫ and take into account the relation between the coupling constants in the two schemes [107] α DR for fermions in the fundamental representation, and for fermions in the adjoint representation. Notice that scalars do not contribute to (6.8) at three loops. This leads to the following relation for the cusp anomalous anomalous dimension in the two schemes In the special case of N = 4 SYM theory, for n s = 6 and n f = 4, we use the relations (6.8) and (6.9) together with (5.6) to find the three-loop cusp anomalous dimension in DR scheme This confirms a conjecture made in our previous paper [19].
We use (5.5) to find at three loops in QCD We verify that this expression is in perfect agreement with the known result [7,44].
In a similar manner, we obtain an analogous expression in a supersymmetric Yang-Mills theory, with n f fermions and n s scalars in the adjoint representation, and, then, convert the result into the DR scheme with a help of (6.8) to get To obtain from this expression the three-loop light-like cusp anomalous dimension in N = 4 SYM, we adjust the parameters following (5.7), in agreement with [87].

Universal scaling function
We can use the large angle asymptotics of the cusp anomalous dimension (6.11) to introduce a new effective coupling constant a: 19 Inverting this relation we can expand the cusp anomalous dimension in powers of a and define the following function where Γ(φ, α s ) and K(α s ) are evaluated in the same scheme. The expansion coefficients of the two functions are related to each other as

17)
19 It is also known in QCD literature as physical coupling constant [108].
with Ω (i) (φ) being the coefficients of the expansion of Ω(φ, a) in powers of a/π. According to (6.9), the change of the renormalization scheme (from MS to DR) amounts to a finite renormalization of the coupling constant. An immediate consequence of (6.9) is that the coefficients Ω (i) are the same in the two renormalization schemes. This is not the case however for the expansion coefficients of the cusp anomalous dimension, Γ (i) and K (i) . Let us first compute the function Ω(φ, a) in N = 4 SYM. Using (6.10) and (6.14), we obtain from (6.17) By construction, this function takes the same form in MS and DR schemes. Similarly, we can apply the relations (5.6) and (6.12) to compute the corresponding function Ω(φ, a) in QCD and in a generic Yang-Mills theory containing fermions and scalars. Since the cusp anomalous dimension depends on the particle content of the theory, we should expect to find different results for Ω(φ, a). Using the obtained results for the cusp anomalous dimension, we found that the function Ω(φ, a) is independent on the number of fermions and scalars! This remarkable property immediately implies that, at least to three loops, the function Ω(φ, a) is the same in any gauge theory, a) . (6.19) Combining this relation with (6.17), we conclude that all n f and n s dependent terms in Γ(φ, α s ) are generated from lower-loop terms through expansion of K(α s ) in powers of α s . It would be interesting to elucidate the origin of the relation (6.19) as well as its validity beyond three loops. We would like to mention that similar phenomenon has been also observed in other supersymmetric Yang-Mills theories. In particular, various quantities in three-dimensional N = 6 supersymmetric ABJM theory [109] and in N = 2 superconformal Yang-Mills theory [110,111] can be obtained from their counter partners in N = 4 SYM by replacing the coupling constant by the universal 'effective' coupling.
Let us examine the properties of the function Ω(φ, a).
In the large angle limit, for φ = −i log x with x → 0, we combine together (6.11) and (6.16) to see that Ω(φ, a) has universal asymptotic behavior 20) where the coefficient in front of the logarithm does not receive corrections and is one-loop exact, that is Ω (i) = O(x 0 ) for i ≥ 2. Matching this relation into (6.18), we find that the linear combinations ofÃ andB functions that appear in the expansion of Ω(φ, a) at two and three loops remain finite in the large angle limit.
In the small angle limit, for φ → 0, the integration contour in figure 1 reduces to the straight line leading to the vanishing of the cusp anomalous dimension. For small cusp angle φ we expect that Ω(φ, a) = −φ 2 B Ω (a) + O(φ 4 ) , (6.21) where B Ω (a) is an analog of the bremsstrahlung function (2.22). We use (6.18) to obtain the three-loop result As before this function takes the same form in any gauge theory (at three loops at least) and does not depend on the choice of the renormalization scheme. Substituting (6.21) into (6.16) we find for the bremsstrahlung function (2.22) Then, we use the obtained three-loop results (6.22) and (6.12) to get in QCD The two-loop correction to B MS QCD (α s ) agrees with [7,12], the three-loop result is new. In N = 4 SYM we find from (6.23) and (6.14)

The relation to the quark-antiquark potential
As another check of our results, let us consider the limit φ = π − δ with δ → 0, or equivalently x = e i(π−δ) → −1. In this limit, the two rays forming the cusp become antiparallel and the one-cusp anomalous dimension (2.14) develops a pole Γ (1) ∼ −C R π/δ. It is expected that the cusp anomalous dimension should have the same behaviour up to three loops, whereas at four loops it receives corrections of the form (log δ)/δ 20 with V cusp (α s ) = 1 + (α s /π)V (1) + (α s /π) 2 V (2) depending on the renormalization scheme. As before, it is convenient to examine the asymptotic behavior of the universal function Ω(φ, a). Indeed, we find from (6.18) that it develops a pole 1/δ at three loops We note that this relation comes about as a result of nontrivial cancellation of more singular contributions coming from various terms in (6.18). See discussion in section 5. We substitute (6.27) into (6.17) and use the three-loop result for the light-like cusp anomalous dimension (6.12) and (6.14) to verify that the cusp anomalous dimension satisfies (6.26) in QCD and in N = 4 SYM. The corresponding functions V cusp (α s ) are given by Let us compare the relation (6.26) with an analogous expression for color-singlet contribution to the static potential of two heavy color sources carrying the SU (N ) charge C R in generic Yang-Mills theory. In the momentum representation, it has the form where the function V QQ depends on the coupling constant normalized at the scale µ 2 = q 2 with the expansion coefficients a 1 and a 2 known both in QCD [114][115][116] and in N = 4 SYM [117]. In the coordinate representation, the potential is given by withᾱ s = α s (µ 2 = e −2γ E /r 2 ) and ∆V (ᾱ s ) is proportional to the beta-function As was observed in [47], the one-loop correction to (6.28) coincides with analogous correction to heavy quark-antiquark static potential (6.31) in QCD, i.e. a MS 1,QCD = 31 36 C A − 5 9 n f T F . Of course, the coincidence is not accidental and can be understood in the conformal limit of QCD.
Namely, for small δ we can define a conformal transformation x → y that maps two almost antiparallel semi-infinite rays, shown in Figure 1 for φ = π − δ, into two (infinite) lines separated by distance δ. To show this, we assume that the cusp point is located at the origin and introduce the radial and angular coordinates x 0 = r cos φ, x = r n sin φ with x 2 µ = r 2 = e 2ρ and φ = π − δ, so that the metric takes the form where in the last relation we took δ → 0 and introduced new coordinates y 0 = ρ and y = n δ. As follows from the last relation, the transformation x → y is conformal at small δ.
If the conformal symmetry were exact, as it happens in N = 4 SYM theory, the conformal transformation x → y would allow us to identify the Wilson loops evaluated in two different configurations, thus leading to the expected relation between the cusp anomalous dimension (6.26) and the static potential (6.32) for r = δ Note that, in virtue of conformal symmetry, the coupling constant does not depend on the renormalization scale,ᾱ s = α s , and, in addition, ∆V = 0 in (6.32).
In the case of QCD, the conformal symmetry is broken by a nonzero beta function. As a consequence, the Wilson loop receives additional, conformal symmetry breaking corrections under the transformation x → y which generate the difference between the cusp anomalous dimension and the static potential in QCD. Since these corrections are necessarily proportional to the beta-function, we expect that the difference between (6.26) and (6.32) (for r = δ) should be also proportional to β(α s ), see e.g. [118]. 21 Notice that the expansion of the cusp anomalous dimension (6.26) and the static potential (6.32) runs in powers of coupling constant normalized at different scales, α s (µ 2 ) and α s (e −2γ E /r 2 ), respectively. In agreement with our expectations, the difference between the two couplings is proportional to beta-function multiplied by logarithms of the ratio of the two scales. Choosing µ 2 = e −2γ E /r 2 we can eliminate such logarithms and arrive at the following relation 22 with β(α s ) = ( 11 3 C A − 4 3 T F n f )α s /(4π) + O(α 2 s ) and C(α s ) being some function of the coupling constant.
The relations (6.35) and (6.36) can be tested using the known two-loop result for the static potential (6.31) in N = 4 SYM [117] and in QCD [114][115][116]. Replacing V cusp (α s ) by its expressions (6.28) and (6.29), we verified the relations (6.35) and (6.36) and identified 21 The situation here is similar to that for the Crewther relation in QCD. The conformal symmetry breaking corrections to this relation have been studied in [119,120] 22 We did not include ∆V (αs) into this relation since, by definition (6.33), this function is proportional to beta-function and, therefore, can be absorbed into C(αs). the lowest order correction to C(α s ) in the MS scheme where O(α 2 s ) term depends on Γ cusp at four loops. 23 It was found in [121][122][123] that the three-loop correction to the static potential V QCD (r) involves higher SU (N ) Casimirs defined in (6.6). As we argued in section 6.1, the same happens for the cusp anomalous dimension (6.5) at four loops. Applying (6.36) we can relate the corresponding terms order-by-order in the coupling. In particular, assuming that the two-loop correction to (6.37) does not involve higher Casimirs, we can use (6.36) to predict the four-loop correction to Γ cusp (π − δ) proportional to higher Casimirs in δ → 0 limit. Together with (6.5) this leads to the following asymptotic behavior of the functions with the numerical coefficients κ A and κ F obtained in [121][122][123]

Nonplanar corrections at four loops
We recall that nonplanar corrections first appear in Γ(φ, α s ) at four loops and have the general form (6.5). Applying (6.17) we can relate them to nonplanar correlations to the function Ω(φ, a) and to the light-like cusp anomalous dimension (6.15) with Ω (1) = C RÃ1 . In general, four-loop nonplanar corrections ∆Ω (4) and C R ∆K (3) have the same form as (6.5) and are given by a sum of two higher Casimirs. Notice that one of the Casimirs is accompanied by the factor of n f . Assuming that (6.19) is valid at four loops, we find that ∆Ω (4) should be n f independent and, therefore, involve only one Casimir leading to with K A and K F independent of the cusp angle as well as of the number of flavours n f . Substituting these relations into (6.40) and matching the resulting expression into (6.5) we obtain Since K F does not depend on φ, we can fix its value by examining the asymptotic behavior of the both sides of the last relation for φ → π. Taking into account (6.38) together with with κ F given by (6.39). This leads to the following prediction for the n f dependent part of nonplanar correction (6.5) to the cusp anomalous dimension In distinction with f F (φ), the expression for f A (φ) in (6.42) involves in addition the function f Ω (φ) defined in (6.41). Although the explicit form of the function f Ω (φ) is unknown, we can use (6.41) and (6.42) to deduce some of its properties. Namely, examining the asymptotic behavior of both sides of the first relation in (6.42) for φ = π − δ with δ → 0 we find that this function has to satisfy In addition, in the large angle limit, for φ = −i log x with x → 0, it follows from (6.20) and (6.41) that ∆Ω (4) should stay finite in this limit leading to f Ω (φ) = O(x 0 ). This property excludes the possibility for f Ω (φ) to be proportional toÃ 1 .
To summarize, we demonstrated in this subsection that assuming the validity of (6.19) at four loops leads to a definite prediction (6.44) for n f dependent part of the nonplanar correction to the cusp anomalous dimension (6.5).

Comparison with the supersymmetric cusp anomalous dimension
It is instructive to compare (6.10) with the analogous result for the supersymmetric Wilson loop (2.4) with Γ (1) , Γ (2) and Γ (3) defined in (4.73). In comparison with (6.10), this expression depends on the internal cusp angle θ on S 5 . The θ−dependence enters into (4.73) through ξ 0 given by (4.67). For θ = π/2 we find leading to Γ (1) = A 1 , Γ (2) = A 3 and Γ (3) = A 5 . In this way, we arrive at (1), we observe that this expression involves the same coefficient functions as (6.10). However, in distinction with (6.10), it does not vanish for φ = 0 but for φ = θ = π/2, or equivalently x = i. Then, defining It is interesting to analyze the properties of the terms on the right-hand side of this equation. First of all, in the light-like limit x → 0, they have to give a finite limit, since the scalar coupling to the supersymmetric Wilson loop (2.4) is suppressed in this limit. Indeed, we observe thatÃ 2 ,Ã 4 ,B 3 andB 5 all go to constants or vanish in this limit. Second, by definition, they are also well-behaved in the small angle limit, where they modify the coefficients in the Taylor expansion. Third, the limit of the backtracking Wilson line is more interesting. At two loops, the functionÃ 2 has a term ∝ log δ/δ, which is required to cancel a corresponding term inÃ 3 . Such terms are present in the supersymmetric Wilson line operator at two loops due to certain ultrasoft effects, but not in the case of the bosonic Wilson line operator. Likewise, at three loops, the functionsÃ 4 andB 3 are required to cancel 1/δ 2 , (log δ) 2 /δ, and log δ/δ terms not present in the final result. Finally, the functionB 5 just contributes a term 45π 5 /δ in this limit.

Conclusions
In this paper, we computed the angle-dependent three-loop cusp anomalous dimension in QCD and in a supersymmetric Yang-Mills theories. The obtained expressions are rather compact and are given in terms of harmonic polylogarithmic functions that can be readily evaluated numerically. We discussed in detail special physical limits of the cusp anomalous dimension and, in particular, placed special emphasis on the backtracking Wilson line limit that is related to the quark-antiquark potential. We showed that this relation holds in QCD up to a conformal symmetry breaking corrections proportional to the beta function and identified the leading contribution to the conformal anomaly. It would be interesting to investigate whether the latter can be computed from the first principles.
We found that, unexpectedly, the results for the different theories considered are very similar. In fact, up to three loops at least, they can be written in terms of a single universal function evaluated at an effective charge given by the light-like cusp anomalous dimension. Assuming that this property holds at higher loops, we derived the contribution of the n fdependent term that violates Casimir scaling and produces a nonplanar correction to the cusp anomalous dimension at four loops.

A Definition of Yang-Mills theories
Throughout the paper, we consider two Yang-Mills theories with different particle content.
In the first case, for gauge fields coupled to n f species of Dirac fermions, we have where F µν = F a µν T a and D µ = ∂ µ −igA a µ T a with T a being the generators of the fundamental representation of the SU (N ) normalized as The fermion fields ψ i are defined in the fundamental representation of the SU (N ) and carry the additional flavour index i = 1, . . . , n f . In the second case, for gauge fields coupled to n s scalars and n f fermions, we have The reason for choosing the Lagrangian in the form (A.3) is that, by fine tuning the number of fermions and scalars, we can use it to describe supersymmetric Yang-Mills theories with different number of supercharges. In particular, the maximally supersymmetric Yang-Mills theory corresponds to the special case of (A.3) with n f = 4, n s = 6 and matrices T I (with I = 1, . . . , 6) given by (chiral blocks of) Dirac matrices in six-dimensional Euclidean space, T I T † J + T J T † I = δ IJ (see Appendix B in ref. [124] for details).

B Wilson lines and HQET
Wilson lines can be conveniently studied using the heavy quark effective theory. 24 In fact, the HQET Lagrangian was first introduced as a technical device for this purpose [2,4]. The heavy quark effective fields h v (x) and h † v (x) depend on the unit four-vector v 2 µ = 1 which has the meaning of the heavy quark velocity. The correlation function of two HQET fields ( figure 12(a)) is x ν is the projection of x onto the subspace orthogonal to v. Also, W (t) is the expectation value of the Wilson line evaluated along the segment of length t oriented along v µ .
It is convenient to work in momentum space. Introducing the notation for the Fourier transformed HQET fieldh with ω being the so-called residual energy of heavy quark, we find from (B.1) Expanding S v (ω) in powers of the coupling constant yields Feynman diagrams shown in figure 3 for ω = −1/2. Within the HQET framework, the cusp anomalous dimension can be identified as anomalous dimension of the local gauge-invariant operator To see this, we consider the correlation function of two HQET fields and the current (B.4) where x i⊥ stands for the component of x i orthogonal to v i and W (t 1 , t 2 ; φ) is the expectation value of Wilson line evaluated along two segments of lengths t 1 and t 2 separated by a cusp angle φ (see figure 12(b)). In particular, for v 1 = v 2 , or equivalently φ = 0, we have Going to the momentum space, we obtain where ω 1 and ω 2 are the residual energies, S v (ω) is the propagator of the HQET field (B.3) and V (ω 1 , ω 2 ; φ) is the one-particle irreducible vertex function (without the external-leg propagators). For ω 1 = ω 2 = −δ the vertex function V (ω 1 , ω 2 ; φ) coincides with V (φ) in (2.16) and is given by Feynman diagrams shown in figure 2.
It is convenient to extract the renormalization Z−factor from Note that Z(φ) should be gauge invariant and independent on ω 1 and ω 2 . In order to avoid infrared divergences, ω 1 and ω 2 should be different from zero. It is convenient to choose ω 1 = ω 2 = −δ. Then, the dependence of the HQET integrals on δ can be trivially obtained by dimension counting. To simplify the calculation, we can set δ = 1/2 and evaluate the resulting dimensionless integrals in D = 4 − 2ǫ dimensions. The cusp anomalous dimension can be found by matching log Z(φ) into the expected result (2.11). At φ = 0 we find from (B.6) that the vertex function V (ω 1 , ω 2 ; 0) satisfies the Ward identities As a consequence, UV divergences of V (ω, ω; 0) match those of the heavy quark propagator (B.3) leading to log V (ω, ω; 0) = − log Z h + finite , (B.10) where Z 1/2 h is the renormalization factor for the HQET field h v (x). Note that Z h is not gauge invariant; at three loops it has been calculated in [51,52]. We reproduced this result from our three-loop calculation of V (ω, ω; φ) by setting φ = 0 and using (B.10).

C Abelian large-n f terms
In this appendix, we compute the special class of QCD corrections to the cusp anomalous dimension (2.20) of the form (T F n f ) L−1 α L s and C F (T F n f ) L−2 α L s . They originate from QED like diagrams which have the form of the one-loop diagram shown in figure 3(a) with a free gluon propagator dressed by fermion loop corrections (see figures 13(a) and (b)). 25 24 Methods of calculation of multiloop Feynman diagrams in HQET are reviewed, e.g., in [10,96,125,126]. 25 Starting from (TF n f ) L−3 α L s order, we also have to take into account the additional abelian diagrams shown in figure 13(c). They involve the light-by-light scattering and their calculation is more involved To compute the contribution of such diagrams it is sufficient to consider QED with n f massless lepton flavors. In this case, we put C F = T F = 1, C A = 0 and treat the one-loop beta-function as a large parameter. Then, the above mentioned corrections take the form β L−1 0 α L s and β L−2 0 α L s . We shall refer to them as the leading (Lβ 0 ) and next-to-leading (NLβ 0 ) large-β 0 corrections, respectively.
For our purposes we need the expression for the photon self-energy (g µν k 2 −k µ k ν )Π(k 2 ) with the NLβ 0 accuracy Π(k 2 ) = Π 0 (k 2 ) +Π (k 2 ) the next-to-leading correction to (C.2). It can be written in the form [127,128] Π(k 2 ) = 3ǫ where the function F (ǫ, u) is given by with the Euclidean integral (with p 2 = 1) that can be expressed via a hypergeometric 3 F 2 −function of unit argument [129,130]. The function F (ǫ, u) is regular at the origin and admits a double series expansion F nm ǫ n u m , (C. 6) with the coefficients F nm that can be calculated to any order in terms of multiple ζ values. For u = 0, the function F (ǫ, 0) reduces to Euler gamma functions [128] (the same holds for F (ǫ, 2ǫ)). To save space, we do not presentan explicit expression for F (ǫ, u).
It is convenient to introduce the renormalized coupling constant b = β 0 α(µ) 4π . (C.7) In the large β 0 limit, we keep b fixed and use 1/β 0 as an expansion parameter. In the MS scheme the renormalized coupling is defined as where e 2 0 is a bare coupling constant and the charge renormalization constant is given with the NLβ 0 accuracy by Here expansion ofZ α,1 (b) starts from b 2 , that ofZ α,2 (b) from b 3 , etc. The charge renormalization constant satisfies the renormalization group equation that allows us to express Z α (b) in terms of beta-function.
In the abelian theory, Z α (b) is related to the photon self-energy (C.2) expressed in terms of renormalized coupling constant log(1 − Π(k 2 )) = log Z α (b) + O(ǫ 0 ) = −ˆb Then, we use (C.4) to obtain −(1 + b/ǫ)Π(k 2 ) = −3b ∞ L=2 F (ǫ, Lǫ)(b/(ǫ + b)) L−1 /L. Expanding (b/(ǫ + b)) L−1 in powers of b and replacing F (ǫ, Lǫ) with (C.6), we find that all coefficients but F n0 cancel leading tõ . (C.12) We can use this relation to find the β function with NLβ 0 accuracy [127,128] Finally, we substitute β(g) into the last relation in (C.10) and obtain O(1/β 0 ) correction to the charge renormalization constant (C.9) (4 + F 10 ǫ) b 3 ǫ 2 − 1 4 9 + 3F 10 ǫ + F 20 ǫ 2 b 4 ǫ 3 + · · · (C.14) We are now ready to determine the cusp anomalous dimension at the NLβ 0 order. To this end, we have to repeat the one-loop calculation of the vertex function V (ω, ω; φ) (see (2.8) for ω = δ), with a free photon propagator modified by self-energy corrections and, then, express the result in terms of the renormalized coupling constant (C.7). Performing the calculation we obtain where the coupling constant b is defined at the scale µ 2 = e − 5 3 (2ω) 2 , the function F (ǫ, L ′ ǫ) is given by (C.6) and the notation was introduced for The function f (ǫ, u; φ) is regular at the origin: f nm (φ)ǫ n u m . (C. 18) In particular, for u = 0 we have The cusp anomalous dimension is related to the residue at the pole Z 1 (b; φ)/ǫ in the expression (C.16) (C.20) Replacing f (ǫ, Lǫ; φ) in (C.16) with its general expression (C.18), we find that the coefficient in front of 1/ǫ on the right-hand side of (C.16) only depends on the coefficients f n0 (φ) = −2(φ cot φ − 1)f n . Then, at the NLβ 0 order the cusp anomalous dimension is given by coefficient function C(α s ) defined in (6.36) and (6.37). As before, it is sufficient to perform calculations in QED with n f lepton flavors. With the NLβ 0 accuracy, the potential (6.32) is determined by the full photon propagator in the Coulomb gauge, where Π(−q 2 ) is given by (C.2) for q µ = (0, q). Beyond the NLβ 0 order, this relation is modified by corrections due to light-by-light scattering. At the Lβ 0 order, we have from (C.3), (C.8) and (C.9) where the coupling constant b is defined at the scale µ 2 = q 2 . At the NLβ 0 order we use (C. 15) and (C.4) to get Replacing F (ǫ, L ′ ǫ) with (C.6) we find after some algebra where v n = (5/3) n /n! are the expansion coefficients of [D(ǫ)] u/ǫ = n v n u n + O(ǫ). We use the known results [128] for the coefficients F n0 and F 0n to obtain V (q) are in agreement with the known C F (T F n f ) 2 α 3 s and C 2 F T F n f α 3 s terms in the two-loop potential [116], as well as the C F (T F n f ) 3 α 4 s and C 2 F (T F n f ) 2 α 4 s terms at three loops [121]. We can use (C.22) together with (6.26) to determine the function V cusp in the large n f -limit. Comparing this function with the potential (D.1) we verify the anomaly relation (6.36) and compute the corresponding coefficient function Here the expansion can be extended to any desired order. We verify that the n f −dependent term in (6.37) matches the first term in the expression for C 0 (b). Notice that 1/b term cancels in C 1 (b) leading to the absence of the b/β 0 term in C(b). This explains why (6.37) does not contain an abelian color factor C F .