Soft gluon evolution and non-global logarithms

We consider soft-gluon evolution at the amplitude level. Our evolution includes Coulomb exchanges and applies to generic hard-scattering processes involving any number of coloured partons. We emphasise the special role played by a Lorentz-invariant evolution variable, which coincides with the transverse momentum of the latest emission in a suitably defined dipole zero-momentum frame. We also relate the evolution algorithm, which was used originally in the derivation of super-leading logarithms, to renormalization group evolution equations that have been encountered recently. Handling large colour matrices presents the most significant challenge to numerical implementations and we present a means to expand systematically about the leading colour approximation.


Introduction
Precise predictions for observables at particle colliders often cannot be achieved without the resummation of logarithmically enhanced contributions to all orders in perturbation theory. Leading logarithms of soft or collinear origin are accounted for in general purpose event generators such as Herwig [1,2], Pythia [3] and Sherpa [4]. Resummation can be based on the direct analysis of contributing Feynman graphs in QCD, and this is the approach we have taken in the past, though effective field theories have also been recognized as powerful tools to organize resummed calculations through a renormalization group evolution [5]. For a large class of observables, which are fully inclusive below some resolution scale in all phase-space regions, no contributions originate from unresolved parton emission due to a perfect cancellation of real and virtual corrections. This eases the resummation procedure, as typically only one or very few emissions need to be taken into account. Observables of this kind are referred to as global observables and include, for example, many of the event shape variables measured at LEP.
On the other hand, observables that are insensitive to emissions into a certain patch of phase space are called non-global and they are subject to contributions from an arbitrary number of emissions, effectively probing QCD dynamics in the blinded phase-space region [6]. This effect makes the all-orders resummation (even of the leading contributions) a much more complicated endeavour, mainly because non-trivial colour correlations become unavoidable. Fortunately, these colour correlations simplify dramatically in the leading colour approximation and they can therefore be approximately accounted for in the general purpose event generators.
In recent years there has been a good deal of progress in developing the technology to tackle non-global observables and, in many cases, go beyond the leading colour approximation [7][8][9][10][11][12][13][14][15][16][17]. Subleading colour contributions have also been addressed in the context of parton shower algorithms [18][19][20][21]. Often, attention has focussed on processes with no coloured particles in the initial state, not least because the simpler colour structure eliminates the need to consider Coulomb (a.k.a. Glauber) gluons. Coulomb gluon interactions are particularly interesting since they have been shown to induce a breakdown in the factorization of wide-angle soft gluon emission from hard-collinear emission [22][23][24][25][26].
In spite of the existing progress, it remains to develop an automated approach to resummation beyond the leading colour approximation for general hard processes. Progress in this direction has been made by Nagy and Soper [19][20][21]. In this paper, our aim is to present a general framework that can be used as a basis for future automated resummations.
To be more precise, we consider algorithmic, recursive definitions of QCD amplitudes for the radiation of many soft gluons and including leading virtual corrections to all orders. Such an approach is at the heart of direct QCD analyses of observables involving many coloured legs, and it was used to identify the aforementioned violations of strict collinear factorisation that occur at hadron colliders.
The present work consists of two main parts, and a number of appendices devoted to more technical details. In Section 2 of the paper we lay down the general evolution algorithm, in a form that is suited to the calculation of multiple soft-gluon contributions to any observable in a fully-differential way. We show how this approach connects to earlier work. We then proceed to reformulate the algorithm in such a way as to make the cancellation of infrared divergences explicit. We also highlight the role of the specific ordering variable first identified in [27].
In Section 3, we focus on the colour structures encountered when solving the evolution equations. We present a systematic procedure to calculate the resulting colour traces, which is based on the colour flow basis and the work presented in [28]. Identifying the leading contribution leads us to re-derive the Banfi-Marchesini-Smye equation [29]. However our formalism is more general and can systematically perform resummation of contributions enhanced by the t'Hooft coupling α s N ∼ 1, along with successive perturbations that are parametrically suppressed by powers of 1/N .
In Appendix A we show how our approach connects to the work presented by Becher et al [14], as well as Weigert and Caron-Huot [7,11,30]. In Appendix B we make the cancellation of infrared divergences explicit for observables that are inclusive below a resolution scale, and in Appendix C we explicitly calculate contributions in a fixed-order expansion. Finally, Appendix D sets up the machinery to deal with the fact that most colour bases are not orthogonal (see [31] for a notable exception).

