The role of colour flows in matrix element computations and Monte Carlo simulations

We discuss how colour flows can be used to simplify the computation of matrix elements, and in the context of parton shower Monte Carlos with accuracy beyond leading-colour. We show that, by systematically employing them, the results for tree-level matrix elements and their soft limits can be given in a closed form that does not require any colour algebra. The colour flows that we define are a natural generalization of those exploited by existing Monte Carlos; we construct their representations in terms of different but conceptually equivalent quantities, namely colour loops and dipole graphs, and examine how these objects may help to extend the accuracy of Monte Carlos through the inclusion of subleading-colour effects. We show how the results that we obtain can be used, with trivial modifications, in the context of QCD+QED simulations, since we are able to put the gluon and photon soft-radiation patterns on the same footing. We also comment on some peculiar properties of gluon-only colour flows, and their relationships with established results in the mathematics of permutations.


Introduction
Barring an unlikely turn of events, searches for new physics at the LHC will increasingly rely upon high-precision methods, on both the experimental (with the collection of largestatistics samples and the reduction of systematics) and the theoretical sides. An obvious benefit of high-precision results is that they also allow one to carry out progressively more stringent checks on the (supposedly) well-known Standard Model physics, which in turn help stress-test, and ultimately hone, theory predictions.
From the theoretical viewpoint, meeting the precision targets that current and future experimental data demand almost always requires performing computations that are perturbative in the series expansion in terms of a coupling constant (the so-called fixed-order predictions) at the highest possible accuracy. At the LHC, one is obviously dominated by QCD effects, but electroweak contributions must nowadays be considered as well, especially in the large transverse momentum regions where some characteristic new-physics phenomena are expected to show up. In the past few years, the progress in fixed-order technologies has been extremely significant, for example with next-to-next-to-leading order (NNLO) results becoming routinely available. That being said, fixed-order predictions cannot generally be directly compared to experimental data: even when they are fully differential (a feature that is more difficult to achieve the higher the perturbative accuracy of the computation), hadronization and some acceptance corrections generally need be applied.
It is therefore desirable that such corrections also be carried out at the highest possible accuracy, so as not to spoil that of the underlying fixed-order calculations. This operation more often than not involves the matching and/or merging of fixed-order results with parton-shower Monte Carlo (PSMC) simulations. Unfortunately, the systematic application of matching and merging techniques is currently limited to next-to-leading order (NLO) at most. Furthermore, the PSMCs employed in the context of these techniques are formally of leading-logarithmic (LL) and leading-colour (LC) accuracy, although in practice they perform (in comparison to analytically-resummed computations) better than that.
Another issue that emerges when one considers large-statistics data samples is that many-jet observables become accessible; loosely speaking, this may also be viewed theoretically as a high-precision problem, although of a different nature from those posed by higher-order computations. In particular, here one must be able to calculate the matrix JHEP11(2021)045 elements 1 associated with a large number of partons in an efficient way. The complexity of this task grows rapidly, and sooner or later one must consider approximated, as opposed to exact, matrix-element results.
The primary aim of this paper is to illustrate how several of the problems outlined thus far can be addressed by using colour flows. This applies to both tree-level matrix element computations, and to the inclusion of beyond-LC effects in PSMCs. We argue that by employing colour flows at both the matrix-element and the PSMC level one has the advantage of a language common to the two cases, which should facilitate the interplay between them in the attempt at increasing the overall precision of the final predictions. Here, we focus mainly on the use of colour flows in tree-level matrix elements to provide a full-colour treatment of the initial conditions for parton showers. The same approach can in principle be applied in the subsequent evolution of the showers. We remark that the extension of PSMCs beyond leading (logarithmic and/or colour) accuracy is the subject of much ongoing activity (see e.g. refs. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]) in directions that are possibly orthogonal to the methods discussed here; a comparison between all these developments is beyond the scope of the present work.
We point out that while the concept of colour flow is well known, our work exploits it at the level of matrix elements (i.e. of amplitudes squared), in contrast to what is typically done in the literature (see e.g. refs. [17][18][19][20][21][22] and references therein), where the emphasis is at the amplitude level. By doing so, we are naturally led to introduce quantities (such as closed flows and secondary flows) that help to express the matrix elements in a very compact way that does not require any further colour algebra, and is independent of whether the original amplitudes are written in terms of a fundamental or a colour-flow basis. An interesting question that thus emerges is whether similar expressions could be obtained for matrix elements beyond the tree-level (which, at one loop, would be a natural complement of the real-emission subtraction procedure presented in ref. [23]). While the general ideas on colour flows we employ here are certainly valid beyond tree-level, there are specific aspects of the representations of loop amplitudes in terms of colour bases (see e.g. refs. [24][25][26][27][28][29]) that require an explicit analysis, and whose implications are therefore left to future work. This paper is organized as follows. In view of the fact that we often rely on the methods of ref. [23], in section 2 we give a brief summary of the results of that paper that are directly relevant here. In section 3, we define the colour loops that stem from a given choice of colour flows. In section 4 we discuss the so-called flow representation of matrix elements in geometrical terms. Sections 5 and 6 present the results for tree-level matrix elements and their soft limits, respectively, expressed in terms of elementary quantities associated with colour flows. In section 7 we show how the results obtained for the matrix elements could be employed in the context of PSMCs that include beyond-LC effects, and we propose two different but equivalent ideas that generalize the concept of colour connections used in existing PSMCs. In section 8 we illustrate how colour flows provide a natural language to deal with the simultaneous perturbative expansion in two parameters, JHEP11(2021)045 using the QCD+QED case as a definite example. Finally, in section 9 we summarize our findings. The appendices collect some technical material. In appendix A miscellaneous results are presented that are relevant to gluon-only colour flows, with the reinterpretation of some of their properties in the context of the mathematics of permutations, whereas appendix B gives some additional information on an idea introduced in section 7 that is particularly relevant to dipole-based showers.

Summary of basic concepts
The present work makes extensive use of the results of ref. [23]. For the reader who is not familiar with that paper, in this section we give a very brief summary of the ideas that underpin the derivations that will follow later. Henceforth, eq. (x.y) of ref. [23] will be denoted by eq. (I.x.y).
All partons 2 that participate in a hard scattering are regarded as outgoing. For a process with q quark-antiquark pairs and n gluons 0 −→ 2q + n (2.1) the labelling conventions for the i th parton are as follows: Tree-level scattering amplitudes are written as sums of products of dual (or colour-ordered) amplitudes times colour factors. For any given choice of the representation of the SU(N ) colour algebra, colour factors are in one-to-one correspondence with colour flows, which are ordered sets of parton labels, and can be obtained from the latter by means of simple rules. For the process of eq. (2.1), the colour flows have the following form:

JHEP11(2021)045
The set of m-integer permutations is denoted by P m ; the set of the partitions of eqs. (2.7) and (2.8) is denoted by T n|q . Thus, the set of colour flows relevant to the process of eq. (2.1) is given by eq. (I. 3.27) F 2q;n = P n , P q , T n|q , (2.10) whose number of elements is equal to |F 2q;n | = |P n | |P q | T n|q = n! q! n + q − 1 q − 1 = q (n + q − 1)! . (2.11) Each of the ordered sets in γ p in eq. (2.5) is called an open colour line, that connects quark −p with antiquark µ(−p − q) and features the emission of gluon σ(t p−1 + 1) (which is the closest to the quark) up to gluon σ(t p ) (which is the closest to the antiquark). Denoting by the colours of gluons, and of quarks and antiquarks, respectively, the colour flow of eq. (2.5) is associated with the following colour factor in the fundamental representation 5 The integer ρ(γ) is equal to the number of open colour lines that coincide with flavour lines, 6 minus one in the case of maximal coincidence. Given a colour configuration, i.e. a specific choice for the colours of quarks, antiquarks, and gluons, the scattering amplitude is written as was anticipated, i.e.

JHEP11(2021)045
denoted by S {a} . The set of all colour configurations corresponds to a complete basis of orthonormal vectors in S {a} , thus 7 |a −2q , . . . a n ∈ S {a} , (2.17) b −2q , . . . b n |a −2q , . . . a n = n i=−2q |a −2q , . . . a n a −2q , . . . a n | = I . (2.19) The colour space is separable into the subspaces S The dimensions of these vector spaces are With this, one can introduce vectors associated with amplitudes at given colour configurations or colour flows, thus |A (2q;n) (a −2q , . . . a n ) = A (2q;n) (a −2q , . . . a n ) |a −2q , . . . a n , The r.h.s. of eq. (2.29) is expressed in terms of scalar quantities, namely the dual amplitudes already introduced, and the colour-flow matrix element . (2.30) Equations (2.27)- (2.30) give us the opportunity to stress a fact left implicit so far. Namely, a colour flow is a quantity defined at the amplitude level. In the computation of cross sections, what is relevant is always a pair of two flows, one on each side of the cut. Such a pair, i.e.
γ , γ (2.31) in the expressions above, is called a closed colour flow. In eq. (2.31), it is understood that γ (γ) is the flow entering the amplitude on the left-hand side (right-hand side) of the cut. Because of this, when it is necessary to distinguish between the colour flows on either side of the cut they will be called L-flows and R-flows. The formulae written so far are relevant to the process of eq. (2.1). Although in principle that could be used to describe purely gluonic processes as well (i.e. when q = 0), in practice such a case entails significant simplifications that justify the use of a slightly different notation. In particular, for a purely-gluonic n-body process 0 −→ n , (2.32) the analogue of the colour flow of eq. (2.5) features a single colour line, denoted as follows: σ ≡ (σ(1), . . . σ(n)) . (2.33) Furthermore, in view of the fact that its associated colour factor, i.e. the analogue of eq. (2.14), reads in the fundamental representation the information embedded in the colour flow must be invariant under cyclic permutations. Therefore, the set of the n-gluon colour flows is the set of the non-cyclic n-integer permutations, which we shall denote by and which constitutes the analogue of the set introduced in eq. (2.10). These, and the formal replacements 37) are all that is needed in order to convert the expressions for the quark-gluon matrix elements and related quantities into those relevant to the gluon-only case.

Colour loops
In the context of current parton shower Monte Carlos, colour flows stemming from matrix element computations are responsible for giving crucial information on the initial conditions for the showers (and for hadronization); therefore, ultimately they have a major role in the determination of the radiation pattern. In retrospect, this might appear strange for two reasons. Firstly, as was discussed in section 2, the quantities that are relevant to matrix elements are not flows, but closed flows. Secondly, such flows are Born-level (i.e. m-body) objects, while even the firstemission radiation pattern is determined, in cross section computations, by (m + 1)-body objects. As far as the first item is concerned, the key point is that existing PSMCs only consider closed flows for which the L-flow coincides with the R-flow; in other words, only the diagonal of the colour-flow matrix is taken into account, which is consistent with the leading-N nature of the MCs. Obviously, then, the information on the closed flow is redundant, and that of the flow is used instead. In turn, this immediately suggests that if one is interested in going beyond the leading-N approximation in MCs, retaining the information on closed flows (e.g. in determining colour connections) will be essential. As far as the second observation above is concerned, this has again to do with the leading-N character of PSMCs, which is such that the radiation pattern at the (m + 1)-body level can be described by m-body quantities. Beyond the leading-N this is not true any longer. However, we shall show how one can circumvent this difficulty, and thus to keep on working with m-body objects.
Before proceeding, we need to define precisely colour loops, which will be used throughout this paper. Firstly, colour loops are properties of a given closed colour flow. The set of all colour loops relevant to the closed flow (γ , γ) will be denoted by L(γ , γ) . (3.1) Each loop is defined iteratively, as follows: 0. Write the L-and R-flow as two vertical sets of indices (the former to the left of the latter). Each open colour line that belongs to these flows is a separate subset of such sets. The relative order of these subsets is irrelevant.

