Cwebs beyond three loops in multiparton amplitudes

Correlators of Wilson-line operators in non-abelian gauge theories are known to exponentiate, and their logarithms can be organised in terms of collections of Feynman diagrams called webs. In [1] we introduced the concept of Cweb, or correlator web, which is a set of skeleton diagrams built with connected gluon correlators, and we computed the mixing matrices for all Cwebs connecting four or five Wilson lines at four loops. Here we complete the evaluation of four-loop mixing matrices, presenting the results for all Cwebs connecting two and three Wilson lines. We observe that the conjuctured column sum rule is obeyed by all the mixing matrices that appear at four-loops. We also show how low-dimensional mixing matrices can be uniquely determined from their known combinatorial properties, and provide some all-order results for selected classes of mixing matrices. Our results complete the required colour building blocks for the calculation of the soft anomalous dimension matrix at four-loop order.


Introduction
The infrared (IR) structure of scattering amplitudes in gauge field theories is universal. It is this universality -the independence from the specific nature of the hard process under consideration, that makes it an important object of study. Remarkable insights on how IR singularities organise themselves in perturbation theory have been accumulated over a long history spanning almost a century [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. While these studies are of course interesting in their own right, uncovering subtle properties of gauge field theories, it is important to note that they also have practical applications to high-energy scattering. First of all, even though IR singularities cancel when one constructs an experimentally measurable IR-safe observable, they often leave behind potentially large logarithms of kinematic variables that jeopardise the applicability of fixed-order perturbation theory. The universal IR structure of scattering amplitudes is one of the key ingredients allowing for the all-order summation of these large logarithms, leading to more precise and controlled physical predictions [19][20][21]. Furthermore, even at fixed orders, implementing the cancellation of IR singularities for complex collider observables at high perturbative orders is a difficult and important undertaking: indeed, the construction of general and efficient subtraction procedures beyond next-to-leading order (NLO) is a broad ongoing effort (see, for example, [22][23][24][25][26][27][28][29][30][31][32]).
To be somewhat more precise, gauge theory scattering amplitudes in the IR limit factorise into universal soft and collinear functions [15][16][17][18][19]33], multiplying finite matching coefficients. These soft and collinear functions, in turn, can be expressed as matrix elements of field operators and Wilson lines, which are the object of the present study. We note in passing that such matrix elements play a ubiquitous role, not only for the factorisation of scattering amplitudes, but also in many effective theories based on QCD [34][35][36]. The methods developed in [1] and the objects of interest -Cwebs -concern the evaluation of Wilson-line correlators of the general form S n (γ i ) ≡ 0| n k=1 Φ (γ k ) |0 , (1.1) where Φ(γ) are Wilson-line operators evaluated on smooth space-time contours γ, where A µ (x) = A µ a (x) T a is a non-abelian gauge field, and T a is a generator of the gauge algebra, which can be taken to belong to any desired representation. If the smooth contours γ are closed, then the correlator is gauge-invariant, while, if they are open, the correlator is a colour tensor with open colour indices in the chosen representations, attached to the endpoints of each Wilson line. In the present paper, we will restrict ourselves to soft colour operators associated with multi-particle scattering amplitudes in gauge theories, which encode all their soft singularities. These operators are of the form of Eq. (1.1), with the contours γ k given by rays extending from the origin along directions β k , corresponding to the four-velocities of the particles participating in the scattering. In this case, we write the soft operator as where each Wilson line is taken in the representation of the gauge algebra corresponding to the high-energy particle it represents. We note that soft operators of the form of Eq. (1.3) are affected by ultraviolet, soft, and, for β 2 i = 0, collinear singularities: as a consequence, special care is required to evaluate them [37][38][39][40]. In what follows, we will generally assume β 2 i = 0, so that collinear divergences will not arise; for soft divergences, we will make use of the smooth exponential suppression of gluon interactions at large distances, discussed in Refs. [39,40]; finally, we will retain dimensional regularisation for ultraviolet singularities. With this scheme, the bare correlator is finite, and can be evaluated and renormalised. The physically relevant quantity to be extracted is then the set of ultraviolet counterterms for the correlator: indeed, in the absence of an IR regulator, all radiative corrections to Eq. (1.3) vanish, since all relevant Feynman diagrams are scale-less; the renormalised correlator is therefore given precisely (in a minimal scheme) by the set of its UV counterterms.
Renormalised correlators of the form of Eq. (1.3) obey renormalisation group equations which lead to exact exponentiation. This means that they can be written as S n β i · β j , α s (µ 2 ), = P exp − 1 2 where Γ n is the soft anomalous dimension matrix, which is the focus of our studies. Γ n was computed at one loop in [41] (see also [42]); at two loops in the massless case in [43,44], and in the massive case in [45][46][47][48][49]; finally, at three loops in the massless case in [50,51]. Pushing these calculations to higher perturbative orders is of great interest, since the resulting structures are highly constrained by symmetries, and display non-trivial mathematical properties, tied to interesting physical configurations: preliminary results on the four-loop structure of Γ n have been given in [52,53].
In practice, an approach alternative to the above renormalisation-group based methods is often useful for studying general Wilson-line correlators of the form of Eq. (1.1). All such correlators are known to obey a non-trivial form of diagrammatic exponentiation, so that one can write S n (γ i ) = exp W n (γ i ) , (1.5) where the exponent, W n (γ i ), can be directly computed in terms of a subset of the Feynman diagrams contributing to S n (γ i ). For non-abelian gauge theories, this was first pointed out in Refs. [54][55][56], in the case of two straight, semi-infinite Wilson lines. For general configurations, it was proven in Refs. [57,58]. These studies show that Feynman diagrams contributing to W n (γ i ) are organised in sets that are called webs: in particular, • For an abelian theory, webs are given by connected photon diagrams attaching to the Wilson lines.
• For a non-abelian theory, if only two Wilson lines are present, webs are given by two-eikonal irreducible diagrams, i.e. diagrams that do not become disconnected upon cutting only the two Wilson lines.
• For general, multi-line correlators, webs are sets of diagrams that differ among themselves by the ordering of their gluon attachments to the Wilson lines.
Clearly, by means of webs, one can directly compute the soft anomalous dimension matrix Γ n . It is important, however, to note that, for general multi-line correlators, multiplicative renormalisability and exponentiation combine non-trivially, due to the non abelian nature of webs, and the calculation of renormalised correlators includes commutators of counterterms and bare webs, as required [59]. Not surprisingly, non-trivial simplifications of the final result for Γ n appear only when the commutator terms are properly taken into account. We also recall the special properties arising in the massless case, when the countours γ i are light-like, straight, semi-infinite Wilson lines: for such amplitudes, scale invariance imposes strong constraints on the functional dependence of the soft anomalous dimension matrix. Up to two loops, Γ n can only involve dipole colour correlations between Wilson lines [16,17,33,43,44,60]; beyond two loops, scale-invariant conformal cross ratios of the form ρ ijkl ≡ (β i · β j β k · β l )/(β i · β k β j · β l ) correlating four partons can also appear: the first such correlations arise at three loops, with at least four Wilson lines, and were computed in Ref. [50]. At four loops, one finds the first occurence of higher-order Casimir operators of the gauge algebra, which have recently been computed in the case of the cusp anomalous dimension [61,62]: the interplay of this class of contributions to the cusp with similar contributions to multi-particle correlators is a very interesting open problem, with subtle connections to the factorisation of collinear poles. We note that an interesting alternative approach to IR exponentiation, focusing not on diagrammatics but on the symmetries and renormalisation properties of Wilson-line correlators, was developed in Refs. [63,64]. The present paper focuses on the colour structure of the perturbative exponent W n for the soft operator in Eq. (1.3), mostly at the four-loop order. Specifically, we complete the task of computing the mixing matrices determining the exponentiated colour factors for all four-loop Cwebs, initiated in Ref. [1], and we present some general properties of mixing matrices, valid to all orders. We begin by reviewing existing results on diagrammatic exponentiation in Section 2, where we discuss the concept of Cweb and the main properties of webs and Cwebs, and we list the four-loop Cwebs connecting two and three Wilson lines. Section 3 describes the replica method [58] for generating the Cweb mixing matrices, and its implementation in our code. In Section 4, we give two detailed examples, presenting the mixing matrices for one two-line Cweb and for one three-line Cweb. In Section 5, we present some observations on the structure of mixing matrices, including a derivation of two-and three-dimensional mixing matrices directly based on their combinatorial properties, bypassing the replica algorithm. Finally, we summarise our results and perspectives in Section 6. In the Appendix, we list in detail the mixing matrices for all the Cwebs discussed in the main text.