The general algorithm
Our starting point is the cross section for emitting n soft gluons. At this stage we will assume that it is ok to order successive emissions in energy. This assumption is ok for processes that are insensitive to Coulomb gluon exchanges but appears not to be valid otherwise [29,32]. We have that where H(Q) is the hard scattering matrix, H = |M M| and, in the eikonal approximation, If partons i and j are both in the initial state or they are both in the final state then δ ij = 1 otherwise δ ij = 0. Note that the sum over partons j in the definition D µ i is contextspecific, i.e. it runs over any prior soft gluon emissions in addition to the partons in the hard scattering. Likewise, the colour charge operators, T j , and the Sudakov operators, V a,b , are in a context-specific representation of SU(3) c . The operators A n satisfy the recurrence equation where Θ(E ≤ E n ) is the Heaviside function. A general observable, Σ, can be computed using where the u n are the observable dependent measurement functions and the k i are soft gluon momenta. We suppress dependence on the hard partons and integration over their phase space.
In the above, we should take the limit µ → 0, though we will consider non-zero values in what follows. The carat onk reminds us that ω ij (k) is dependent only upon the direction of the vector k in the ij rest frame. The path-ordering, P, in the definition of V a,b is not actually needed here, because the expression in curly brackets in Eq. (2.2) is independent of the ordering variable, E k . The cross sections in Eq. (2.1) are the general building blocks for any observable and they can be used as the basis for a Monte Carlo computer code to generate partonic events.
We can also write an evolution equation 2 : When the measurement function factorizes, i.e. u n (k 1 , k 2 , · · · k n ) = u(E 1 ,k 1 ) · · · u(E n ,k n ), we can define and Eq. (2.6) becomes This can be re-written as In Appendix A, we show that Eq. (2.10) is the same as the leading-logarithmic accuracy RG equations considered in [7,11,14,30].
Note that, if we are interested in a specific observable that we know is fully inclusive of real emissions with E < Q 0 then we can use the Bloch-Nordsieck cancellation in order to fix µ = Q 0 and integrate the real-emission phase space over E > Q 0 . This is proved in Appendix B.

An infra-red finite reformulation
As it stands, there are fixed-order, infra-red divergences in A n that cancel in the sum over n after integration over the real emissions. In this section we present a reformulation of the algorithm in which the cancellation of infra-red divergences (both soft and soft-collinear) arising in the eikonal approximation is manifest.
General observables may be defined by dividing the angular phase-space into two complementary sub-regions, which we refer to as the "in" and "out" regions, such that the observable is fully inclusive over emissions in the "out" region. If the "out" region is of zero extent then the observable is referred to as a global observable, otherwise it is known as a non-global observable. Phrased this way, we see that all observables are non-global to some extent (since 4π-detectors do not exist).
In order to expose the infra-red cancellation, it is useful to break apart the virtual loop-integral in Γ (see Eq. (2.7)) so that we expose the part which is destined to cancel against a corresponding real emission contribution. To this end, it is useful to consider the measurement function in the soft gluon limit: (2.11) This is quite general, the important thing is that, u(q j , {q 1 , · · · , q j−1 , q j+1 , · · · , q m }) → 1 in the limit that gluon j has zero energy. Generally, we can write (2.12) The set {q} corresponds to all other real emissions and Θ in/out (k) is defined to be unity if k is in the "in"/"out" region and zero otherwise. For global observables, the "out" region is of zero extent, in which case u(k, {q}) = u in (k, {q}). Armed with this we define (2.14) Notice that the virtual gluons are summed to all orders only if they are in the "in" region, i.e. V involves virtual gluons integrated over the "in" region. Since 1 − u(k, {q}) → 0 when E k → 0 it follows that there are no soft singularities in V except those arising from Coulomb gluon exchange. The poles from Coulomb gluon exchange cancel though because they always appear in terms ∼ Tr(V 0,a · · · V † 0,a ). The cyclicity of the trace ensures that the iπ terms can be combined into the unit matrix, since the real part of V 0,a vanishes in the limit E → 0.
We can now re-write the observable as where the operators B n satisfy the recurrence relation (i.e. the analogue of Eq. (2.4)): It should be understood that ,Q , and µ should be set equal to zero. In Eq. (2.15), Φ n encodes the measurement functions for any number of real emissions: The P a n set is indexed by a indicating which of the n gluons are real, and m is the cardinality of the set. For example, if n = 2 then P 1 2 = {1, 2} indicates that gluons q 1 and q 2 are both real and m = 2, P 2 2 = {1} indicates that q 1 is real and q 2 is virtual (m = 1), P 3 2 = {2} indicates that q 2 is real and q 1 is virtual (m = 1), and P 4 2 = {} indicates that both gluons are virtual (m = 0). In this case, and, recalling that u(q, {}) = u 1 (q) and u(q, {q 1 })u 1 (q 1 ) = u 2 (q 1 , q),

etc.
This is the "out of gap" expansion used in [22,23] to derive the super-leading logarithmic contribution to gaps-between-jets.