1.
Start from either a quark or a gluon that belongs to the R-flow; this is the first element of the loop.

2.
Go down one place in the R-flow; a gluon or an antiquark is found.

3.
Jump to the L-flow, landing on the element whose label is identical to that of the gluon or the antiquark of the previous step.

4.
Go up one place in the L-flow; a gluon or a quark is found.

5.
Jump to the R-flow, landing on the element whose label is identical to that of the gluon or the quark of the previous step.

JHEP11(2021)045
6. When the landing element coincides with the first element of the loop then go to step 7, otherwise iterate steps 2-5.

7.
A loop is the ordered set composed of the starting element and of all of the landing elements subsequently found, except for the last one (because it coincides with the starting element).
Once a loop has been constructed, the next one (if it exists) can be found by following again steps 1-7, with the condition that the starting element in step 1 must not be already included as an R-element of one of the loops previously constructed. We call L-and Relements of a given loop ∈ L the landing elements when jumping from right to left and from left to right, respectively; on top of those, the starting element of the loop is included in the set of R-elements as well. By construction, colour loops have the following properties: i. The sets of R-elements of any two loops are disjoint (i.e. have an empty intersection); the same is true for the L-elements.
ii. All quarks are R-elements, and all antiquarks are L-elements.
iii. A quark label will appear only once, i.e. it belongs to a single loop; the same is true for antiquarks.

iv.
A gluon label will appear twice, either in the same loop or in two different loops.
In the former case, one instance corresponds to an R-element, and the other to an L-element (i.e. the gluon label cannot correspond to two L-or two R-elements).
These properties will turn out to be important in deriving results for the soft limit of matrix elements. We point out that the idea that underpins the construction of colour loops as given above is based on the usual procedure of following colour (as opposed to anticolour) indices. 8 This is straightforward in the case of quarks and antiquarks; for gluons, one needs to consider the pair of λ matrices associated with any gluon index, and to replace them with the first term of the Fierz identity. Roughly speaking, this corresponds to "splitting" colour lines by replacing each gluon there with a qq pair (of a fictitious flavour, different from all of the other flavours); this idea will play an important role starting from section 5.2. We can also introduce the idea of distance within a colour loop. More precisely, for any two parton labels k and l that belong to a loop k , l ∈ ∈ L(γ , γ) , (3.2) the distance between k and l, denoted by is defined as the number of partons in between k and l; thus, if k and l are contiguous, the distance is equal to zero. Furthermore, according to item ii. above and taking into account Figure 1. Colour loop sets for 0 → qqgg: L(γ 1 , γ 1 ) (left) and L(γ 1 , γ 2 ) (right).
that colour loops alternate L-and R-elements, the distance between any quark and any antiquark is an even number, whereas the distance between two quarks or two antiquarks is an odd number.

JHEP11(2021)045
the number of colour loops of the closed flow (γ , γ); this quantity will appear throughout this paper. The construction of colour loops works also in the case of gluon-only amplitudes. Given a closed flow (σ , σ), steps 0-7 previously described are unchanged (with "quark" or "antiquark" there to be ignored), except for the necessity of including a cyclicity condition, which is simply achieved in the following way: • Going down one place from the last element in the R-flow means arriving at the first element of that flow.
• Going up one place from the first element in the L-flow means arriving at the last element of that flow.

Extended colour space, and the flow representation
The colour-anticolour structure of gluons suggests the possibility of working in a vector space larger than that introduced in eq. (2.20). Let p be any gluon; we consider an The relationship between S {a} and S {ij} is established as follows. Let be an orthonormal basis of S The vectors |a that constitute an orthonormal basis of S (p) {a} are then Therefore (see footnote 5 for the normalization) which shows that eq. (4.4) behaves as expected. Furthermore

JHEP11(2021)045
Equation (4.6) suggests defining the vector so that, given eq. (4.3), from eq. (4.6) one obtains It also follows that Therefore, P (p) and P  {ij} is in one-to-one correspondence with what is usually called the flow representation for amplitudes, which in the present context will acquire a natural geometrical meaning.

Tree-level matrix elements
In this section, we present results for tree-level matrix elements, regarded as quantities where all partons are hard and well separated from each other.

Quark-only matrix elements
The processes of interest are those of eq. (2.1), with n = 0 there. Because of this, the colour factors are simply products of Kronecker δ's (see footnote 5) which, when employed in the computation of the colour-flow matrix element of eq. (2.30), are contracted with each other, alternating one Kronecker δ from the R-flow with one Kronecker δ from the L-flow. Eventually, one finds a δ whose right index is equal to the left index of the δ from which the contraction had started, thus resulting in the trace of the identity matrix, equal to N . It is easy to convince oneself that the indices encountered when contracting the δ's are nothing but the colour labels of the partons that enter the colour loop to which any of these partons belong (bear in mind that any quark or antiquark appears in a single loop). Therefore, each loop is responsible for a contribution equal to N to the colour-flow matrix element which implies, from eqs. (2.14) and (2.29)

Quark-gluon matrix elements
The remarkable simplicity of eq. (5.1) stems from the trivial nature of the colour factor relevant to quark-only processes. The presence of gluons complicates the picture; however, we shall show in this section that results for quark-gluon matrix elements can essentially be cast in the same form as eq. (5.1), by introducing quantities that we shall call secondary colour flows, and that can be derived from the given colour flows. The starting point is the vector |A (2q;n) (γ) ∈ S {a} that represents the amplitude for a given flow γ, eq. (2.24). Consider one of the gluons that enter the process; it is not restrictive to assume its label to be 1. There exists a quark with label −p 1 such that In other words, gluon 1 belongs to the colour line that begins at quark −p 1 and ends at antiquark µ(−p 1 − q). Write the colour factor associated with γ as follows: with Λ (A) and Λ (B) suitable products of λ matrices; the definition ofΛ(γ) can be derived from eq. (2.14), but will not be needed here. Equation (2.24) can then be re-written more compactly as follows: and we have used eq. (4.4), understanding |ij ∈ S (1) {ij} . By means of the Fierz identity one obtains having used eq. (4.11).

JHEP11(2021)045
The vector on the l.h.s. of eq. (5.7) belongs to S (1) {a} . It is written in the r.h.s. as the difference of two vectors, the second of which belongs to S (1) {ij} \ S (1) {a} . In other words, a physical vector is written as the difference of two unphysical vectors, constructed so that their unphysical components cancel each other. The colour factors associated with these two vectors are also interesting. In that of the first term on the rightmost side of eq. (5.7) the gluon is replaced by a QQ pair with colour and anticolour indices equal to i and j, respectively, whereas in that of the second term the gluon disappears altogether. Thus, these colour factors are those one would obtain by employing eq. (2.14) with flows 9 potentially up to an overall power of N , which can however be easily fixed. Namely, the original colour flow γ has a pre-factor equal to N −ρ(γ) , which here is part of the definition of |v in eq. (5.6). For ρ(γ) to be equal to ρ(γ a ), it is sufficient that the flavour of the fictitious pair QQ be different from all of the other quark flavours relevant to the current process, which is a condition that we can impose by construction. Conversely, since the flavour structure ofγ b is identical to that of γ, then ρ(γ) = ρ(γ b ). This however does not take into account the factor −1/N that appears in eq. (5.7); indeed, it will turn out to be convenient (in particular, in view of its negative sign) to account for that factor separately from ρ(γ).
The procedure that leads one from eq. (5.5) to eq. (5.7) can be iterated, by considering all of the other gluons one after another. At each step, the number of contributions will grow by a factor of two (two being the number of terms on the r.h.s. of the Fierz identity), with each gluon either replaced by a fictitious quark-antiquark pair, or eliminated; with n gluons, the overall number of contributions will be thus equal to 2 n . In order to enumerate them, and to construct the analogues of the flows in eqs. (5.8) and (5.9), we introduce the following notation. Let Note that the definition of S The set S (n) k encompasses all possible choices of k gluon labels out of those of the given n gluons; an element s k ∈ S (n) k is one specific such choice. Then, for a given colour flow γ consider the 2 n -element set n k=0 s k ∈S & & s k element is called a secondary colour flow, and is constructed by taking γ and removing from it the k gluons associated with the labels in s k . As far as the other n − k gluons are concerned, eq. (5.8) tells us that they need be replaced by fictitious quarkantiquark pairs. Since such pairs are essentially a bookkeeping device, it is easier to continue using the same gluon labels, with a bar on top. For future reference, we point out that according to the procedure of section 3 exactly the same colour loops will emerge from closed secondary flows regardless of whether quark-antiquark labels, or barred gluon labels, are employed.
It is a matter of simple algebra to see that the iteration of eq. (5.7) leads to the following result: is understood. Note that, apart from a factor N −ρ(γ) , the colour factor Λ in eq. (5.16) is a string of n − k + q Kronecker δs. Concerning the −1/N factors mentioned above, they collectively give rise to the (−1) k /N k factor on the r.h.s. of eq. (5.16). As one can see, it is extremely easy to keep track of them, by means of the summation index k, which essentially counts the number of gluons whose labels are removed (thanks to the second term in the Fierz identity) from the secondary flows.

JHEP11(2021)045
subspace, which one projects onto by means of P (s k ) ⊥ . Therefore, U(1)-gluon vectors are orthogonal to physical-gluon vectors, which belong to the S (p) {a} subspaces; it is therefore natural to expect that the interferences between such terms vanish.
The k = 0 term in eq. (5. 16) is what is usually called the "flow representation" of the given amplitude. That equation clarifies the extent to which this is an improper terminology (strictly, two representations are two different ways for writing the same object). As can be seen there, the fundamental and flow representations differ by non-null vectors, which have at least one component in the U(1)-gluon subspaces.
We now turn to computing the analogue of eq. (5.1), i.e. the interference of two amplitudes Equation (5.16) implies an immediate simplification. By using it to write symbolically one obtains (see eq. (5.18)) owing to the fact that and eq. (4.10). From eq. (5.16) one can read the explicit form of w(s 0 )| w(s 0 )| = A (2q;n) (γ ) 1 2 n/2 In the computation of eq. (5.22), one therefore finds and