Webs and Cwebs for correlators of two and three Wilson lines
In this section we introduce mixing matrices and their main properties, using the language of webs for multi-particle amplitudes, developed in [57][58][59][65][66][67][68][69]. We then proceed to introduce correlator webs, or Cwebs, following [1], exploiting the fact that the combinatorial problem associated with gluon attachments to Wilson lines is structurally factored from the colour algebra associated with gluon subdiagrams. Finally, as an example, we enumerate Cwebs with lowest-order contributions at three and four loops and connecting two or three Wilson lines.
A web is defined as a set of fixed-order Feynman diagrams contributing to a Wilson-line correlator, which differ only by the permutation of their gluon attachments to each Wilson line. If we write every Feynman diagram D as the product of its kinematic factor K(D) and its colour factor C(D), then the logarithm of the correlator, W n , can be written as where the Exponentiated Colour Factor (ECF) for diagram D, denoted by C(D), is a linear combination of the colour factors of the diagrams which are present in the web. Thus we can write where the sum runs over all the diagrams present in web w and R w is the web mixing matrix. Using Eq. (2.1) and Eq. (2.2) we obtain an expression for the correlator in Eq. (1.5), where the block-diagonal matrix R is built using the web mixing matrices R w as blocks. In turn, each web w can be written as In this language, W n = w n , where the sum extends to all webs arising in the presence of n Wilson lines, order by order in perturbation theory.
3 (2, 2, 2) It is evident that a computation of the soft anomalous dimension matrix crucially involves these mixing matrices: they were studied systematically up to three loops, and several of their all-order properties were uncovered, in [57][58][59][65][66][67][68]. The first important point is that they are idempotent i.e. R 2 w = R w , and thus they are projection operators; in particular, their eigenvalues can either be 1 or 0. Acting on the right, according to Eq. (2.4), they select a subset of the possible colour factors, while, acting on the left, they identify the set of linear combinations of kinematic factors that will contribute to the exponent. The rank r w of the mixing matrix R w equals the number of independent colour factors onto which it projects. Web mixing matrices, further, enforce the non-abelian exponentiation theorem [67], stating that all exponentiated colour factors correspond to Feynman diagrams where the gluon sub-diagram is completely connected. Two further combinatorial properties have important physical consequences: first, web mixing matrices of dimension d w > 1 obey a row sum rule, stating that the sum of the elements of each row vanishes, D R DD = 0; this ensures that diagrams contributing to the exponent are a proper subset of all diagrams contributing to the correlator. In addition, mixing matrices are conjectured [59] to obey a column sum rule, which can be written as D s(D)R(D, D ) = 0, where s(D) is called the column weight vector. For a diagram D, s(D) denotes the number of different ways in which the connected gluon sub-diagrams can be sequentially shrunk to the common origin of the Wilson lines. This second combinatorial property has important consequences for the renormalisation of the correlator: it ensures that web mixing matrices project on kinematic factors which are free from ultraviolet sub-divergences; this, in turn, simplifies the evaluation of the exponent, implying that the UV counterterms are independent of the selected IR regulator.
A close examination of the properties of mixing matrices confirms the intuitive notion that the internal structure of connected gluon sub-diagrams does not affect the combinatorics of their attachments to the Wilson lines: a connected multi-gluon correlator may generate a number of independent 'internal' colour factors, but each one of them can be treated separately when considering the connection to Wilson lines. This observation motivated the introduction, in Ref. [1], of the concept of Cwebs, building Wilson-line correlators in terms of connected gluon correlators instead of individual diagrams.
More precisely, a Cweb is a set of skeleton diagrams, built out of connected gluon correlators attached to Wilson lines, and closed under permutations of the gluon attachments to each Wilson line. Clearly, the main difference between webs and Cwebs is the fact that Cwebs are not fixedorder quantities, but admit their own perturbative expansion in powers of g. Cwebs strongly simplify the counting and organisation of contributions to the exponent, especially at high orders, where radiative corrections to gluon sub-diagrams become important and proliferate. They may also provide ingredients to generalise to multi-line correlators the fascinating arguments given in Ref. [70] for the two line case (see also Ref. [71]). Cweb mixing matrices are identical to those of their constituent webs, as the properties of mixing matrices are blind to gluon interactions away from the Wilson lines, and only depend on ordering of their attachments on the Wilson lines, which is same in both the definitions.   To identify Cwebs, we adopt the notation of Ref. [1], writing W (c2,...,cp) n (k 1 , . . . , k n ), where n denotes the number of Wilson lines, c m denotes the number of m-point correlators forming the Cweb, and k l denotes the number of attachments on each Wilson line, where, without loss of generality, we take k 1 ≤ k 2 ≤ . . . ≤ k n . When a further degeneracy arises, we distinguish the Cwebs sharing the same notation by roman numerals.
A systematic way of generating all Cwebs recursively at any perturbative order was presented in [1]. Assuming one has already enumerated all Cwebs at O(g 2n−2 ), the steps to generate all Cwebs at O(g 2n ) are listed below, keeping in mind the possibility of adding a Wilson line with no attachments at lower orders.
4. Discard Cwebs that are given by the product of two or more disconnected lower-order webs.
5. In a massless theory, discard all self-energy Cwebs, where all gluon lines attach to the same Wilson line, as they vanish as a consequence of the eikonal Feynman rules. 6. Discard Cwebs that have been generated by the procedure more than once.
Following the above steps, and focusing on the results for two and three Wilson lines, we find four Cwebs connecting two Wilson lines at three loops, shown in Fig. 1, and six Cwebs connecting three Wilson lines at three loops, shown in Fig. 2. At four loops, we find eight Cwebs connecting two lines, shown in Fig. 3, and twenty-two Cwebs connecting threee lines, shown in Fig. 4. Cwebs connecting four and five lines were studied in [1]. Note that Cwebs that differ only by the permutation of the Wilson lines are identical in structure and we do not include them in our counting.
3 Implementing a replica algorithm to generate Cweb mixing matrices In this section, we will briefly describe our implementation of the replica method, introduced in [72], and used in this context by [58,73], for the calculation of the web mixing matrices. The starting point is the path integral expression for the Wilson line correlator where S(A a µ ) is the classical action. The method is based on the construction of a replicated theory, where we replace the gluon field A a µ with N r identical copies, A a, i µ (i = 1, . . . , N r ), which do not interact with each other. Further, we associate a copy of each Wilson line to each replica, effectively replacing each Wilson line in Eq. (3.1) with the product of N r Wilson lines. The correlator in the replicated theory can then be written as Using Eq. 3.2, one can calculate W n (γ i ) by computing the coefficient of the O(N r ) terms of the replicated Wilson line correlator. As there are no interaction vertices involving gluons belonging to different replicas, each connected gluon correlator in a Cweb is naturally assigned a unique replica number: this means that the replica method designed for webs in [58,73] immediately generalizes to Cwebs. The general structure of the algorithm can be summarised as follows.
• Assign a replica number i, 1 ≤ i ≤ N r , to each connected gluon correlator in a Cweb.
• Introduce a replica ordering operator R, which acts on the colour generators on each Wilson line ordering them according to their replica numbers. More precisely, if T (i) k denotes a generator associated with the emission of a gluon belonging to replica i from Wilson line k, then R acts on a product of T k by preserving the given order if i ≤ j, and by reversing it if i > j. The action of R effectively replaces a skeleton diagram with another one belonging to the same Cweb.
• In order to compute the colour factors in the replicated theory, one then needs to determine the number of possible hierarchies of replica numbers occurring in a Cweb with m connected pieces, which we call h(m), and the number of occurrences of a particular hierarchy in the presence of N r replicas, which we call M Nr (h). h(m) counts the number of weak orderings on a set of m elements, and is known in number theory and enumerative combinatorics as ordered Bell number or Fubini number [74]. The first few ordered Bell numbers are given by h(m) = {1, 1, 3, 13, 75, 541} for m = 0, 1, 2, 3, 4, 5, where we included the case m = 0 to follow the mathematical convention. On the other hand, the multiplicity of a given hierarchy h, containing n r (h) distinct replicas, is easily determined, and it is given by • With these ingredients, one can compute the colour factor of a skeleton diagram D in the replicated theory, which is given by where R C(D) h is the replica-ordered colour factor of diagram D for a given hierarchy h. The exponentiated colour factor for diagram D is finally given by the O(N r ) terms in Eq. (3.4).
In order to compute the Cweb mixing matrices by applying the replica trick algorithm, we developed an in-house Mathematica code which was used to calculate the mixing matrices for Cwebs connecting four and five Wilson lines at four loops in Ref. [1]. Interestingly, lowering the number of Wilson lines increases the combinatorial complexity of the problem, which may lead to a critical slowing down of the algorithm for Cwebs involving many connected correlators attached to few Wilson lines. The basic reason is the fact that, with few Wilson lines, at high orders the typical number of gluon attachments to each Wilson line grows, and so does the size of a typical mixing matrix. To keep the runtime under control, one has to take full advantage of the symmetries of the problem, including Bose symmetry for each connected gluon correlator, and the symmetry under the exchange of identical gluon correlators connecting the same set of Wilson lines. Below, we briefly present the main steps in our current implementation of the algorithm.
• The starting point is the building blocks for four-loop Cwebs, which are the sets of connected gluon correlators shown in Fig. 5. The code generates all four-loop Cwebs by combining the building blocks and attaching them to the Wilson lines in all possible ways. This is done starting from the most intricate case, which is the attachment of a set of four two-point correlators in all possible ways to the Wilson lines. Proceeding with other combinations of correlators, the combinatorial problem simplifies, and finally the five-point correlator shown in Fig. 5 only produces completely connected Cwebs with a mixing matrix R = 1.
• The above step generates Cwebs which are related to one another by the permutation of the Wilson lines. As far as the colour structure is concerned, they are all identical, so the duplicates are removed.
• The code then assigns a distinct replica number to each connected correlator.
• Starting from a single diagram, the code generates all diagrams for a Cweb by identifying the gluon attachments on each Wilson line, and placing the attachments originating from different gluon correlators in different packs. These packs are then shuffled.
• A subroutine generates all possible hierarchies h for a Cweb. As an example, for a Cweb with two connected gluon correlators, the subroutine generates the three hierarchies h = {i = j, i > j, i < j}, where i and j denote the replica numbers. The code then determines the number of distinct replicas n r for each hierarchy. For example, for the hierarchy i = j one has n r = 1, whereas for hierarchies i > j and i < j one has n r = 2. Using Eq. 3.3, one then computes the multiplicity M Nr (h).
• At this stage, one can construct, for each Cweb, a table of the form of Table 1 of Ref. [58], and finally determine the mixing matrix R w for Cweb w. The code then diagonalizes where d w is the dimension of the mixing matrix, and r w its rank. The matrix Y w , acting from the left on the vector of colour factors for the diagrams in w, computes the exponentiated colour factors.
• In its present form, the code does not generate self-energy Cwebs, where a connected gluon correlator is attached only to a single Wilson line, since such Cwebs vanish in the massless theory. The code however generates disconnected Cwebs (where the set of Wilson lines can be partitioned in subsets whose elements are not connected by any gluons): the vanishing of the mixing matrices for disconnected Cwebs works as a check on the code. Further checks are provided by verifying two known properties: the idempotence of mixing matrices, and the row sum rule. The conjectured column sum rule, on the other hand, is verified a posteriori. Figure 5: Combinations of connected correlators that can form Cwebs at four loops.
We note that, for a number of Cwebs with non-trivial symmetries, the procedure of generating all diagrams by shuffling the packs on each of the Wilson lines will produce multiple duplicates of each diagram, which can be obtained from each other by simply relabelling some of the correlators.
In the case at hand, this happens for the Cwebs: W 2 (4,4), shown, with all other Cwebs connecting three and four lines at four loops, in the Appendix. To contain the computation time, it is crucial to remove these duplicates before applying the replica algorithm. For example, the Cweb W (4) 2 (4, 4) has a total of 24 diagrams, but each diagram is generated in 24 duplicate copies: the present version of the code automatically identifies and removes all duplicates.
To further illustrate this issue consider the simple case of the two-loop Cweb W 2 (2, 2), which has two attachments on both Wilson lines, and is shown in Fig. 6. The shuffle (which is just a single permutation in this case) of the correlators on each Wilson line generates the four diagrams shown in Fig. 6. As the colour indices a and b in Fig. 6 are summed over, diagrams C 1 and C 3 are duplicates, and similarly C 2 and C 4 . The code identifies the duplicates and deletes C 3 and C 4 .
The exponentiated colour factors generated by the code have open colour indices for each gluon attachment to the Wilson lines: these open indices must then be contracted with the colour structures arising from the connected gluon correlators. These colour structures, in turn, are constrained by Bose symmetry and gauge invariance. For example, a two-gluon correlator joining any two lines i and j will generate only the dipole structure T i · T j , while a three-point gluon correlator is conjectured to be proportional to the structure constants f abc to any perturbative order. For four-gluon correlators, three colour structures are available at tree level, while more complicated ones will arise when loop corrections are considered, including contributions proportional to quartic Casimir operators of the gauge algebra [75][76][77].