The ordering variable
So far we have presumed energy ordering in the virtual gluon operators V a,b . However, this is known not to generate the correct super-leading logarithms and instead transverse momentum ordering should be used [32]. Interestingly, and working in the eikonal approximation, but only for gluons that couple to the original hard partons, explicit calculation of all relevant Feynman diagrams (at one loop) reveals that the corrections associated with the exact triple and four-gluon vertices can be largely subsumed into a Lorentz invariant ordering variable, in a potentially simple extension of the algorithm that we described in [27,32]. This intriguing physics may not be so clearly visible in an effective field theory treatment, where the evolution is in an arbitrary renormalization scale. Figure 1 illustrates the key finding of [27]. The upper graph illustrates that the virtual gluon loop integral has its transverse momentum limited by the transverse momenta of the two nearest real emissions. The relevant transverse momentum of a gluon, q (ij) , is defined by its Sudakov decomposition over the momenta p i and p j involved in dipole (ij), and is given by When one of the gluons to which the virtual gluon couples also happens to be one of the nearest emissions, the relevant dipole transverse momentum vanishes. In this case however, dk T k T Figure 1: Illustrating how the 'dipole transverse momentum' serves to limit the virtual loop integration. The index [n] refers to the gluon with momentum q n .
the explicit calculation reveals that the relevant dipole momentum is that of the parent of the parton which couples to the virtual gluon, i.e. parton i in the lower graph of Figure 1. The corresponding differential cross section has a similar structure to Eq. (2.1) but with where it should be understood that q = q (ij) when placed in the argument of I ij . The Lorentz invariant operator I ij (a, b) is given by restricts the region of the angular integration to be the same as in the phase space integral for a real gluon with the same transverse momentum. This can be written Up to the limits on the cos θ integral, I ij is equal to the exponent in the V operator defined in Eq. (2.2). In [27] we presented the form of this operator in 4 − 2 dimensions after integration over the solid angle, i.e.
where c Γ = 1 − γ E , and γ E is the Euler-Mascheroni constant. This expression is accurate up to non-logarithmic terms of order 0 in the real part and order 1 in the imaginary part. When both scales a and b are non-zero, I ij is finite and given by Eq. (2.29) would be identical to the result of Catani and Grazzini [33] (see also [34,35]) if we replaced the factor 1+iπ δ ij in Eq. (2.29) by cos(π )+i sin(π δ ij ). It should be stressed that the calculations in [27,32] were performed only at one loop and it remains to be seen how the improved resummation proceeds to all orders.

Large-N structures
In practical calculations, the colour algebra rapidly becomes intractable after only a few real gluon emissions. To simplify matters, we shall now identify the leading contributions in an expansion in the number of colours. We will work in the colour flow basis [36], which is closely related to the way colour is treated in parton shower algorithms [37][38][39][40].
In the following subsection, we will show that, to leading order in the number of colours, our algorithm gives rise to a dipole-type parton shower and that it reproduces the Banfi-Marchesini-Smye equation [41]. We then turn our attention to setting up a framework to calculate the first subleading-colour corrections.