JHEP11(2021)045
having employed the leftmost part of eq. (4.11). The colour factor Λ({aij},γ & & s k ) in eq. (5.16) does not depend on any of the indices i p and j p if p ∈ s k . One can therefore carry out the sums over these indices, exploiting the second Kronecker δ in eq. (5.27): this results in a factor N k , that cancels an identical factor in the denominator of eq. (5.27). Conversely, the dependence on the indices that appear in the leftmost Kronecker δ in eq. (5.27) is non trivial, owing to the colour factor in eq. (5.25). However, it can be simplified by observing that (5.28) By putting everything back together, and by exploiting the rightmost Kronecker δ in eq. (5.27) one finally obtains The colour factor in eq. (5.29) has the same form as any colour factor in a quark-only process. By using the result of section 5.1, one obtains , (5.30) and therefore the sought final result The similarity of this result to that of eq. (5.1) is striking. It implies that quark-gluon matrix elements have essentially the same structure as quark-only ones, provided one expresses them by means of secondary colour flows, thereby introducing U(1) gluons into the picture. The role of the latter is quite prominent in eq. (5.31): their number is equal to the summation index k, and they induce alternating signs ((−1) k ) in the summands. We stress that this is true also on the diagonal of the colour-flow matrix (i.e. when γ = γ), which is a real number squared. Furthermore, U(1) gluons give eq. (5.31) a hierarchy in powers of 1/N . Unfortunately, this hierarchy is not absolute for off-diagonal elements. This is due to the fact that while L(γ JHEP11(2021)045

Summing over flows
The matrix elements relevant to a given closed flow given in eq. (5.31) need be summed, in most physical applications, over flows. Owing to the summations over secondary flows, this operation is less trivial than its quark-only counterpart, which had led immediately to the result of eq. (5.2). We perform it in this section, showing in the process that it leads naturally to the definition of linear combinations of dual amplitudes where U(1) gluons play a special role. The relevant summations can be carried out directly on the scalar quantities on the r.h.s. of eq. (5.31); this implies two sums, as both the L-and R-flows must be summed. It is therefore more convenient to arrive at the sought result by summing amplitudes, as opposed to interferences of amplitudes, since this entails carrying out a single sum; in other words, we exploit eq. (2.26) rather than the rightmost side of eq. (2.27). Thus, we consider |A (2q;n) = γ∈F 2q;n |A (2q;n) (γ) , (5.33) where the summands on the r.h.s. are given in eq. (5.16). The key idea is the following: so far, secondary colour flows are quantities derived from some given colour flows. However, they have a structure which is perfectly well-defined in its own right; therefore, when a sum over flows has to be carried out, it is possible to exchange the order in which the sums over flows and secondary flows are performed, thus inverting the logic, since the former are now seen as quantities derived from the latter. In order to do this, one starts by introducing a notation for the sets of q-quark, (n − k)gluon colour flows as follows: 2q;n−k ≡ F 2q;n−k gluon labels in 1, . . . n \ j 1 , . . . j k . (5.34) 2q;n−k is identical to F 2q;n−k up to the labels of the gluons that enter its flows, which must be different from the indices j i that belong to the given set s k . This implies With this, one can generalize the operator I + of ref. [23] (see appendix B there). In particular for any k ≥ 1; when k = 0, I + is defined to be the identity. The single-label I + operator employed k times on the rightmost side of eq. (5.36) coincides with that ref. [23], except for a minor notational detail: here, the upper index j α indicates the label of the gluon to be inserted to the right of position i α ; this index was not necessary in ref. [23], because j α = n + 1 there. Implicit in eq. (5.36) is the fact that the I + operators are to be used

JHEP11(2021)045
in a right-to-left order: first the gluon with label j k is inserted, then that with label j k−1 , and so forth up to j 1 . As far as the values that the position indices i α can assume are concerned: for the first gluon to be inserted (j k ) one has q + n − k possibilities (i.e. to the right of each quark and of each gluon in the given flowγ); for the second one, the possibilities are q + n − k + 1, since that gluon is now inserted in the flow I (j k ) + (i k )γ which has one additional gluon w.r.t. the original oneγ. The procedure is then iterated till the last gluon to be inserted (j 1 ). Since the I + operator of eq. (5.36) is single-valued given the position parameters i α , the dimension of the set that contains all possible flows generated by I + acting on a givenγ is equal to the dimension of the codomain of the acceptable values of the i α parameters. We have therefore for any given flowγ ∈ F Taken together, the above facts imply that and I (s k ) Therefore, for any function F that is summed over all flows and depends on a q-quark, n-gluon flow and on one of its secondary flows, the following identity holds: for any k and s k ∈ S (n) k ; the number of summands on the two sides of eq. (5.42) is the same, owing to eqs. (2.11), (5.37), and (5.40).
We now consider the sum on the r.h.s. of eq. (5.33), by replacing the amplitude there with its expression given in eq. (5.16). There are no cross-constraints between the sum over γ and those over k and s k , and this implies the latter two can be used as the outermost ones, which leads us to the following result:

JHEP11(2021)045
In this expression, the part to the right of the sum over F 2q;n (including it) is in the same form as the l.h.s. of eq. (5.42). Thus, by employing the r.h.s. of that equation instead we arrive at where we have defined the following linear combinations of the standard (because of eq. (5.39)) dual amplitudes: In words: for a given q-quark-antiquark, (n−k)-gluon colour flowγ, an amplitudeĀ (2q;n) (γ) is associated with it that is obtained by summing the dual amplitudes relevant to the colour flows constructed by inserting k gluons inγ in all possible manners. Those extra gluons are colourless as far as the total amplitude is concerned, which is consistent with their (implicit) role inγ as U(1) gluons. Starting from eq. (5.44), we can repeat the procedure that has brought us from eq. (5.20) to eq. (5.31), to obtain As was already the case for the matrix elements at given closed colour flow (eq. (5.31)), eq. (5.46) is very similar to its quark-only counterpart, eq. (5.2). The same comments made about eq. (5.31) vs eq. (5.1) in this regard apply here, with the additional observations on theĀ (2q;n) (γ) amplitudes of eq. (5.45). Note that eqs. (5.31) and (5.46) share an important feature, namely: not only do quantities associated with different numbers of U(1) gluons not interfere with each other but also, for the interference to have a non-zero result, the identities of such gluons must be the same on both sides of the cut (technically, the set s k is relevant to both the L-and R-flow). We remark that expressions whose content is analogous to that of eq. (5.44) appear in ref. [19], where use is made of the same linear combinations of dual amplitudes as in eq. (5.45). What we have shown here is that: a) at the level of amplitudes, the structure of eq. (5.44) is already present at fixed colour flow (eq. (5.16)); b) the number of summands associated with different sets of U(1) gluons (the differences being in their numbers and/or identities) is equal to 2 n (see eq. (5.14)) at the level of both the amplitudes (eqs. (5.16) and (5.44)) and the matrix elements (eqs. (5.31) and (5.46)). This fact is non-trivial: one

JHEP11(2021)045
can show that by naively (i.e. without using the properties of projectors) squaring eq. (5.44) to obtain eq. (5.46), an individual term with k U(1) gluons in the latter expression results from the sum of 3 k terms, 3 k − 1 of which mutually cancel.

Gluon-only matrix elements
The results of section 5.2 show that ultimately gluons can be made to behave as extra quark-antiquark pairs in quark-gluon amplitudes. This is a strong hint that the same property will hold true also in the case of gluon-only amplitudes, i.e. when there are no quarks or antiquarks to start with. In this section, we shall prove that this is indeed the case. At the end of the day this will not matter when summing over colour flows, since in that case Ward identities will ensure that U(1)-gluon contributions vanish. However, at fixed flows U(1) gluons do play a non-trivial role, which puts gluon-only amplitudes on the same footing as quark-only ones, as has already happened with quark-gluon amplitudes. Let be a generic n-gluon colour flow. We define two operators, qR and qL, that by acting on σ result in the following 2-quark-antiquark, n-gluon colour flows: which among other things rely on the fact that for any σ. With these definitions, for any closed colour flow (σ , σ) we have

JHEP11(2021)045
Furthermore, the definitions of eqs. (5.48) and (5.49) respect the colour-loop structure. In other words, by using the results of section 3, it is straightforward to prove the following properties: • The closed colour flow (σ , σ) has the same number of colour loops as (σ qL , σ qR ) • The colour loops in L(σ , σ) and L(σ qL , σ qR ) are pairwise identical, except for that which contains σ(n) as an R-element, and that which contains σ (1) as an L-element; 10 these differ only because the loops in L(σ qL , σ qR ) include the fictitious quarks and antiquarks labelled −1, . . . − 4.
• Given any loop 1 ∈ L(σ , σ), let 2 ∈ L(σ qL , σ qR ) be the loop corresponding to it according to the previous item, so that for any two gluon labels k and l we have Then the distance introduced in eq. (3.3) fulfills the property with i an integer. In other words, the distance between k and l has the same parity in the 1 and 2 loops.
Equations (5.53)-(5.57) are sufficient to allow one to repeat, with only notational changes, what has been done in section 5.2. The analogue of eq. (5.16) reads as follows: Note a technical peculiarity: because of the necessity to use the qR and qL operators, vector space. However, the fictitious-quark degrees of freedom only play the role of giving the correct normalization and index-contraction conditions when an interference of amplitudes is considered. More significantly, the contributions to eq. (5.58) with k ≥ 1 (i.e. with at least one U(1) gluon) do not vanish. This is contrary to what one normally expects, but in fact that expectation is based on results fully summed over colour degrees of freedom. In

JHEP11(2021)045
other words, not summing over those, and fixing e.g. the colour flow, the contributions due to U(1) gluons are non-trivial. The implication of this is that, at fixed flows, the flow and the fundamental representations do not coincide for gluon-only amplitudes, thus making the latter behave exactly as their quark-gluon counterparts. As we shall see later, this symmetry is lifted upon summing over flows, and this thanks to the dual Ward identities, which have no analogue in the quark-gluon case. By interfering the amplitude of eq. (5.58) with an analogous one on the left-hand side of the cut (relevant to an L-flow equal to σ ), one obtains the analogue of eq. (5.31), namely having used eq. (5.52). As was the case at the amplitude level, U(1) gluons do contribute non-trivially to eq. (5.60).

Summing over flows
The sum over colour flows of the results of eqs. (5.58) and (5.60) follows the same steps as in section 5.2.1. Therefore, one readily arrives at the analogue of eq. (5.44), namely where, similarly to what was done in eq. (5.45), we have been led to definē 2q;n−k , i.e. the set of (n − k)-gluon colour flows with indices not equal to any of the integers that belong to the set s k (the presence of the fictitious quark labels is understood, and ignored notationwise). Furthermore, the I + operators that appear on the r.h.s. of eq. (5.62) are defined in a manner formally identical to that of eq. (5.36), with the condition that the indices i α can only assume positive values here; this implies that no gluon will be inserted to the right of either quark 11 inσ. This constraint is necessary in order to preserve the cyclicity properties of gluon-only flows. Because of it, one observes that, for any given values of i 2 . . . i k , the sum over i 1 in eq. (5.62) features n − 1 terms, which are exactly those that enter into the dual Ward identity. Thereforē This is the anticipated, and expected, result: when the sum over flows is carried out, all of the contributions due to U(1) gluons (which are non-zero vectors that live in the S subspaces, with p = 1 . . . n) add to zero. One must bear this fact in mind when claiming that, for gluon-only amplitudes, the fundamental and the flow representation are strictly equivalent.
One can finally use eq. (5.65) to arrive at the result, analogous to that of eq. (5.46), for the matrix elements which again is as expected. Although the notation of eq. (5.66) bears memory of the usage of secondary flows, in the absence of U(1) gluons, and in view of the definition of colour loops adopted here, there is no difference between the original flows and the secondary ones in eq. (5.66).

Soft limits
We consider here the tree-level matrix elements relevant to processes that feature at least one gluon, in the limit in which one of such gluons is soft, while all of the other partons are hard and well separated from each other. These soft limits have been studied in detail in ref. [23]; we shall show below how the formulae of that paper assume an extremely simple form with a geometrical character by exploiting what has been done in section 5. The results of ref. [23] were derived with the goal of employing them in the FKS subtraction formalism [30,31] for the computation of NLO-accurate cross sections and, to that end, were given in three different forms, namely: at fixed colour configurations, and at fixed colour flows either at the real-emission or at the Born level. Since the emphasis of the current paper is on PSMCs, working at fixed colour configurations is of no interest. Fixing the colour flows at the real-emission level is also problematic. Such kinematics configurations are those one obtains after the first emission in a PSMC; therefore, information gathered from the matrix elements at that level cannot easily be used by the PSMC to generate those very same emissions. In keeping with its Markovian structure, a shower will generate emissions in its k th step by using information available at the (k − 1) th step. The most convenient option, then, is to work at fixed Born-level flows, which is what we are going to do in the present section.
Consistently with the previous statement, when referring to processes that feature n gluons we understand those gluon to be at the Born level. In particular, a quark-only case is that of the soft limit of a matrix element relevant to those quarks plus an additional gluon that becomes soft.