Examples of four-loop Cwebs connecting two and three Wilson lines
In this section, we present in some detail the calculation of two four-loop Cwebs connecting respectively two and three Wilson lines. We present their mixing matrices, the diagonalizing matrices Y, and we give explicit results for the exponentiated colour factors. The results for all other four-loop Cwebs connecting two and three Wilson lines are presented in the Appendix.

A two-line Cweb at four loops
As a first example, we consider the two-line Cweb W (1,0,1) 2 (2, 4), which contains one two-gluon correlator, no three-gluon correlators, and one four-gluon correlator; furthermore, there are two gluon attachments on line 1 and four gluon attachments on line 2. Clearly, the perturbative expansion for this Cweb starts at O(g 8 ), i.e. at four loops. A representative skeleton diagram for the Cweb is shown in Fig. 7. The available shuffles on line 1 and line 2 produce a total of eight skeleton diagrams for this Cweb which we denote by the order of the gluon attachments on the two lines, and are given by   (2, 4).
Note that, only C 4 and C 5 have nonzero s values. We find that the mixing matrix R, the diagonalising matrix Y and the diagonal matrix D for this Cweb are given by One may easily verify the properties of the mixing matrix: R is idempotent, the matrix elements in each row sums to zero, and the column sum rule is obeyed. After diagonalising the mixing matrix, we find that the rank of the mixing matrix is r w = 7, which means that there will be 7 independent exponentiated colour factors for this Cweb, which are given by We observe that all the exponentiated colour factors correspond to completely connected Feynman diagrams, which verifies, as expected, the non-abelian exponentiation theorem [67].