Colour flow basis
To start with, we collect together some of the key results concerning the colour flow basis. We label the set of basis tensors as {|σ }, and we assign a colour or anti-colour index, c i orc i , to each external leg i of any scattering amplitude. Gluons carry both colour and anti-colour and incoming quarks carry anti-colour. We start to count colour index labels from 1, and choose c i = 0 (c i = 0) if i only carries anti-colour (only carries colour). The basis tensors are labelled by permutations σ of the colour indices and are given by products of Kronecker δ's as  Table 1). The right hand part of the figure illustrates how inner products, i.e. elements of the scalar product matrix, relate to powers of N depending on how many loops are formed after contraction (see Eq. (3.2)). n = n q + n g = nq + n g possible colour lines and n! colour flows (i.e. there are n! basis tensors). Inner products of colour flow basis tensors are given by where #transpositions(σ, τ ) is the number of transpositions by which the permutations σ and τ differ. This is equal to n minus the number of loops obtained after contracting the Kronecker symbols, see the right-hand part of Figure 2. In Figure 2, we show three of the six colour flows that represent the four-parton state on the left, and in Table 1 we specify the corresponding colour and anti-colour indices for each of the four partons (labelled by i).
We also include the binary variables λ i andλ i , where λ i = √ T R ,λ i = 0 for a quark, λ i = 0,λ i = √ T R for an antiquark and λ i =λ i = √ T R for a gluon (T R = 1/2 in QCD). Note that in the figure we use the more compact notation: We express amplitudes as |A = σ A σ |σ , where σ labels the individual basis tensors, and the evolution and traces in colour space can be performed in terms of ordinary complex matrices with elements A τ σ , which relate to the basis independent objects via The coefficients A τ σ are not matrix elements of the operator A since the colour flow basis is not orthonormal. Consequently, we will introduce a dual basis in which Table 1: The colour index specifications for the four partons in Figure 2, which correspond to n = 3. .
We refer to Appendix D for more details on the properties of the dual basis vectors (see also [28]). The scalar product matrix S τ σ = τ |σ has to be considered when evaluating the traces of operators in colour space: The colour charge (or emission) operator associated to each leg i can be decomposed as where the colour-line operators t,t and s are defined through their action on the basis states, i.e.
Colour-line operators and their products, such as t α · t β = t β · t α , are referred to as colour reconnectors in [28]. We note that s · t α = t α · s = 1 and s · s = N 1. Matrix elements involving colour reconnectors are straightforward to compute because of the important property that where R is a general reconnector (see Appendix D). Note that there is no sum over τ since reconnectors constitute a unique map from one colour flow to another. The matrix elements of the colour correlators are for i = j, where (ab) denotes an ordered pair ((ab) = ba if a > b) and τ (ab) denotes swapping the elements a and b in the permutation τ . δ τ σ is zero if the permutations τ and σ are not equal and unity otherwise. The sum over (ab) is rather cumbersome, since each of the four terms can be written without any summation after implementing the colour-flow Kronecker delta. However, this way of writing things ensures that the second and third line in Eq. (3.13) do not contribute if i and j are colour connected in σ. By 'colour connected' we mean 3 Note also that the off-diagonal elements in the matrix representation of T i · T j are non-vanishing only if the permutations labeling the two basis tensors in question differ by at most one transposition. A similar expression can be obtained for colour charges multiplied to the left and right of a colour matrix, A, which corresponds to real emission: The colour lines associated with the emitted particle, n, are labelled by c n andc n , and σ\n denotes the permutation with the entries associated with c n andc n merged and removed, i.e.
. Our aim is to organize contributions to the cross section in terms of a series of leading powers in N , to extract both the large-N limit as well as corrections to it. To this end we introduce the operation where the notation indicates to pick those terms in A τ σ which are suppressed by a factor of 1/N k with respect to the leading power present in A τ σ . Contributions to the trace of A then all yield an enhancement or a supression by the same power of N by virtue of either an explicit suppression in A τ σ or by picking up a subleading element in the scalar product matrix. In other words, if A ∼ A 0 + A 1 /N + ... is an operator in the space of n colour lines, then We are specifically interested in traces originating from soft-gluon evolution: and where λ refers to all other quantum numbers of the emission, and it is understood that

Leading contributions
For the leading contributions of the virtual evolution operators we find where and c.c. here means colour connected, i.e.
i,j c.c. in σ Also, we have defined W For single gluon exchange, with r And the emission contribution is We are now able to compute traces in the leading-N limit. We must sum over diagonal colour flow contributions, i.e. A σσ . For each colour flow, σ, we multiply by N raised to the number of colour lines present in the colour flow σ. Notice that the number of possible colour flows at this level of approximation is equal to the number present at the level of the hard process plus the number of real emissions. Each of the contributions A σσ can then be computed by a set of recursive rules that correspond to working inwards from the outer matrices (multiplied from the left and right) towards the hard matrix in between. The rules are as follows: • For a pair of evolution operators, V n A n V † n , multiply by for each colour connected dipole i, j in σ.
• For a virtual gluon insertion, γ n A n + A n γ † n , multiply by and sum over the dipoles (i, j) that are colour connected in σ.
• For a pair of emission operators, D n A n−1 D † n , combine the dipoles (i, n) and (n, j) in the colour flow σ, leaving behind a dipole (i, j) in the colour flow σ\n, and include a factor (3.32) This procedure is illustrated in Figure 3, where we take the opportunity to show a specific contribution at leading-colour. This graph would contribute to the soft gluon evolution of qq → qq scattering with one emission and one virtual correction.
• If the hard process has been reached, multiply by the square of the corresponding amplitude, |M σ | 2 .