Quark-only matrix elements
In this case, the quantity of interest is that in eq. (I. 3.86) (with n = 0), which we report here for convenience The r.h.s. of this equation features the usual eikonal factors, where p is the four-momentum of the emitted soft gluon, p k and p l are those of the corresponding quarks or antiquarks, and M (2q;0) kl (γ , γ) is the colour-linked matrix element for the given closed flow. The latter have been presented in eq. (I.3.115); again we report here their expressions Fuller details can be found in ref. [23]. It turns out that the explicit computation of eq. (6.1) leads to the following simple result: The proof of eq. (6.4) proceeds as follows. Start by considering eq. (6.3): of the four terms on the r.h.s., only one may not be equal to zero, owing to the fact that both k and l are either a quark or an antiquark, and the conclusion follows from the presence of the Kronecker δ's. Thus, a single colour factor is relevant, which we write as follows (the notation is symbolic: the appropriate choices of the signs where ± appear are given in eq. (6.3)) The general form of eq. (6.5) can be obtained from eq. (2.14). Taking footnote 5 into account, and denoting by a the colour of the only gluon present in the process, there are only two possibilities for the colour factor of eq. (6.5), namely

JHEP11(2021)045
The case of eq. (6.6) (eq. (6.7)) occurs when the operators I ± insert the gluon onto open colour lines which do not (do) belong to the same closed colour line. 12 It is easy to see that this happens when the partons k and l do not (do) belong to the same Born-level colour loop. This explains the two summations that appear on the r.h.s. of eq. (6.4). By considering the only non-trivial case, eq. (6.7), one sees that the rightmost side is equivalent to removing the gluon from the closed colour line where it appears (with the factor C F accounting for it); the remaining trace of the identity is the contribution of that closed colour line stripped of the gluon (which now contains only quarks and antiquarks, and thus coincides with a colour loop) to the overall colour factor. Then, since the other closed colour lines have not been affected by the I ± operators, and are therefore identical to their Born-level colour-loop counterparts, eq. (6.7) can be written more precisely as follows having repeated the procedure that led us to eq. (5.1). Equation (6.8) essentially completes the proof of eq. (6.4), bar the sign in front of each eikonal. In order to establish that, one observes that the sign in front of each colour-flow matrix element on the r.h.s. of eq. (6.3) (taking the overall minus sign into account) is equal to −1 if k and l are both quarks or both antiquarks, and equal to +1 if one of k and l is a quark and the other an antiquark. According to item ii. in section 3 the former case is that of k and l being both R-elements or L-elements, while the latter case is that of one R-element and one L-element. As such, in the former (latter) case k and l are separated by an odd (even) number of partons. By putting all this together, one obtains the factor (−1) δ (k,l) on the r.h.s. of eq. (6.4). Finally, the factor 1/2 in eq. (6.1) is necessary to avoid double counting, the summands on the r.h.s. being symmetric when k ↔ l; such a symmetry is removed in the sum on the r.h.s. of eq. (6.4) (owing to the k < l constraint), and therefore no factor 1/2 is present there. This concludes the proof. There are two remarkable things about eq. (6.4). Firstly, it is an entirely geometrical formula; no colour algebra is required, and its expression is fully determined by the structure of the colour loops. Secondly, such loops are Born-level quantities, whereas one would expect that soft radiation pattern be determined by (n + 1)-body matrix elements. In fact, there is no contradiction between these two statements: colour-wise, the present case is so tightly constrained (by the Kronecker δ's in eq. (6.3), and by the structure that leads to eqs. (6.6) and (6.7)) as to give a one-to-one correspondence between Born-and real-emission-level quantities.

Quark-gluon matrix elements
The calculation of the soft limit of interest here is trivial if one exploits the results of sections 5.2 and 6.1. The essence of the procedure of section 5.2 is to write the matrix 12 A closed colour line is obtained by merging open colour lines, after reversing those that belong to the L-flow: starting from one given such line in the R-flow, one attaches to its end the colour line in the L-flow whose antiquark label is identical to that of the line in the R-flow; the procedure is iterated until the quark label in the colour line of the last L-flow coincides with that of the first line considered in the R-flow.

JHEP11(2021)045
elements for a quark-gluon process in terms of secondary (i.e., quark-only) flows. However, the procedure is iterative, and can be stopped after dealing with any number of gluons. This is not of particular interest when all partons are hard and well separated from each other, but in the present context it is helpful, since one can get rid of all gluons except that which becomes soft. By doing so, one is left with a sum of q-quark-antiquark, one-gluon contributions, each of which has a soft limit given by eq. (6.4). By putting all together one obtains Equation (6.9) allows one to sum over flows, by following again the same procedure as in section 5.2. Defining from eqs. (5.46) and (6.9) one obtains For future use, we denote the sum in the last line of eq. (6.11) as follows: and call it the radiation pattern associated with the (secondary) closed flow (γ ,γ). The similarity between eqs. (6.9) and (6.11) on the one hand, and eq. (6.4) on the other, is remarkable, in view of the fact that, at variance with the case of a single finalstate gluon, matrix elements that feature at least two gluons have a more involved colour structure, which prevents one from establishing in a straightforward manner a one-to-one correspondence between Born-and real-emission-level quantities. As is shown here, this can be done by employing secondary flows.

JHEP11(2021)045
Needless to say, eq. (6.11) can be cast in the standard 13 form for soft matrix elements summed over flows (see e.g. eq. (I.3.33)), namely where use is made of the colour-linked Borns, M (2q;n) kl , whose definition can be found in eq. (I.3.35). These quantities satisfy the completeness relations (see e.g. eq. (3.14) of ref. [32]) which, as is shown in ref. [23] using the language of the present paper, stem from the conservation of colour. In eq. (6.14) by C(I k ) we have denoted the Casimir factor (equal to C F for quarks and to C A for gluons) of a parton with identity I k .

Gluon-only matrix elements
The case of gluon-only matrix elements can be dealt with as has been done in section 6.
Note that, at variance with what happens in eq. (6.9), the innermost sum on the r.h.s. of eq. (6.15) constrains the indices k and l to be larger than zero, which in turn implies that no quarks or antiquarks will appear in the eikonals. This is an obvious requirement, since the quarks and antiquarks introduced in eqs. (5.48) and (5.49) are unphysical. With a procedure analogous to that which leads to eq. (5.66), from eq. (6.15) one can also obtain the soft matrix element summed over flows

General considerations
Equations (6.4), (6.9), and (6.15), or their flow-summed counterparts (eqs. (6.11) and (6.16), the quark-only case being trivially obtained from eq. (6.4)) tell one that processes that feature only quarks at the Born level are essentially maximally complicated as far as softradiation patterns are concerned. Gluons are either treated as equivalent to a quarkantiquark pair, or they just decouple completely, consistently with their interpretation as U(1) colourless objects.
These results also suggest a natural generalization of the concept of colour partner we are used to in the context of current PSMCs: given a parton (understood here as either a quark, or an antiquark, or one of the two colour indices of an SU(N ) gluon), its colour partners are all of the partons that belong to the colour loop to which that parton also belongs. This idea indeed encompasses the standard one, which is relevant to leading-N contributions, where γ = γ or σ = σ, and no U(1) gluons in the case of quark-gluon and gluon-only processes. In such a situation, loops are trivial, and it is easy to see that one ends up with quark and antiquarks having one colour partner, and gluons having two of them (one for the colour, and one for the anticolour), as usual. In general, however, a parton can have up to (n − 1) colour partners (i.e. up to 2(n − 1) for a gluon, if its colour and anticolour are simultaneously accounted for), corresponding to a more involved radiation pattern than that relevant to a leading-N picture.
It should be clear that, having chosen a parton, its colour partners are those that enter all of the eikonals that also feature the former one. From the results presented before, one sees that eikonals contribute to the cross section with either sign, and therefore cancellations must be expected. Therefore, in a refinement of the previous paragraph, one is entitled to consider as colour partners only those partons that enter the eikonals that survive algebraic simplifications.
Some such simplifications can be understood in a fairly general manner. Consider in particular the case of a gluon, and its role in colour loops, as detailed in item iv. in section 3. Suppose that the two colour indices of a given gluon, which we denote by g , appear in the same loop. It follows that all of the eikonals that feature g will cancel in the physical cross section. This is because those two indices, per item iv., will correspond to one L-element and to one R-element. Any other parton k in the loop will therefore result in two eikonals, both equal to [k, g ] but with opposite signs in front -this is because the two distances δ (k, g ) when computed with g regarded as an L-or an R-element will have opposite parity. A corollary of this is that when a closed colour flow in a quark-gluon process corresponds to a single colour loop, then no eikonal that contains at least one gluon contributes to the cross section.

6.5.1qq → gg
The set of primary flows was given in eq.    Table 1 reports the various quantities that appear in eq. (6.11), for each non-trivial secondary closed flow -number of U(1) gluons (p, first column), ρ parameters for the L-and R-flows (second and third columns), set of colour loops (fourth column, eq. (3.1)), number of colour loops (fifth column), and colour factors (sixth column; these are defined as the products of all the N monomials in eq. (6.11), times the factor (−) p ).
The radiation patterns (see eq. (6.12)) associated with the secondary closed flows of table 1 read as follows: Notice that (6.23) illustrates the fact, mentioned in section 6.4, that gluons appearing twice in any loop do not contribute to the associated radiation pattern. It trivially follows from eq. (3.4) that there are two independent dual amplitudes Thus, the amplitudes associated with secondary flows are the following (from eq. (5.45), and the definition of the I + operator, eq. (5.36)): This leads to four quantities that enter eqs. (5.46) (tree level) and (6.11) (soft limit). Note that amplitudes associated with different sets of U(1) gluons do not interfere with each other. By direct computation, defining s = ( With this, eq. (5.46) leads to (note that ρ(γ α ) = 0 for any α) with the first two, the third, and the fourth term on the r.h.s. stemming from the contributions associated with zero, one, and two U(1) gluons, respectively. Likewise, from JHEP11(2021)045 eq. (6.11) we obtain with the first, second, and third line on the r.h.s. stemming from the contributions associated with zero, one, and two U(1) gluons, respectively. In terms of the colour-linked Borns introduced in eq. (6.13), in this case denoted by M As expected, these expressions satisfy the completeness relations (6.14)  The secondary flows areγ The sets of associated colour loops are given in table 2, together with the computation of their colour factors according to eq. (6.11).
Defining we find by explicit calculation, summing over helicities, 64) Note that only six of these are independent, owing to the property of the dual amplitudes The products of secondary flow amplitudes follow from eq. (5.45), which gives in this casē

Monte Carlo simulation
In a parton shower Monte Carlo, each parton of a primary hard subprocess has the potential to generate a shower, with initial conditions to be determined by the kinematics and colour structure of the hard process. In the case of a dipole shower, each parton is paired with another to form a dipole which sets the initial scale of the shower according to its invariant mass. In a single-parton shower, the partner is again involved in setting the initial scale of evolution. For example, in an angular-ordered shower it determines the maximum angle of emission in the frame of shower evolution. In the leading-colour approximation, the parton and its partner form a colour singlet, so that the end products of the shower from a dipole, or from the angular-ordered showers of a parton and its partner, can hadronize without involving the products of showers initiated by other hard partons. In order to reproduce the leading-order soft gluon emission pattern from an m-parton hard process, the dipoles or parton pairs have to be associated with eikonal terms in the JHEP11(2021)045 soft limit of the (m + 1)-parton process. In the leading-colour approximation, considering all partons as outgoing, every quark can be associated with an antiquark or the anticolour of a gluon, every antiquark with a quark or the colour of a gluon, and the remaining colours and anticolours of gluons with other gluons.
As we have seen, at subleading colour there are eikonal terms that pair quarks with quarks and antiquarks with antiquarks. There may also be a term where a given gluon is paired with a quark and a term where it is paired with another quark, and similarly for pairing with two antiquarks. We can still associate dipoles or parton pairs with such terms, and use them to generate showers, but the final products cannot form colour singlets without communicating with other showers. However, since colour is conserved overall, it is always possible to hadronize after some colour reconnection between showers.
In this section we discuss some strategies that could be adopted to include subleading colour contributions from the hard subprocess. Which of these strategies may prove to be the most practical or efficient probably depends on the process under consideration and remains for further investigation. A related but separate issue is the inclusion of subleading colour contributions inside parton showers, studied for example in refs. [1][2][3][4][5][6]. We envisage that the use of colour flows in matrix elements may prove useful in this case as well, since it has the advantage that it leads to an easy accounting of the number of colours, and thus renders it straightforward to see which contributions need to be taken into account to achieve any desired level of accuracy w.r.t. N .

Single-parton showers
In a single-parton shower, the eikonal pairings are used to set the evolution scales but the showering of partons is otherwise independent. To achieve this the dipoles associated with eikonal terms have to be split in some way between the two partons involved. As an example, we consider here the splitting that leads to angular-ordered shower evolution. For simplicity we adopt the original formulation used in the Herwig6 event generator [33], rather than the more sophisticated treatments in Herwig+ + [34] and Herwig7 [35,36], which are equivalent for our purposes.

Angular-ordered shower
For an angular-ordered parton shower, one makes the following substitution for the eikonal (6.2): where ω is the energy of the emitted soft gluon, ξ kl = 1 − cos θ kl where θ kl is the angle between the three-momenta of partons k and l; ξ k = 1−cos θ k where θ k is the angle between the three-momenta of the soft gluon and parton k, and similarly for ξ l . This approximation follows from separating the collinear singularities of the eikonal and azimuthally averaging each term. The averaging destroys Lorentz invariance and so one has to perform it in the same frame as that used by the parton shower. The subsequent evolution of each shower is then ordered with respect to the relevant angular variable, ξ k or ξ l . This ordering takes

JHEP11(2021)045
into account the coherence of emission from the ensemble of partons in the shower at angles smaller than the current value. Making the substitution (7.1) in eq. (6.13), we have For any given k, using the phase-space element for soft gluon emission, and integrating over φ k , the azimuthal angle of emission with respect to parton k, the differential distribution of soft gluon emission becomes The emission of non-soft collinear gluons can be included by replacing dω/ω by P (z) dz, where P (z) is the appropriate collinear splitting function of z = 1 − ω/p 0 k . In a leading-colour treatment [37], one first generates a kinematic configuration of the Born process using its full matrix element and then chooses among the leading-N colour flows (if there are more than one) according to their relative contributions. In the radiation pattern of a leading-N flow, all eikonal terms have coefficient +1, each quark or antiquark appears in only one eikonal term, and each gluon appears in two. Therefore, if k is a quark or antiquark, there is a unique colour partner l in eq. (7.4), which specifies the angular region of radiation from k, ξ k < ξ kl . On the other hand, if k is a gluon it has two partners, say l and m such that ξ km < ξ kl , and it radiates with half intensity when ξ km < ξ k < ξ kl and full intensity when ξ k < ξ km .
These radiation patterns, together with subsequent multiple emissions, can be generated by a parton showering process controlled by a Sudakov form factor where C(I k ) is the relevant Casimir factor (see eq. (6.14)) and ξ 0 is a cutoff below which emissions are not resolved. When parton k has a single colour partner l, the probability of no emission between the maximum angular scale ξ kl and some lower scale ξ k is given by ∆ k (ξ kl )/∆ k (ξ k ). In the above case of a gluon k with partners l and m such that ξ kl > ξ km , the probability of no emission between ξ kl and ξ k > ξ km is ∆ k (ξ kl )/∆ k (ξ k ) while that between ξ km and ξ k < ξ km is ∆ k (ξ km )/∆ k (ξ k ). In other words, the leading-colour noemission probability for a given colour flow is controlled by

JHEP11(2021)045
the sums being over all colour partners l of k in the chosen leading-N colour flow. Equations (5.46) and (6.11) suggest an obvious extension of (7.7) to include subleading colour. After selecting a closed flow (γ ,γ) according to its contribution to the Born matrix element in (5.46), the evolution of parton k in that flow is controlled by Consider for example the processqq → gg. We first choose one of the six flows in table 1 according to their contributions to the Born matrix element in eq. (6.38). Note that the flows (γ 1,1 ,γ 1,1 ) and (γ 1,2 ,γ 1,2 ) make negative contributions and would have to be treated as counter-events or with a resampling method such as that of refs. [38][39][40]. Showering of partons that contribute to the radiation pattern of the chosen flow, according to eqs. (6.22)-(6.27), then proceeds as in leading colour. However, if any of the subleading colour flows is chosen, then one or both of the gluons do not contribute and the power (7.9) is not defined for them. These are either gluons that appear twice in the same colour loop (in (γ 0,1 ,γ 0,2 )) or else (in (γ i,j ,γ i,j ) with i = 1, 2) the U(1) gluons discussed in section 5.2. They do not shower and can be treated in a manner similar to the self-connected gluons to be discussed in section 7.2.2. For our other example processqq → q q g, the situation is more complicated. The radiation patterns of the 13 contributing closed flows in table 2 are given in (6.74)-(6.86). Here some of the eikonal terms in subleading flows have negative coefficients, leading to negative contributions in eq. (7.9). Note however that the denominator is always +1 for quarks and antiquarks, and either 2 or 0 for gluons. In the latter case the gluon is nonshowering as before. In the former, the gluon may have either two partners with coefficients +1 or four partners, three with +1 and one with −1. Similarly a quark or antiquark may have one partner or three. Consider for example the flow (γ 0,1 ,γ 0,4 ) (which is maximally colour-suppressed). According to eq. (6.77) the gluon is non-showering, while each quark or antiquark has two partners with coefficient +1 and one with coefficient −1. For example, quark −1 has antiquark partners −3 and −4 with coefficients +1 and quark partner −2 with coefficient −1 (this shows in a specific case the implications of the generalization beyond leading colour of colour connections/colour partners, which we discussed in section 6.4). There are six possible angular orderings

JHEP11(2021)045
The last entry in each case shows the corresponding sequence of values of the power (7.9), as ξ −1 evolves from the highest to the lowest values. Note that the power is always 1 at the lowest angles, but can be 2 (higher emission), 0 (no emission) or even -1 at higher values. Negative powers, corresponding to a non-negative-definite Sudakov exponent, have already been considered in the literature (see e.g. refs. [41][42][43]), and can be dealt with for example by using a veto and weighting technique. The incidence of negative flow contributions and/or negative powers (7.9) could lead to low Monte Carlo efficiency, which might be improved by not preselecting a flow, thus allowing all flows to contribute to shower evolution. This corresponds to using (7.4) directly to define the evolution power in place of (7.9), where we have used the completeness relation (6.14). The colour-linked Borns M (2q;n) kl are given explicitly in (6.41)-(6.46) forqq → gg and in (6.89)-(6.102) for qq → q q g. Then since the leading-colour flows contribute to all stages of showering there will be on average a lower incidence of negative powers. Statistically, at the parton level the results should be equivalent to evolving flows separately. On the other hand, all colour information is lost by the end of the shower, which could be problematic for hadronization.

Dipole showers
A dipole Monte Carlo implementation including subleading colour will in general involve contributions from dipole configurations of all possible types, connecting quark-antiquark, quark-quark, quark-gluon, gluon-gluon and their conjugates. With each configuration we can associate a graph specifying the dipole connections. The theory of such graphs is linked to that of flows discussed above, as elaborated in the following subsection.

Dipole graphs
For a given process that features q quark-antiquark pairs (all considered as outgoing) and n gluons let be the sets of the quark/antiquark and gluon labels, respectively, in keeping with eqs. (2.2)-(2.4). We define a dipole graph to be a set of q + n pairs, as follows: In other words: a dipole graph is a set of pairs of non-identical parton labels; each quark and antiquark (gluon) label will appear in it exactly once (twice). Each pair that belongs

JHEP11(2021)045
to a dipole graph is a dipole, and its two elements are the dipole ends. The order in which the pairs enter a dipole graph is irrelevant from a physics viewpoint: they can be re-shuffled as is convenient. We denote by G 2q;n = Γ Γ eq. ( 7.19) , (7.20) the set of all dipole graphs relevant to q quark lines and n gluons. By analogy with eq. (6.12), we call the quantity the radiation pattern associated with the dipole graph Γ. Any dipole graph admits a representation by means of a square (2q + n) × (2q + n) matrix, 14 defined as follows: (7.22) That is, each dipole in Γ gives a unit contribution to the matrix element whose indices coincides with the dipole ends, both in the given and in the reverse order. Owing to eq. (7.19), the matrix defined in eq. (7.22) has the following properties: a) it is symmetric; b) it is traceless; c) its elements are equal to either 0, 1, or 2; d) each of 2q of its rows sums to one, and each of the remaining n rows sums to two. By construction We stress that eqs. (7.21)-(7.23) are invariant under reordering of the dipoles that enter the given dipole graph. The counting of distinct matrices of the type defined in eq. (7.22) tells us the number of possible dipole configurations for q quark-antiquark pairs and n gluons, |G 2q;n |. It is convenient to define the generating function F (x, y) = q,n |G 2q;n | x q q! y n n! , (7.24) so that For a process involving only q quark-antiquark pairs and no gluons, the dipole matrix is a 2q × 2q symmetric, traceless matrix with elements 0 or 1 and each row summing to 1. It is readily seen that the number of distinct matrices of this type is (2q − 1)!!, so that the generating function is 14 In keeping with the labeling of ref. [23] (see section 2), parton labels equal to zero must be ignored. 0  1  2  3  4  5  6  0  1  0  1  1  6  22  130  1  1  1  3  10  46  252  1642  2  3  6  21  93  510  3306  24762  3  15  45  195  1050  6750  50280  425490  4 105 420 2205 13965 103110 867510 8183490 Table 3. Number of distinct dipole graphs for q quark-antiquark pairs and n gluons.