A three-line Cweb at four loops
As a second example, we select a three-line Cweb at four loops, labelled as W (1,0,1) 3,I (2, 1, 3), with a two-gluon correlator and a four-gluon correlator. Here the roman numeral I appears to distinguish this Cweb from a second one, shown in the Appendix, and denoted by W , which shares the same number of attachments to the Wilson lines and the same number of two-, three-and four-gluon correlators. A representative skeleton diagram is displayed in Fig. 8.  With this ordering, the skeleton diagram depicted in Fig. 8 is diagram C 1 . The weight factors s(C i ) for the diagrams in this Cwebs are given by We find the mixing matrix R, the diagonalising matrix Y , and the diagonal matrix D, as so that the rank in this case is r w = 5. Once again, this mixing matrix satisfies all the desired properties: idempotence, row sum rule, and column sum rule. The exponentiated colour factors are As expected, they all correspond to connected Feynman diagrams.

Observations on mixing matrices and their direct construction
The present work, together with the results presented in [1], as well as earlier work in Refs. [58,59,65,67], provides a considerable amount of empirical data about mixing matrices, partly summarised in Table 1 and in Table 2. In a massless theory, the total number of non-vanishing Cwebs having lowest-order contributions at O(g 8 ) or below is 79. There is just one Cweb at O(g 2 ), four at O(g 4 ), fourteen at O(g 6 ) and sixty at O(g 8 ). Of these Cwebs, thirteen are fully connected, so they consist of only one skeleton diagram and have a trivial mixing matrix, R = 1. Setting those aside, the dimensions of the remaining sixty-six mixing matrices range between d w = 2 and d w = 36, while their ranks range between r w = 1 and r w = 29. To give a flavour of the distribution In Table 2 we present the ranks of all mixing matrices appearing up to four loops, for different matrix sizes. It is clear that the combinatorial problems associated with mixing matrices are nontrivial, and indeed interesting connections to abstract combinatorics have already been uncovered and exploited in Refs. [65,66,68,69]. In particular, while the dimensions of mixing matrices are easily computable from the diagrammatic structure of their Cwebs, their ranks are not in general predictable with current knowledge, and Table 2 does not display a discernible pattern.
In what follows, we provide some simple results emerging from the empirical data up to four loops, which in some cases allow for all-order generalisations. In particular, low-dimensional mixing matrices are highly constrained and can be uniquely determined from their general properties up to dimension d w = 3, and partially at d w = 4.  A first natural question to ask is how many different mixing matrices can be generated with a given dimensionality d w . This requires establishing when two mixing matrices should be considered different: given their definition, we take the viewpoint that two matrices are in the same equivalence class if they are connected by a permutation of their rows or of their columns: this corresponds to permuting the diagrams in the web and the components of the vector of colour factors; row and column sum rules are preserved by these permutations, provided one permutes the components of the s vectors appropriately. On the other hand, taking general linear combinations of rows or columns of mixing matrices does not have a diagrammatic interpretation and we do not consider it. With this definition of equivalence, we note that matrices of dimension d w = 2 occur ten times up to four loops, but there are only two different matrices, one of which occurs only once. This is easily understood and is discussed in Section 5.1. Matrices of dimension d w = 3 occur six times up to four loops, but actually these are six occurrences of the same matrix: this is discussed in Section 5.2. At d w = 4, for the first time we have matrices of different rank (r w = 1, 3), which occur both at three and at four loops: it turns out that the rank-three matrix is uniquely determined, while in principle different d w = 4 matrices of lower rank could appear at higher orders. For d w = {6, 8,9,12,16,18,24, 36}, the numbers of inequivalent mixing matrices arising up to four loops is n i = {6, 4, 1, 6, 1, 3, 8, 2}, and in several cases inequivalent matrices of the same rank are found. Given these data, we now examine in more detail the low-dimensional cases.