Dipole evolution and the BMS equation
We will now show how the rules of the preceding section give rise to the BMS equation [41].
Our rules apply to a general process with any number of outgoing partons. The algorithmic incarnation of the generalized BMS equation that we present here corresponds to a dipole shower algorithm. The evolution of dipoles is universal, i.e. at this level of approximation the process dependence solely enters through selecting an initial colour flow weighted by the modulus squared of the corresponding amplitude |M σ | 2 . To illustrate how things work out, we will consider the same example as in Section 2.1.1. In this case and ω The evolution with the in-region anomalous dimension contributes a factor per colour flow. Figure 4: The leading-N graphs corresponding to Eq. (3.37) .
The above expressions have a very simple diagrammatic interpretation, illustrated in Figure 4. To simplify the discussion, we consider the case of e + e − scattering, i.e. we take H = 1 N 1. Each double line in the figure corresponds to a Sudakov factor, V E 1 ,E 2 ij , where i and j label the directions associated with the corresponding colour and anti-colour lines. The shaded circles correspond to a factor ω ij (k), and the vertical dashed line indicates the associated energy. The arguments of the Sudakov are also determined by these vertical dashed lines. We can immediately see how the algorithm maps onto a classical dipole shower at leading N . These diagrammatic rules can be used to compute the leading colour contribution to the non-global logarithms: where the hard partons have momenta p a and p b , t i = (N α s /π) ln(E i /ρ) and we used the notation ω i ab = ω ab (q i ). The Σ n can also be obtained by iteratively solving the BMS equation, as we will now illustrate. The BMS equation can be written as follows, and our observable corresponds to with G ab (0) = 1 and t Q = (N α s /π) ln(Q/ρ). To solve the BMS equation iteratively, we will first rewrite it by replacing G ij (t) = V ρ,E ij g ij (t), which gives Putting g which gives the desired result after integrating over 0 < t < t Q . The next iteration gives Σ 2 , i.e. we substitute g ij (t) on the RHS of the BMS equation by g ak (t)g    where we left g (0) ak = 1 explicit for clarity. It is easy to show that V Q,ρ ab g So we see that, at leading N , our algorithm generates the iterative solution to the BMS equation.