JHEP11(2021)045
In the case of purely gluonic processes, q = 0, the corresponding matrices have all rows and columns summing to 2. The counting of such matrices is a known problem, whose solution is given by the generating function [44] F (0, y) = e −y/2+y 2 /4 √ 1 − y . (7.27) For processes that involve both quarks and gluons, we can suppose that n − k of the gluons are connected amongst themselves while k are inserted onto one or more of the quark lines. All the permutations of the latter are distinct, and the number of ways of distributing them is the number of ways of locating q − 1 subset boundaries amongst k + q − 1 objects, so the number of distinct dipole graphs is 15 which yields the generating function This gives the numbers shown in table 3.

Self-connected gluons
At leading order in colour, and also at subleading order for processes involving only quarks or only gluons, the eikonal combinations that appear in the differential cross section for emission of an additional soft gluon can be represented by a weighted sum of dipole graphs involving all the partons of the primary process. In processes involving both quarks and gluons, however, at subleading order in colour, it is in general necessary to include graphs with fewer participating gluons than in the primary process. In essence, this is due to the presence of gluons both of whose colour labels enter the same colour loop; as is shown in section 6.4, such gluons (denoted by g there) cannot contribute to radiation patterns. Further discussions on this point will be given in section 7.4 and appendix B. 15 Note therefore that this is analogous to the counting problem that leads one to eq. (2.11). 0  1  2  3  4  5  6  0  1  1  2  5  17  73  388  1  1  2  6  23  109  618  4096  2  3  9  36  177  1035  7029  54462  3  15  60  300  1785  12315  96720  852630  4 105 525 3150 21945 173985 1546965 15250200 Table 4. Number of distinct dipole graphs for q quark-antiquark pairs and n gluons, allowing any number of self-connected gluons (see eq. (7.32)).