On two-dimensional mixing matrices
All matrices of dimension two appearing at three and at four loops are the same and have the form (5.1) At two loops, however, the matrix in Eq. (5.1) appears once, but, in addition, one more two dimensional matrix appears, in the case of Cweb W 2 (2, 2). It is It is easy to show that indeed these are the only two possible two-dimensional mixing matrices at any perturbative order, and furthermore Eq. (5.2) occurs only once, at two loops. With n ≥ 3 Wilson lines, a two-dimensional mixing matrix requires a Cweb with precisely two gluon attachments from distinct correlators to a specific Wilson line, while all other Wilson lines must allow no shuffles, i.e. they must each be attached to gluons from a single correlator. The s vector for such a configuration is always s = {1, 1}, as seen from the example in Fig. 9, since in each contributing diagram one of the two correlators involved can always be shrunk to the origin without affecting the other.
Armed with the knowledge of s, we can impose the row and column sum rules on a generic 2 × 2 matrix, obtaining 3) The normalisation can be fixed by observing that the matrix has rank r = 1, and, since it is a projection operator, the trace of the matrix must equal the rank. This leads uniquely to Eq. (5.1). The only exception to the above reasoning occurs at two loops, when there are only two Wilson lines, and two dressed gluon propagators attach to both lines, as depicted in Fig. 10. The dimension of the mixing matrix is reduced from 4 to 2 by the symmetry of the correlators, but in this case the s vector is given by s = {0, 1}, since for diagram (a) there are no possibilities to shrink sequentially the two correlators. Once again, imposing the row and column sum rules, with the new s vector, and normalising the trace of the matrix to its rank, we find uniquely Eq. (5.2). We conclude that, at any perturbative order, all two-dimensional mixing matrices will have the form of Eq. (5.1), with the sole exception Eq. (5.2) in the two-loop, two-line case we just examined.