Subleading contributions
Subleading colour contributions are substantially more difficult to compute. In this section we present some initial steps towards a systematic approach to including 1/N k corrections to the leading result. Figure 5 illustrates the general structure of the calculation (of which Figure 3 is a specific example). Figure 6 provides an overview of the power counting we use to define successive orders -we hope its interpretation will become clear after the following paragraphs.
There are subleading colour contributions arising from the hard scattering matrix, from the 1/N and 1/N 2 suppressed terms in the real emission operator (see Eq. (3.14)) and the virtual correction operator (see Eq. (3.13)), and from off-diagonal contributions to the scalar product matrix. In the following, we will use the general form of the anomalous dimension resulting from Eq. (3.13), i.e.
Each of Γ, Σ and ρ are of order α s . To compute a correction of order 1/N k we need to consider states σ and τ in Figure 5 that differ by k − l permutations, where 0 ≤ l ≤ k. Then we must determine the 1/N l corrections arising from the soft gluon evolution and from the hard scattering matrix. The leading colour contributions from the virtual evolution operator come from Γ and so are all enhanced by powers of α s N which, owing to the fact that the leading contribution is diagonal, can easily be accounted for to all orders in a simple exponential. This evolution does not result in any difference between the colour structure in the amplitude and that in its conjugate, and it corresponds to the blue boxes in Figure 6. If this evolution is then supplemented by those pieces of the real emission operator that also preserve the identity of the colour structure in the amplitude and its conjugate (such as the example in Figure 3) then we recover the leading-N picture of the last two sections. Subleading colour contributions may result in differences between the colour in the amplitude and that in the conjugate amplitude. To keep track of this, we will count the number of colour reconnections (or transpositions or flips or swings) by which the two colour structures differ. It turns out that pure 1/N corrections can only originate from interference contributions in the hard process matrix. We will ignore subleading colour contributions from this source in what follows, though they could easily be included. The most important subleading colour contributions due to real emission are suppressed by a power of 1/N 2 relative to the leading contribution and they originate as a result of the following three possibilities: (i) two colour flips accompanied by no explicit factor of 1/N (coming from contributions of the type t[· · · ]t); (ii) one flip and a factor of 1/N (coming from contributions of the type t[· · · ]s and s[· · · ]t); (iii) zero flips and a factor of 1/N 2 (coming from contributions of the type s[· · · ]s). See Eq. (3.14) to appreciate the factors of 1/N . We note that real emissions never reduce the number of flips by which the amplitude and its conjugate differ. We will present the explicit rules corresponding to these real emission contributions below but first we consider subleading virtual corrections.
A single insertion of a perturbation Σ σ σ comes with a factor of (α s N )/N relative to the leading contribution and it results in a single flip, and hence an additional 1/N suppression via the scalar product matrix. This flip can undo that induced by a previous real emission of the type s[· · · ]t or t[· · · ]s. However, since this will require the action of a single s operator, it re-introduces the additional factor of 1/N . Thus, both of these fixed-order contributions, when combined with the all-order summation of contributions from Γ, are suppressed by (α s N )/N 2 ∼ 1/N 2 relative to the leading contributions. These contributions are illustrated by the dark orange boxes in Figure 6. A similar reasoning applies to the contribution of a single ρ perturbation (light orange boxes), which contributes at the same order α s /N ∼ 1/N 2 since it generates zero flips.
We finally need to consider two insertions of Σ σ σ combined in such a way that the net number of flips is zero or two. The zero-flip case is clearly proportional to α 2 s and hence contributes a (α s N ) 2 /N 2 correction (green boxes). The two-flip case can also contribute at this order provided it compensates a t[· · · ]t two-flip real emission, i.e. so the net result is that the amplitude and its conjugate differ by zero flips and there is no suppression from the scalar product matrix. These two contributions should therefore be included along with the contributions discussed above. However, contributions from the diagonal below the green boxes lead to a factor of (α s N ) 2 /N 4 , which means they are beyond the next-to-leading colour approximation.
These corrections to soft-gluon evolution can be considered as fixed-order corrections to the leading-N rules, though these need to be extended to include the possibility that the permutations σ and τ need no longer be equal. This means we should update the rules at the end of Section 3.2 as follows.
• For a pair of evolution operators V n A n V † n multiply by V , i.e. include a factor of the virtual amplitude exponentiated for each colour connected pair of legs in σ, and each colour connected pair in τ .
• For a virtual gluon insertion γ n A n + A n γ † n , include a factor w (n) (3.44) • Real emission operators, D n A n−1 D † n , only contribute if the emitted gluon is connected to the same, identically connected, colour line in σ and τ , as otherwise the operation of merging the emitting dipoles would alter the number of flips by which σ and τ differ.
We will now present the rules to compute the first 1/N 2 corrections to the leading colour trace. For the virtual contributions we need to consider the next-to-leading colour approximation for the evolution operator, which (using the notation in [28]) is where the colour flow transition matrix elements can be expressed as and This source of subleading correction contributes when σ and τ are identical or differ by a single flip. In the latter case (one flip), we include a factor of 2 Re Σ (n) στ and the other factor of 1/N comes from the scalar product matrix. In the former case there are three possibilities. Specifically, we may include either one or two factors of Σ (in either the amplitude or the conjugate amplitude) in such a way as to undo the effect of a one or two-flip real emission (see below for the rules for including real emissions), these contribute to the dark orange and green boxes in the figure. Or else we may include a factor of |Σ (n) τ σ | 2 in the case that the two flips (one from each Σ) cancel each other out (green boxes).
At this order we also need to include corrections which are suppressed by 1/N 2 and are proportional to s · s (light orange boxes) in the virtual evolution operator. This term is diagonal in colour and As discussed above, we must also consider subleading corrections to real gluon emission. Recall that the 1/N 2 corrections arise when τ and σ differ by at most two flips. In this case, we are to include contributions where the gluon is emitted off either colour line c i or c k in the amplitude (c i andc k are colour connected in σ\n), and either colour line c l orc j in the conjugate amplitude (c l andc j are colour connected in τ \n). Evolving towards the hard process, we are to combine the dipoles (i, n) and (n, k) in σ, and (l, n) and (n, j) in τ , leaving behind dipoles (i, k) and (l, j). The corresponding factor is which comes from the first two lines in Eq. (3.14). Notice that a potential 1/N contribution arising when σ and τ differ by only one flip vanishes because of the first bullet point in the list above Eq. (3.13), i.e. contributions of the type t[· · · ]t require σ and τ to differ by two flips 4 . If σ and τ differ by one flip and the gluon connects to itself in τ but not in σ, then we should combine the dipoles (i, n) and (n, k) in σ and include a factor of or the corresponding conjugate. This corresponds to the third line in Eq. (3.14). Again a possible 1/N correction arising when σ and τ are equal vanishes because contributions of the type s[· · · ]t require σ and τ to differ by one flip (see the second bullet point in the list above Eq. (3.13)). Finally, if the gluon is connected to itself in both σ and τ , we include a factor of corresponding to the fourth line in Eq. (3.14). Armed with these rules it is possible to go ahead and compute the first subleading colour contributions to the BMS equation. We leave such a phenomenological study to future work.

Conclusions
Accounting systematically for partonic radiation in short-distance scattering processes is of practical importance and theoretical interest. Progress in accurately accounting for this physics has been dominated by coherence-improved parton/dipole shower Monte Carlo programs [2][3][4] though to date these are all limited to leading N , with some subleading improvements [18]. Probably the main challenge in going beyond leading N arises because of the need to include quantum interference effects, which would seem to necessitate an amplitude-level approach. This paper represents our first steps towards the implementation of a general algorithmic approach to amplitude-level parton evolution, which has also been advocated in [19]. We anticipate that numerical results using the technology outlined in this work will be available soon, and we postpone a detailed discussion of the computational methods to a follow-up work.

