Relating amplitude and PDF factorisation through Wilson-line geometries

We study long-distance singularities governing different physical quantities involving massless partons in perturbative QCD by using factorisation in terms of Wilson-line correlators. By isolating the process-independent hard-collinear singularities from quark and gluon form factors, and identifying these with the ones governing the elastic limit of the perturbative Parton Distribution Functions (PDFs) -- $\delta(1-x)$ in the large-$x$ limit of DGLAP splitting functions -- we extract the anomalous dimension controlling soft singularities of the PDFs, verifying that it admits Casimir scaling. We then perform an independent diagrammatic computation of the latter using its definition in terms of Wilson lines, confirming explicitly the above result through two loops. By comparing our eikonal PDF calculation to that of the eikonal form factor by Erdogan and Sterman and the classical computation of the closed parallelogram by Korchemsky and Korchemskaya, a consistent picture emerges whereby all singularities emerge in diagrammatic configurations localised at the cusps or along lightlike lines, but where distinct contributions to the anomalous dimensions are associated with finite (closed) lightlike segments as compared to infinite (open) ones. Both are relevant for resumming large logarithms in physical quantities, notably the anomalous dimension controlling Drell-Yan or Higgs production near threshold on the one hand, and the gluon Regge trajectory controlling the high-energy limit of partonic scattering on the other.