The three-dimensional mixing matrix
There are a total of six three-dimensional mixing matrices up to O(g 8 ): two arising at three loops and four at four loops. It turns out that these are six occurrences of the same matrix, which is It is not difficult to prove that this is a general result, and Eq.  The column weight vector for this Cweb is s = {1, 0, 1}. Importantly, we notice that diagram (b), which has s = 0, cannot be generated from diagrams (a) and (c), which have s = 1, by the action of the replica ordering operator: under that action, all attachments belonging to a single gluon correlator (and thus to a single replica) move together. In other words, diagram (b) is a singlet under the mixing operation, whereas diagrams (a) and (c) form a doublet. Taking this into account, and applying the row and column sum rules, one sees that the mixing matrix must be of the form Finally, using the idempotence property R 2 = R and imposing that the trace of the matrix must equal its rank (which is r w = 2 in this case) one can fix the two remaining parameters and obtain precisely Eq. (5.4). We conclude that all three-dimensional mixing matrices at any perturbative order will be given by Eq. (5.4).

Results for higher-dimensional mixing matrices
The reasoning of the previous two sections generalises in two different directions. On the one hand, one may consider mixing matrices whose dimension is a prime number, d w = p. As in the case p = 3, such matrices can only arise from Cwebs where one Wilson line has a single attachment from a gluon correlator G 1 , together with p − 1 attachments from a second gluon correlator G 2 , with no other Wilson lines providing further shuffles. The column weight vector for such a web is again of the form s = {1, 0, . . . , 0, 1}, with p − 2 zeroes. Furthermore, once again, the p − 2 diagrams where the G 1 gluon is inserted between attachments of G 2 gluons cannot be reached by the action of the replica ordering operator, and are singlets under the mixing operation, while the two diagrams where the G 1 gluon is on one or the other side of the G 2 cluster form a doublet. This information, together with the row sum rule, is sufficient to conclude that the rank of the mixing matrix is p − 1. Imposing idempotency, and normalising by setting the trace of the matrix equal to the rank, one finds that mixing matrices of prime dimension p are unique, and they can be written as This reasoning explains the fact that the largest prime number appearing as the dimension of a mixing matrix at four loops is p = 3. The first mixing matrix of dimension d w = p = 5 will appear at O(g 10 ), arising when four gluons from a five-point gluon correlator are shuffled with a single gluon from a two-point correlator attaching on the same Wilson line. Conversely, we can consider any Cweb where a single gluon from a correlator G 1 is shuffled on a Wilson line with a set of n − 1 gluons from a second correlator G 2 , and no other shuffles are available on other lines, even when n is not a prime number. The reasoning above still holds, so that the mixing matrix will have dimension n and rank n − 1, and it will be of the form of Eq. (5.6). In this case, however, we cannot argue that the matrix is unique, since other sets of correlators and attachments will likely lead to other matrices with the same dimension. A case in point is d w = 4, where we find two distinct mixing matrices, occurring in two three-loop Cwebs and in seven four-loop Cwebs. One of these two matrices, as expected, has rank r = 3 and is of the form of Eq. (5.6): it arises, among others, from the Cweb W 1,0,1 3 (1, 1, 4), depicted in Fig. 4(p), which is of the form just described. The second type of d w = 4 mixing matrix emerging from the algorithm has rank r w = 1. It arises for example in the case of Cweb W (2,1) 4, I (1, 2, 2, 2), discussed in Ref. [1], and it has a column weight vector given by s = {1, 2, 2, 1}. In this case, there are no 'singlet' skeleton diagrams in the set, so the available constraints are not sufficient to determine this second matrix uniquely. Clearly, we are barely scratching the surface of this combinatorial problem, and further dedicated studies are likely to bring about a much deeper understanding.