JHEP11(2021)045
Although such non-participating gluons are not associated with eikonal terms in the soft gluon cross section, they cannot be altogether ignored in a Monte Carlo implementation, since to do so would violate momentum conservation. Instead they can be treated as having self-connected colour and anticolour. Such a self-connected gluon i is associated with a null eikonal [i, i], which cannot radiate: it propagates without showering and must be hadronized through momentum transfer with other showers in the same event.
To make connection with the theory of dipole graphs, for any 0 ≤ p ≤ n and s p ∈ S (n) p as in eq. (5.11), we define G 2q;n−p,sp = Γ Γ eq. ( 7.19) , g s −→ g s s p (7.30) to be the set of q-quark, (n − p)-gluon dipole graphs, where the gluon labels do not assume the values in the set s p . The gluons whose labels belong to the set s p are understood to be self-connected gluons. Note that p = 0 =⇒ s 0 = =⇒ G 2q;n,s 0 ≡ G 2q;n . If we denote the generating function for this by F (x, y), then F (x, y) = e y F (x, y) = e y/2+y 2 /4 √ 1 − 2x − y , (7.33) which gives the numbers shown in table 4.

7.3.1qq → gg
Consider again the processqq → gg. According to tables 3 and 4 there are six relevant dipole graphs, Γ 1 . . . Note that these are not independent, since The closed-flow radiation patterns (6.22)-(6.27) can be expressed in terms of the above dipole radiation patterns as follows: Thus, for the soft radiation matrix element (6.39) we now have Using (7.40) we may, for example, eliminate D (Γ 6 ) which involves two self-connected gluons.

Figure 3.
Dipole graphs for qq → q q g. In the cases of Γ 7 , Γ 8 , and Γ 9 the self-connected gluon is also depicted.

JHEP11(2021)045
We note the following identities: which allow one to eliminate two (but not all) of the dipole graphs with a self-connected gluon.

Closed colour flows versus dipole graphs
In this section, we shall show how the radiation patterns that stem from closed colour flows, eq. (6.12), can be represented in terms of analogous quantities, defined as properties of suitably-chosen dipole graphs, which involve their radiation patterns of eq. (7.21). In order to do so, we start by observing that the radiation pattern of any closed flow (γ ,γ) can be decomposed into the radiation patterns of the individual colour loops associated with that flow, thus Furthermore, the reader must bear in mind that the number of elements in a colour loop is an even number; half of the elements are L-elements, the other half are R-elements (see section 3); they are alternating, beginning with an R-element. We write this as follows: ∈ L(γ ,γ) , = (e 1 , e 2 , . . . e 2m ) , (7.73) so that 2m is the number of elements in the loop, and summands in E + (γ ,γ) and E − (γ ,γ), respectively. Note that is indeed the total number of summands in E (γ ,γ) (before any algebraic simplifications).
Since the two sums in eq. (7.75) are linear combinations of eikonals with coefficients all equal to one, we see the possible emergence of radiation patterns of dipole graphs JHEP11(2021)045 in eq. (7.72). More precisely, the argument goes as follows. Equation (7.76) suggests organizing the parton labels associated with the summands in E + (γ ,γ) in m sets, each composed of m pairs of labels. Likewise, from eq. (7.77) one obtains m − 1 sets, each composed of m pairs of labels, relevant to E − (γ ,γ). These sets of pairs of labels are candidates to be subsets of suitable dipole graphs. For this to happen, all quark and antiquark labels that belong to the colour loop must appear exactly once in each of such sets, whereas gluon labels must appear at least once and at most twice. These properties must be fulfilled for each of the loops in L(γ ,γ), so that when the candidate subsets are combined one forms properly-defined dipole graphs. This can happen since m = q + n , (7.79) that is, the number of elements in the joined candidate subsets is equal to the number of dipoles in the dipole graphs. The fact that candidate subsets emerging from different loops have disjoint sets of labels is trivially true for quarks and antiquarks, but not necessarily so for gluons (see items i.-iv. in section 3). We point out that there are m , (m − 1) , (7.80) possible ways of joining the candidate subsets that emerge from different loops, for E + (γ ,γ) and E − (γ ,γ), respectively. In appendix B we shall sketch a proof of the equivalence of colour flows and dipole graphs. Here, we limit ourselves to stating that the consequence of this equivalence is the existence, for any given closed colour flow (γ ,γ), of two sets of dipole graphs for suitable integer coefficients c ± (Γ) (which, as we shall show in appendix B, are typically equal to one). We point out that the sets in eq. (7.81) are not necessarily disjoint, and that their choice is not necessarily unique. While the proof of eqs. (7.82) and (7.83) is complicated in general, an exception is that of a diagonal closed flow (i.e. one for which the L-and R-flow coincide, (γ,γ)). In that case, there are q + n colour loops, 16 each with two elements; that is, m = 1 ∀ , and E − (γ,γ) = 0 =⇒ E (γ,γ) = E + (γ,γ) .
The set on the r.h.s. of eq. (7.85) is a dipole graph that belongs to G 2q;n ; therefore, the leftmost set in eq. (7.81) contains a single element, and the rightmost one is empty. When the L-and R-flows do not coincide the sums on the r.h.s. of eqs. (7.82) and (7.83) are non-trivial. Without knowing their explicit forms, one can still establish a rather general property, namely that the leftmost set in eq. (7.81) is not a subset of G 2q;n , but rather of ∪ p,sp G 2q;n−p,sp , for suitable values of p and s p . The argument is the following: a necessary condition for the following relations to hold true is the absence of self-connected gluons. Indeed, when this is not the case, of the m 2 summands in the leftmost sum in eq. (7.74) there are n g , which are equal to zero, being self-eikonals. Here, by n g , we have denoted the number of gluon labels that appear twice in the loop . Therefore, in general out of m 2 − n g , pairs one cannot identify an integer number of non-overlapping sets (the candidate subsets of suitable dipole graphs), each with m pairs of labels. The fact that in general n g , = 0 for some implies that eq. (7.86) does not hold true, and must be modified as follows: We point out that this is necessary in spite of the absence of U(1) gluons in the flow (or, more precisely, for any given number of U(1) gluons -see footnote 16), which further clarifies the fact that the p gluons missing from the sets on the r.h.s. of eq. (7.87) are not U(1) gluons, but rather the self-connected ones introduced in section 7.2.2. In order to avoid any confusion, when possible we shall number the self-connected and U(1) gluons by means of labels k and p, respectively. Thus, for each 0 ≤ k ≤ n and s k ∈ S (n) k , we construct the subsets By inserting eqs. (7.82) and (7.83) into eq. (6.11) we arrive at an alternative expression for the soft limit of a generic quark-gluon matrix element .
We point out that the matrix elements defined in eq. (7.91), at variance with all of their counterparts used so far, generally contain contributions associated with different powers of N . Because of this, it is necessary to include the monomials in N into their definitions, again contrary to what was done so far. There is an ultimate physical reason for that. Namely, the structure of dipole graphs, although it extends that of the colour connections employed in current leading-colour PSMCs, is still based on a picture which takes a single side of the Cutkosky cut into account. However, the results of section 6 show that beyond leading colour it is necessary to include simultaneously the information on both sides of the cut, which can be done by means of closed colour loops. If the more complicated structure of such loops is traded for the simpler one of the dipole graphs, this entails a set of more complicated matrix elements. By construction, the sets on the l.h.s. of eq. (7.88) may not be unique. However, different choices of G q,n−k,s k , while changing the individual summands on the r.h.s. of eq. (7.90), will leave the sum invariant -in other words, there are subsets of dipole graphs whose radiation patterns are not linearly independent of one another, as we have seen in the examples of section 7.3. With increasing numbers of quarks and gluons, the possibilities of choosing different G q,n−k,s k also increase. This flexibility may turn out to be useful, since a suitable choice could lead to a reduction of negative weights.

Subevents
In order to implement a dipole shower simulation including subleading colour, one needs to take account of all relevant dipole graphs for a given hard process configuration, and then combine the results. This is necessary because, as we have seen, contributions from graphs with self-connected gluons cannot in general be avoided. Simply choosing a single one out of the set of dipole graphs according to their relative weights would lead to final states where hard hadrons appear that stem from non-showering self-connected gluons, and are mostly isolated. Instead, those hadrons can be more conveniently represented as components of a jet that includes showering of other contributing graphs. We call each dipole graph contribution a subevent of that hard configuration. The final state of that event is then the result of summing all subevents.
As an example, consider the processqq → gg studied in section 7.3.1. For any given Born configuration, one would generate six subevents by showering and hadronization starting from the dipole configurations in graphs Γ 1 . . . eq. (7.47), then combined to produce a complete Monte Carlo event from that Born configuration. Each final-state object would be entered into histograms with the weight of the subevent that generated it. In this example, all the products of dual amplitudes M 1 . . . M 4 , given in eqs. (6.34)-(6.37), are positive, so (7.47) implies that only the subevents from graphs Γ 4 and Γ 5 have negative weights, while Γ 3 does not contribute. Alternatively, we may use eq. (7.40) to eliminate Γ 4 and Γ 5 (or Γ 6 ), but then other subevents will get negative weights for certain Born configurations. This flexibility in the dipole graph approach may be useful for optimizing Monte Carlo efficiency. A possible way to deal with the remaining negative weights could be a resampling method such as that of refs. [38][39][40].
A similar approach would appear necessary in single-parton shower implementations that treat individual colour flows separately, as in eq. (7.9). Flows that do not involve one or more gluons from the Born process will again give rise to isolated hard hadrons coming from such non-showering gluons. If the result of showering each flow is treated as a subevent from the same Born configuration, these hadrons become a subleading correction to the hard component of the associated gluon jet.

The mixed QCD-QED case
An interesting aspect of the formulation presented in this paper is that it applies, essentially unchanged, to a mixed-coupling scenario, i.e. one in which perturbative expansions are obtained as power series in two parameters, both of which are regarded as "small". The foremost candidate is that of QCD+QED, where the two parameters alluded to before are the coupling constants, α S and α, of the respective theories; in order to be definite, in the following we shall understand this case. We point out that while α α S , the existence of a clear hierarchy between the expansion parameters is not a necessary condition for this formulation to work. We also remind the reader that matrix elements are proportional to α n S α m ; at the LO, n and m are constrained to obey n+m = k 0 , with k 0 a process-dependent integer, while at the NLO n + m = k 0 + 1, and so forth. The interested reader can find in ref. [32] a classification scheme for the various contributions to a given N k LO accuracy.
In the QCD+QED case, several kinds of results become relevant at the LO and NLO accuracy for any given combination α n S α m : a) tree-level matrix elements; b) soft-gluon emission patterns; c) soft-photon emission patterns.
The statement we make is that the formulae of section 5 can be used, with minimal or no changes, to deal with item a), and those of section 6 to deal with item b); as far as item c) is concerned, its treatment entails a significantly simplified version of its gluon-emission counterpart.
A detailed discussion of the items above is beyond the scope of this paper. Therefore, here we shall limit ourselves to presenting a simple example, which should be sufficient to give the reader an idea of the procedures to follow in a general case.

JHEP11(2021)045
• Example. Consider the following process: Incidentally, the fact that there are two quark-antiquark pairs of the same flavour allows us to show how to handle in practice such a case, which has been neglected so far in the interest of a simpler notation. Following ref. [23], we start by distinguishing flavour-line configurations i.e. f A (f B ) corresponds to the graphs where there is a fermion line that connects quark −1 with antiquark −3 (−4), and quark −2 with antiquark −4 (−3). Conventionally, we identify f A -and f B -type graphs with s-and t-channel exchanges, respectively. Although the two flavour-line configurations of eqs. (8.2) and (8.3) induce the same colour flows, such flows need to be distinguished, because their contributions to the matrix elements must be accounted for separately (also in view of the fact that the associated colour factors might differ). In the case of gluon exchange we have whereas in the case of a photon exchange