Introduction
It is well known that perturbative QCD at fixed order in α s , which is highly successful in describing hard processes at colliders, loses its predictive power in kinematic regions where there is a large hierarchy of scales. Familiar examples are Drell-Yan or Higgs production near threshold, see e.g. [1][2][3][4][5][6], or at small transverse momentum, which are dominated by soft-gluon radiation. Another example is the high-energy limit of QCD scattering, where the centre-of-mass energy is much larger than the momentum transfer [7][8][9][10][11][12][13][14]. In each of these cases, and many others, factorisation techniques allow us to derive all-order resummation formulae, which extend the predictive power of QCD, leading to highly successful phenomenology in many cases.
The theory underlying factorisation relies on identifying the origin of any parametrically-enhanced corrections through operators, which capture the relevant divergences. Independently of whether one uses QCD fields [15,16], or Soft-Collinear Effective Theory [17] ones, the relevant operators involve Wilson lines, which follow the trajectory of fast-moving partons, and capture their interactions with soft gluons. These operators obey evolution equations, governed by corresponding anomalous dimensions, which are computable order by order in QCD perturbation theory. The most familiar amongst these is the (lightlike) cusp anomalous dimension [18], γ cusp , which in particular describes double poles in the Sudakov form factor, originating in overlapping soft and collinear singularities. While the cusp anomalous dimension occurs universally, governing the leading singularities in any kinematic limit, single-logarithmic contributions characterising separately large-angle soft or hard-collinear or rapidity divergences, are somewhat less universal, and yet -as we shall see -recur in a variety of physical quantities that are not a priori related. Resummation formulae are obtained upon solving the aforementioned evolution equations, leading to exponentiation. The anomalous dimensions therefore have a central role in the predictive power of QCD, and in certain cases their computation has been recently pushed to threeloop order, e.g. [19][20][21][22][23][24], with very recent progress towards four loops [25][26][27][28][29][30][31][32][33] (even more is known in maximally supersymmetric Yang-Mills theory, see e.g. [34][35][36][37][38][39]). Despite this impressive progress, there remains several unresolved questions regarding the anomalous dimensions governing single-logarithmic corrections and their universality, some of which we address below.
In the present paper we study two fundamental physical quantities, which are recurrent ingredients in the factorisation of amplitudes and cross sections [16,40]. The first is the massless on-shell form factor, associated e.g. with an electromagnetic vector current in the case of quarks, or effective Higgs production vertex, gg → H, in the case of gluons. The second is parton distribution functions (PDFs), or more precisely, the large-x limit of diagonal qq and gg Altarelli-Parisi splitting functions, governing the scale dependence of PDFs according to the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation [41][42][43]. Each of these physical quantities is important in its own sake, and their infrared factorisation will be discussed in some detail in sections 2 and 3, respectively. The main motivation to our study comes from the relation between the two, namely a particular combination of single-pole anomalous dimensions, which respectively capture collinear singularities in these two quantities. The relation holds separately for quarks and for gluons: where γ q G (γ g G ) is defined by the function G (see eq. (2.6)), which along with the cusp anomalous dimension, governs the infrared structure of the quark (gluon) form factor in eq. (2.1) below; and B q δ (B g δ ) is the coefficient of the δ(1 − x) term, in the large-x limit of the quark-quark (gluon-gluon) splitting function, see eq. (3.6) below. It was observed long ago [44,45] that while the separate perturbative results for γ G and B δ are very different between quarks and gluons (this is expected: collinear singularities are known to depend on the parton's spin), the combination (1.1) vanishes at one loop in both cases, and admits a Casimir-scaling relation 1 at two loops, namely The same Casimir-scaling property persists at three loops [45]. This is a clear indication that f eik has an interpretation purely in terms of Wilson lines -hence the name, an eikonal function. A Wilson-line-based definition would explain why the result does not depend on the parton's spin, while it depends on its colour representation in proportion to the relevant quadratic Casimir through three loops. The question we would like to address is what is the Wilson-loop correlator corresponding to f eik .
Before describing our approach to answer this question, let us note that the combination in (1.1) has a direct physical interpretation as the soft anomalous dimension associated with Drell-Yan production near partonic threshold [1][2][3][4][5][6], namely γ q G − 2 B q δ = 1 2 Γ DY . Similarly γ g G − 2 B g δ is associated with Higgs production through gluon-gluon fusion near threshold. The corresponding soft function is defined at cross-section level, by replacing the energetic partons, which move in opposite lightlike directions (before annihilating at the hard interaction vertex), by Wilson lines that follow the same trajectory, in both the amplitude and its complex conjugate. The cusp where the complex-conjugate amplitude Wilson lines meet is displaced by a timelike distance with respect to the amplitude: this distance is the Fourier conjugate variable to the energy fraction carried by soft partons. 2 Final-state radiation, namely the set of soft particles connecting the amplitude side to the complexconjugate amplitude side, are described by cut propagators. This soft function admits an evolution equation governed by γ cusp and Γ DY (see e.g. eq. (9) in ref. [3], or eqs. (43)(44) in ref. [6]). The latter was computed through three loops directly based on the aforementioned Wilson-line definition [20,52], and the results agree with the combination of anomalous dimensions in (1.1), which were extracted from independent QCD computations of the form factor [44,45,53] and DGLAP splitting functions [54][55][56][57][58][59][60][61]. Thus, from this perspective, this physical quantity is well understood, and its Casimir-scaling property simply follows from the above-mentioned Wilson-line definition.
Our own investigation starts with the simple observation that the two-loop result for γ G −2 B δ in (1.2) also agrees, up to an overall factor of 4, with the result for the parallelogram Wilson loop made of four lightlike segments (see figure 1c), which was computed in 1992 by Korchemsky and Korchemskaya [62]. It is a highly appealing proposition that 3 holds to all orders. 4 The parallelogram Wilson loop, is a very simple object: being compact it has no infrared divergences, so the singularities arise here from short distances, and the calculation can be done directly in dimensional regularisation. Importantly, in contrast to the Drell-Yan soft function described above, real corrections and cut propagators do not arise here. The natural questions to ask then are first, does the relation in (1.3) indeed hold to all orders, and second, can we see how a parallelogram Wilson loop arises from the definitions of the objects on the left-hand side of eq. (1.3), the form factor and the PDF. Establishing this relation is one of the main goals of the present paper. The infrared factorisation of the form factor is well understood [16,40,63], and has been used as the starting point for the factorisation of massless amplitudes with any number of legs in general kinematics [64][65][66][67][68][69][70][71]. The form-factor factorisation gives rise to a different Wilson-line configuration, namely a couple of semi-infinite lightlike Wilson lines (with different 4-velocities) meeting at the hard-interaction vertex, see figure 1a. We shall refer to this contour as the ∧ geometry. We emphasise that in contrast with the Drell-Yan soft function described above, where the cross section was considered [20,52], here the Wilson-line configuration is defined at amplitude level. At a difference with the parallelogram of [62], the ∧ geometry is non-compact, and thus gives rise to infrared divergences, in addition to ultraviolet ones. We shall return to the ∧ geometry and its properties below. At this point it suffices to say that considering the infrared factorisation of the form factor, the origin of the relation between γ G − 2B δ and the parallelogram geometry remains obscure: the ∧ geometry has no finite segments while the parallelogram consists exclusively of such.
An important step in explaining the eikonal nature of f eik in (1.1), based on the infrared factorisation properties of the form factor and the PDF, was taken in 2008 in a paper by Dixon, Magnea and Sterman [63]. The fundamental explanation is that spin-dependent hard-collinear contributions are common to both γ G and 2 B δ and drop in the difference, leaving behind a purely eikonal component. This is the premise we shall follow here as well. However, ref. [63] relied on the assumption that B δ , as the coefficient of δ(1 − x), is a purely virtual quantity and hence the factorisation of the PDF could be done at "amplitude level". According to the factorisation outlined in [63] the eikonal component of B δ should correspond to Wilson lines with a ∧−geometry, much like the form factor. Taking this at face value, if the eikonal components of γ G and B δ on the right-hand side of (1.3) indeed both correspond to the ∧−geometry, one concludes that the ∧ and the anomalous dimensions must be proportional to each other, at least through two loops, or, put differently, one may deduce the anomalous dimension of the ∧−geometry from (1.2).
The first direct two-loop computation of ∧−geometry Wilson loop was performed only in 2015, by Erdogan and Sterman [72]. This calculation is an important step forward also in the sense that it presents a new method for dealing directly with (semi)-infinite lightlike Wilson lines in configuration space (which a priori lead to scaleless integrals) without resorting to an extra regulator. This is done by cleverly using the exponentiation properties and isolating a well-defined integrand, before renormalising ultraviolet divergences by means of a suitable cutoff. We shall adopt and generalise this method in section 4 below. The result of ref. [72] is that the anomalous dimension corresponding to the ∧−geometry Wilson loop is given by where C i = C F for Wilson lines in the fundamental representation and C A for the adjoint. As with f eik and Γ ∧ above, we omit the superscript q/g for Γ ∧ wherever it is not necessary. While the result in (1.4) bears a striking resemblance to f eik in (1.2), it is evidently not identical; the coefficient of the ζ 3 term is entirely different. The authors of ref. [72] further provided a detailed diagrammatic analysis, comparing their calculation to that of the parallelogram in ref. [62], and explaining the origin of the difference in the coefficient of ζ 3 as emanating from endpoint contributions that are present in finite lightlike segments, but are absent in infinite ones. This conclusion can be confirmed by a momentum-space computation.
It is useful to bear in mind that infinite and semi-infinite Wilson-line configurations (but not finite ones!) are of direct relevance to partonic scattering amplitudes in the high-energy limit (the Regge limit) [11][12][13][14]. Also, the explicit two-loop combination in (1.4) appeared in the literature in that context long before the computation of the ∧ configuration in ref. [72]. Specifically, considering gg → gg, qq → qq or qg → qg scattering in the limit where the centre-of-mass energy is much larger than the momentum transfer, s ≫ −t, the leading and next-to-leading logarithms in s/(−t) in the (real part of the) amplitude exponentiate according to a simple replacement of the t-channel gluon propagator (dubbed gluon Reggeisation): where α(t, ǫ) is the gluon Regge trajectory 5 [74][75][76][77][78] given by: where α s = α s (−t, ǫ), with ǫ = (4 − d)/2 the dimensional regularisation parameter,b 0 is the one-loop QCD beta function of (2.3a), γ g (n) cusp are the coefficient of the cusp anomalous dimension of eq. (2.4) for the gluon, and Γ g (2) ∧ is the two-loop coefficient in eq. (1.4), again with C i = C A . We further recall that the overall similarity between the parallelogram Wilson loop in [62] and the gluon Regge trajectory in (1.6), as well as the peculiar difference between them in the coefficient of ζ 3 , were already observed early on, in ref. [79], where an evolution equation for the Regge trajectory was derived, considering the forward limit of crossed Wilson lines. However, this raises no difficulty: as stressed above, it is the infinite Wilson-line geometry which is expected to be relevant for the factorisation of partonic scattering amplitudes, not the parallelogram.
A real puzzle arises, however, upon considering the explicit result for the ∧−geometry anomalous dimension in eq. (1.4) in view of eq. (1.2), if the conclusion of ref. [63] is taken at face value. Given that the factorisation of the form factor is well understood, and the eikonal component of γ G is determined by the ∧−geometry, we are compelled to revisit the assumption of ref. [63] that B δ is a purely virtual quantity, systematically establish the infrared factorisation of the PDFs at large x, and identify the eikonal component of B δ , which clearly must not be proportional to Γ ∧ .
We proceed as follows. In section 2 we review the factorisation of long-distance singularities of the QCD form factor and identify the process-independent spin-dependent hard-collinear component of γ G . In turn, in section 3 we discuss the factorisation of PDFs in the limit x → 1. We show explicitly that the calculation of B δ requires both real and virtual corrections. To this end we perform an explicit two-loop calculation of the splitting functions at large x (the details are presented in appendix A). Next we identify the eikonal component of B δ as the anomalous dimension associated with a ⊓-shaped Wilson-line geometry, see figure 1b. By using the known value of B δ along with the hard-collinear anomalous dimension extracted from the form factor, we then predict the Γ ⊓ anomalous dimension at two loops. Then, in section 4 we compute Γ ⊓ directly to this order, finding agreement with the extracted result of section 3. In section 4 we also derive an evolution equation for the ⊓-shaped Wilson-line and show that while in the ultraviolet it is characterised by double poles, as any other cusped Wilson loop, its infrared properties are different, displaying strictly single poles, in agreement with single-pole nature of PDFs themselves. In section 5 we put together our results for the factorisation of the form factor and the PDF, and establish the relation of (1.3) with the parallelogram to all orders. We further summarise the 5 See also a more recent observation in ref. [73] that the two-loop coefficient Γ g (2) ∧ occurs also in the QCD impact factor. state-of-the-art knowledge of higher-order corrections to Γ in view of its relations with other physical quantities. We briefly summarise our conclusions in section 6.

Infrared factorisation of the on-shell form factor
Let us review the well-known factorisation of a colour-singlet on-shell form factor of coloured massless particles (quarks or gluons) in QCD [15,16,40,45,63,80]. We label the external momenta by p 1 (incoming) and p 2 (outgoing) with the momentum transfer Q 2 ≡ −(p 1 −p 2 ) 2 , and, as usual, we renormalise all ultraviolet singularities in the MS scheme, denoting the renormalisation scale by µ 2 .
The quark form factor is defined in terms of the electromagnetic vector current, proportional toψγ µ ψ, which does not renormalise. The gluon form factor in turn is defined using an effective local interaction vertex with the Higgs field, HG a µν G µνa , and it does renormalise, proportionally to the QCD beta function [44]. The distinct ultraviolet properties of the quark and gluon form factors will be of little relevance for us: we focus instead on the infrared singularities of the form factor, which have a rather similar structure for massless quarks and gluons.
For large Q 2 the form factor F Q 2 /µ 2 , α s (µ 2 ), ǫ features large logarithms in the ratio Q 2 /µ 2 , and fixed-order perturbation theory breaks down. These large logarithms can be resummed using a renormalisation-group equation (see e.g. [80]), giving the following allorder formula for the form factor, where we set the renormalisation scale µ 2 = Q 2 for simplicity. Note that we have absorbed into the function G any operator renormalisation terms -see [44,45] for more details. Infrared singularities are generated in eq. (2.1) through an integration, from λ 2 = 0, over the d = 4 − 2ǫ dimensional running coupling α s (µ 2 , ǫ), which obeys We report the coefficientsb 0 ,b 1 andb 2 of the QCD beta function respectively at one [81][82][83][84], two [85][86][87][88] and three loops [89,90], because we will use them in the rest of this paper with the quadratic Casimir defined by T a T a = C R 1, with T being the SU(N c ) generator in the representation R and C A corresponds to the Casimir in the adjoint representation; n f is the number of light quarks and the normalisation of the generators t a in the fundamental representation, Tr(t a t b ) = T f δ ab , is conventionally set to T f = 1/2. Equation (2.1) applies for both quarks and gluons, but with distinct functions γ cusp (α s ) and G(Q 2 /µ 2 , α s , ǫ), which do depend on the type of particles (although this is suppressed in our notation). The former, which captures all double poles, depends solely on the colour representation of the particles (fundamental and adjoint for quarks and gluons, respectively) while the latter, which controls single poles, depends also on their spin. This distinction will be crucial in what follows and it is a direct consequence of the fact that γ cusp is an eikonal quantity, namely one that can be defined exclusively in terms of Wilson lines, while G(Q 2 /µ 2 , α s , ǫ) instead, contains hard-collinear effects, which cannot fully be described by Wilson lines. Specifically, γ cusp is the lightlike cusp anomalous dimension [18], defined as the coefficient of the leading ultraviolet divergence occurring in a cusped Wilson loop, which evaluates to where C i , defined above, is the quadratic Casimir in the fundamental or the adjoint representation for quarks and gluons, respectively. The three-loop value of γ cusp was computed in [91], and recently there has been significant progress towards a four-loop determination [29][30][31][32][33]. Through three loops, the cusp anomalous dimension, much like other quantities that are defined exclusively in terms of Wilson lines, depends on the colour representation proportionally to the quadratic Casimir C i , as in (2.4), adhering to the so-called Casimir scaling property. Starting at four loops quartic Casimirs, d ij ≡ d abcd i d abcd j , appear as well, making the dependence of the colour representation more involved. Differently from γ cusp , the function G(1, α s (λ 2 , ǫ), ǫ) has an expansion both in α s and ǫ, as follows therefore it generates both infrared poles and non-negative powers of ǫ upon integrating over the scale λ 2 of the running coupling as in eq. (2.1). We isolate the divergent contribution order-by-order in α s , by defining the anomalous dimension γ G such that where γ G depends on ǫ only through the coupling. The coefficients γ G for the quark and for the gluon are well known in the literature; they are referred to sometimes as "collinear anomalous dimensions" and were denoted by G in [92], by G 0 in [39] and by γ q and γ g in appendix I of [17]. The latter has a conventional factor of −2. In practice, we derive here γ G to four loops by substituting the d−dimensional running coupling of eq. (2.2) into eq. (2.6) and then identifying the singularities arising on the two sides of equation (2.6), getting where G(l, n) are defined in eq. (2.5) and their values can be extracted from refs. [44,45,53] where the form factors have been computed to three loops. For the purpose of this paper we only need explicit results for the collinear anomalous dimensions through two loops, which read where we added superscripts i = q, g to distinguish between quarks and gluons.

Infrared factorisation
At high energy (Q 2 → ∞) the infrared behaviour decouples from the hard scattering where the jet function J i , one for each external leg, captures the collinear singularities, the soft function S contains the contribution of any long-wavelength gluons exchanged between the external particles and the eikonal jet function J i captures all the singularities that are present both in S and in the jet function J i , which are associated with exchanges that are both soft and collinear to the massless external particles. Therefore, the ratio S J 1 J 2 in eq. (2.9) includes only the divergences associated with soft wide-angle emissions. H is the hard function, found from matching to the full form factor. Each other factor in eq. (2.9) has an operator definition which dictates their functional dependences in eq. (2.9), involving the momenta p i of the external particles and their lightlike velocities β i , defined by where Q 0 is an arbitrary normalisation and would typically be of the order of the hard scale of the process, Q. The operator definitions of S, J i and J i are written in terms of the expectation values of Wilson lines, defined as where v is the direction of the line and x and y are its endpoints. In general, the vector v can be either lightlike v 2 = 0 or non-lightlike v 2 = 0. In the context of the on-shell massless form factor, lightlike kinematics for the external legs, β 2 i = 0, is dictated by eq. (2.10), and we define the functions entering the factorisation formula (2.9) by: where n i is an auxiliary non-lightlike vector and the dependence on its choice must cancel in eq. (2.9). The contour defining S is shown in figure 1a. In eq. (2.14) we presented the jet function J i for fermion fields; for a definition of the gluon jet function see refs. [93][94][95].
The representation of the Wilson lines in eq. (2.13) is the representation of the corresponding external particle. Any function built solely from Wilson lines, such as S and J i , is called eikonal. As mentioned in the context of the cusp anomalous dimension, one of the properties of eikonal quantities is that they admit Casimir scaling up to three loops; this is a consequence of non-Abelian exponentiation. Beyond three loops there are quartic (and eventually higher order) Casimir contributions, but given that the same Wilson-line diagrams contribute for quarks and gluons, differing just by the representations of the Wilson lines, one expects a relation between these quantities. Indeed, a conjectural relation was proposed in [29] based on partial four-loop computations; we shall return to this in section 5.2 below. The individual functions in eqs. (2.12)-(2.14) are heavily constrained by kinematic considerations, such as the dependence on the auxiliary vectors n i , and by renormalisation group evolution. These constraints can be solved to give explicit formulae [70,71], where Γ J and Γ ∧ are constants to be determined by direct calculation. Note that Γ ∧ was denoted in the literature [63,72] as −G eik . As in eq. (2.1), the infrared singularities of J i and S are generated by integrating over the d dimensional running coupling α s (λ 2 , ǫ) from λ 2 = 0. We notice that the soft function and the product of the eikonal jets share the same dependence on γ cusp ln µ 2 /λ 2 , which is associated with the overlapping soft-collinear singularities of these two quantities. This fact ensures that the ratio S J 1 J 2 is free of overlapping divergences and depends only on the logarithm of the kinematic variable which is insensitive to the normalisation of the vectors β i in the definition eq. (2.10). Using the factorisation equation eq. (2.9), we determine the partonic jet function by dividing the form factor in eq. (2.1) by the ratio S J 1 J 2 , yielding is a matching coefficient that captures the finite parts of the jet function and γ i , with i = q for the quark and i = g for the gluon, is the anomalous dimension of the field i in axial gauge. The latter is only concerned with the ultraviolet behaviour of the jet function and indeed it is not associated with any IR pole, because the contribution from the IR region λ 2 ≃ 0 is absent in the second term of eq. (2.18). All the IR poles of the form factor are generated by the second integral in the equation above, involving the anomalous dimensions γ cusp , Γ ∧ , Γ J and the resummation function G(1, α s , ǫ). The dependence on γ cusp is such that the combination with S J 1 J 2 reconstructs the kinematic dependence of the form factor eq. (2.1) through definition in eq. (2.6), and we get the ratio where on the last line we have defined the anomalous dimension γ J/J As mentioned above, the collinear anomalous dimension γ G is known to three loops [44,45,53] for both quarks and gluons, and we quoted the corresponding expressions through two loops in eq. (2.8). The anomalous dimension Γ ∧ , in turn, is derived from the renormalisation of the soft function S, that can be read off eq. (2.16) The equation above clarifies the meaning of the subscript ∧, which symbolises the contour of the lightlike Wilson loop in the definition of the soft function in eq. (2.13) that defines Γ ∧ . This notation will be used throughout this paper and it will be generalised for different contours. Γ ∧ is known to two loops [72] by direct computation of the equation above where C i is the quadratic Casimir dependent on the representation of the Wilson lines in eq. (2.13). Using the results in eqs. (2.8) and (2.24) we determine γ J/J to two loops. First for quarks we have Then for gluons, We have thus isolated the hard-collinear singularities of the form factor and found the quantity γ J/J that governs this behaviour for quark and for gluons according to eq. (2.21). We emphasises that in contrast to the conventional collinear anomalous dimension γ G given in eq. (2.8), which is specific to the form factor (recall eqs. (2.6) and (2.1)), the hard-collinear anomalous dimension γ J/J defined here is process independent. This universality will now be put to use. In the next section we will consider the factorisation of parton distribution functions (PDFs) at large x where we will use the above two-loop results for γ q J/J and γ g J/J given in eqs. (2.25) and (2.26) respectively, and ultimately identify the eikonal anomalous dimension relevant to the PDF evolution.

Parton distribution functions at large x
Parton distribution functions, f AB (x), describe the probability of finding parton A with momentum fraction x inside hadron (or parton) B. We will be interested here in PDF evolution, which is the same for the partonic and for the hadronic quantities, and will therefore consider partonic PDFs. PDFs are inherently defined at cross-section level with the need to combine real and virtual radiation to cancel soft singularities such that only pure collinear singularities associated with the massless initial-state parton are kept. We will see that in the elastic limit, x → 1, the contributions from different regions factorise and claim that the hard-collinear behaviour of the initial-state partons is described by γ J/J , the same anomalous dimension we identified in the factorisation of the form factor.

Definition
The light-cone PDF for a quark (gluon) in a parton P of momentum p with longitudinal momentum fraction x is given by [96] f bare The Wilson-line operator W u is defined in eq. (2.11) and |P is either an on-shell quark or gluon, P = q, g. We take the lightlike momentum p to be in the (+) direction and then the velocity four-vector u is in the (−) direction. It is worthwhile noting here that the bare PDFs f bare j ′ j (x, ǫ) are scaleless. This will be used later in the context of factorisation. They are renormalised through a convolution, where Z jj ′ is a renormalisation factor, removing the UV divergences from the bare PDF in the MS scheme and f jk is the renormalised PDF. From Z jj ′ (x, α s , ǫ) we can get the splitting functions, The RG evolution of the PDFs is governed by the DGLAP equations [41][42][43]: The DGLAP splitting kernels P jk are known to three loops [43, 54-61, 91, 97-100] with some recent results at four loops [25,27,29].

Perturbative calculation at large x
In the limit x → 1 the diagonal terms in the splitting functions, P qq and P gg , feature divergent contributions [46,[101][102][103], namely where the label i = q, g indicates quarks and gluons, respectively, and the plus distribution is defined as usual, see e.g. [43]. The splitting functions are determined from the UV singularities of the PDFs defined in eqs. (3.1) and (3.2), which can be computed perturbatively. We can relate these definitions to time-ordered products by the discontinuity in x, This relation, which is illustrated diagrammatically in figure 2, can be derived as follows.
One first splits the Wilson line in eq. (3.1) into two Wilson lines that extend to infinity, W u (y, ∞)W u (∞, 0), one then inserts a complete set of states between them and finally identifies the result as the discontinuity of the time-ordered product. This relies on the fact that the condition x ≤ 1 selects the cuts with positive energy [16,104]. One can think that the coefficient B δ in eq. (3.6) is entirely determined by the contribution of the virtual diagrams, such as the second term in the left-handside in figure 2, however the explicit calculation will lead to a different conclusion. At one loop, the relevant diagram is shown in the right-hand side of figure 2, which in Feynman gauge reads (3.8) where we used p and k respectively to denote the incoming and outgoing quark momenta, and q the gluon momentum. For brevity, we also drop the superscript bare. It is straightforward to compute the integral over q − by complex analysis. This places a bound on q + i.e. p + > q + > 0. The q T integral is scaleless but as we are interested only in the UV divergence it is simply a matter of replacing, We then scale out p + by defining q + = p + w to produce an elegant integral representation, where we have absorbed the (e γ E 4π) ǫ factors in the MS coupling. The representation in eq. (3.10) has the advantage of compactly displaying the sum over cuts: individual cuts can be isolated by computing the residues corresponding to each of the propagator poles. Using partial fractioning, so the full discontinuity of the integrand equals, and we find Here the first term is a real emission cut, while the second, a virtual correction. As usual, the endpoint divergence in the first term is combined with the divergence as w → 0 in the second to give, (3.14) We emphasise that it is ambiguous to determine which cuts have contributed to the δ(1− x) term, as its coefficient is only finite after the cancellation of the soft divergences between the real and the virtual cuts. We combine eq. (3.14) with the mirror diagram representing the correction of the right vertex, which yields an identical result, and with the box-type diagram, which does not contribute to divergent terms at large x. We complete the calculation of the (bare) PDF by including the two diagrams featuring radiative corrections on the external legs where we used the wavefunction renormalisation Z 2 at one loop. The expression of the UV singularities of the bare PDF at one loop reads Following eq. (3.3), we derive the renormalisation factor Z qq that cancels the ultraviolet divergence in the equation above Finally, we obtain the splitting function by computing the derivative with respect to the renormalisation scale eq. (3.5), which yields the well-known result for the qq splitting function The one-loop calculation with on-shell states is straightforward but at two loops and beyond it becomes complicated to disentangle the UV from the IR in the transverse integrals. To regularise the IR we can take the initial states to be off-shell p 2 = 0. The intermediate expressions become more verbose but introduce no major conceptual issues. As the states are now unphysical the correlators become gauge dependent. It means that the running of the gauge parameter, ξ → Z A ξ has to be taken into account in O(ǫ 0 ) finite terms. A similar observation was made in [32]. Using this method we are able to arrive at the integral representation similar to eq. (3.10) for each two-loop diagram. For two loops it is a two parameter integral with integrals over the plus component of the two loop momenta. As an example, the diagram in figure 3 can be represented as .
The three denominators correspond to the three Wilson-line propagators after integration over the (−) and transverse components of the two loop momenta. We distinguish the contribution of the real emission and the ones of the virtual corrections by applying partial fractioning as in eq. (3.11). The discontinuity of the first propagator in eq. (3.19) is proportional to δ(1 − x) and it determines the virtual contribution. The other two propagators in eq. (3.19) correspond to real emissions. Each term features infrared divergences, which cancel in the sum of all cuts. Furthermore, we notice that the real emission cuts yield UV poles that are proportional to δ(1−x) and therefore contribute to the P qq splitting function. This particular calculation is detailed in section A, where we also present the full two-loop results for quarks and gluons, diagram by diagram. Our final result for the splitting functions eqs. (A. 19) and (A.22) reproduces the known results [43,[54][55][56][57][58][59][60][61][97][98][99]. These previous splitting function calculations have been performed using different methods, including extracting them from corresponding deep inelastic structure-function calculations [61], by means of the operator product expansion [54-57, 60, 97-99], by means of light-cone axial gauge [105,106], or by relating them to splitting amplitudes [107]. To our knowledge, our direct calculation is the first of its kind. This method has the advantage to show that not all the diagrams contribute to the singular behaviour of the splitting functions in eq. (3.6) and that the coefficient B δ includes both the virtual and the real corrections.

Factorisation
As x → 1 the momentum of the final-state parton tends to the initial-state one, meaning that the contribution from soft gluon radiation dominates. It then implies a factorisation of the renormalised PDFs at large x, allowing us to separate the hard-collinear divergences from the soft divergences [46,108]. In the following we shall only consider diagonal splitting functions and since the formulae apply to both quarks and gluons we shall drop the subscript jj on the partonic PDF and related quantities and only specialise when needed. To factorise the PDFs we shall transform into Mellin space, where convolutions become products. In this space the divergent terms become, The large-x limit corresponds to the large-N limit. The factorisation works in much the same way as the form factor by defining two jet functions and two corresponding eikonal jet functions along with a soft function [46,108], where the four-velocity β is in the p direction and L and R indicate which side of the cut the jet functions are (see figure 2). The renormalised parton distribution functions are defined as pure counterterms in minimal subtraction schemes, because they can only depend on the factorisation scale. Since the hard function H and the jet functions J i are the only functions with finite terms it must mean that their non-divergent terms cancel such that eq. (3.22) contains only poles, where J| pole has the same meaning as in eq. (2.20), that it is only the poles of the jet function. As in the case of the form factor, the soft functionS ⊓ resums the emission of gluons with vanishing momenta in the eikonal approximation. We shall shortly see however that while its ultraviolet behaviour is qualitatively the same as that of the form-factor soft function in eq. (2.13), its infrared behaviour is qualitatively different, as it presents only single poles.
The functionS ⊓ is defined by the Mellin transform of the x−space soft function where W ⊓ is the Wilson loop with ⊓−shaped contour, see figure 1b (in ref. [46] it is defined in the axial gauge), Note that the time-ordering operation here acts on the product of the three Wilson lines together. The soft function can be written in this way, despite coming from a cross-section definition because of the particular relation between path-ordering and time-ordering [101]. The definition (3.24) determines two important properties concerning the analytic structure of W ⊓ , as argued in [101]. First of all, the soft function has support in the physical region with x ≤ 1 only if the singularities of W ⊓ are located on the positive imaginary axis in the complex y−plane. Indeed, if this is the case, for x > 1 we can close the integration contour in y in eq. (3.24) through the lower half-plane getting a vanishing result. Furthermore, the reality of soft function implies that W ⊓ is unchanged by the transformation y → −y followed by complex conjugation. Both these conditions are satisfied if W ⊓ is a holomorphic function in the variable In section 4 we show that we can write the renormalised W ⊓ as, where the factor √ 2 was introduced in order to identify µ as the MS renormalisation scale. The quantity Γ ⊓ will admit Casimir scaling to three loops and the scaling is determined by the representation of the Wilson lines in eq. (3.25). Following ref. [101], the soft functionS in the limit of large N , which is conjugate to the behaviour of W ⊓ at large y through the Fourier transform in eq. (3.24), is obtained to leading power in N by replacing y → −iN in eq. (3.27), which leads to so thatS ⊓ admits the following evolution equation Note that the UV behaviour ofS ⊓ is double logarithmic: the right-hand side of eq. (3.29) is dominated by γ cusp (α s (µ 2 )) log µ 2 and therefore it has the same UV behaviour as the one of the form-factor soft function S in eq. As before with the form factor, we seek to isolate the hard-collinear and the purely soft contributions from the Mellin transform (3.20) of the splitting functions in eq. (3.4), P (N, α s ). The following argument is in the spirit of [63]. As mentioned earlier, the bare PDFsf bare (N, ǫ) formally vanish because they are scaleless in dimensional regularisation [109]. They feature UV divergences which are renormalised by the splitting func-tionsP (N, α s ) throughZ(N, α s , ǫ), see eq. (3.3). They are also infrared divergent because there are massless on-shell incoming partons. The IR divergences are the same as in the renormalised PDFs described by eq. (3.23). In perturbation theory it must mean that iñ f bare (N, ǫ) the IR poles match the UV poles. In a minimal subtraction scheme the factor Z in eq. (3.3) consists of only poles. We are then able to constructf bare (N, ǫ) in a way that separates the UV from the IR, The kinematic dependence in the argument of the logarithm cancels upon identifying p·n β·n = p·u β·u = p + β + . We now require the Mellin transform of eq. (3.6) at large N [46,[101][102][103], Substituting this into eq. (3.32) the dependence on γ cusp drops. This shows that the factor √ 2 present in eq. (3.27) is indeed necessary for µ to be identified as MS scale. Comparing the non-logarithmic terms in eqs. (3.32) and (3.33) we finally arrive at the relation, (3.34) The above equation mirrors the form factor equation for γ G in eq. (2.22). In both equations the same hard-collinear anomalous dimension γ J/J is present. We now proceed to use its universality to extract Γ ⊓ at two loops from the above equation. As in the form factor case to specialise to quarks or gluons we simply add a superscript i = q, g. Up to two loops the expressions for B δ may be read off the results in eq. (A.19) and eq. (A.22) of the calculation in the appendix, in agreement with refs. [54][55][56][57][58][59][60][61]. They read Substituting these results into eq. (3.34) along with the values of γ i J/J calculated in eq. (2.25) and eq. (2.26) we arrive at the same quantity for Γ ⊓ for quarks and gluons up to an overall Casimir scaling: The fact that Casimir scaling is recovered is expected of course, as this quantity is defined by Wilson lines. Nevertheless, recovering it by subtracting non-eikonal quantities is a non-trivial consistency check. It is worthwhile noting that only the ζ 3 term is different between Γ⊓ 2 and Γ ∧ in eq. (2.24). The factor of two is present because there are two cusp contributions for the ⊓ contour as opposed to one for the ∧ contour. The different coefficient in front of ζ 3 will be discussed further in section 5.
We have found the anomalous dimension that controls the non-collinear soft divergences of the diagonal DGLAP kernels by separating it from the hard-collinear behaviour that is identical to that in the form factor. We shall now verify the above result, eq. The derivation of eq. (3.27) consists of two parts: firstly we will compute the bare diagrams and the UV counterterms related to the renormalisation of the QCD coupling constant, then we will subtract the short-distance singularities associated with the Wilsonline operators, thus completing the renormalisation of log W ⊓ . The non-Abelian exponentiation theorem [110][111][112][113] allows us to determine directly log W ⊓ by computing only the webs that capture the maximally non-Abelian colour factors of each Feynman diagram, as defined in [113]. Moreover, log W ⊓ has a simpler singularity structure compared to W ⊓ , which allows us to setup the renormalisation procedure directly at the level of the webs.
We introduce the following parameterisation for the contour of the Wilson loop We use the following Feynman rules in configuration space for the gluon propagator in Feynman gauge and for the gluon emission from the eikonal lines, respectively 2) where T a is the SU(N ) generator in the appropriate representation and N = − Γ(1−ǫ) 4π 2−ǫ . In section 4.1 we consider the one-loop calculation of log (W ⊓ ) and then establish its general form before and after renormalisation. In section 4.3 we perform the calculation at two loops, verifying the general structure and obtaining an explicit result for Γ ⊓ consistent with eq. (3.36).

One-loop calculation
As a direct consequence of the Feynman rules given above, all the diagrams that feature a gluon exchange between two lines with the same lightlike velocity v are proportional to v 2 and therefore they are automatically zero. At one loop there will be only two non-vanishing webs contributing to log (W ⊓ ) which differ only by a translation and therefore yield the same result where C i , with i = A, F is the quadratic Casimir in the adjoint or in the fundamental representation. We notice that the integral over the parameter t diverges both in the UV limit t → 0 and in the IR regime t → ∞. This fact is a consequence of the absence of any scale associated with the integration over an infinite Wilson line and it implies that the bare diagram in eq. (4.4) yields a vanishing contribution. Nevertheless, the diagram is non-trivial after the renormalisation procedure, which subtracts the divergence for t → 0 and allows us to define uniquely the integrand in eq. (4.4). In order to expose the analytic structure of eq. (4.4) in terms of the variable ρ defined in eq. (3.26), we rotate the path along the negative imaginary axis in the complex t−plane. Then we change variables t = −i √ 2λ, The complete result for log (W ⊓ ) at one loop is given by twice the contribution of eq. (4.5).
It is convenient to write it with the factor (4πe γ E ) ǫ absorbed into the MS running coupling as follows The label "bare" reminds us that eq. (4.6) still has the UV divergences associated to the cusps of the Wilson loop in eq. (4.1), which must be subtracted before IR singularities can be identified. Indeed, it is convenient to show explicitly that eq. (4.6) is independent on the renormalisation scale, by writing the running coupling as which leads to the expression

Exponentiation and renormalisation
The integrand in eq. (4.8) is finite in the limit ǫ → 0 and the singularities of log W ⊓ arise only after the integration over λ, σ. In particular, following the coordinate-space analysis of refs. [72,114,115], we distinguish three possible types of singular behaviour: cusp singularities, which are associated to the limit λ ≃ σ → 0 in which all the vertices approach a cusp of the Wilson loop; collinear singularities, which arise if either λ or σ approaches the cusp, while the other parameter stays finite; finally, the large-distance region with λ → ∞, which determines the IR pole. At higher perturbative orders, individual diagrams feature soft and collinear subdivergences when a subset of the vertices approaches one of these limits, which give rise to poles of higher order compared to those in eq. (4.8). However, owing to its exponentiation property, upon considering the logarithm of the Wilson-line correlator, all the subdivergences cancel in the sum of webs at each perturbative order [72,108,[115][116][117]. It is always possible to organise the calculation of log (W ⊓ ) such that the integral over the position of the vertex that is located at the largest distance along the infinite Wilson line is performed last. Thus, the single infrared pole will be generated only in the final integration, while all the subdivergences of individual diagrams cancel in the sum of webs. This procedure, which follows the prescriptions of ref. [72], allows us to generalise the representation of eq. (4.8) to all orders where the integrand w has an expansion in ǫ that involves only non-negative powers The representation eq. (4.9) is analogous to the one derived in [72] for the soft function of the form factor, defined in eq. (2.13), with the difference that in the latter case the integrals over both the parameters are unbounded. This is consistent with the presence of a double pole of long-distance origin in the form factor, as compared to the single pole of this type arising in eq. (3.27). We now proceed with the renormalisation of the singularities of short-distance origin that are present in the bare expression of eq. (4.9). Following [72], we notice that the integral of w 0 in eq. (4.9) generates double UV poles, which are subtracted by cutting the integration domain with λ < 1 µ , σ < 1 µ in eq. (4.9), where µ defines the subtraction point. The contributions of w i with i ≥ 1 generate at most one UV singularity, which we subtract in the last integration. In conclusion we derive the representation for the sum of renormalised webs in configuration space where we performed the integral over σ by expanding the coupling constant α s 1 λσ at the scale 1 λ 2 , as in eq. (4.7) Eq. (4.11) directly leads to the result eq. (3.27) from the web integrals in coordinate space and it allows us to extract the coefficients γ cusp and Γ ⊓ . At one-loop order, we expand the web in eq. (4.8) and we get Applying the renormalisation procedure described above we find where we have used the fact that W ⊓ consists of pure poles. The pole is infrared and is exactly the one that replicates the soft divergence of the PDF. We compare eq. (4.14) with the poles of eq. (3.27) getting (4.15)

Two-loop calculation
We now apply the renormalisation procedure to the two-loop webs. Only a few diagrams contribute to this order and they are represented below 3s = , (4.16) where we omit the configurations that are simply obtained by mirror symmetry. The diagrams in the first row of eq. (4.16) are computed following the same steps as the oneloop case. As in the one-loop case, we write the bare webs using the representation in eq. (4.9) and we define the integrand w From now on we drop the arguments on d i and w i which are understood to have the above arguments unless otherwise stated. The first diagram, d SE is obtained from eq. (4.4) by replacing the gluon propagator eq. (4.2) with the one-loop expression where we discarded the longitudinal components of the propagator, that are proportional to ∂ µ ∂ ν , because they decouple from the amplitude via Ward identities [72]. The result is (4.19) in agreement with the results of [62,72]. We notice immediately that, at two-loop level, the representation eq. (4.17) of the individual webs has subdivergences, which are manifest as explicit poles in the integrand w i . In this case the subdivergence is cancelled by the coupling renormalisation in the QCD Lagrangian that will be taken into account later in this section. The double gluon exchange diagrams give w (2) Both results are in agreement with the maximally non-Abelian contributions of the diagrams W c and W d reported in [62]. The integrand of the diagram d These subdivergences are not related to QCD renormalisation and they will cancel in the sum of all webs. Eq. (4.21) is finite when ǫ → 0 as we discuss more in detail in appendix B.2. The diagrams in the second row of eq. (4.16) involve the three-gluon vertex, whose Feynman rule in configuration space reads We notice that the diagrams d Ys and d Y L are not related by symmetry transformations, because the former has two gluon attachments on the segment of finite length y, while the latter has two emissions from the semi-infinite line. We begin with the calculation of d where we introduced the normalisation factor K Y = ig 4 We write the differential operators in eq. (4.23) in terms of total derivatives as follows which allows us to perform immediately the integrals over s 2 and over s 1 , respectively in the first and in the second term in curly brackets, by evaluating the appropriate propagator at the endpoints of the integration interval. Eq. (4.23) becomes in terms of the functions where the prescription +i0 is understood in every factor appearing in the integrals. Each function has a clear diagrammatic interpretation, because the integrands are products of scalar propagators in coordinate space. Thus, d Ys is decomposed in a sum of diagrams, as discussed in [72], giving in one-to-one correspondence with the three terms in eq.
Let us discuss the singularity structure of the separate integrals above. The only one which is separately finite is eq. (4.29a), which corresponds to the integrand of the first diagram in eq. (4.28). Eq. (4.29b) has single and double poles that will cancel the corresponding singularities in eq. (4.20). Indeed, the second diagram in eq. (4.28), which is associated to the integrand in eq. (4.29b), has subdivergences of short distance origin when the threegluon vertex approaches the cusp, similarly to the behaviour shown by diagram d (2) X 2 . The single pole in eq. (4.29c) is entirely due to the presence of a one-particle-irreducible UV divergent subgraph in the last diagram in eq. (4.28). Therefore, this singularity is removed by the counterterms of the QCD Lagrangian. Using these results, the total contribution of the diagram d (2) Ys in eq. (4.25) agrees with the corresponding expression for diagram W e in [62] and in the notation of eq. (4.17) it reads w (2) The next diagram, d Y L , differs from eq. (4.23) only by the presence of the two gluon attachments on the semi-infinite Wilson line rather than on the finite one. Once again, we write the three gluon vertex in terms of total derivative and we decompose the diagram as Y L , the Wilson line is infinite and this term is absent. This result was shown in [72], by introducing a cutoff on the infinite line and carefully taking the limit to infinity, which does not commute with the integration over z. The same conclusion is found by computing d (2) Y L in momentum space, as shown in appendix B.1. Using eqs. (4.29a), (4.29b) and (4.29c) we get By construction, the expression above has the same singularities as w (2) Ys , because the integrand of the diagram d (2) Y L differs from d (2) Ys only by the function in eq. (4.29a), which is finite.
We compute the diagram d (2) 3s using the same procedure which is finite because it involves only the function in eq. (4.29a). We renormalise the UV divergences associated with the QCD vertices and propagators by means of the one-loop counterterm where d (1) is the result of the one-loop diagram eq. (4.5). Finally, we sum all the diagrams depicted in eq. (4.16), including the symmetric configurations which are not shown there, getting 36) where to get to the second line we used the identity 2d (2) Ys = 2d (2) Y L + d (2) 3s that is obtained by comparing eqs. (4.28), (4.31) and eq. (4.33). The terms in curly brackets in the final expression are the same that appear in the calculation of the cusped Wilson loop with two semi-infinite lightlike lines, discussed in [72]. The last two contributions in eq. (4.36) are special to the configuration of W ⊓ , where the semi-infinite lines are connected by a finite lightlike segment. The final expression in eq. (4.36) follows the decomposition of polygonshaped Wilson loops presented in [72]. The distinction between the terms inside and outside the curly brackets in eq. (4.36) stems from the structure of their singularities. The former ones give rise to cusp configurations characterised by double UV poles and therefore they can be written in terms of the representation in eq. (4.9) with a finite integrand. The last contributions in eq (4.36) generate at most a single pole, associated to the configurations where all the vertices simultaneously approach a lightlike segment, and therefore their combination will give rise to an integrand of order ǫ in eq. (4.9), as we verify by expanding eqs. (4.21) and (4.34) (4.38) By expanding the equation above in ǫ we get where Γ (2) ⊓ is the two-loop contribution to Γ ⊓ eq. (3.36), Now we renormalise using the procedure outlined in the one-loop case, see eq. (4.11). For the terms of O(ǫ 0 ), the cusp terms, we integrate from 1/µ on both integrals. For all the subsequent terms, the σ integral is performed first, integrating from 0 to ρ √ 2 . Then the parameter λ integrated from 1 µ . By doing this we get, Again it is reminded that we have used the fact that W ⊓ consists of pure poles. The above is reproduced by eq. (3.27) By this point we have determined the anomalous dimension Γ ⊓ in two different ways, first by extracting it from the evolution of PDFs using the universality of the hard-collinear poles J/J and now by a direct computation of the renormalisation of the corresponding Wilson-line correlator.

Relating Wilson-line geometries to physical quantities
In this section we establish a set of relations between different physical quantities, based on the properties of the Wilson loops discussed in section 4. In section 5.1 we will show that the single infrared poles in the quark and in the gluon form factors are related to the corresponding diagonal term in the DGLAP kernels by a precise eikonal quantity that is associated to the geometry of the Wilson loops with lightlike lines. The latter emerges as the difference between the anomalous dimensions associated with a wedged Wilson loop with two semi-infinite lines and a ⊓−shaped Wilson loop. This difference, in turn, can be expressed as the anomalous dimension associated with a parallelogram (or more generally) a polygon with lightlike segments. In section 5.2 we use this relation to extract the anomalous dimensions associated with a polygon Wilson loop to three loops, which is related to the soft anomalous dimension appearing in the resummation of threshold logarithms in the Drell-Yan process. Finally we extract the fermionic components of the four-loop result in the planar limit.

Relating the form factor with the DGLAP kernels
The direct calculation of the anomalous dimension Γ ⊓ in section 4 confirms the identity in eq. (3.34) which follows from the factorisation of the parton distribution functions for large x. This identity is interpreted as a decomposition of B δ , which was defined in eq. (3.33) as the coefficient of the delta distribution in the splitting functions in the limit x → 1, into the contribution of the hard-collinear radiation, γ J/J , and the purely soft one, which is encoded by Γ ⊓ . In eq. (5.1) we suppressed the dependence on the external parton: the relation holds for both quarks and gluons. The hard-collinear contribution γ J/J is process independent, as discussed in section 2.2 in the context of the infrared factorisation of the form factor. Indeed, eq. (2.22) provides the analog of eq. (5.1) where γ G is the anomalous dimension that determines the single poles of the form factor. By comparing eq. (5.1) and eq. (5.2) we derive the relation which connects the single poles in the form factor with the diagonal DGLAP kernels. The two quantities appearing on the left-hand side of eq. (5.3) depend on both the spin and the colour representation of the external particles in a non-trivial way. In contrast, the righthand side involves the anomalous dimensions of two eikonal quantities, which depend only on the colour representation of the particles and obey Casimir scaling up to three loops. Therefore, eq. (5.3) allows us to interpret the function f eik eq. (1.1), which was introduced in [44,45] as the difference f eik ≡ γ G − 2B δ , in terms of the anomalous dimensions of Wilson-line correlators. By substituting the two-loop expressions of Γ ⊓ and Γ ∧ from direct calculations, respectively in eqs. (3.36) and (2.24), into the right-hand side of eq. (5.3) we reproduced the two-loop result obtained from the difference of γ G and B δ in ref. [45], namely thus verifying eq. (5.3) through two loops. The difference of anomalous dimensions appearing on the right-hand side of eq. (5.3) has also a geometric interpretation, which suggests to define it as universal quantity. Following the analysis of the singularities of the Wilson loops with lightlike lines detailed in ref. [72] and the calculation in section 4 above, the anomalous dimensions Γ ⊓ and Γ ∧ receive contributions only from the singular configurations, in which all the vertices approach one lightlike line. In this sense, these anomalous dimensions depend only on the features of each lightlike line separately and they are insensitive to the global shape of the Wilson loop. Both Γ ⊓ and Γ ∧ encode the collinear singularities associated with the two semi-infinite lightlike lines, but the former receives an additional contribution from the configurations that are collinear to the finite segment. Such singularities differ from the ones originating from infinite lines by the presence of endpoint contributions, as we showed by computing the diagrams d (2) Ys and d (2) Y L in eqs. (4.28) and (4.31). It is therefore useful to define the difference of Γ ⊓ and Γ ∧ as the anomalous dimension that captures the collinear singularities of a finite lightlike segment. Similarly, we define also the collinear anomalous dimension associated to infinite lines in terms of Γ ∧ only The two-loop expression of Γ fin co coincides with the right-hand side of eq. (5.4), while Γ inf co to the same order is obtained from eq. (2.24). Comparing the two expressions we get The factor of two multiplying Γ inf co is consistent with the fact that the finite Wilson line is obtained as a contour involving two semi-infinite lines. The remaining discrepancy proportional to ζ 3 is related to the endpoint contributions in eq. (4.37).
The geometric interpretation of Γ fin co and Γ inf co allows one to derive the anomalous dimensions of Wilson loops with the contour consisting of arbitrary, possibly open, polygons with lightlike lines. The first example is the parallelogram-shaped Wilson loop W that features four lightlike segments of finite length (see figure 1c), whose renormalisation was given in [62] where x and y are the four-vectors that define the sides of the parallelogram. Γ receives contributions from the collinear divergences of four finite segments in lightlike directions, therefore it is By replacing in the equation above the two-loop value of Γ fin co in eq. (5.4), we reproduce the results of Γ in ref. [62]. In the case of a generic polygonal Wilson loop W i with lightlike lines the evolution equation in eq. (5.8) generalises [62,72] where the sum is extended over all the cusps in the contour and x a , x a−1 define the sides adjacent to the cusp a. The anomalous dimension Γ i collects all the collinear contributions and it can be derived by summing the appropriate multiples of Γ fin co and Γ inf co , given by the number of finite and infinite sides, respectively.
Finally, having identified the difference Γ ⊓ − Γ ∧ = Γ 4 in eq. (5.9), we may notice that eq. (5.3) provides yet another identity relating the form factor, the DGLAP kernel and the Wilson loop W computed in ref. [62], namely thus explaining the numerical agreement of these two quantities computed respectively in ref. [45] and in ref. [62].

The Drell-Yan soft function and Γ beyond two loops
We move to relate the abstract W to a physical quantity relevant for soft-gluon resummation. It is known that the Drell-Yan cross-section factorises near threshold [1][2][3][4][5] (see also the more recent literature in the Soft Collinear Effective Theory [6,17]). The hard-collinear region is described by the PDFs, the hard function by a squared timelike form factor and the soft region by Wilson lines in the DY configuration [52]. This leads to the all-order relation γ G − 2B δ = Γ DY /2, where Γ DY is the anomalous dimension associated to the DY configuration of Wilson lines (see e.g. [5,6]). Using eq. (5.11) we have, The ideas in section 5.1 will allow us to test the identification in eq. (5.12). The three-loop value for γ G − 2B δ was first extracted in [45] using the three-loop results for γ G and B δ . If we expand Γ as, using the values in [45] we can then state Γ (1) = 0 (5.14) As mentioned in section 5.1, the two-loop Γ (2) was calculated explicitly using Wilson lines in [62] and agrees with the extracted value in eq. (5.15). The three-loop Γ (3) displayed in eq. (5.16) should be regarded as a prediction to be verified by direct calculation. For the Drell-Yan configuration of Wilson lines, Γ DY was computed at two loops in [52] and three loops in [20]. The three-loop Γ DY coincides with eq. (5.16). This is a non-trivial three-loop test of the identification in eq. (5.12), we arrive at the same value for Γ At four loops the complete picture in QCD for Γ is unknown but in planar N = 4 super Yang-Mills the four-loop result for the difference γ G − 2B δ was found in [39]. We identify this as Γ planar N =4 and quote the result here, It is well-known that to reach the above result one can simply take the QCD result and take the limit N c → ∞ and the maximal transcendental weight term at each order in α s . We can do this at two and three loops by looking at eqs. (5.4) and (5.16) respectively.
Above three loops, strict Casimir scaling has been proven to fail [27,28,35]. As such, we need to distinguish between quarks and gluons or rather particles in the fundamental and adjoint representation. We focus on the quark case. To compute Γ (4),q we need B q δ and γ q G at four loops. The state of the art is that some colour structures are known for both B q δ [25,27] and γ q G [26] in the planar limit, N c → ∞. Using the values in [26,27] we can extract the following terms in that limit for Γ (4),q , where we have used T f = 1 2 . We are unable to deduce the N 4 c n 0 f term as it is unknown for γ G but it is known for B δ [27]. In the planar limit N c → ∞ the quartic Casimirs d (4) F F ≡ d abcd F 2 contribute to the colour factor in eq. (5.18) since, It means that we are unable to fully construct the Casimir scaling of N 3 c n f . The full (planar and non-planar) contribution of the quartic Casimir colour factor d (4) F F to γ G is known [30] but not to B δ . Only the low-N values of the splitting functions or γ cusp is known [29,32]. In [29] it was also found that, within numerical errors, quartic Casimir contribution to the cusp anomalous dimension did not depend on the representation, i.e. it is the same for gluons and quarks. It was conjectured that although Casimir scaling is violated there is a generalised version where quartic factors are simply exchanged depending on gluon or quarks, where N F/A denote the dimensions of the corresponding representations, namely N F = N c = C A and N A = N 2 c − 1 = 2N c C F . The relation in eq. (5.3) may be used as an interesting test for a generalised Casimir scaling extension to the anomalous dimension Γ .
However, the quartic Casimirs do not appear in the n 2 f or n 3 f terms of eqs. (5.19) and (5.20). We are then able to use these terms for a leading-N c , Casimir-scaling part of Γ (4) . We put these terms together with the conjectured generalised scaling to create an ansatz for Γ (4) , where the ellipsis represents all terms subleading in N c , including the n 1 f and n 0 f terms, which are not found from the quartic Casimirs when they are expanded in N c .

Conclusions
We have presented a detailed study of the infrared factorisation of form factors and PDFs at large x using a common formalism. By identifying the universal contributions from the hard-collinear region in both quantities, those controlled by the anomalous dimension γ J/J , we were able to derive the relation in eq. (5.3), That is, the difference between anomalous dimension describing single poles in the on-shell form factor of quarks (gluons) and that associated with the δ(1 − x) term in the largex limit of the quark (gluon) diagonal DGLAP splitting function, reduces to a difference between corresponding eikonal quantities, Γ ∧ and Γ ⊓ defined directly in terms of Wilson loops. Furthermore, based on the configuration-space origin of the contributions to these two eikonal quantities we concluded that their difference simply corresponds to the anomalous dimension associated with a closed polygonal Wilson loop, such as the parallelogram analysed first in ref. [62]. The contributions of the semi-infinite Wilson lines in W ⊓ and W ∧ cancel in the difference. We emphasise that while each of the quantities on the left-hand side of eq. (6.1) depends in a non-trivial way on the spin of the partons, in addition to their colour representations, yielding very different results for quarks and for gluons, the eikonal quantities, by definition, depend only on the colour representation of these partons, and in particular admit Casimir scaling through three loops. We stress that the relation in eq. (6.1) is expected to hold to all orders in perturbation theory. An obvious next step is to compute Γ to three loops in order to check it explicitly to this order. In establishing the relation between Γ and Γ ⊓ − Γ ∧ we used the fact that singularities arise only from configurations where all the vertices approach a cusp or one where they all approach a particular lightlike segment [72]. This underlies the cancellation of the two infinite segments, isolating a remaining finite segment. The very same logic may be applied to other, more complicated Wilson-line geometries involving both finite and semi-infinite lightlike segments. Specifically, the double pole is always governed by γ cusp while the singlepole anomalous dimension is written as a sum of terms, building blocks, each corresponding to either a finite or semi-infinite segment, which contribute Γ fin co and Γ inf co , respectively. An example of such a construction with only finite segments can be found in refs. [118][119][120], where polygons of up to six sides were computed to two loops. Following our discussion in section 5.1 it may be interesting to explicitly compute other Wilson-line configurations involving both finite and infinite segments. A simple example of direct relevance to physics is the non-forward amplitude, generalising the ⊓ configuration.
One interesting aspect that we have encountered is that W ⊓ behaves very differently in the ultraviolet as compared to the infrared, as can be seen explicitly in eq. (4.42). In the ultraviolet, one encounters a double logarithmic dependence on the scale µ 2 , originating from the cusp singularity, while in the infrared there is just a single pole. This stands in sharp contrast to the W ∧ , corresponding to the soft function of the form factor (2.16) (or more generally, in soft functions corresponding to multi-leg amplitudes) where the infrared behaviour entails a double pole, mirroring the ultraviolet. The absence of any distance scale in the relevant Wilson-line contour implies such mirroring. Indeed the symmetry between the ultraviolet and the infrared is broken in W ⊓ due to the presence of the scale β · y. The single-pole character of W ⊓ can be seen as intermediate in comparing W , which, lacking infinite rays, is infrared finite, to W ∧ , which is double logarithmic.
The relation in eq. (5.12) between the soft anomalous dimension in Drell-Yan production and the parallelogram W is interesting in its own right. The Drell-Yan soft function involves real gluon emission diagrams where the propagators connecting the amplitude side to the complex-conjugate one are cut, while in W there are no cut propagators. A possible way to explain 6 this is to recall that a parallelogram made of four lightlike segments features two cusps where the exchanged gluons span timelike distances and two others where gluons span spacelike distances. The latter correspond to diagrams that feature in virtual corrections to the Drell-Yan process (these propagators are not cut). In turn, the former are naturally time-ordered, because there path-ordering coincides with time-ordering (just as in the case of the W ⊓ , discussed below eq. (3.25)) and could be computed using either cut propagators or ordinary ones, giving the same answer. This way the calculation of the parallelogram can be mapped into that of the Drell-Yan soft function. It would be interesting to turn this argument into a proof relating the two Wilson line configurations directly. It would also be interesting to explore in this context the conformal mapping techniques of refs. [50,51].
Another interesting direction to explore is the connection between partonic amplitudes in the Regge limit and anomalous dimensions of Wilson lines. In particular, one would like to derive the relation between the Regge trajectory and W ∧ in eq. (1.6) and understand its generalisation to higher orders.

A Direct calculation of the splitting functions at large x
In this appendix we present a calculation for parton distribution splitting functions directly using the definitions (3.1) and (3.2). As explained in the main text we take incoming partons to be off shell p 2 = 0 but with zero transverse momentum p = (p + , p 2 2p + , 0 d−2 ). This regulates the infrared such that we are only exposed to UV poles.
To calculate a single diagram there is a general strategy talked through in the main text with a slight change for the off shell case: • Write down the integral using Feynman rules • Integrate over the transverse component of all loop momenta. As we are only interested in the UV divergent terms, this can be simplified to just calculating iterated bubbles at two loops. Although we do need to calculate the finite terms of one loop graphs to perform the renormalisation.
• Rescale the plus component to arrive at a general form such as, The denominators correspond to the Wilson line propagators.
• Now we take the discontinuity in x and perform the final integrations. Often these integrals evaluate to 2 F 1 functions at two loops.
• Finally we expand in ǫ using, For brevity of results we shall define, Also, all the following expressions are valid up to but not including terms that diverge as slowly as log(1 − x). For PDFs we expand in powers of αs 4π , Example. Let us illustrate the above steps. For this we choose the two loop diagram in figure 4e. The Feynman rules for the diagram, in Feynman gauge, give, where k · u = xp + and the +iε prescription is implied. It is reminded that the Wilson line direction u is in the (-) direction, u = (0, 1, 0 d−2 ).
We shall define f (2),(e) qq to be the contribution to f (2) qq of diagram 4e. When integrating the q 1− and q 2− components using the residue theorem, constraints are placed on the plus momenta, The integrations over the transverse components are just iterated bubbles. Rescaling the plus component we arrive at, .
Taking the discontinuity we use, .
There are three separate terms now to calculate. The first is the virtual cut and evaluates to, The second and third are real cuts, In the sum we find, Above we see a salient feature of two loop diagrams: individual cuts are ǫ −4 and ǫ −3 divergent. These are poles from when the emitted gluon goes soft and cancel in the sum of real and virtual cuts. The remaining divergences are UV, whose renormalisation gives the splitting functions. The L and L 2 terms are present in individual diagrams but cancel in the combination such that the splitting functions diverge as in eq. (3.33).
Another feature is that we found that the real cuts, eqs. (A.10) and (A.11), contribute to B δ . Rather than inferring the coefficient of δ from sum or momentum conservation rules, we are able to state that for the off-shell extraction of the splitting functions, real cuts contribute to δ(1 − x).
We now calculate the two loop diagonal splitting functions at large x for quarks and gluons.

A.1 Calculating P qq
The one loop contributions are figures 4a and 4b and the self energy on each external leg. They sum to, where ξ is the gauge parameter in a general covariant gauge. The two loop contributions are shown in figures 4c-4m. They exclude self energies on external legs. Their calculation was performed in Feynman gauge ξ = 1 and the results are,

Figure 4:
Large-x divergent contributions to the quark-quark parton distribution up to two loops. The grey blob represents a self energy insertion. Each diagram has a multiple factor displayed. Insertions on external legs are excluded.
Summing the two loop contributions with the factors shown in figure 4 we find, At two loops we need to take into account the running from the one loop contribution, . This is found by replacing ξ → 1 + αs 4πǫ 10 6 C A − 4 3 T f n f ξ and α s → 1 + αsb 0 πǫ α s . We then specialise to Feynman gauge ξ = 1. We then find the Z qq that minimally subtracts the divergences in δ + αs 4π f qq . As the renormalisation is multiplicative, convolutions need to be taken into account for one loop squared terms. For example, Equivalently the renormalisation can be transformed to Mellin space, eq. 3.20, where the convolutions become products ensuring that, is finite in ǫ. We can then extract the splitting functions to two loops from, where Z q is the wavefunction renormalisation in MS for the quark. Up to two loops, Converting back to x space we find, Notice that we find that all L n terms cancel. This reproduces B q δ , the coefficient of δ in eq. (3.35), and shows that the coefficient of P is γ cusp as in eq. (2.4).

A.2 Calculating P gg
The one loop contributions for the gluon gluon distribution function are shown in figures 5a and 5b. The total one loop contributions are, The two loop contributions are shown in figures 5c-5p. The two loop contributions are, f (2),(j) The total two loop contribution is, The extraction of the splitting function from above is the same as in the quark case. Instead of Z q we use the gluon field renormalisation in MS, Performing those steps we find, Again this aligns with B g δ in eq. (3.35) and γ cusp in eq. (2.4). We have replicated previous splitting function calculations at large x directly from the definitions (3.1) and (3.2) in a covariant gauge. By taking the incoming partons off shell, p 2 = 0, we regulate the infrared divergences allowing the extraction of the UV poles of the PDFs. Although the divergent terms remain gauge independent the finite terms become gauge dependent. It means that we need to take into account the running of the gauge parameter ξ → Z A ξ in finite terms, even when working in Feynman gauge.

B Particular two-loop diagrams contributing to W ⊓
In this appendix we elaborate on aspects of the calculation of W ⊓ presented in section 4. We consider two specific diagrams where some subtle points arise. In section B.1 we discuss the endpoint contributions in diagram d  In section 4 we revisited the analysis of non-Abelian contributions to the correlators of finite and semi-infinite Wilson lines [72]. Specifically, we derived the representations of the twoloop diagrams that contain a three-gluon vertex and made a clear distinction between ones where two gluons are emitted from a finite Wilson-line segment as compared to the case where two emissions emerge from a semi-infinite line, corresponding respectively to diagrams d (2) Ys and d (2) Y L in (4.16). The difference is that in the former case both endpoint contributions appear, as in (4.28), while in the latter case there is no endpoint contribution from infinity, so the representation of d (2) Y L simplifies to (4.31). Let us now present this calculation in detail using momentum space and show explicitly that this endpoint contribution is indeed absent.
Using the Feynman rules given in section 4, diagram d which is analogous to eq. (4.23). In the equation above, we integrate over z using the momentum-space representation of the propagators After taking the derivatives with respect to s 1 β, s 2 β and integrating over the infinite line we get where the prescription +i0 in the denominators ensures the convergence of the integrals for s 1 → −∞. The expression above may be conveniently rewritten as This directly leads to the representation of eq. (4.31), as we now show. Upon introducing an auxiliary integration constrained by momentum conservation we obtain: The representation of the delta function (2π) d δ d (k 1 + k 2 + k 3 ) = d d z e i(k 1 +k 2 +k 3 )·z is interpreted as an integral over the position of the scalar "three gluon" vertex in eq. (4.31). Using eq. (B.2) we recover the expression of the three gluon propagators in coordinate space, carrying momenta k 1 , k 2 and k 3 , obtaining Substituting the definitions in eqs. (4.29a), (4.29b) and (4.29c) we verify the result in eq. (4.31).

B.2 The diagram d (2)
X3 connecting three Wilson lines In this section we derive the representation of eq. (4.21) of the diagram d (2) X 3 that connects two cusps with a lightlike segment of finite length. Following the discussion of ref. [72], the singularities of the webs of this kind are associated with the configuration where all the vertices approach the lightlike segment of finite length. These webs do not contribute to the cusp singularities because there is not any region of configuration space where all the vertices are in proximity of the cusp. By using the Feynman rules in eq. (4.2), the diagram d where the factor − C A C F 2 corresponds to the maximally non-Abelian part of the colour factor of the diagram, which is exponentiated [110][111][112][113]. We expose the overall infrared singularity in the last integration by rewriting the integration domain using θ(t 1 − t 3 ) + θ(t 3 − t 1 ) = 1 and changing the order of integrations. Thus we obtain d (2) We stress that the expression above still has infrared singularities from the limit t → ∞ in the upper bound of the t ′ integral. Therefore we decouple the infrared contributions by applying the changes of variables The parameters a 1 and a 2 are integrated immediately, leading to d (2) Thus we apply the change of variables introduced before eq. (4.5) and we get d (2)