Summary and Outlook
The study of diagrammatic exponentiation of Wilson-line correlators has a long history, and has provided many important insights concerning the infrared structure of perturbative gauge amplitudes. In this paper, we developed the idea of correlator webs, or Cwebs, introduced in Ref. [1], and we completed the study of all Cweb mixing matrices appearing up to four loops. Specifically, we have listed all four-loop Cwebs connecting two and three Wilson lines, following the recursive algorithm developed in [1]. Using an improved version of our in-house Mathematica code, we have computed the mixing matrices for these Cwebs using the replica method [58,73]. In all cases, the mixing matrices we have found verified both the proven and the conjectured properties listed in the literature [58,59,[65][66][67][68], in particular the general form of the non-abelian exponentiation theorem: all the exponentiated colour factors correspond to fully connected gluon sub-diagrams [67].
As shown in the examples reported in Section 4.1 and in Section 4.2, once the mixing matrix is determined, exponentiated colour factors are readily computed by multiplying the diagonalising matrix Y w for the given Cweb times the column vector formed out of the colour factors of individual skeleton diagrams. One must then use the commutation relations of the gauge algebra to simplify products of colour generators attached to each Wilson line. All mixing matrices for four-loop Cwebs connecting two and three Wilson lines are presented in the Appendix. When the number of Wilson lines is small, as in this case, at high orders many attachments are present on each Wilson line, and the commutation relations must be applied repeatedly, generating a large number of ECFs, many of which have lengthy expressions. As a consequence, we refrain from presenting explictly the ECFs in the Appendix: given the mixing matrices, the steps required to derive the ECFs are straightforward. In any case, ECFs can be obtained from the authors upon request in the form of Mathematica code.
Given the considerable accumulated data about mixing matrices, from the knowledge of the 66 non-trivial Cwebs arising up to four loops, we have taken the opportunity to examine their emerging properties. We noted that, while the dimension d w of mixing matrices is directly computable, given the set of gluon correlators and attachments in the selected Cweb, the rank r w of the resulting matrices is not easy to predict in general. On the other hand, we have observed that low-dimensional mixing matrices are highly constrained, and we have derived some general results. We have shown that two-and three-dimensional mixing matrices are uniquely determined to all orders in perturbation theory by their general properties, given the limited possibilities available for their column weight vectors. We have also shown that mixing matrices whose dimension is a prime number are unique, and we have given their form; furthermore, we have uniquely determined the mixing matrices for an infinite series of simple Cwebs, consistently with the low-order examples that we have explicitly computed.
Together with the results presented in [1], the results for three-and two-line Cwebs at four loops presented here complete the computation of the mixing matrices for all four-loop Cwebs. A direct calculation of mixing matrices using posets (partial ordered sets) was also explored in [66,68,69] for a class of webs: our results provide additional useful data for the direct construction of the mixing matrices, and for the study of their interesting algebraic and combinatorial properties.
We emphasize that exponentiated colour factors are not all independent, and it is necessary to reduce them to a basis, implementing the constraint of colour conservation, which implies that, when acting on a physical amplitude, the colour-insertion operators must satisfy i T i = 0, where the sum runs over all hard particles, represented here by the Wilson lines. The consequences of colour conservation were implemented at three loops in [50], and studied at four loops in [52,78]. Together with Bose symmetry, they allow to reduce the form of the soft anomalous dimension matrix to a relatively simple parametrisation in terms of a few scalar functions. Of particular interest is the appearance, at four loops, of contributions proportional to quartic Casimir operators of the gauge algebra, which are fully known for the simple case of two Wilson lines [61,62,79,80], but, as yet, undetermined in the general case, although partial results have begun to emerge [53,81]. Our results provide all the necessary colour ingredients for the complete four-loop calculation, and other tools are available for the study of colour structures at high orders, such as the effective vertex analysis of Ref. [67], and the generating functional approach of Ref. [82]. That being said, of course by far the most difficult challenge remains the calculation of the kinematic contributions, in particular for four-loop connected diagrams. The interplay of colour and kinematics which emerges with striking power and simplicity in the soft anomalous dimension matrix remains one of the most challenging and interesting topics in the study of perturbative non-abelian gauge theories.

Acknowledgments
AT and LM and SP would like to thank MHRD Govt. of India for the GIAN grant (171008M01), "The Infrared Structure of Perturbative Gauge Theories" and for the SPARC grant (P578) "Perturbative QCD for Precision Physics at the LHC", which were crucial to the completion of the present research. SP would also like to thank MHRD Govt. of India for an SRF fellowship, and the University of Turin and INFN Turin for warm hospitality during the course of this work.

A All four-loop mixing matrices connecting two and three Wilson lines
In this appendix we give results for all the Cwebs that appear at 4 loops in the scattering amplitude, that can connect two or three Wilson lines. Throughout the list, R and D denote the mixing matrix and the diagonalized matrix respectively. Clearly, D is just a representation of the rank of R, and we write it as D = (1 r , 0), where r is the rank. We display only one skeleton diagram per web, and we explicitly give the order of the shuffles that generate the other diagrams, which is tied to the order the columns of the mixing matrix in the chosen basis. The tables of shuffles also give the components of the column weight vector for each Cweb. In Tables 3 and 4   of Cwebs, the total number of skeleton diagrams for each Cweb, the number of replica hierarchies generated in the application of the replica method, and the rank of the resulting mixing matrices. We omit from the list the Cwebs that are composed of a single skeleton diagram, whose mixing matrix is just a number, R = 1.
A.1 Cwebs connecting three Wilson lines (1, 2, 3) This Cweb has six diagrams, one of which is displayed below. The table gives the chosen order of the six shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences s-factors The R matrix is given by This Cweb is a second kind of Cweb with notation W (1,0,1) 3 (1, 2, 3) and has two diagrams, one of which is displayed below. The table gives the chosen order of the two shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences s-factors The R matrix is given by 3,I (1, 2, 3) This Cweb has three diagrams, one of which is displayed below. The table gives the chosen order of the three shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences s-factors The R matrix is given by 3,II (1, 2, 3) This Cweb has six diagrams, one of which is displayed below. The table gives the chosen order of the six shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences s-factors The R matrix is given by

W
(2,1) 3,I (2, 2, 3) This Cweb has 24 diagrams, one of which is displayed below. The table gives the chosen order of the 24 shuffles of the gluon attachments, and the corresponding s factors.

Sequences
s-factors The R matrix is given by 6. W

Diagrams
Sequences The R matrix is given by 3,III (2, 2, 3) This Cweb has twelve diagrams, one of which is displayed below. The table gives the chosen order of the twelve shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences The R matrix is given by 8. W (1,0,1) 3 (1, 2, 3) This Cweb has three diagrams, one of which is displayed below. The table gives the cho-sen order of the three shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences s-factors The R matrix is given by 3,I (2, 2, 2) This Cweb has two diagrams, one of which is displayed below. The table gives the chosen order of the two shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences The R matrix is given by R = (2, 2, 2) This Cweb has four diagrams, one of which is displayed below. The table gives the chosen order of the four shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences The R matrix is given by 3,II (2, 2, 2) This Cweb has eight diagrams, one of which is displayed below. The table gives the chosen order of the eight shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences s-factors The R matrix is given by

Sequences
s-factors The R matrix is given by 13. W

Diagrams
Sequences s-factors The R matrix is given by

W
(1,0,1) 3 (1, 1, 4) This Cweb has four diagrams, one of which is displayed below. The table gives the cho-sen order of the four shuffles of the gluon attachments, and the corresponding s factors.

Diagrams Sequences s-factors
The R matrix is given by 15. W (0,2) 3 (1, 1, 4) This Cweb has six diagrams, one of which is displayed below. The table gives the chosen order of the six shuffles of the gluon attachments, and the corresponding s factors.

Diagrams Sequences s-factors
The R matrix is given by

W
(2,1) 3 (1, 2, 4) This Cweb has 24 diagrams, one of which is displayed below. The table gives the chosen order of the 24 shuffles of the gluon attachments, and the corresponding s factors.

Sequences
s-factors The R matrix is given by (1, 2, 4) This Cweb has 12 diagrams, one of which is displayed below. The table gives the chosen order of the 12 shuffles of the gluon attachments, and the corresponding s factors.

Sequences
s-factors The R matrix is given by  The R matrix is given by The R matrix is given by 20. W

Sequences
s-factors  (2,4) This Cweb has eight diagrams, one of which is displayed below. The table gives the chosen order of the eight shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences s-factors The R matrix is given by This Cweb has nine diagrams, one of which is displayed below. The table gives the chosen order of the nine shuffles of the gluon attachments, and the corresponding s factors.  This Cweb has six diagrams, one of which is displayed below. The table gives the chosen order of the six shuffles of the gluon attachments, and the corresponding s factors.

Diagrams Sequences s-factors
The R matrix is given by This Cweb has nine diagrams, one of which is displayed below. The table gives the chosen order of the nine shuffles of the gluon attachments, and the corresponding s factors.

Diagrams
Sequences s-factors The R matrix is given by

W
(2,1) 2 (3,4) This Cweb has thirty six diagrams, one of which is displayed below. The The R matrix is given by 2 (4, 4) This Cweb has 24 diagrams, one of which is displayed below. The table gives the chosen order of the 24 shuffles of the gluon attachments, and the corresponding s factors.

Sequences
s-factors The R matrix is given by This completes our listing of all Cwebs with a perturbative expansion starting at O(g 8 ), and connecting two and three Wilson lines.