Acknowledgments
This work has received funding from the UK Science and Technology Facilities Council (grant no. ST/P000800/1), the European Union's Horizon 2020 research and innovation programme as part of the Marie Skłodowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104). JRF thanks the Institute for Particle Physics Phenomenology in Durham for the award of an Associateship. MDA thanks the UK Science and Technology Facilities Council for the award of a studentship. SP acknowledges partial support by the COST action CA16201 PARTICLEFACE, and is grateful for the kind hospitality of ESI at Vienna, AEC at Bern and MIAPP at Munich, where part of this work has been addressed.
We are indebted to Thomas Becher for valuable discussions. Figures have been prepared using JaxoDraw [42].

A The connection with other approaches
In this appendix we show how the colour evolution algorithm defined by Eq. (2.1) relates to the previous work of Becher et al [14] and Caron-Huot [11,30].

A.1 Becher et al.
In [14], the hard process e + e − → 2 jets is considered with the requirement that the total energy emitted outside of cones centred on the two (back-to-back) jets should satisfy 2E out < βQ where Q = √ s. This observable is of the type described by our Eq. (2.10) and, because there are no coloured particles in the initial state, the Coulomb terms can be neglected. Accordingly, in the leading logarithmic approximation they find (see Section 5.2 of [14]).
where δ = tan(α/2) (α is the opening angle of the jets), µ h = Q and µ s = Qβ. Formally, the evolution operator is given by where {n} is the set of l light-like vectors that fix the directions of the final-state partons, and The Θ in (q m+1 ) restricts the emitted gluon with momentum q m+1 to lie inside either the quark or anti-quark jet (defined by the cones around their directions). There is potential for confusion here, because this "in" region corresponds to what we called the "out" region: in both cases we are speaking of the region where there is no veto on real emissions. The subscripts "L" and "R" on the colour charge operators denote that they sit to the left or right of the object upon which they operate. Expanding the exponential in Eq. (A.2) gives rise to exactly the same series as using Eq. (2.4) (see Eq. (5.17) of [14]).

A.2 Weigert and Caron-Huot
We can also translate Eq. (2.10) into the notation and language of [11]. The starting point is to introduce a rotation matrix, U i , for each parton in the hard subprocess and each soft gluon. Operators L i and R i are also defined such that Their commutation relations are inherited from the colour algebra: Now one defines the one-loop kernel where η → 0 and θ ij → 1 gives rise to ordering in energy. Ordering in dipole transverse momentum is obtained with η → − and The corresponding equation for a general observable (i.e. Eq. (2.5) for energy ordering and Eq.(A.14) for dipole transverse momentum ordering) is which is fully differential. Here N is the number of partons in the hard subprocess, the path ordering acts over λ and the colour matrices U i should be independent of this parameter.
In proving the equivalence of Eq. (A.9) with Eq. (2.5) and Eq. (A.14) it is useful to note that here λ n = E n for energy ordering and λ n = q (ij) n for dipole ordering (see below). Up to the Coulomb gluon term, this is equal to the lowest order resummation contained in equations (2.7) and (2.14) of [11]. This is also very closely related to the work of Weigert [7].
Another way to write Eq. (A.9) is via a generating functional: The case λ n = q (ij) n is reminiscent of (but not the same as) the dipole transverse momentum ordering we discussed in Section 2.2 (see Eq. (2.25)). Indeed, Eq. (A.9) can be re-written as a recurrence relation: where A 0 (µ) = V µ,Q H V † µ,Q and V is defined analogously to Eq. (2.2). Formally, each of the gluons should have an energy E < Q and this is imposed via θ injn (q n ). In [27], direct calculation led to θ ij (q) = Θ(p i · (p j − q) > 0)Θ(p j · (p i − q) > 0) and we introduce it here to cut-off arbitrarily high momentum modes. We cannot avoid the long chain of indices because the observable is obtained by integrating over the multi-gluon phase space subject to Θ q where i m , j m ≤ N + m − 1 and N is the number of hard partons. As pointed out in [11,30], choosing the dipole transverse momentum to order the emissions is, ultimately, a renormalization scheme choice in the effective theory, albeit one that has the virtue of making Lorentz invariance manifest. The dipole transverse momenta of successive real emissions are ordered and this set of ordered momenta acts to limit the virtual gluon loop integrals in Eq. (A.14). We note that dipole ordering avoids all collinear poles except for those associated with the very last emission. This is a very attractive feature. The proof proceeds along the following lines: for a given dipole chain, poles come from (p i · q) (p j · q) = 0. But this quantity is proportional to the ordering variable, so the only possibility of it equalling zero is the case of the final emission (when µ = 0).