JHEP11(2021)045
All of the other combinations are identical to those above, according to the identifications of eqs. (8.4)-(8.7). Conversely, as far as the dual amplitudes associated with the colour flows are concerned, we have the following: Here, we have denoted by A s and A t the s-and t-channel amplitudes, stripped of charge and colour factors (contrary to the conventions used so far), so that There are three tree-level LO matrix elements, understood to be multiplied by the factor α n S α m with n + m = 2, which we denote by M (4;0) (n,m) . (8.24) In order to compute them, it is sufficient to apply eq. (5.2), with the sums over flows determined by the specific contribution we need to derive. Explicitly  .7)). Such contributions would constitute the building blocks for the determination of the shower initial conditions in a leading-N PSMC, e.g. according to the prescription of ref. [37]. But while for eqs. (8.28) and (8.30) the relevant kinematical quantities are "squared" matrix elements (either M s or M t ), for eq. (8.29) what is relevant is the interference between the s-and t-channel amplitudes (M st ). This shows that, even at leading N , when the matrix elements that underpin the hard events that initiate the QCD showers do not have the maximal power of α S one must consider objects that are not necessarily positive definite. What is important is that the nature of such objects is determined by closed colour flows, and not by their kinematical characteristics. Direct evidence of this fact can be seen in the comparison of eqs. the latter equations is reflected in that of the contributions to the former ones, and the correspondence is driven by the underlying closed colour flows.
We conclude this section by considering soft-photon emission patterns. This case is significantly simpler than its soft-gluon counterpart, owing to the fact that the charge-linked Borns (i.e. the matrix elements that multiply suitable linear combinations of eikonals) are all proportional to the LO tree-level matrix elements. As far as the linear combinations of eikonals are concerned, they are entirely determined by the charges of external partons that participate in the hard process. Thus, from eq. (3.12) of ref. [32], one defines the photon radiation pattern, i.e. the QED analogue of eq. (6.12) The fact that on the r.h.s. of eq. (8.36) one sums over all pairs of external partons implies that this radiation pattern is quite similar to a soft-gluon one relevant to a closed flow with a single colour loop that includes all the partons of the process. This is remarkable, because it puts QED and QCD radiation on the same footing, and allows one to treat the former in the same way as (one of the cases of) the latter. However, in general the two patterns are not identical, since the eikonals in E QED have coefficients that, in absolute value, are not all equal to each other, which does not happen in QCD. While this entails only minor modifications to the angular-ordered showering approach outlined in section 7.1.1, it is potentially more problematic for dipole showering, since the equivalence of the radiation pattern of eq. (8.36) and dipole graphs must be re-considered.

Conclusions
In this paper we have investigated the role that colour flows play in both the calculation of tree-level matrix elements and of their soft limits, and in the inclusion of effects beyond leading colour in parton shower Monte Carlos.
For these applications, and in general for any kind of non-approximate approach, it is appropriate to work with what we have called closed colour flows (section 2), which generalize the colour flows usually employed in leading-N computations (which are therefore more appropriately referred to as open colour flows). In essence, closed colour flows are based on a picture that features a double copy of the colour/anticolour labels of the partons that participate in the process (one copy for each side of the Cutkosky cut; each such copy is an open colour flow). Leading-N contributions are (possibly a subset of) those for which JHEP11(2021)045 the orderings of the labels in the two open flows that form a closed flow are identical to each other.
Colour loops, which we have defined algorithmically in section 3, give a representation of closed colour flows that can be conveniently used in matrix element computations and PSMCs alike. By means of colour loops, one generalizes the concept of colour connection that is used at the leading N (see in particular section 6.4). We have discussed how colour loops can in turn be equivalently described in terms of dipole graphs (section 7), which are an alternative to the former in PSMC implementations. Interestingly, the counting of the numbers of dipole graphs gives rise to integer sequences that, to the best of our knowledge, are presented here for the first time (sections 7.2.1 and 7.2.2), and for which we give generating functions.
In the case of purely gluonic processes, the colour loops as defined in this paper are seen (appendix A) to correspond to quantities emerging from the mathematics of permutations (in particular, the alternating cycles of a cycle graph). This is interesting, because it could pave the way to exploiting results that are increasingly important in other branches of science (mainly genome studies -the curious reader can find a list of classical mathematics papers, relevant to biology, for example in chapter 9 of ref. [45]).
While when working with colour loops each colour or anticolour label may have as many as m − 1 colour partners (m being the number of partons of the process), dipole graphs are more similar to the colour connections one is used to at the leading N , in that they feature a single partner per (anti)quark, and two partners per gluon. Furthermore, while colour loops require a double copy of the partons (as the closed colour flows from which they stem), dipole graphs need a single copy, and thus are topologically identical to any set of colour connections employed in existing PSMCs. At variance with the latter, however, dipole graphs feature colour-colour and anticolour-anticolour connections, in addition to the usual colour-anticolour ones. We stress again that, in spite of their topological differences, colour loops and dipole graphs give equivalent representations of the same thing. The simpler topological structure of the latter is compensated by the necessity of working with more complicated objects at the level of amplitudes squared.
Regardless of their possible use in the context of PSMCs, colour flows give an organising principle for the computation of tree-level matrix elements (section 5). We have shown that the treatment of any process is in fact no more complicated than that of one featuring only quarks and antiquarks, thanks to the introduction of the idea of secondary flows (section 5.2); these allow one to express the results for processes that include gluons in terms of the formulae relevant to the no-gluon case.
Although from tree-level matrix elements results it might seem that secondary flows are a computational device, the results for the soft limits (section 6) show that they are in fact more fundamental objects, in that they control the soft-radiation patterns with a clear hierarchy in N . Furthermore, through secondary flows one recovers a fundamental property of the leading-N PSMCs, namely the fact that the radiation pattern of the first emission is entirely controlled by Born-level quantities. This property, which is trivially true at leading N but not beyond it, allows one to incorporate subleading-N effects in PSMCs in a Markovian way.

JHEP11(2021)045
For matrix element computations we often make use of vector spaces (that represent colour degrees of freedom), which are separable into spaces relevant to individual partons (section 4). In the presence of gluons, we have shown how the so-called fundamental and flow representations are associated with vectors in an (N 2 − 1)-and N 2 -dimensional space respectively for each gluon, the former of which can be seen as a subspace of the latter, and whose complement in the latter is interpreted as representative of the degrees of freedom of U(1) gluons.
Matrix elements for processes that feature both quarks and gluons can be written as the incoherent sum of non-null contributions associated with different sets of U(1) gluons, in the same way as for amplitude-level quantities (see in particular eqs. (5.46) and (6.11)). In turn, these sets are associated with non-overlapping sets of secondary flows, thereby confirming the physical nature of the latter. U(1) gluons do appear also in purely-gluonic processes but, as is well known, they do not contribute to the cross section. However, this is true only after summing over all possible flows; conversely, if one restricts oneself to a subset of the flows, U(1)-gluon contributions are generally different from zero.
In section 7 we have discussed possible applications of colour flows to parton shower Monte Carlo simulations. In particular, we have proposed ways in which the subleading colour structure of Born matrix elements could be included in the matching of hard processes to parton showers. One possible approach, based on eq. (6.11), is to select closed flows according to their relative contributions to the Born process and then use the radiation pattern of each flow to influence the evolution of the associated showers. This could be done by using a modified Sudakov factor that reflects the sequence of scales associated with the colour connections of the evolving partons in the selected flow. In the case of an angular-ordered shower, the eikonal terms in the radiation pattern of the flow define a sequence of angular regions for the evolution of each participating parton, the rate of evolution in each region being controlled by the relevant colour-linked Born matrix elements.
Another possible Monte Carlo approach is more suited to the evolution of dipole showers, being based on dipole graphs and eq. (7.90) rather than directly on colour flows. We showed that each flow can be represented by a linear combination of dipoles, which can serve as initiators of dipole showers.
An issue that will affect the implementation of any scheme for including subleading colour contributions is that these can have either sign, leading to negative Monte Carlo weights. There are established methods for handling such weights [38][39][40][41][42][43]; which of them is most suitable in this case will probably depend on the process under study. A less familiar issue that we have highlighted is the appearance of self-connected gluons in both the singleparton and dipole implementations. Rather than treating these as non-showering objects that would give rise to isolated hard hadrons, we have suggested interpreting them as corrections to gluon jets generated by more leading contributions from the same momentum configuration of the Born process. In this approach, each Born configuration would be used to generate several "subevents" corresponding to different flows or dipole graphs.
The inclusion of subleading colour contributions via subevents, based on either colour flows or dipole graphs, will require some modification of hadronization models, since, as JHEP11(2021)045 we have seen, the concept of colour partners is generalized to include quark-quark and antiquark-antiquark pairs, as well as gluons that connect to two quarks or two antiquarks. However, since colour is conserved overall, it is always possible to hadronize after allowing some colour reconnection among showers within each subevent.
We have also briefly discussed (section 8) how the techniques developed here in QCD can be applied with minor or no changes to the case where two perturbative couplings are relevant. Presently, the most important application of this is in the context of QCD+QED computations and, for matrix elements that do not have a maximal power of α S , it affects PSMC simulations even at the leading N .
We have usually presented our results both at fixed (Born-level) flows and after summing over flows. While only flow-summed predictions are physical, it is not a foregone conclusion that they must be obtained by means of analytically-summed formulae. In particular, when the number of partons involved becomes large, a random sum might be the only viable option, and this can be carried out only by using fixed-flow expressions. We point out that, even in the context of an MC@NLO matching [46], where a local connection between the short-distance cross sections and the PSMC must be established, there is no loss of efficiency, if suitable flow-based NLO-subtraction formulae are employed, e.g. as proposed in ref. [23]. 17 Conversely, a random sum behaves in any case poorly from a statistics viewpoint unless one is capable of understanding a priori which flows give the dominant contributions in N (or, more in general, some selected N k contributions), and of generating them with high efficiency. While the generation of flows is a topic not directly relevant to this paper, since flows are assumed to be given here, it may become important for the application of our results to multi-parton processes, and it is well worth pursuing.

Acknowledgments
SF is grateful to Fabio Maltoni for collaborating during the early stages of this work, and for the endless and very entertaining brainstorming sessions during one of the Covid-19 lockdowns. Conversations with Stefan Prestel are also gratefully acknowledged. BRW thanks members of the Cambridge Pheno Working Group for valuable comments. This work was partially supported by STFC Consolidated HEP grants ST/P000681/1 and ST/T000694/1.

A On gluon-only colour flows and colour factors
In this appendix we discuss some properties of colour flows relevant to gluon-only processes, and their associated colour factors. Where necessary, we shall exploit the following property: of the colour-flow matrix element (see eqs. (2.34) and (5.53)) Tr λ a σ (1) . . . λ a σ (n) Tr λ a σ (1) . . . λ a σ(n) . (A.2) 17 A flow-based representation must also be chosen for the finite virtual contributions, the discussion of which was not within the scope of ref. [23].

JHEP11(2021)045
In eq. (A.1), by the symbol • we have denoted the composition of permutations, to be applied from right to left, i.e.: for any two permutations σ A and σ B .

