Multiparton webs beyond three loops

Correlators of Wilson-line operators are fundamental ingredients for the study of the infrared properties of non-abelian gauge theories. In perturbation theory, they are known to exponentiate, and their logarithm can be organised in terms of collections of Feynman diagrams called webs. We study the classification of webs to high perturbative orders, proposing a set of tools to generate them recursively: in particular, we introduce the concept of Cweb, or correlator web, which is a set of skeleton diagrams built with connected gluon correlators, instead of individual Feynman diagrams. As an application, we enumerate all Cwebs entering the soft anomalous dimension matrix for multi-parton scattering amplitudes at four loops, and we compute the mixing matrices for all Cwebs connecting four or five Wilson lines at that loop order, verifying that they obey sum rules that were derived or conjectured in the literature. Our results provide the colour building blocks for the calculation of the soft anomalous dimension matrix at four-loop order.


Introduction
Studies of the structure of infrared (IR) singularities that appear in scattering amplitudes in gauge field theories have a long and rich history, and have led to remarkable all-order insights into the organisation of the perturbative expansion [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. In the computation of loop corrections to scattering amplitudes, IR singularities arise when a virtual particle flowing in a loop becomes soft or collinear to one of the external particles. Upon constructing, from amplitudes, a well-defined physical observable, these singularities are cancelled by the contribution of real emission diagrams, which must be integrated over the phase space for undetected real radiation. The singularities, however, often leave their imprints on the perturbative expansion, in the form of potentially large logarithms of kinematic variables, which may need to be resummed in order to obtain precise predictions. Such resummations 1 are made possible by the universal nature of infrared radiation, which results in factorised expressions for scattering amplitudes, where soft and collinear effects are organised by universal functions, depending only on the quantum numbers of the external particles involved in the scattering, but not on the specific nature of the hard process being considered [14,17,18]. 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 JHEP05(2020)128 study. We note 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 [19][20][21]; furthermore, we emphasise that a detailed knowledge of infrared singularities is also important for collider phenomenology at finite orders: indeed, the cancellation of singularities between squared matrix elements with different numbers of external particles is difficult to implement at high orders for complicated collider observables, which must be evaluated numerically. In practice, the cancellation must be performed analytically, and the construction of general and efficient subtraction procedures beyond next-to-leading order (NLO) is an ongoing effort (see, for example, [22][23][24][25][26][27][28][29][30][31]).
The methods developed in this paper concern the evaluation of Wilson-line correlators of the general form where Φ r (γ) is a Wilson-line operator evaluated on a (smooth) space-time contour γ, defined by where A µ r (x) = A µ a (x) T a r is a non-abelian gauge field, with T a r the generators of the representation r of the gauge group. The smooth contours γ k can be closed (in which case the correlator is gauge-invariant), or open: in this case, the correlator is a colour tensor with open colour indices in representations r k attached to the ends of each Wilson line. While most of our considerations will apply in general to any correlator of the form (1.1), we will be especially interested in the soft colour operator associated with multi-particle scattering amplitudes in gauge theories, which encodes all their soft singularities. This soft operator is of the form of eq. (1.1), with the contours γ k given by semi-infinite straight lines extending from the origin along directions β k , corresponding to the four-velocities of the particles participating in the scattering. In this case we write more explicitly S n β i · β j , α s (µ 2 ), ≡ 0| n k=1 Φ β k (∞, 0) |0 , Φ β (∞, 0) ≡ P exp ig ∞ 0 dλ β · A(λβ) , (1.3) where α s = g 2 /(4π), for simplicity we did not display the representations to which the Wilson lines belong, and we introduced the dimensional regularisation parameter , setting the space-time dimension to d = 4 − 2 .
Soft operators of the form (1.3) are highly singular, being affected by ultraviolet, soft, and, in case β 2 i = 0, collinear divergences: as a consequence, special care is required to evaluate them [32][33][34][35]. In the massless case, this can be done by introducing auxiliary regulators: for collinear divergences, one may set β 2 i = 0, for soft divergences one may for example introduce a smooth exponential long-distance cutoff on gluon interactions, as discussed in refs. [34,35], while retaining dimensional regularisation for ultraviolet singularities. The bare correlator can then be evaluated and renormalised, yielding the desired answer.

JHEP05(2020)128
Wilson-line correlators of the form of eq. (1.1) have the following basic properties.
• After renormalisation, n-line correlators obey renormalisation group equations which lead, in dimensional regularisation, to exact exponentiation in terms of a soft anomalous dimension matrix Γ n . For straight-line correlators, of the form of eq. (1.3), one may write S n β i · β j , α s (µ 2 ), = P exp − 1 2 It is important to note that, in the massless case, the soft anomalous dimension Γ n is affected by collinear singularities, which must be properly organised in terms of appropriate jet functions [14]; collinear-finite contributions can be computed in the massless case by considering Wilson lines slightly tilted off the light cone, and taking the light-cone limit in the intermediate stages of the calculation [36]. The soft anomalous dimension matrix Γ n is a central quantity for the study of perturbative non-abelian gauge theories, and has been the focus of much theoretical work. It was computed at one loop in [37] (see also [38]); at two loops in the massless case in [39,40], and in the massive case in [41][42][43][44][45]; finally, at three loops in the massless case in [36,46].
• General Wilson-line correlators of the form of eq. (1.1) obey a non-trivial form of diagrammatic exponentiation, so that one can write S n (γ i ) = exp W n (γ i ) , (1.5) where the logarithm of the correlator, 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. [47][48][49], for the case of two straight, semi-infinite Wilson lines. For general configurations, it was proven in refs. [50,51]. Feynman diagrams contributing to W n (γ i ) are collectively called webs. For an abelian theory, webs are connected diagrams; for a non-abelian theory, if only two Wilson lines are present, webs are 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. The properties of webs will be further discussed in section 2, and a useful generalisation of the concept of web will be proposed in section 3. Clearly, by means of webs, one can directly compute the soft anomalous dimension matrix Γ n .
• For the case when the contours γ i are semi-infinite, or infinite straight lines, all loop corrections to the bare correlator, in the absence of auxiliary IR regulators, vanish in dimensional regularisation, since they are given by scale-less integrals. Bare correlators are then exactly given by the unit matrix in the tensor product of the colour representations of the n Wilson lines. The renormalised correlator is therefore a 'pure JHEP05(2020)128 counterterm': in order to compute it, one must first construct an IR-regulated version of the bare correlator, whose loop corrections will not vanish, but will in general be regulator-dependent; one proceeds then to renormalise the regulated correlator, extracting the relevant UV counterterms; the set of these counterterms constitutes the renormalized version of the original correlator. It is important 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 [52].
• When the countours γ i are light-like, straight, semi-infinite Wilson lines, scale invariance imposes strong constraints on the functional dependence of the soft anomalous dimension matrix Γ n . Up to two loops, Γ n can only involve dipole correlations between Wilson lines [15,16,53,54]; beyond two loops, only quadrupole correlations can arise, which must depend on scale-invariant conformal cross ratios of the form : the first such correlations arise at three loops, with at least four Wilson lines, and were computed in ref. [36]; further correlations may arise only in association with higher-order Casimir operators, which may start contributing only at four loops, as discussed in section 5.3.
In this paper, we study the properties of the logarithms of Wilson-line correlators, W n (γ i ), with emphasis on their colour structures, extending earlier studies to higher orders in perturbation theory. After reviewing existing results on diagrammatic exponentiation in section 2, in section 3 we propose a natural extension of the concept of web, which we believe will prove useful for classification purposes and high-order studies. Subsequently, in section 4, we discuss how our definition leads to a simple recursive method to generate webs at (p+1) loops from those arising at p loops, and we discuss how this recursion can be implemented in a code based on the 'replica trick' used in ref. [51]. In section 5, we apply our method to construct the web mixing matrices at four loops for all webs connecting four or five Wilson lines, verifying that these matrices satisfy the expected properties, including conjectures that were proposed at lower orders; furthermore, we briefly discuss the colour structures that arise at four loops, verifying the compatibility of our results with the discussion in refs. [55,56]. We conclude, in section 6, with an outlook on possible future developments, while an appendix lists in detail the web mixing matrices for all the webs discussed in the main text.

Diagrammatic exponentiation for multiple Wilson-line correlators
Diagrammatic exponentiation in the eikonal approximation was first observed in QED [6], where W n (γ i ) is given by the sum of all connected photon subdiagrams (obtained by remov-   interpretation. For the case of two Wilson lines [47][48][49], W 2 is constructed from the set of diagrams which remain connected after cutting the two Wilson lines: these two-eikonalirreducible diagrams are called webs (see figure 2 for examples). It is important to note that, even in the simple case of two Wilson lines in a colour-singlet configuration, webs appear in W 2 with modified colour factors: more precisely, the only colour factors appearing in W 2 are those that correspond to connected gluon subdiagrams, such as the first one portrayed in figure 2. This provides an interesting generalisation of the abelian exponentiation, and points to further extensions to the multi-Wilson-line case. The problem of constructing the colour operator W n (γ i ) for the general case of n Wilson lines was solved in the remarkable series of papers [50][51][52][57][58][59][60][61], while an interesting alternative approach was developed in refs. [62,63]. Here we briefly summarise the main results, with special attention to the colour structures, which will be the main focus of our paper.
Let D be a Feynman diagram contributing to the correlator S n . Each such diagram can be written as the product of a kinematic factor K(D), depending on the four velocities β i (or more generally on the contours γ i in the case of non-straight Wilson lines) and a colour factor C(D). In full generality, the logarithm of the correlator, W n , can be written as linear combination of the same diagrams, with modified colour factors. We write and therefore do not contribute to the logarithm of the correlator: for example, for n = 2, one can show that all two-eikonal reducible diagrams have C(D) = 0. In general, ECFs are linear combinations of the ordinary color factors of sets of diagrams that differ only by the order of their gluon attachments to the Wilson lines. This naturally leads to the general definition of a web not as a single diagram, but as a set of diagrams that can be obtained from any representative element by permuting the gluon attachments to the Wilson lines: a simple example of a two-loop, three-line web involving two diagrams is presented in figure 3. As a consequence of these considerations, the sum in eq. (2.1) is naturally organised as a sum over webs, and for each diagram belonging to a given web w the ECF is expressed as where the sum extends to the diagrams belonging to web w, and R w (D, D ) is called web mixing matrix. Combining eq. (2.1) and eq. (2.2) one can express the Wilson-line correlator as where the sum extends to all diagrams appearing in the correlator, order by order in perturbation theory, and the matrix R is block-diagonal, with blocks corresponding to individual webs. 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.
Clearly, web mixing matrices are crucial quantities for the purpose of computing Wilson-line correlators, and therefore the soft anomalous dimension matrix. Their properties were extensively studied in refs. [50][51][52][57][58][59][60], and can be summarised as follows: • Web mixing matrices are idempotent, i.e. ∀w, R 2 w = R w , as a consequence of which their eigenvalues can only be either 1 or 0.

JHEP05(2020)128
• Operating on the vector of color factors for the diagrams of web w, the web mixing matrix R w acts as a projection operator, selecting a subset of the possible colour factors, in number equal to its rank.
• To all orders in perturbation theory, it can be shown [59] that the set of colour factors surviving the projection is the set of colour factors of connected gluon diagrams: this is the general form of the non-abelian exponentiation theorem.
• Acting on the vector of kinematic factors for the diagrams of web w, in the presence of an infrared regulator, the web mixing matrix R w projects onto kinematic factors that do not contain ultraviolet sub-divergences. This is crucial, in order to be able to isolate UV counterterms with no dependence on the infrared regulator.
• The elements of web mixing matrices obey the row sum rule D R w (D, D ) = 0, implying that at least one of the eigenvalues of R w must vanish.
In addition to these properties, which were proved in refs. [52,57,59], in the case of semiinfinite Wilson lines radiating out of a common origin, the matrix elements of web mixing matrices are also conjectured to obey a column sum rule which can be stated as follows. Given a diagram D, consider the set {D i c } of subdiagrams that remain connected after the Wilson lines are removed; we say that a connected subdiagram D i c can be shrunk to the common origin of the Wilson lines if all the vertices connecting the subdiagram to the Wilson lines can be moved to the origin without encountering vertices associated with other subdiagrams. For a given diagram D, we define the column weight of diagram D, s(D), as the number of different ways in which the connected subdiagrams D i c can be sequentially shrunk, so that all gluon attachments to Wilson lines in D are moved to the common origin. This means, in practice, that if all gluon attachments are entangled, so that no subdiagram can be shrunk without shrinking the whole diagram, then s(D) = 0. On the other hand, if, for example, a single subdiagram D i c can be shrunk without affecting the others, this provides a non-trivial sequence for the shrinking of the whole diagram, so that s(D) = 1: this is the situation for the two diagrams portrayed in figure 3. With this definition, it is conjectured that [52]. In what follows, row and column sum rules will be illustrated in a number of examples, and the conjectured column sum rule will be verified to hold for all four-loop webs connecting four or five Wilson lines.
We conclude this section by noting that the characterisation of web mixing matrices as projectors can be made explicit by introducing, for each web, a diagonalising matrix Y w , such that where p w is the number of diagrams for web w, and thus the dimension of the matrix R w , while r w < p w is the rank of R w . Without loss of generality, we have ordered the JHEP05(2020)128 eigenvalues of R w , labelled by λ i , so that null eigenvalues appear in the last positions. With these conventions, eq. (2.4) can be rewritten in matrix notation as where K is the vector of kinematic factors and C is the vector of colour factors for web w. It is clear that only r w ECFs will contribute to web w: the non-abelian exponentiation theorem tells us that they will be colour factors which, by the Feynman rules, would be associated to connected gluon subdiagrams.

From gluon webs to correlator webs
As discussed in section 2, webs for multi-parton scattering amplitudes are defined as sets of diagrams connecting gluons to Wilson lines, containing diagrams which are related to one another by permutations of gluon attachments to the Wilson lines. In this section, we present a definition for a closely related structure, where fixed-order Feynman diagrams are replaced by 'skeleton' diagrams, which are built out of connected gluon correlators, instead of gluon propagators and vertices. The basic reason to introduce these structures is that they strongly simplify the counting and organisation of contributions to the logarithm of Wilson-line correlators, especially at high orders, where radiative corrections to gluon subdiagrams become important and proliferate; furthermore, as we will see, using connected correlators does not affect the definition and structure of web mixing matrices, which are derived exclusively from the ordering of gluon attachments to the Wilson lines, and are not affected by gluon interactions away from the Wilson lines.
With this in mind, we define a correlator web, or Cweb as 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 fixed-order quantities, but admit their own perturbative expansion in powers of g. In a manner similar to webs, it is not easy to find a non-ambiguous notation to uniquely identify Cwebs: indeed, to some extent, they reflect the full complexity of the perturbative expansion. Below, we will use the notation W n (k 1 , . . . , k n ) for a Cweb involving n Wilson lines, 2 with k i gluon attachments to the i-th Wilson line. It is clear that beyond the lowest orders several different n-line Cwebs will share the same number of attachments to the Wilson lines: one can then refine the notation, using W (c 2 ,...,cp) n (k 1 , . . . , k n ) for a Cweb with the prescribed attachments, constructed out of c m m-point connected gluon correlators: one should keep in mind, however, that, at sufficiently high-orders, also this notation becomes ambiguous. We note also that there is an obvious degeneracy in the counting of Cwebs, since Cwebs that differ only by a permutation of their Wilson lines JHEP05(2020)128 2 (1, 1) are structurally identical, and it is trivial to include them in the calculation of the full correlator, simply by summing over Wilson-line labels. As a consequence, for classification purposes, we will break this degeneracy by picking a specific Wilson-line ordering, choosing Taking into account the fact that the perturbative expansion for an m-point connected gluon correlator starts at O(g m−2 ), we can write the perturbative expansion for a Cweb, with a prescribed correlator content and number of attachments, as which defines the perturbative coefficients of the Cweb, W n, j . Based on eq. (3.1), it is natural to classify Cwebs based on the perturbative order where they receive their lowestorder contribution, which is given by the power of g in the prefactor of eq. (3.1); one may then easily design a recursive procedure to construct all Cwebs order by order. We begin by noticing that only two Cwebs have lowest-order contributions at order g 2 : a self-energy insertion with a two-point connected gluon correlator attached to a single Wilson line, which we denote by W 1 (2), and the configuration with a two-gluon correlator joining two Wilson lines, which we denote by W In the massless case (with Wilson lines on the light cone), the self-energy Cweb vanishes identically, since, by the eikonal Feynman rules, it is proportional to the square of the Wilson-line four-velocity vector, β 2 , so one is left with a single non-vanishing O(g 2 ) Cweb. Starting from this initial condition, one may systematically construct Cwebs starting at higher perturbative orders: keeping in mind that one must imagine having an unlimited supply of yet-uncoupled Wilson lines, one may proceed by performing the following moves.
• Add a two-gluon connected correlator connecting any two Wilson lines (including Wilson lines that had no attachments at lower orders).
• Connect an existing m-point correlator to any Wilson line (again, including Wilson lines with no attachments at lower orders), turning it into an (m+1)-point correlator.
• Connect an existing m-point correlator to an existing n-point correlator, resulting in an (n + m)-point correlator. JHEP05(2020)128 3 (1, 2, 1) Executing all moves is clearly redundant, since the same Cweb is generated more than once through different sequences of moves: upon performing all moves, one must therefore remove multiply-counted Cwebs. This procedure can be considerably streamlined by taking into account known properties of webs, which naturally generalise to Cwebs. Specifically, • webs (and thus Cwebs) that are given by the product of two or more disconnected lower-order webs (so that there are subsets of Wilson lines not joined by any gluons) do not contribute to the logarithm of the correlator, W n ; • as mentioned above, in a massless theory all self-energy webs (Cwebs), where all gluon lines attach to the same Wilson line, vanish as a consequence of the eikonal Feynman rules. Thus, any Cweb containing a connected gluon correlator attaching to a single Wilson line will vanish.
It turns out that both these rules can be applied to trim the recursive procedure: more precisely, moves that lead to a disconnected Cweb, or (in the massless theory) to a selfenergy Cweb can be immediately discarded. This is less than obvious, since a disconnected Cweb can become connected at the next recursive step, and similarly a self-energy Cweb can become connected to other Wilson lines upon adding more gluons. It is however easy to convince oneself that all non-vanishing Cwebs that are reached by the recursion through intermediate stages involving either self-energy or disconnected Cwebs, are also reached through sequences of steps involving only non-vanishing Cwebs. The recursion can therefore be stopped whenever a vanishing contribution of these two kinds is reached. Using these recursion criteria, it is easy to enumerate inequivalent Cwebs at low orders. In the massless theory, we find four inequivalent Cwebs starting at O(g 4 ), which we label W (1, 1, 1) and W 3 (1, 2, 3) and W 3 (2, 2, 2). Notice that here we find the first occurrence of two Cwebs JHEP05(2020)128 Figure 6. Representative skeleton diagrams for the four non-vanishing Cwebs in a massless theory connecting two Wilson lines, and whose perturbative expansion starts at O(g 6 ). Cwebs (a) and (b) have a single skeleton diagram, while Cwebs (c) and (d) have six. with the same correlator content and attachments: they are distinguished by different attachments of the three-gluon correlator to the Wilson lines, and we tag them by different roman numeral indices. Finally, still at O(g 6 ) we find four Cwebs connecting four Wilson lines: 4 (1, 1, 2, 2) and W we find a total of 60 new Cwebs: 8 connecting two lines, 22 connecting three lines, 21 connecting 4 lines and 9 connecting five lines. Four-loop Cwebs connecting four and five lines will be discussed in more detail in section 5, and are listed in the appendix.
We conclude this section by discussing briefly the symmetry properties of Cwebs and the counting of skeleton diagrams entering each Cweb. First, we note that Cwebs are built out of connected gluon correlators, each of which is constrained by Bose symmetry and gauge invariance. In simple cases, this poses strong limitations on the colour structures entering the Cweb: for example, colour conservation forces the gluon two-point function to be diagonal in colour, and thus proportional to the unit matrix in the adjoint representation JHEP05(2020)128 4 (1, 3, 1, 1) of the gauge group; thus, to any order, a two-point function joining Wilson lines i and j will always generate the dipole structure T i · T j ≡ T a i T j, a . Similarly the three-point (off-shell) correlator is conjectured to be proportional to the structure constants f abc to all orders in perturbation theory, a conjecture which has been verified to high orders (see for example refs. [64,65]), and which implies the complete antisymmetry of the corresponding kinematic factor. For correlators with n > 3 gluons, several colour structures are possible (including building blocks for higher-order Casimir operators) and decoding the constraints imposed by colour conservation and Bose symmetry becomes more intricate: these constraints are however crucial for the construction of the full soft anomalous dimension matrix (see for example [36,46]), and they are efficiently summarised to all orders in the language of Cwebs.
Concerning the counting of skeleton diagrams contributing to a given Cweb, which gives the dimension of the corresponding web mixing matrix, we note that if a connected correlator is attached to a given Wilson line with p gluons, the permutations of those gluons should not be taken into account in the enumeration of contributing diagrams, since Bose symmetry is embedded in the structure of the correlator and all p gluons are equivalent. The dimension of the web mixing matrix should therefore be computed not by counting permutations of gluons on each Wilson line, but rather by counting shuffles of the subsets of gluons emerging from each correlator. Thus, for example, if a Wilson line has a total of five attachments, three of them emerging from one correlator, and the remaining two from another one, that line will contribute to the dimension of the web mixing matrix a factor of 10 (the number of possible shuffles of two sets of three and two cards), rather than a factor of 120 (the number of permutations of five cards). We note also that, in the language of Cwebs, the bulk of the growth of the number of diagrams at high orders in perturbation theory is hidden in the internal structure of the contributing gluon correlators: thus, for example, in our language, the four Feynman diagrams contributing to W (0,0,1) 4 (1, 1, 1, 1) at O(g 6 ) are understood as four contributions to that same Cweb at that order, rather than four distinct webs; similarly W (0,0,0,1) 5 (1, 1, 1, 1, 1) receives contributions from 25 Feynman diagrams at O(g 8 ).
We now move on to the discussion of the calculation of web mixing matrices with the method of replicas [51], before discussing the four-loop case in greater detail.

JHEP05(2020)128 4 A replica-trick algorithm to generate Cweb mixing matrices
As in other combinatorial problems involving exponentiation, such as the construction of effective actions, Wilson-line correlators can be studied efficiently by means of algorithms constructed with the replica method [66]. Here we will briefly discuss the application of the replica method to infrared exponentiation, following refs. [51,67], and then outline our algorithm for the construction of web mixing matrices.
In order to introduce the replica method, consider the path integral expression for the Wilson-line correlator in eq. (1.1) where S is the classical action. 3 In order to compute W n (γ i ), we may imagine building a replicated theory, replacing the single gluon field A a µ with N r identical copies A a, i µ (i = 1, . . . , N r ), which do not interact with each other. Under this replacement, one has if, furthermore, we replace each Wilson line in eq. (4.1) by the product of N r Wilson lines, each belonging to one replica of the theory, one readily realises that the replicated correlator is given by As a consequence, in order to compute W n (γ i ) order by order in perturbation theory, one may compute the replicated correlator, and then extract from the result the term of order N r . Importantly, while gluon fields belonging to different replicas do not interact, they all belong to the same gauge group: therefore, the colour matrices associated with their attachments to the Wilson lines do not commute, and their ordering is relevant. On the other hand, in a Cweb, each connected gluon correlator can be assigned a unique replica number, since there are no interaction vertices connecting different replicas. It is clear then that the contributions of different skeleton diagrams to the replicated correlator in eq. (4.2) will be simply related to those of the same diagrams in the unreplicated theory, by means of combinatorial factors, counting the multiplicities associated with the presence of different replicas. The computation of W n (γ i ) in terms of the skeleton diagrams contributing to S n (γ i ) is thus reduced to the computation of these combinatorial factors. The necessary steps, listed below, were identified in refs. [51,67], and can be naturally adapted to the language of Cwebs.
• Given a Cweb, assign a replica number i (1 ≤ i ≤ N r ) to each connected gluon correlator present in the web. Note that if only one connected correlator is present, only one replica can contribute to any diagram. It is then easy to show that all JHEP05(2020)128 diagrams of the Cweb contribute to W n (γ i ) with unit weight: in other words, there is no mixing matrix.
• Define a replica ordering operator R which acts by ordering the color generators on each Wilson line according to their replica numbers. Note that the Wilson lines are naturally oriented, so this operation is unambiguous. If we denote by T (i) k the group generator associated with the emission of a gluon in replica i from the k-th Wilson line, the replica-ordering operator acts on a product T k by relabelling the generators, exchanging their replica numbers i and j, if i > j, while leaving them unchanged if i ≤ j.
• In order to compute the colour factor of a skeleton diagram in the replicated theory, it is now sufficient to note that any non-trivial action by the replica ordering operator R effectively replaces the selected skeleton diagram with another one drawn from the same Cweb. The colour factor in the replicated theory is thus a linear combination of the colour factors of all skeleton diagrams in the Cweb, with multiplicities given by the number of possible replica orderings of the gluon attachments on every Wilson line.
• Algorithmically, for a Cweb W (c 2 ,...,cp) n (k 1 , . . . , k n ), built out of m = p r=2 c r connected correlators, one needs to determine two relevant numbers: the number of possible hierarchies between the m replica numbers of the correlators, h(m), and, for every fixed hierarchy h, the multiplicity with which that hierarchy can occur in the presence of N r replicas M Nr (h). The determination of h(m) is made non-trivial by the fact that the case of equal replica numbers must be treated separately: the sequence h(m) is however well known [68], and given by the so-called Fubini numbers. 4 In the first instances, for m = {1, 2, 3, 4, 5, 6} one finds h(m) = {1, 3, 13, 75, 541, 4683}. The multiplicity of a given hierarchy, on the other hand, is easily seen to be given by where n r (h) is the number of distinct replicas present for hierarchy h. To give a concrete example, for m = 5, labelling the replica numbers of the 5 correlators with i k , (k = 1, . . . , 5), and picking the hierarchy • Finally, the exponentiated color factors for a skeleton diagram D is given by The Fubini numbers admit a generating function, and they can be defined by

JHEP05(2020)128
where R C(D) h is the color factor of the skeleton diagram obtained from D through the action of the replica-ordering operator R for the case of hierarchy h. The Cweb mixing matrix is built out of the coefficients M Nr (h), which are polynomials in N r , picking terms that are linear in N r : many examples will be described in section 5 and in the appendix. Note that in the presence of m-point correlators, with m ≥ 4, each such correlator can contribute different 'internal' colour factors, for example different permutations of products of structure constants. Since, however, the information on the mixing matrix is encoded in the coefficients M Nr (h), the different colour factors arising from the internal structure of the correlators can be treated one by one, without affecting the mixing of the diagrams.
• We note that, given a Cweb W (c 2 ,...,cp) n (k 1 , . . . , k n ), the dimension d W of its mixing matrix R W can be expressed as follows. Let s W (k) be the number of shuffles that can be performed with the attachments to the k-th Wilson line, given the arrangement of connected correlators for the Cweb under consideration. One would naïvely conclude that the dimension of R W is the product of the factors s W (k) over all Wilson lines. This however fails to take into account an important degeneracy, which already appears in the simple two-line Cweb W (2) 2 (2, 2) at two loops: counting shuffles separately on each line yields d W = 4, which is wrong, because the two shuffles on the second Wilson line can be obtained from the shuffles on the first line by exchanging the two gluons, which is manifestly a symmetry of the Cweb. In order to take into account this degeneracy, one must divide by the number of available permutations of subsets of m-point correlators that have the property of being attached to the same sets of m Wilson lines.
In order to compute Cweb mixing matrices at four loops, extending the results refs. [51,52,59], we have developed an in-house Mathematica code which we describe very briefly below.
• The first step is to generate all the Cwebs that appear at four loops (O(g 8 )), in particular those involving four and five Wilson lines. To do so, we note that at four loops (O(g 8 )), all possible Cwebs can be obtained by combining the connected correlators shown in figure 9, and attaching them to 2 ≤ n ≤ 5 Wilson lines. One may begin by attaching four two-point correlators in all possible ways to the Wilson lines. Next one proceeds to Cwebs generated by attaching two two-point correlators and one three-point correlator, in all possible ways, and similary with the other types, obtained by using the other building blocks in figure 9. We note that the five-point correlator at this order can only give a trivial Cweb, since it contains only a single skeleton diagram.
• The code assigns a distinct replica variable to each of the correlators present in a given Cweb: for example, four two-gluon correlators will be assigned replica indices i, j, k and l. Then the correlators are sequentially attached to the Wilson lines in lexicographic order, beginning with the one with index i, attached between Wilson line 1 and 2, proceeding to the one with index j, attached in all possible ways, so as to generate a set of 'partial skeleton diagrams'. This procedure continues till all JHEP05(2020)128 correlators are exhausted. Clearly, this procedure will generate several Cwebs that are identical (as far as color is concerned), since they are related to each other by a mere renaming of the Wilson lines: duplicates are discarded. At this stage, we have only one diagram per Cweb.
• The code then takes one of the above diagrams and generates all the other diagrams of the corresponding Cweb by permuting (or more precisely shuffling) the gluon attachments on each of the Wilson lines.
• The next step is to generate all the multiplicities associated to the possible hierachies, corresponding to the entries of the table 1 of ref. [51]. A subroutine creates all possible different hierarchies h for each diagram; the code then proceeds to obtain all other colums of the table, and finally leads to the exponentiated colour factorsC(D), from which we obtain the mixing matrix R. Finally, R is diagonalised and the diagonalising matrix Y is recorded.
• The code automatically discards self-energy Cwebs, where a connected correlator is attached only to a single Wilson line, but, in principle, keeps disconnected Cwebs, so that the vanishing of the corresponding R matrix works as a test.
We note that the run time of the code increases steeply with the increase in the number of connected correlators, which causes a rapid increase in the number of available replica hierarchies. The code has been checked by reproducing all lower-order results available in the literature, and by verifying the two known properties of mixing matrices: their idempotence, and the row sum rule, which hold true for all the R matrices that we have computed. Furthermore, one can verify that all different Cwebs at O(g 8 ) have been constructed, by applying the recursive construction described in section 3. Finally, as shown explicitly in section 5 and in the appendix, all R matrices we have computed satisfy the conjectured column sum rule discussed in section 2.

A selection of four-loop Cwebs
In this section we present in some detail the calculation of two four-loop Cwebs, involving respectively four and five Wilson lines. This will allow us to introduce some notation which will be useful to simplify the full listing of results for similar webs, which is presented in the appendix. For each Cweb, we will present the mixing matrix R, the diagonalising matrix Y , and the exponentiated colour factors.

A four-loop, four-line Cweb
As an example, we have selected the Cweb W (1, 1, 2, 2), which comprises four skeleton diagrams depicted in figure 10. We note that in this case the Cweb label includes a roman numeral, to distinguish it from a second Cweb with identical correlator and attachment content, W  1, 2, 2), where however the four-gluon correlator has two attachments to the same Wilson line. That Cweb, discussed in the appendix, involves only two skeleton diagrams, and therefore has a 2 × 2 mixing matrix. In this case, it is evident by looking JHEP05(2020)128 Figure 9. Combinations of connected correlators that can form Cwebs at four loops.
at any one of the diagrams that there are four shuffles to be considered, so that the mixing matrix will have dimension four. In order to write down an explicit expression for the matrices R and Y , it is necessary to introduce an ordering among the diagrams: in this case, the ordering is displayed in figure 10, but in general it can be identified by labelling the gluon attachments to be shuffled, and giving the sequence of the shuffles to be considered. As an example, in figure  We find that the mixing matrix R and the diagonalising matrix Y are given by The expected properties of the mixing matrix are easily verified: R is idempotent, i.e. R 2 = R, the matrix elements in each row sum to zero, and furthermore the column sum rule is obeyed. Indeed, in this case the vector built out of the indices s(D) for the four diagrams in figure 10 is given by S ≡ {s(C 1 ), s(C 2 ), s(C 3 ), s(C 4 )} = {1, 0, 0, 1}: in the first and last diagrams one can move the gluon attachments of one of the two correlators to the origin without affecting the other correlator, which is not possible for the second and the third diagram. One then readily verifies that the column vector S · R is a null vector.
We observe that these exponentiated color factors correspond to fully connected Feynman diagrams, so that the non-abelian exponentiation theorem of ref. [59] is, as expected, satisfied.

A four-loop, five-line Cweb
As an explicit example of a four-loop, five-line Cweb, we select the one labelled by W (2,1) 5, I (1, 1, 1, 2, 2), one of two Cwebs with this particular set of attachments and correlators. A sample skeleton diagram contributing this Cweb is shown in figure 11: it is immediate to note that there are four possible shuffles of the labels, so that the mixing matrix will have dimension four. One observes that both the row and the column sum rules are satisfied. Furthermore, diagonalising the mixing matrix, one finds that it has rank r = 1, so that there is only one exponentiated colour factor. It is

JHEP05(2020)128
In the appendix, we will similarly treat all the other four-loop Cwebs connecting four and five Wilson lines.

A note on colour structures at four loops
Four-loop colour structures are especially interesting for the study of gauge-theory scattering amplitudes, since at this order for the first time the soft anomalous dimension matrix can receive contributions from quartic Casimir operators of the gauge algebra. The relevance of higher-order Casimir invariants was first noted in this context in [15], but it was in fact known from the early days of QCD [69]. The presence of quartic Casimir contributions at four loops in the cusp anomalous dimension is one of only two possible sources for the violation of the dipole structure of the soft anomalous dimension matrix for massless theories, which holds up to two loops. Remarkably, the cusp anomalous dimension for QCD was JHEP05(2020)128 recently computed analytically at four loops [70,71], following the precise numerical predictions of refs. [72,73], and the presence of a non-vanishing contribution from quartic Casimir operators was confirmed. Studying the implications of such contributions for the multiparticle soft anomalous dimension matrix is therefore now an open and very interesting problem, also in light of the results of refs. [74,75]. Of course, the ultimate understanding of the role played by these contributions, and their implications for collinear factorisation, will depend upon a full calculation of the kinematic factors for the relevant Cwebs: indeed, we note that collinear factorisation remains compatible with the three-loop expression of the soft anomalous dimension matrix only thanks to a remarkable connection between matrices with different number of partons, and after enforcing the constraints of colour conservation [36]. We may, in any case, make a few remarks in this issue, based purely on the analysis of the colour structure of four-loops Cwebs. First, we observe that all the exponentiated colour structures that arise in four-loop Cwebs connecting four and five Wilson lines (listed in the appendix) correspond to the connected gluon diagrams depicted in figure 12, with the open ends of the gluon lines attaching to generic permutations of Wilson lines. The gluon-loop structure on the right of figure 12 is of course built out of products of structure constants f abc , however, upon symmetrization, could in principle yield a quartic Casimir contribution. We note, however, that, for all Cwebs whose perturbative expansion starts at O(g 8 ), the gluon-loop structure cancels in the exponentiated colour factors (Y C) i , before their recombination with kinematic factors. It would appear that the only possible source of quartic Casimir contributions from four-and five-leg Cwebs at four loops is in the oneloop correction to the O(g 6 ) Cweb W (0,0,1) 4 (1, 1, 1, 1), depicted in figure 8(a). Indeed, at one loop, the four-gluon correlator featuring in that Cweb (which must be Bose symmetric), can develop a symmetric colour structure, yielding a contribution of the form where K is a kinematic factor, and d (r) abcd is the completely symmetrized trace of four generators of the gauge algebra in representation r. We note, however, that the colour structures directly arising from the diagrams are not all independent, and, upon enforcing colour conservation, which implies the operator constraint i T i = 0, they must be reduced to a basis, following for example the analysis of ref. [59]. Since antisymmetric combination of generators on the same Wilson line can be eliminated using the gauge algebra, this step JHEP05(2020)128 leads to the appearance of symmetric combination of generators, which can again yield contributions of quartic Casimir type. These results and considerations are in agreement with the analysis of refs. [55,56], with the effective vertex analysis in ref. [59], and with the generating functional approach of ref. [76], where the role played by symmetric colour structures in the exponentiated soft function is emphasised.

Summary and outlook
The study of the infrared structure of perturbative gauge amplitudes is an active developing field, with implications for both theoretical studies, concerning the mathematical properties of gauge theories, and for high-energy phenomenology at colliders. These studies are very advanced, and we now have a well-developed understanding of infrared factorisation and exponentiation, with the soft anomalous dimension matrix fully known at three loops for massless particles [36], the cusp anomalous dimension computed to four loops [70,71], and several high-order radiative amplitudes available (see for example [74,75]). The current frontier is the calculation of the soft anomalous dimension matrix for multiparticle scattering at three loops including massive particles, and the exploration of the four-loop domain; in general, all these studies bring to the fore the relevance, for gauge theory calculations, of correlators involving Wilson lines, possibly together with gauge and matter fields, which provide leading-power approximations to scattering amplitudes and cross sections in soft and collinear limits.
In the present paper, we have developed a set of tools for the analysis of soft anomalous dimensions at high orders, and we have studied the exponentiated colour structures arising at the four-loop level in multiparticle amplitudes. We have introduced the concept of a correlator web, or Cweb, which, we believe, will be useful for the classification and study of exponentiated correlators at high orders: Cwebs include their own radiative correction, are easily generated and enumerated, since their number grows only moderately as a function of the perturbative order, and may help to clarify and implement symmetry properties of gluon correlators and the consequences of colour conservation.
We have enumerated all Cwebs for a massless theory up to four loops, and we have computed their mixing matrices and exponentiated colour factors for all cases involving four and five Wilson lines. We observe that all exponentiated colour factors correspond to completely connected gluon diagrams, as expected from the non-abelian exponentiation theorem [59]. Furthermore, we note that structures compatible with the presence of quartic Casimir contributions are present, as expected, but they cancel in Cwebs arising at O(g 8 ), while they survive in radiative corrections to Cwebs arising at lower orders. Finally, we verify the properties of mixing matrices that were proved or conjectured in earlier studies [51,52,[57][58][59][60].
The combinatorial complexity of the calculation of web mixing matrices grows rapidly with the perturbative order, most notably because of the proliferation of hierarchies that arise in the application of the replica method, which are counted by the Fubini numbers. We have developed a code which can handle this complexity at four loops with minimal computing resources, but we believe that going to yet higher orders is likely to require JHEP05(2020)128 significant optimisation, the development of new tools, or the deployment of considerably larger computing power. That not withstanding, our results in this paper provide a number of needed ingredients for the calculation of infrared divergences at four loops, and we believe that the tools that we have developed will be useful to further our understanding at higher orders as well.

Acknowledgments
We thank Einan Gardi and Niamh Maher for useful discussions concerning the choice of colour bases and the role of colour conservation. AT and LM 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. AT and SP would also like to thank the University of Turin and INFN Turin for warm hospitality during the course of this work. AD would like to thank CSIR, Govt. of India for a JRF fellowship. SP would also like to thank MHRD Govt. of India for an SRF fellowship. LM would like to thank IIT Hyderabad for their warm hospitality, and acknowledges the contribution of the Italian Ministry of University and Research (MIUR), grant PRIN 2017LNEEZ.

A Mixing matrices and exponentiated colour factors
In this appendix we give results for all the webs that appear at 4 loops in the scattering amplitude, that can connect 4 or 5 Wilson lines. Throughout the list, R, Y and D denote the mixing matrix, the diagonalizing matrix and the diagonal matrix respectively. D is represented as D = (1 r , 0), where r is the rank of the mixing matrix R. 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. Finally, we give the expressions for the exponentiated colour factors, which, in all cases, correspond to fully connected gluon diagrams, as expected. For Cwebs involving four-point correlators, the colour factors that we present correspond to one of three possible permutations of structure constants arising from the internal structure of the correlator, as in section 5.1. We omit from the list the Cwebs that are composed of a single skeleton diagram, such as W (0,0,0,1) 5 (1, 1, 1, 1, 1), whose mixing matrix is just a number, R = 1. 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 S-factors

JHEP05(2020)128
The R, Y and D matrices are given by Finally, the exponentiated colour factors are

W
(1,0,1) 4,II (1, 1, 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 S-factors
The R, Y and D matrices are given by Finally, the only exponentiated colour factor is 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.

JHEP05(2020)128
The R, Y and D matrices are given by Finally, the exponentiated colour factors are 4,I (1, 1, 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 S-factors
The R, Y and D matrices are given by Finally, the only exponentiated colour factor is This is a three-diagram Cweb, represented by Diagrams Sequences S-factors The R, Y and D matrices are given by

JHEP05(2020)128
Finally, the exponentiated colour factors are 1, 1, 4) This is a more complicated Cweb, containing twelve diagrams. We find Diagrams Sequences S-factors The R, Y and D matrices are given by

JHEP05(2020)128
Therefore, there are six exponentiated colour structures, given by This is a six-diagram Cweb, the first of three with same correlator and attachment content, for which we find Diagrams Sequences S-factors The R, Y and D matrices are given by This yields two colour structures, (A.14) JHEP05(2020)128 8. W (2,1) 4, II (1, 1, 2, 3) Another six-diagram Cweb, the second of three with same correlator and attachment content, for which we find Diagrams Sequences S-factors The R, Y and D matrices are given by The exponentiated colour factor are (A.16)

W
(2,1) 4, III (1, 1, 2, 3) Yet another six-diagram Cweb, the third of three with same correlator and attachment content. We find Diagrams Sequences S-factors The R, Y and D matrices are given by

JHEP05(2020)128
which leads to the colour factors

W
(2,1) 4, I (1, 2, 2, 2) This is the first of five Cwebs with the same correlator and attachment content. Its four diagrams are

Diagrams
Sequences S-factors The R, Y and D matrices are given by There is therefore only one colour structure, (A.20)

W
(2,1) 4, II (1, 2, 2, 2) This is the second of five Cwebs with the same correlator and attachment content. It has eight diagrams, which can be organised as follows.

Diagrams
Sequences S-factors The R, Y and D matrices are given by The three colour factors are

Diagrams
Sequences S-factors The R, Y and D matrices are given by

JHEP05(2020)128
The four colour structures are 13. W The fourth Cweb of this set has twelve diagrams.

Diagrams Sequences
S-factors The R, Y and D matrices are given by There are thus six colour structures 1, 2, 3) The last Cweb of this set has six diagrams.

Diagrams
Sequences S-factors The R, Y and D matrices are given by

JHEP05(2020)128
The three colour structures are 4 (2,2,2,2) This Cweb, comprising sixteen diagrams, generalises a memebr of an infinite series of highly symmetrical 'Multiple Gluon Exchange Webs', studied in ref. [35]. Explicit calculations confirm the structure of the mixing matrix predicted there. We find The R, Y and D matrices are given by and the five colour structure are 4 (1, 1, 2, 4) We now come to one of the largest Cwebs of the set, with twenty-four diagrams.
The R, Y and D matrices are given by The R, Y and D matrices are given by   The R, Y and D matrices are given by

JHEP05(2020)128
The R, Y and D matrices are given by   and there is a single colour factor, 3. W (2,1) 5 (1, 1, 1, 1, 3) A six-diagram Cweb, Diagrams Sequences S-factors where the R, Y and D matrices are given by and there are two colour structures, The single colour factor is (A.48)

Diagrams
Sequences S-factors The R, Y and D matrices are given by The single colour factor is (A.50) 6. W  While large, this Cweb has only one exponentiated colour factor,  Diagrams Sequences S-factor

JHEP05(2020)128
with R, Y and D matrices given by