B On the cancellation of infrared divergences below the inclusivity scale
The aim here is to show that, for observables fully inclusive for E < ρ, we can simply impose E > ρ in the algorithm. This fact follows because (q 1 , . . . , q n+1 ) − u n (q 1 , . . . , q n )) , hence if u n+1 = u n for E n+1 < ρ then one can set ρ as the lower bound on the energy integrals for both real emissions and virtual exchanges.
To prove Eq. (B.1) we make use of the identity where we used the shorthand D 2 (q) ≡ D µ (q) D µ (q) and the Sudakov operator is given in Eq. (2.2). Using Eq. (B.2), we can rewrite the contribution to the observable from n real emissions as where it should be understood that E 0 = Q, i.e. A 0 (E 0 ) = H. Eq. (B.1) trivially follows from this expression by grouping terms that depend on the same trace. In this appendix, we compute the fixed-order expansion of the non-global logarithmic contributions to the hemisphere mass. Apart from checking the correctness of the algorithm, this allows us to confirm that the expansion proposed in Section 2.1 is indeed free from infra-red divergences at each order, i.e. the Σ n are separately finite. We will be very explicit in the hope that it will be useful to see how a calculation proceeds in detail. First we compute the hemisphere jet mass to fixed order, as in Dasgupta-Salam [6]. As they do, we start by computing the lowest order non-global correction to the cumulative event shape where the jet mass is required to be less than ρ. We can do this by using the algorithm with µ = ρ and taking the "out" region to be the region of phase space that does not contribute to the hemisphere jet mass, i.e. it is the wrong-side hemisphere. The "in" region is the complement of this. Note it is only to leading ln(Q/ρ) accuracy that the observable is fully inclusive over gluon emissions with E < ρ. In which case we may write (see Eq. (2.24) with V → V in and − 1 2 D 2 → γ): We have set the Born matrix element equal to the identity (since we are considering a two-jet e + e − event shape) and the factor of 1/N removes the colour factor for the lowest order cross section. To lowest order, Expanding out gives (note the lower case notation, t a , which (only in this appendix) indicates that these operators act on 3-parton objects): which reduces nicely to To compare to [6] we write p a = Q 2 (1, 0, 0, 1) (C.5) (1, 0, 0, −1) (1, sin θ 2 sin φ, sin θ 2 cos φ, cos θ 2 ) and (1 + cos θ 1 ) (1 + cos θ 2 ) (1 − sin θ 1 sin θ 2 cos φ − cos θ 1 cos θ 2 ) ω ab (q 1 ) = 2 sin 2 θ 1 ω ab (k) = 2 sin 2 θ 2 . (C.6) We can do the azimuthal integral using × Tr(t a · t q 1 T a · T b ) (1 − cos θ 1 )(1 + cos θ 2 ) cos θ 2 − cos θ 1 + Tr(t b · t q 1 T a · T b ) (1 − cos θ 2 )(1 + cos θ 1 ) cos θ 2 − cos θ 1 Now do the colour traces, i.e.
C.2 Calculation of Σ 1 and Σ 2 at order α 3 s The same methodology as in the previous subsection can be used to compute Σ 2 at order α 3 s , and Σ 1 at the same order. The sum Σ 1 + Σ 2 then gives the non-global contribution at order α 3 s . The result for Σ 2 is A C F ω ab (q 1 )(A q 2 k a1 + A q 2 k q 1 b − A q 2 k ab ) + In order to facilitate comparison with the work of Delenda & Khelifa-Kerfa [8], we have used the notation A ij ab = ω ab (q i )(ω aq i (q j ) + ω q i b (q j ) − ω ab (q j )) , (C.14) see Eq. (2.2b) of [8]. In addition, Σ 1 at order α 3 s is + ω ab (q 1 ) ω ab (k) (ω aq 1 (q 2 ) + ω bq 1 (q 2 ) − ω ab (q 2 ))) = 2 α s π 3 ln(Q/ρ) 3 3! C 2 A C F ζ(3) , (C.15) whereĀ q i q j ab = A q i q j ab /ω ab (q i ) (see Eq. (3.8) and Eq. (3.11) in [8]). This is in agreement with the result in [8,9], which is written as where O = V|H H|V † , and |H represents the hard scattering process while V accounts for the subsequent evolution (real and virtual). Since the basis is non-orthonormal it is useful to introduce dual basis vectors, |α] defined so that Our interest is to compute the matrix elements [σ|O|τ ] for a specified pair of external states, σ and τ . This we do by evolving inwards from the external states, stripping off soft-gluon operators as we head towards the hard scattering (which lies at the heart of O). The key result in allowing us to accomplish this is the fact that we can write O = LO R where L and R are colour reconnectors, which means Note that there is no sum over β on the right-hand side, i.e. reconnectors constitute a unique map from one basis vector into another. A similar relation holds for L. To make the equations slightly simpler, we will put L = 1 in what follows. We want to calculate where in the final line the state α satisfies R|α = C α R |τ . This state α is unique since the state τ is fixed. In this way, we see that it is possible to recursively strip off evolution operators leaving behind c-number factors and reduced matrix elements in the dual basis.