A.1 Colour flows
A gluon-only closed flow (σ , σ) has an interesting property. Define the permutation where the subscripts → and ← indicate that the corresponding permutation has to be rotated by one position to the right and to the left, respectively, i.e.
The definition in eq. (A. That is: for a given closed flow, each colour loop is in a one-to-one correspondence with a cycle of the permutation defined in eq. (A.4). Furthermore, the labels of such a cycle are mapped onto the labels of the R-elements of the corresponding colour loop by means of the R-flow permutation. By establishing a direct connection between colour loops and ordinary cycles of permutations (which are a much better studied subject) these results are helpful in that they give an alternative viewpoint on colour loops. Note that eq. (A.8) is manifestly invariant under the operations on the r.h.s. of eq. (A.1).

JHEP11(2021)045
There is a further remarkable fact: the idea of colour loops, and in particular their construction given in steps 0-7 of section 3, is essentially identical to that of alternating cycles, which (as far as we know) has been introduced in ref. [47] in the context of the mathematics of permutations, and which is nowadays a powerful tool in the study of genome re-arrangements; a good summary can be found e.g. in chap. 9 of ref. [45]. The essence is that, given a permutation p, one constructs a double-edged graph, denoted by G(p) and called the cycle graph of p, which connects the various elements of p. Such connections are represented as ordered sets of elements, called the alternating cycles of p, whose number is usually denoted by c (G(p)). The interested reader can easily verify that the rules to construct a cycle graph and its alternating cycles are very closed related to the rules given in section 3. An immediate consequence is therefore 19 It is intuitively clear that this equation leads to eq. (A.8); in fact, we shall show below that the analogy is even stricter. The advantage of being able to establish a direct connection between colour loops and alternating cycles lies in the ability to exploit the growing (thanks to the interest in biology) number of mathematics results relevant to the latter. A first example is that of counting. From definition 9.12 of ref. [45]: • The number of n-element permutations with k alternating cycles is the Hultman number [48] S H (n, k).
Therefore, thanks to eq. (A.11), we have In other words, the number of closed colour flows that correspond to a given number (k) of n-gluon colour loops is equal to the Hultman number, computed at (n − 1, k). Fortunately, Hultman numbers are known in a closed form [49] S H (n, k) = n + 2 k n + 2 2 , n − k odd , (A.13) S H (n, k) = 0 , n − k even , (A.14) where [n k] is the unsigned Stirling number of the first kind, which counts the number of n-element permutations with k ordinary cycles. Note that from eq. (A.12) one must have 19 Theorems on permutations usually include cyclic ones. Since any non-cyclic n-element permutation (say, σ) is in one-to-one correspondence with an (n − 1)-element permutations (say, p), a relabelling (e.g. σ(i) = p(i − 1), with 2 ≤ i ≤ n), plus the insertion of a conventional fixed point (e.g. σ(1) = 1), allows one to straightforwardly apply the results of such theorems to colour flows. This is understood in eq. (A.11), and is the reason for the n − 1 argument on the r.h.s. of eq. (A.12).

JHEP11(2021)045
which can indeed be found with eqs. (A.13) and (A.14). Note, also, that eqs. (A.12) and (A.14) confirm that the number of colour loops for any given flow has the same parity as the number of gluons. A theorem by Doignon and Labarre [49] (specifically, theorem 8 there) is directly relevant to the counting of colour loops. It reads as follows: • Let 20 α = (0 . . . n − 1). The mapping is bijective.
As is customary, for the map F one introduces the domain and codomain (eq. (A.16)) and a symbol ( • π≡ F (π)) for its result when acting on a specific element (π) of the domain (eq. (A.17)). As before, c(G(π)) is the number of alternating cycles of π, while c(Γ(ρ)) is the number of ordinary cycles of ρ. S n−1 is the set of n-object permutations; S(n) is the set of factorizations of a given n-cycle (in this case, α) into some n-cycle (ρ) and some permutation (σ). By taking into account the different labelling conventions of ref. [49] and of the present paper, and by observing that The bijective nature of the map F and eq. (A.11) give a formal proof of eq. (A.8) (taking the property of eq. (A.1) into account). We point out that ref. [49] need only introduce ρ, and not the combination on the l.h.s. of eq. (A.21), which in our case stems directly from the construction of colour loops. The consistency of the two procedures follows from the following fact: • For any given n-element permutation σ, the quantity σ → • σ −1 is an n-cycle.
Proof: first observe that cycles must be invariant under a general index relabelling; in other words, label indices by σ(i) rather than by i. Then 21 20 As is standard in the mathematics of permutations, (. . .) is the cycle notation. Furthermore, ref. [49] labels elements starting from zero, as opposed to one as is done here. 21 The ⊕ and symbols are normal addition and subtraction operators, augmented by a cyclicity condition: n ⊕ 1 = 1 and 1 1 = n.

JHEP11(2021)045
whence which concludes the proof. We finally quote two theorems relevant to permutations that could be useful in dealing with closed colour flows. The first is reported as theorem 9.9 of ref. [45], and stems from ref. [50]. It uses the concept of block interchange distance bid(σ) for a given permutation σ, defined as the smallest number of block interchanges that sort σ (i.e. transform σ into the identity). That theorem states that (in the language of this paper) This gives an interesting geometric interpretation to the number of colour loops which, as is explicit in eq. (A.25), decreases linearly with the block interchange distance. 22 The second theorem alluded to before is reported as theorem 9.22 of ref. [45]. Again using the language of this paper, it states that |L| = n + 1 2 Since max(|L|) = n, eq. (A.28) tells one that the larger n, the more unlikely that a flat random generation results in a flow with a large number of colour loops. Conversely, a flat random generation should be relatively efficient in producing highly-suppressed colour contributions.

A.2 Colour factors
We now turn to the discussion of some features of the colour factors. Firstly, the colourflow matrix elements of eq. (A.2) on the diagonal (i.e. σ = σ) can be easily computed by means of secondary flows, by employing eqs. (5.54) and (5.60) since, for any n-gluon colour flow σ, we have The block interchange distance can be defined between two permutations, exploiting the same property as in eq. (A.1). Therefore, eq. (A.25) holds in the case of a generic closed colour flow, as it should for consistency.

JHEP11(2021)045
Thus The rightmost side of eq. (A.30) exhibits the usual factor On top of that, it is easy to show analytically for any n that for some polynomial q(N ). We stress that while the factorization of f N occurs also for off-diagonal elements, this is not the case for the factor N 2 − 2.
As is the case for the diagonal colour-flow matrix elements, for the off-diagonal ones the emergence of the factor f N is not immediately apparent, the r.h.s. of eq. (5.60) being a linear combinations of monomials N p . It is the contributions of the U(1) gluons that add and subtract suitable terms to the leading N n monomial, so that f N eventually appears. This poses an interesting question, namely: the matrix elements summed over flows, with the colour factors computed with traces of λ matrices, manifestly factor out f N because each of the products of two traces of λ does so. On the other hand, as has been shown in eq. (5.66), such matrix elements can be computed by simply counting colour loops, and by ignoring U(1)-gluon contributions. Therefore, in view of the fact that it is those contributions that cause the emergence of the f N factor prior to summing over flows, where does f N come from if eq. (5.66) is employed?
In order to answer that question, we re-write eq. (5.66) more compactly as follows: 34) i.e. we temporarily denote the flows by integer indices. Next, consider one of the contributions to eq. (A.33) stemming from summing over columns after choosing a row -to be definite, take the i = 1 row generally after a relabelling of the flows, so that those with the same number of colour loops are contiguous. The symbol restricts the sum over k to run over values whose parity is that same as that of n. In eq. (A.36) we have used the sequence T H (n, k) that is monotonically increasing for decreasing k T H (n, k) = It is also easy to see that eq. (A.43) is not only a sufficient condition for the emergence of f N , but also a necessary one. Indeed, if the polynomial X of N factors out f N , then X(N = 1) = 0. Thus, from eq. (A.35) which proves the statement above.
These arguments about f N allow us to point out that, while the fundamental and the flow representations are equivalent for flow-summed gluon-only quantities, this is not true in the intermediate steps of any computation. Whether this has any practical consequences depends on which role such intermediate steps play in attaining the final result. An example relevant to the case discussed here is the following: suppose one performs the sum over flows by means of Monte Carlo methods. The difference between the results obtained with the fundamental (where U(1)-gluon contributions are kept non-null) and the flow representation will vanish only in the infinite-statistics limit. At finite statistics, it is likely that the flow representation will induce a larger number of cancellations among different terms, because this is the only manner in which the factor f N can arise. In order to increase the efficiency, then, it is wise to set up a flow-representation computation in a way that incorporates the dual Ward identities exactly and in each step of the procedure.

B Equivalence of closed colour flows and dipole graphs
In this appendix, we sketch the proof of the equivalence between closed colour flows and dipole graphs that underpins eq. (7.90). Although we are unable to prove formally this equivalence for any number of quarks and gluons, we show a possible way to do this. We stress that, in all the cases we have explicitly dealt with, the equivalence has been found to hold true. We also show that, should the equivalence break down for some values of q and n, there is an alternative whereby one still employs eq. (7.90), by means of a minor extension of the definition of dipole graphs.
We start by observing that the radiation patterns of eq. (7.74) can be represented as a matrix, whose row and column indices are the positions of the elements in the colour loop, and whose entries are equal to either +1 or −1, according to the sign in front of the eikonal associated with those elements; the diagonal elements are set equal to zero (since they correspond to self eikonals). This matrix being symmetric, we limit ourselves to writing the block matrices as was done thus far. However, one may treat the colour and anticolour labels of a gluon as a quark and an antiquark. One then proceeds as above and, at the very end, merges the two rows and columns relevant to each gluon (by summing their contents). The resulting rows and columns will sum to two, rather than to one (to four rather than to two in the case of the M (E − ,k ) matrices), with the symmetry and traceless properties obviously unspoilt. One thus obtains again matrices that represent dipole graphs (or generalized dipole graphs, in the case of E − ) that include gluons, except for the following important detail. By treating, during the construction, the colour and anticolour indices of a gluon as independent from each other, one obtains self-eikonals which are not on the diagonal (see e.g. eq. (B.23)). However, such a case is easy to diagnose, because after merging the rows and columns it results in an entry equal to two on the diagonal. When this happens, that gluon has to be removed, and thus interpreted as a self-connected gluon. After the removal, one is left with a dipole graph of dimension smaller than the original one, which is what has already been discussed in section 7.2.2.
Example. We consider the case of the closed flow (γ 0,1 ,γ 0,2 ) in table 1. Equation (B.1) reads (there is a single loop, so the index is omitted) Since n g = 0, we expect the radiation pattern E + to be given in terms of dipoles graphs at least one of which will feature a number of gluons smaller than two (the "missing" gluons being self-interacting ones). Furthermore, we shall see that the radiation pattern E − can also be expressed in terms of dipole graphs, notwithstanding the lack of formal proof that this is possible. However, before going into that, we observe that by applying the mergingline procedure advocated above to the each of the gluons 25 that appear in the matrix in eq. (B.22) we obtain the following: As has been anticipated in the discussion given before, in the matrix on the r.h.s. of eq. (B.27) we see the presence of two elements equal to two on the diagonal, which is equivalent to saying that there are two self-connected gluons. Indeed, it is immediate to see that this matrix is representative of the dipole graph Γ 6 of figure 2, with the corresponding radiation pattern of eq. (6.23).
We now consider the matrices of eq. (B.11), which we construct as explained before, namely by first treating the colour and anticolour of each gluon as a quark and an antiquark, and by merging the relevant lines afterwards. Thus so that with the matrices on the r.h.s. understood to be those obtained by line-merging. The present case being derived from that of a 3-quark-line process, we know (see the comment before eq. (B.16)) that the equivalent dipole graphs must include quark-antiquark eikonals. These are the elements of the matrices on the l.h.s. of eqs. (B.31)-(B.33) whose row and column indices have opposite parities. Since these eikonals cannot stem from E − , they are subtracted by the matrix on the l.h.s. of eq. (B.34), whose elements are all of this kind. We point out that while this implies that the c − coefficients on the r.h.s. of eq. (7.83) do not have the same sign for the 3-quark-line process, this is not the case after line-merging. which coincides with eq. (7.40). This is a direct manifestation of the non-uniqueness of the choices in eq. (7.88), although of limited practical interest (since the present example is relevant to a single closed flow, rather than to the set of all closed flows that one needs to consider in physical applications).

JHEP11(2021)045
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.