Parton branching at amplitude level

We present an algorithm that evolves hard processes at the amplitude level by dressing them iteratively with (massless) quarks and gluons. The algorithm interleaves collinear emissions with soft emissions and includes Coulomb/Glauber exchanges. It includes all orders in $N_{\mathrm{c}}$, is spin dependent and is able to accommodate kinematic recoils. Although it is specified at leading logarithmic accuracy, the framework should be sufficient to go beyond. Coulomb exchanges make the factorisation of collinear and soft emissions highly non-trivial. In the absence of Coulomb exchanges, we show how factorisation works out and how a partial factorisation is manifest in the presence of Coulomb exchanges. Finally, we illustrate the use of the algorithm by deriving DGLAP evolution and computing the resummed thrust, hemisphere jet mass and gaps-between-jets distributions in $e^+ e^-$.

1 Introduction Modern day experimental particle physics is often performed at hadron colliders. As an unavoidable consequence, QCD corrections play a large role. Contributions from coloured radiation, when evaluated Feynman diagrammatically, diverge at multiple points in the phase space. When regularised and cancelled, the divergences may leave behind large -1 -logarithms. The accurate inclusion of logarithmically enhanced corrections is of importance to both the theoretical and experimental communities. Historically there have been two main approaches to dealing with QCD radiative corrections: resummations and parton showers.
Resummations look to re-organise the perturbative expansion by classifying the large logarithms and then summing the perturbation series such that the most dominant logarithmically enhanced terms are included. Towers of logarithms may be further simplified by making the leading colour (LC) approximation. The re-organised expansions are referred to by their logarithmic accuracy; leading log (LL), next-to-leading log (NLL), etc. This procedure has recently been further formalised by work in soft-collinear effective field theories [1][2][3][4]. From this perspective, resummations are renormalisation group flows that evolve 'safe' perturbative predictions into regions of phase space where perturbative expansions would be otherwise 'unsafe'.
In contrast, parton showers may be thought of as providing an all-purpose approximation to the resummation procedure. Modern parton showers generate an evolving, classical system of partons whilst cleverly encoding quantum interference effects (made possible by working in the LC approximation). The majority of currently available parton showers claim LL accuracy using the LC approximation [5][6][7][8][9][10][11]. The quest to better understand the data from the LHC is a major driver for increasingly precise parton showers. At present, there is a growing list of phenomena that parton showers do not encapsulate. This includes effects sub-leading in colour, some non-global logarithms [12], Coulomb/Glauber exchanges, super-leading logarithms [13][14][15] and the violation of QCD coherence (or collinear factorisation) [16,17]. Moreover, recent fixed-order studies have cast further doubt on the accuracy of modern parton showers. It has been shown in [18] that the Pythia [10,11] and Dire [8] showers suffer from both incorrect next-to-leading logarithms at leading colour and incorrect contributions from sub-leading colour (NLC) at LL. Although these showers never claim NLL or NLC accuracy, the findings of Dasgupta et al questions the fruitfulness of attempts to extend conventional parton showers beyond LL and LC in general. In recent years, there has been movement towards finding new constructions for partons showers; constructions more suited to including NLC or NLL corrections [9,[19][20][21][22][23][24][25][26]. However, as of yet, success has been limited.
The algorithm we present here aims to provide a framework for the development of future parton showers, enabling them to be systematically improved. We hope it will also help make more rigorous the link between resummations and parton showers. Our starting point is the soft-gluon evolution algorithm explored in [25], which we refer to as the FKS algorithm. The evolution generated by the FKS algorithm is systematic to all orders in colour and it accounts for the leading soft logarithms. The FKS algorithm was originally used to derive the super-leading logarithms that may occur in hadron-hadron collisions [13,14]. It has been analytically verified for a general hard process dressed with up to two soft real emissions and one loop [27,28]. It has also been shown to generate the BMS equation [29] (it presumably also includes the NLC corrections to it) and it correctly accounts for the leading non-global logarithms for various observables [25]. The main goal of this paper is to improve the FKS algorithm by including collinear emissions, -2 -spin dependence and kinematic recoil. The algorithm we present is Markovian and can be solved iteratively, making it well suited for use as a parton shower.
The remainder of the paper is organised as follows. In the next section, we introduce the algorithm in a form we refer to as variant A, in which we interleave soft and collinear emissions. Variant A has the virtue of being a simple extension of the FKS algorithm, though it suffers from unnecessarily complex colour evolution in the soft-collinear sector. It also suffers from the fact that we cannot uniquely identify a parent parton in the case of soft-gluon emission, which complicates the issue of longitudinal momentum conservation. We are thus motivated to re-cast the algorithm in a more convenient form, which we refer to as variant B. Specifically, in variant B we manipulate the colour structures of variant A to isolate the full collinear splitting functions, after which we are able to implement longitudinal momentum conservation in a simple way. We also spend some time illustrating how recoils may be included in both variants, though this will only be relevant beyond the LL approximation. As it stands, either variant A or B could be used to create a fully functioning parton shower, though B will be computationally more efficient. In Section 2.5 we present a manifestly infra-red finite version of the algorithm. This reformulation is particularly useful for the resummation of specific observables, though it is not so well suited for use as a general purpose parton shower. This is because the infra-red singularities are regularised by the explicit inclusion (and exponentiation) of a measurement function.
The second half of the paper is devoted to issues of collinear factorisation and to providing examples to illustrate how the algorithm is used. In Section 3 we discuss the factorisation of collinear physics from soft physics. We start by considering the case when Coulomb/Glauber exchanges are turned off (such as would be the case in e + e − collisions). After this we discuss how Coulomb exchanges can be introduced one-by-one. We will see that collinear factorisation occurs below the scale of the last Coulomb exchange. This discussion shows consistency between our approach and the proofs of collinear factorisation by Collins, Soper and Sterman [30,31]. After this, we show how DGLAP evolution for the parton distribution functions emerges [32][33][34]. We finish the paper by illustrating the use of the algorithm; by calculating the thrust, hemisphere jet mass, and gaps-between-jets distributions in e + e − . We leave an extensive discussion of spin correlations to an appendix.

The algorithm
In this section we present the algorithm. It is Markovian and interleaves soft emissions and virtual corrections with collinear emissions and virtual corrections, see Figure 1. Successive real emissions are strongly ordered in an appropriately defined transverse momentum. We will present two variants of the algorithm, which we refer to as A and B. The two differ only in where we put the soft-collinear emissions: in A they are in the soft sector and in B they are in the collinear sector. The second approach allows us to exploit the colourdiagonal nature of collinear emissions and it makes kinematic recoil more straightforward to implement. Variant A has the virtue that it is an almost trivial extension of the purely soft evolution presented in [25]. We present both A and B with the momentum mappings after each real emission parametrised into two initially unspecified functions. This is so the -3 -M(Q)| Figure 1. A term contributing to evolution of the conjugate amplitude (it contributes to A 9 ). Red dashed lines represent the emission of soft gluons and collinear partons are represented by blue dotted lines. Loops (Sudakov factors) have been neglected to avoid clutter. We draw all particles heading to the right, away from the hard process, including incoming particles. In contrast, evolution of the amplitude will have all particles drawn heading to the left and away from the hard process.
algorithm is able to accommodate partonic recoil. Later, in Section 2.4, we discuss specific examples of recoil in action. For processes with coloured incoming partons, the algorithm should be convoluted with parton distributions functions. We leave a full description of how to do this to Section 4.
Before plunging in, we should explain the theoretical basis for what follows. Our algorithm is based on Feynman diagram calculations [14,25,27,32,[35][36][37] and, in its present form, captures all of the logarithms associated with the leading amplitude-level singularities. Therefore the algorithm captures leading logarithms from wide-angle soft emissions, hard-collinear emissions and simultaneously soft and collinear emissions. This means the algorithm is guaranteed to capture only the most leading logarithms at crosssection level. That said, it is also able to capture the leading single non-global logarithms, even if the global part is double-logarithmic, as is the case with the hemisphere jet mass for example (see Section 4). For any process involving incoming hadrons or measured outgoing hadrons, the single logarithms from DGLAP evolution are recovered as well (i.e. parton distribution function and by a simple extension fragmentation function evolution). We believe our framework to be sufficiently flexible that we can, in the future, extend it beyond the LL approximation.

Parton branching with interleaved soft and collinear evolution (A)
The algorithm evolves a hard-scattering matrix, H(Q; {p}), which is defined at some hard scale Q and is a function of the hard-particle four-momenta, {p}. It does so by dressing with successive soft and/or collinear real emissions and virtual corrections. H(Q; {p}) is a tensor in the product space of colour and helicity 1 , defined as H(Q; {p}) ≡ (|colour ⊗ |spin ) ⊗ ( colour| ⊗ spin|). The hard-scattering matrix is defined so that Tr H(Q; {p}) is the hard matrix element squared, summed over colour and spin 2 . Successive real emissions are added via 'rectangular' operators, which act as a map increasing the dimension of the representation of SU(3) × E (2) in which H(Q; {p}) resides. The virtual evolution operators are 'square' and preserve the representation of H(Q; {p}). Specifically, (2.2) At each step, the emission operators (D n ) add one new particle, of four-momentum q n , to the set {p} n−1 , to produce the set {p} n . We use p j ∈ {p} n = {P 1 , P 2 , · · · P n H , q 1 , · · · q n } to denote the momentum of the j th parton and 1 < j < n H + n, where n H is the number of partons associated with the original hard process and n is the number of emitted partons. Hidden in the emission operators is a map from {p} n−1 to a new set, {p} n−1 . The difference between these two sets is determined by the way we implement energy-momentum conservation (i.e. the recoil prescription) and it is why there is an extra integral over p i (it is not a phase-space integral). The virtual evolution operators V a,b encode the loop corrections. To avoid cumbersome notation we write {p} n = {p} n−1 ∪ {q n } is the set of n momenta including the last emission, q n . We have not yet defined the ordering variable, q i⊥ ; we will do that shortly. A generalised observable Σ, with measurement function u n (q 1 , ..., q n ), is then given by 3 where dΠ i is the phase-space for the i th emission (see below). µ should be taken either to 0 or to the scale below which the observable is inclusive over all radiation. The virtual 2 We may also choose to include averaging factors, a flux factor and the hard process phase-space, so that it is then the hard-process differential cross section. 3 For fixed Born-level kinematics. Generally the measurement function will depend upon the hard process momenta Pj, which we do not show explicitly.
-5 - where i and j run over all external legs (those from the initial hard process and also previous emissions in the evolution).δ ij = 1 if both partons i, j are incoming or both outgoing and δ ij = 0 otherwise. θ ij (k) = Θ(p i · p j − k · (p j + p i )) and ensures that the phase space of the integration corresponds to that of a real gluon. Likewise, the z integral is over the range 5 which can be expressed via a single theta function The vectors p and n = (1, n) will be defined shortly: to LL accuracy p = p i and α = 1. υ i , υ ∈ {q, g} label parton species.ῡ = g in all cases except when υ i = g and υ = q, then υ = q. P • υυ i is the υ i → υ hard-collinear splitting function and it is defined in Appendix A along with the conventions we use for helicity states and antiparticles. R soft ij (k, {p}) and R coll i (k, {p}) are concerned with the recoil prescription and are included to preserve unitarity, they are defined in (2.11) and (2.12) below. k (ij) ⊥ and y are the transverse momentum and rapidity in the ij zero-momentum frame. To make the (unitarity) link to the real emissions more explicit, we choose not to use the substitution The real-emission operator is built using two operators: such that D i acts as The path ordering ensures that the operators are ordered in k with the largest to the right. 5 We specify the range corresponding to emission off a final state particle, for emission off an initial state particle exchange pi →pi.
-6 -j again runs over all external legs and i labels the emitted parton. S i generates soft emissions and C i hard-collinear emissions. The symbol ∆ ij is defined so that ∆ ij ∆ ik = δ jk and δ final j (δ initial j ) is unity when parton j is in the final (initial) state and zero otherwise. P ij are the amplitude-level hard-collinear splitting functions and are defined in Appendix A. The splitting functions encode DGLAP evolution [32][33][34] including the spin-dependence.
T g i is a basis independent colour charge operator. We have indexed each T g i with the leg on which it acts, i, and by whether it corresponds to the emission of a gluon or not (i.e. the index q refers to a g → qq splitting). S s i updates the helicity state by adding the helicity of the emitted parton, s i . The operators S and T are also defined in Appendix A.
In the soft sector we have introduced an auxiliary vector m. It is uniquely determined, but only at cross-section level, since we require q ⊥ ) 2 , which corresponds to choosing m to lie in the direction of j and m (the corresponding vector in the conjugate amplitude) in the direction of i. It is only ever this combination that appears at cross-section level. In the collinear sector, the momentum fraction z i is defined by (see Figure 2): z i =p j · n p · n for final-state emissions and z i = p j · n p · n for initial-state emissions, (2.9) where the light-like four-vector n satisfies n · q (j n) i⊥ = 0. The light-like four-vector p satisfies p · q (j n) i⊥ = 0. Neglecting terms suppressed by the transverse momentum of the emission (which is permissible in the LL approximation) we may take p = p j for final-state emissions and p =p j for initial-state emissions, in which case z i is the light-cone momentum fraction. The precise definition of p is dependent on the recoil prescription, as we illustrate in Section 2.4. Now we can define the ordering variable, i.e. the definition of a and b in the Sudakov operator V a,b . We use transverse momentum ordering, where the transverse momentum should be defined by the parent partons of the emitted parton. Doing this means that we really ought not to sum over partons in (2.7) and we should replace (2.2) by where D jn n is defined by The ordering variable is then q n⊥ = q (jnj n ) n if the emission is soft or q n⊥ = q (jn, n) n if the emission is collinear. At LL, this choice of ordering variable is somewhat arbitrary but being a transverse momentum it is able to generate the super-leading logarithms correctly [27]. That said, it is not equivalent to the ordering indicated by the results in [27,28], which is based on fixed-order Feynman diagram calculations. We have not yet figured out Figure 2. Defining the kinematics: (a) q i is emitted off an incoming leg; (b) q i is emitted off an outgoing leg. a way to implement the latter ordering to all orders. In the remainder of the paper, we will use the simpler (though potentially misleading) notation of equation ( and R coll * ij R coll ij , encode the maps that implement energy-momentum conservation. As the algorithm proceeds, each R soft ij and R coll ij will always collect into the pairs just given. The functions only ever appear singularly to aid book keeping 6 . The recoil functions should be constructed out of delta functions and algebraic pre-factors relating the momentum from the current step in the algorithm, {p}, to the momentum that will be carried forwards to the next step of the algorithm, {p} and In the LL approximation, R soft jj = R coll j = 1. S i generates soft emissions and one might suppose that a suitable choice of recoil (to LL accuracy) is We will shortly see that things are not quite so simple, and that this requires modification. C i generates hard-collinear emissions, however only the longitudinal component of the recoil is hard. Therefore, in the LL approximation, we may implement recoils in the collinear sector via (2.14) Finally, the phase-space is included via 6 Singular definitions would require R functions to contain integrals of delta functions which independently evolve momenta in the amplitude and similarly for R * in the conjugate amplitude. The external momentum integrals, k d 4 p k , would the be used to force the two separately evolving momenta to coincide. Giving these definitions provides an entirely unnecessary extra complexity.

-8 -
The pre-factor has been included to simplify the definitions of S i and C i , as well as to make each term in the algorithm dimensionless and keep explicit dependence on the ordering variable in D i . To simplify the notation, in (2.15) and elsewhere, we will drop the dipole labels on transverse momenta. It should be clear from the context which partons are intended. In the case of (2.15), it means we should use the transverse momentum defined by the parent parton and the vector n in the case of collinear emissions or q (ij) ⊥ in the case of soft emissions. It is often useful to note that where y is the rapidity in the frame defining q i ⊥ . Using these last two relations the link between real emissions and virtual corrections is clear, i.e. the square of the emission operators D † i D i dΠ i is, for e + e − collisions 7 , equal to minus twice the real part of the exponent in (2.4).
Using the naive recoil prescription of (2.13) and (2.14), the array of parton momenta gets modified after a collinear emission, generated by P ij , but not after a soft emission (except to add one new soft gluon of course). Specifically, this means acting with P ij R coll ij maps p j →p j = z i p j + O(q ⊥ ) (for final state partons), and a parton with momentum q i = (1 − z)p j + O(q ⊥ ) is added. As usual, p j is the momentum of parton j prior to the action of P ij R coll ij . A more careful treatment of momenta is not required to reproduce the leading logarithms for many observables. However, any observable dependent upon parton distribution functions or fragmentation functions will be incorrectly calculated because this naive recoil prescription does not reproduce DGLAP evolution. This is because the terms with soft-collinear poles are handled in the 'soft side' of the algorithm and do not conserve longitudinal momentum. This manifests as DGLAP evolution with an incorrect plus prescription, i.e.
as the soft poles have been removed from the hard-collinear splitting functions defining P ij . On the flip side, the algorithm works well for event shape observables in e + e − collisions. We will refer to the framework in this section as variant A of the algorithm. Within variant A, this problem could be solved by implementing a universal recoil for all emissions, soft and collinear, i.e. (2.18) We will not consider universal recoils in this paper and will instead solve this 'plus prescription problem' another way; by putting the soft-collinear emissions in the collinear sector of the algorithm. Doing this will lead us to variant B of the algorithm. In Section 2.4, we will use the insight gained from formulating B to show how to solve the plus-prescription problem within the framework of A.

Parton branching using complete collinear splitting functions (B)
Soft-collinear poles can be exchanged reasonably simply between eikonal currents and collinear splitting functions. We will now define variant B of our algorithm, which restores the soft-collinear poles in the collinear splitting functions and removes them from the eikonal currents. This is a good thing to do for two reasons: 1. Collinear evolution is generated by unit operators in colour space. Making this manifest for the soft-collinear poles simplifies the colour evolution of the algorithm.
2. Putting the soft-collinear poles into the collinear 'side' of the algorithm simplifies the recoil prescription because every collinear emission has a uniquely identifiable parent.
Variant B is very similar in form to variant A: The form of S i is the same as in A and the Sudakov changes as The only changes relative to variant A are the appearance of f jj (q i , {p}, {p}) and F ij (k, {p}), which specify the prescription for the subtraction of soft-collinear poles from the eikonal currents, and the replacement of P ij with P ij . Explicit dependence on P ij in C i means that C i OC † i now contains the full spin-dependent DGLAP splitting functions [38]. Unitarity requires that The functional forms of f jj (q i , {p}, {p}) and F jj (q i , {p}) are uniquely fixed by the choice of R coll ij and R soft ij once we have fixed P ij . Specifically, we can derive variant B from A by adding and subtracting a function: where we have labelled each operator with a superscript indicating which variant it corresponds to and where O is some general operator in colour and spin. The subtraction term was constructed so that and after some manipulation is equal to q ipj and [q ipj ] are Weyl products in the spinor-helicity formalism [39]. Also note, s i Os † i is equal to the collinear limit of To see the equality we express polarisation vectors using spinor products, where σ µ ∓ = (1, ∓σ 1 , ∓σ 2 , ∓σ 3 ) T are vectors of Pauli matrices and n is an auxillary light-like vector (best chosen to be either p j or p j ).
To complete the definition of variant B we must compute f jj (q i , {p}, {p}) and F jj (q i , {p}). Note that s † i s i is proportional to the unit operator in colour and helicity. After taking the trace over helicity space, (2.23) leads to (2.27) We can use colour conservation to factorise the colour operators and simplify the second term on the right-hand side, i.e. (2.29) -11 -For a universal recoil it is possible to employ colour conservation and write which enables us to re-write (2.28) as . (2.31) In the case of a universal recoil prescription, the effects of the recoil can be factorised out of the emission operators and into a redefinition of the phase space measure. Recoil schemes that may be universal include the more 'true to Feynman diagrams' global prescriptions which put momenta of partons higher in the chain of emissions off-shell (e.g. see [40] and references therein) and schemes which democratically share recoil across a jet or every parton in the shower. Depending on their implementation, such schemes can be universal since they globally redistribute momentum across the whole event as a n → n + 1 parton processes. We leave the specification a universal recoil scheme to future work. For now we re-iterate that it is only when considering effects beyond LL that (2.31) and (2.28) differ. Our implementation of recoil is not unique, and it remains to be seen (by performing analytic calculations at NLL and beyond) the extent to which we will eventually be able to capture the salient sub-leading logarithms in the framework of our algorithm. A slightly different approach would be to start with variant B (recall we started from variant A above) and assume (2.31) holds true. Variant A could then be constructed but it would now include subtraction functions akin to f jj . In the case of universal recoils, none of this matters of course.

Collinear subtractions and the ordering variable
Before moving on, we'd like to present a slightly more general approach to subtracting the soft-collinear contribution. This calculation will shed some light on the role played by the ordering variable. We start by writing 8 (2.32) which holds for a general definition of the ordering variable, K 2 (p i , p j , k). In order to isolate the collinear divergence, we should first expose, and factor, the soft divergence. To do this, it is sufficient to consider any scaling which is linear in the emitted gluon's momentum components, such that we can re-write (2.33) 8 We ignore ignore hard-collinear corrections and the effects of recoil in this section.
-12 -where n i = q i /(S · q i ), n = k/(S · k) and S is any time-like four-vector, which we choose to satisfy S 2 = 1. The soft divergence is now isolated from the eikonal term, which is singular only in the collinear limits n i,j · n → 0. The collinear divergences can be subtracted. We want the ordering variable to become independent of the other parton's direction in the collinear limit, such that the entire collinear divergence can be moved into a jet factor that is trivial in colour space. We choose to re-write the virtual evolution as ln V ab = ln W ab + ln K ab , where and colour conservation can now be used to obtain This factor contains the ordering variable in terms of a single emitter direction, which is the limiting case of the dipole-type definition in each collinear limit, i.e.
Given the Lorentz invariance of the virtual evolution and the integration measure we can choose S = (1, 0). In the case of energy ordering, we obtain the following for the subtracted soft evolution: where the angular integral can be performed using the same integral that gives rise to angular ordering. And for the collinearly divergent factor: There is no need for θ ij since this simply enforces that the emitted gluon should have energy smaller than 1 2 p i · p j in the ij zero momentum frame, which is automatically satisfied since a < E < b. Now let us consider the case of transverse momentum ordering. This can be implemented through where the similarity sign refers to any function which approaches unity in the limit p·k → 0. Making the minimal choice, the full evolution becomes . This is the same as the subtraction prescription we introduced in the last section, with the only difference being that the lower limit on the z integral is (approximately) equal to α in that case. Finally, using colour conservation we can compute ln V ab − ln K ab and get where the second line illustrates the equivalence with the subtraction scheme presented in the previous section (recall we are ignoring recoil in this section). The approximatelyequal-to sign is because we neglect terms suppressed by powers of k 2 ⊥ . As expected, this finite term is the same as the energy ordering case in equation (2.36). The form factor exp(ln W) captures all of the truly wide-angle soft-gluon physics and is essentially the same as the fifth form factor introduced by Dokshitzer & Marchesini [41].

A local recoil prescription
Next we will show how a more sophisticated recoil prescription (than (2.13) and (2.14)) can be implemented. The recoil we choose is based on the one in [42,43], but extended to work with colour off-diagonal evolution. The dipole recoil is itself based on Catani-Seymour dipole factorisation and furthers the work in [44] so that recoil can be implemented in a dipole parton shower. As a result, this recoil scheme shares similarities with the spectator recoil prescriptions used in modern dipole showers such as Pythia and Dire [8,11]. The idea is not to present a definitive recoil prescription but rather to illustrate how one can be implemented in our algorithm. To that end, we calculate F ij and R coll ij . We also provide -14 -a short discussion on the successes and limitations of the prescription. We then go on to show that, at LL, the recoil prescription can be reduced to the naive recoil prescription when implemented in variant B (but not variant A).
To keep things as simple as possible, we will consider the dipole recoil scheme in the case of only coloured final-state partons. The extension to coloured initial-state partons is straightforward and can be found by following Section 3.2 of [43]. First we will summarise the dipole recoil for colour-diagonal evolution. It works by adding a spectator particle to the standard description of a 1 → 2 collinear splitting (p j →p j , q i ). This spectator particle absorbs the recoil from the splitting, which would otherwise put p j off-shell. The spectator particle has a second function: to specify the frame in which the transverse momentum of the emission is computed. In [43] it was shown that one can obtain the correct colourdiagonal evolution by choosing the parton that is colour connected to parton j as the spectator. We will denote the momentum of the spectator parton by p j LR (the reason for the LR subscript will become clear). The Sudakov decomposition is This now defines a 2 → 3 splitting (p j , p j LR → q i ,p j ,p j LR ) in which q i is emitted collinear to p j . The prescription is momentum conserving, i.e.
and it ensures that all particles are on-shell at each stage in the evolution. One can check that this Sudakov decomposition does not change the functional form of the leading-order collinear splitting functions. Comparing to (2.9), we see that p and n are now fixed: p = p j and n = p j LR . Working in the LC approximation, the effect of this prescription amounts to a correction to the single-particle emission phase space [43], i.e.
This correction contributes soft-collinear NLLs and hard-collinear NNLLs [43]. The dipole recoil prescription was developed for a leading N c shower and as such is not completely sufficient for our purposes. That is because, beyond the LC approximation, the left evolution (of the amplitude) and right evolution (of the conjugate amplitude) are independent, which means they can evolve to produce colour off-diagonal terms. These are terms for which the parton j is colour connected to different partons in the left and right evolution. In such a case p j LR cannot be defined. Instead, we must introduce parton p j L , which is the colour connected parton to j in the left evolution, and parton p j R , which is -15 -colour connected to j in the right evolution. We will now construct a recoil prescription that extends the dipole recoil to include colour off-diagonal terms but collapses back to the dipole recoil in the LC approximation.
To begin, we define a Sudakov decomposition for a 3 → 4 splitting (p j , p j L , p j R → q i ,p j ,p j L ,p j R ). We aim to construct the decomposition so that recoil is shared equally between p j L and p j R . We also wish to leave the collinear splitting functions unchanged. Finally, we also require all partons involved to be on-shell. These constraints are fulfilled by the decomposition: where δ j L ,j R is the usual Kronecker delta symbol. Note that p j +p j L +p j R = q i +p j +p j L +p j R , and so momentum is conserved in the 3 → 4 splitting. Also note that when j L = j R = j LR , i.e. the emission is colour-diagonal, this reduces to the dipole 2 → 3 scattering with p j + p j LR = q i +p j +p j LR . Finally, note that every new term relative the dipole recoil is accompanied by a factor γ, which is two orders higher than the leading terms in the collinear limit and one order higher in the soft limit. Using this decomposition, the recoil prescription for collinear emissions is where the Jacobian, J ij , can (in principle) be evaluated using One factor of z i J ij ensures that the integral over the delta functions is correctly normalised whilst the additional factor of J ij encodes the recoil corrections. This is the factor that was absorbed into the phase-space in [43], i.e.
We can extend this prescription to the soft sector using Catani-Seymour dipole factorisation [44]. The dipole factorisation provides a unique way to split the parent dipole of a soft emission into two halves, identifiable by their separate collinear poles: This provides the means to implement a local recoil using the parton contributing to the collinear pole in each half of the dipole. Thus we can write From this, R soft ij and R coll i can be evaluated: We can now go ahead and determine the subtraction functions used to define variant B. Using (2.28) and (2.29) we get (2.52) Before moving on we want to comment on the discussion in [18], which shows that the dipole recoil scheme, as implemented in a dipole shower, fails at the level of the NLL even at LC due to incorrectly assigning the longitudinal recoil after multiple emissions. It remains to be seen whether this is also true in the scheme discussed here.

A LL recoil prescription
We can now consider constructing a recoil prescription where we only keep the parts contributing at LL. Firstly note that in the strictly LL soft limit We can use this fact with variant B restricted to LL accuracy and find . (2.53) Using the naive recoil prescription, (2.55) With this, variant A also captures the correct DGLAP evolution.

A manifestly infra-red finite reformulation
It is possible to re-cast both variants of our algorithm such that the IR divergences, other than those renormalised into parton distribution and fragmentation functions, explicitly cancel at each iteration of the algorithm. We will demonstrate this for variant A, though the procedure is pretty much identical for variant B. Our method closely follows that in [25]. We begin by expressing a generalised measurement function in the soft and collinear limits as follows (2.56) and u m (q 1 , ..., q m ) where u(q j , {q}) → 1 as j becomes exactly soft or collinear. For many observables, it is possible to further pull apart the measurement function by defining an 'out' region, where there is total inclusivity over radiation. For such observables we can write where Θ in/out (q j ) = 1 when q j is in the in/out region and zero otherwise. For a global observable, the out region has zero extent and so Θ out (q j ) = 0. First we define (2.59) After a simple path-ordered operator expansion, where We define δ R n = 1 when parton n is real and δ R n = 0 when parton n is virtual, and similarly δ V n = 1 − δ R n . Φ n (q 1 , ..., q n ) is a measurement function on the phase-space of real particles. We refer to [25] for its precise definition and here just present an illustrative example: (2.63) -19 -Written in this form, each B n is explicitly infra-red finite provided the measurement function is infra-red-collinear safe, and that the evolution is not convoluted with parton distribution or fragmentation functions. In this case µ can be safely taken to zero. In the case that the evolution is convoluted with parton distribution or fragmentation functions, collinear poles remain (one for each hadron). These poles are removed by renormalisation of the parton distribution functions or fragmentation functions, generating their µ dependence. Finally, note that for recursively-infra-red-safe continuously-global observables [9] in e + e − collisions B n = 0 for n ≥ 1 at single-log accuracy. If the observable depends on fragmentation functions the n ≥ 1 contributions give rise to the DGLAP evolution of the fragmentation functions (see Section 3).

Collinear factorisation
Up to this point we have been interleaving the emission of soft and collinear partons to build up the complete amplitude. As is well known, it is possible to factorise the collinear emissions into the evolution of parton distribution and fragmentation functions. In this section, our aim is to explore collinear factorisation within the context of our algorithm.
The plan is as follows. First, we will derive the factorisation of collinear physics into jet functions; one for each parton in the initial hard process. At first we do this ignoring the presence of Coulomb exchanges. This result is sufficient to derive DGLAP evolution (which we do in Section 4). After this, we go on to construct the complete factorisation of collinear physics into jet functions on every hard or soft leg. Finally, we use a path-ordered expansion of our Sudakov operators, V a,b , to re-insert Coulomb exchanges one-by-one into the previous results. The result is that collinear evolution below the scale of the last Coulomb exchange can be factorised. The outcome of which is the general factorisation of collinear poles into parton distribution functions, as anticipated after the work of Collins, Soper and Sterman [31]. We provide diagrammatic proofs where possible and only sketch in the text the algebra that is going on behind the scenes. Throughout this section we leave aside the recoil functions R soft , R coll , and the integrals corresponding to the momentum maps between each iteration of the algorithm. This is to reduce the length of the algebra that remains, and it is certainly valid to LL accuracy since tracking longitudinal recoil is sufficiently simple. We also drop the inclusion of measurement functions since they too have no affect on the discussion. That said, in Sections 3.1.1 and 3.1.2, we present a summary of the results with all of these functions re-instated (for both of variants A and B).

Factorisation on hard legs without Coulomb interactions
The main result of this subsection is the factorisation of collinear physics into jet functions; one for each leg emerging from the hard scatter. We do this with Coulomb gluons removed (δ ij = 0) and will discuss their re-introduction in Section 3.3. The following manipulations can equally well be performed using either variant A or B of the algorithm. For concreteness -20 -  we will use variant B whenever an operator needs to be given an explicit definition 9 . Let us begin by simply stating the result: (3.1) Figure 3 illustrates what is going on diagrammatically (it shows a contribution with n = 8, m = 5 and p = 1). The collinear evolution operators for hard legs, which provide an operator description of a jet function, are constructed iteratively according to In the case of variant A, for the most part, all that must be done is exchange Pij and Pυ i υ j with the overlined versions Pij and Pυ i υ j .
In both operators j sums only over hard partons. The circle operation, •, indicates the sharing of partons between Col m (µ) and Col † m (µ), i.e.
n−m−p (µ) are the scattering matrices evolved using the algorithm, modified so that all collinear emissions from hard legs have been removed. Specifically, whereÃ n−m (µ) is computed using (2.2) with the collinear evolution off hard legs removed, i.e. with the replacements The number of collinear emissions not off hard legs is indexed by p and n − p − m is the number of soft emissions (in equation (3.1), m is the number of collinear emissions off hard legs). An example term, contributing to B 3 6 (µ), is presented in Figure 4(a). It may be helpful to contrast B p n (µ) with scattering matrices evolved using the FKS algorithm [25]. We denote the FKS matrices as A soft n (µ) and they can be evaluated using variant A with P υ i υ j = 0; an example is shown in Figure 4(b). Note that B 0 n (µ) = A soft n (µ) for n ≥ 1 since B 0 n (µ) still contains the collinear Sudakov factors 'attached' to soft partonic legs. Also note that B 0 0 (µ) = A soft 0 (µ) and B i 0 (µ) = 0 for all i > 0. In Section 3.2 we will generalise the arguments presented here so that we can factorise collinear physics into jet functions that multiply A soft n (µ). However, in this section we will not make any further use of A soft n (µ). Equation where the traces are over colour, c, and helicity, s. However, we avoid working with collinear factorisation in this form because it does not apply when Coulomb exchanges are present.
Having stated the result, let us now proceed to show how it is derived. The following commutation relations are important: (3.7) -22 - Here denotes equality when the operator acts on a matrix element, which is all we ever encounter.
] 0 are trivial to show as V col a,b is proportional to the identity in colour-helicity space. Diagrammatically, proving and As ever, a red dashed line is used to represent a soft parton and a blue dotted line represents a collinear parton. The black dashed line indicates a cut (cut lines are on shell).
can be shown by factorising kinematic factors from the colour and helicity operators, then carefully tracking the action of the colour operators so that colour conservation can be applied. Proving the commutators also requires noting that both soft real emissions and soft Sudakov factors are identity operators in helicity space, and that helicity states are orthogonal. [V a,b (V col a,b ) −1 , C j ] 0 presents the biggest challenge. The derivation follows closely the discussion in [14]. It is most easily illustrated by expressing the operators diagrammatically. Doing so reduces the problem to showing which trivially follows from (3.10). Using these commutation relations, reconstructing the separate strong orderings of collinear and soft physics in (3.1) is simply a case of careful combinatorics and relabelling of momenta. For instance For the sake of completeness, in the next section we will go ahead and put back R soft , R coll , and the measurement functions. However, we have only proven correctness at LL accuracy. As such, Sections 3.1.1 and 3.1.2 are conjectures. It might be the case that only certain classes of recoil prescription factorise in this way. We will focus on variant A in Section 3.1.1 and turn to variant B in Section 3.1.2, where we show how to re-instate the plus prescription in the collinear splitting functions. We caution that both these versions of the algorithm will not produce super-leading logarithms because Coulomb interactions have been neglected. Therefore they are only suitable for processes with fewer than two coloured particles in either the initial or final state of the hard process, i.e. e + e − , deep-inelastic scattering and Drell-Yan.

The details
Now we summarise the results of the previous section and make a conjecture regarding the inclusion of recoils (recall we left these out of the discussions in the previous section). For concreteness we use variant A. The evolution has two phases. In the first phase soft gluons -24 -are emitted: where l sums only over hard partons. Tr(Col † 0 (µ) • Col 0 (µ) B p n−p (µ)), -25 -whereÃ n n+m obeys the recurrence relatioñ An observable can be calculated using Σ(µ) = n dσ n u n (q 1 , ..., q n ), (3.17) where dσ n = n m=0 dσ(n − m, m).

Recovering the 'plus prescription'
Now let us turn to variant B. Recall that, in this variant of the algorithm collinear evolution proceeds using the full DGLAP splitting functions. Things are precisely as in the last subsection except that we now use the splitting operators without overlines (P il ) and the functions f and F are to be included in the soft terms. We can go a little further and expand out the Sudakov factors in order to recover the familiar DGLAP plus prescription.
To that end, we expand the collinear Sudakov factors (V) that appear in the second phase of the evolution: Once again, the sum over l is only over hard partons. Using this, we can regroup terms in the same way as Section 2.5 to generate the plus prescription: dσ(n, 0) = Tr (dσ(n, 0)) = TrÃ n n+0 , and The sum overl only includes hard partons. P il has been redefined to include the plus prescription and labelled P + il . The plus prescription is defined in Appendix A. Observables are computed using (3.17).

Complete collinear factorisation without Coulomb interactions
Now we are going to go ahead and factorise the collinear physics completely. Once again, the manipulations are essentially the same for either variant of our algorithm. To be concrete, we will use variant B whenever an exact definition must be given. As before, we will begin by stating the final result: where have we defined the following operators: and the index j runs over all partons (hard, collinear and soft). We continue to leave aside the functions R soft , R coll , from emission operators and Sudakov factors. These can readily be re-instated as in Sections 3.1.1 and 3.1.2. A soft n (µ) is as defined in Section 3.1 and evolves the same as A n (µ) except that C i → 0 and V a,b → V a,b (V tcol a,b ) −1 ≡ V soft a,b , i.e. it corresponds to summing over diagrams such as the one in Figure 4(b). Ignoring the effects of recoil (and Coulomb exchanges) and using variant A, A soft n (µ) corresponds to FKS evolution [25]. One of the possible contributions to tCol † 4 is represented in Figure 5. In the figure, we have sacrificed the intuitive picture of a parton cascade in lieu of providing more detail on the evolution. To construct Figure 5, we have employed the Casimir structure of V tcol a,b to split it apart as 25) and the product over j is over all partons. In the figures we will be explicit with the labelling so that it is clear whether U j a,b is associated with a hard parton (labelled by Q), a soft parton or a collinear parton, i.e. j ∈ {Q, 1 soft , 2 soft , ..., (n − m) soft , 1 coll , 2 coll , ..., m coll }. Something more in the style of our previous diagrams is illustrated in Figure 6. We will now prove (3.23) by induction. First, we assume that where A n is computed as usual from (2.2). We can see that this is true for n = 1 by expanding out the tCol operators and A soft : where we have usedC 1 ≡ C 1 ≡ C 1 as it only acts on hard legs. We have also used the , derived in the previous section, and V c,a = V c,b V b,a . Notice in the above expressions the theta functions present inC 1 and V tcol q 1 ⊥ ,Q are always unity on hard legs as the ordering guarantees their argument is satisfied. We will now show that if (3.26) is true for A n , it is also true for A n+1 . We begin by noting that from the Markovian way our algorithm evolves, we can write A n+1 ≡Â n (µ, q 1 ⊥ ) whereÂ n (µ, q 1 ⊥ ) is computed using our algorithm (as described in (2.2)) however with the evolution initiated byĤ(q 1 ⊥ ) = D 1 V q 1 ⊥ ,Q H(Q)V † q 1 ⊥ ,Q D † 1 and with the parton momentum indexed as 2, 3, 4, ... . From this we can use (3.26) to write whereÂ soft n−m (µ, q 1 ⊥ ) are generated by the same algorithm as A soft n−m (µ) however usinĝ H(q 1 ⊥ ) as the initial condition.tCol m (µ, q 1 ⊥ ) are generated using the iterative relation in (3.24) but with an initial conditiontCol 0 (q ⊥ , q 1 ⊥ ) = V tcol q ⊥ ,q 1 ⊥ . Next we split apart H(q 1 ⊥ ) asĤ Using the commutation relations from Section 3.1, we can move the collinear operators in -29 -Ĥ (q 1 ⊥ ) past the soft operators which constructÂ soft n−m (µ, q 1 ⊥ ) to arrive at We can now combine the collinear operators usinĝ , where in the second equality we need to relabel the momenta of collinear partons again so that they are indexed as 1, 2, 3, ... . We have denoted the momentum of the hardest collinear emission generated by the collinear operators, tCol, as q coll 1 and the hardest soft momentum in A soft n−m (µ, Q) as q soft 1 . Combining the two sums and theta functions, we arrive at Thus we have proven that (3.23) holds for n → n + 1. It is important to note the role of the theta functions in the definitions ofC i and V tcol a,b . These ensure that the commutation relations from Section 3.1 can always applied. They do this by squeezing to zero the phase space of any collinear partons generated byC 1 and V tcol q 1 ⊥ ,Q from not-hard legs. To illustrate this point, we will consider the relevant Feynman diagrams: (3.32) -30 -Note that the last term on either side of the equation cannot be manipulated using our commutators and there are no more diagrams we could include which they may cancel against. Nevertheless terms of this form are generated by our algorithm. They represent a collinear parton, emitted from a soft parton, restricted so that its transverse momentum is smaller than the transverse momentum of the soft parton. UsingC i , equation (3.32) reduces to (3.33) Though it is not possible to use the identities in (3.7) to factorise collinear physics past a Coulomb exchange (iπ term), it is possible to perform a partial factorisation. Our approach is to expand each Sudakov operator as a series in the number of Coulomb exchanges -31 -it resums. Consequently, this enables A n to be expanded as a series in the number of Coulomb exchanges. We can then factorise soft physics from collinear physics either side of a Coulomb exchange, using our work in the previous section. The partial factorisation is illustrated in Figure 7. We begin by expanding the Sudakov operator:

Partial collinear factorisation with Coulomb interactions
whereV a,b is equal to V a,b withδ ij = 0. Consider using this expanded Sudakov in a parton cascade. The theta functions describing the integral limits on each iπ term can be used to constrain the limits on the transverse momenta of subsequent emissions (after the iπ). For instance Therefore, we can treat each Coulomb scale as hard relative to the emissions that follow it and soft relative to the emission before it. Thus we can perform a factorised evolution on a hard process up to the scale of the first iπ term (k 1⊥ ). We can take the output from this evolution, 36) and use it as a new hard process H(k 1 ⊥ ) from which a second factorised evolution can be initiated. This process can be iterated for each iπ term in the expansion, as illustrated in Figure 7. To complete the computation of Σ, each k i ⊥ must be integrated over the range [µ, Q]. Interestingly, note that any term in the evolution terminating on a iπ term to the left of the hard process will cancel against an equivalent term terminating with an iπ term to the right. Hence collinear emissions can always be factorized below the scale of the last Coulomb exchange. This is consistent with the collinear factorisation shown by Collins, Soper and Sterman [30,31].

Observations on factorisation
Before we leave our discussion on factorisation a few comments are in order. Firstly, we have not been able to achieve factorisation of collinear emissions past Coulomb exchanges. This is to be expected and there is already extensive literature exploring this subject [14,16,17,31,[45][46][47][48][49]. That said, it should be possible to factorise more completely than we have done, by re-expressing the evolution so that all Coulomb terms are only attached to the initial state partons [17], i.e. so we would have complete factorisation on all final state legs. Secondly, in order to factorise the collinear physics on all legs we had to keep track of intermediate soft scales, from which to initialise the collinear evolution. The number of scales required is equal to the number of soft emissions that occurred prior to factorisation. This means the fully factorised algorithm is no-longer Markovian. We anticipate that our attempts to factorise the collinear physics should bring us in to contact with exact resummations and soft-collinear effective theory (SCET).
It also should be noted that by factorising collinear emissions from the soft evolution, the soft evolution can be explicitly seen to be independent of spin. This is less evident in the interleaved variants of the algorithm. Soft gluons, and subsequent collinear partons, trapped between Coulomb exchanges might conceivably contribute non-trivial spin correlations. This is because, despite equal probabilities for the probability of emission of positive and negative helicity gluons, a collinear emission originating from a soft gluon may depend on its helicity (specifically g → qq splitting). This has also been explored in the literature, where it has been noted that soft gluons in the presence of Coulomb/Glauber exchanges can generate spin asymmetries [48]. Further discussions on the spin evolution of the algorithm after factorisation can be found in Appendix B.
It is also interesting to consider the consequences of factorisation in the case of variant B with a universal recoil. A universal recoil allows B to be partitioned in terms of colourdiagonal evolution generated by C i and colour off-diagonal evolution generated by S i . Hence, provided the recoil prescription does not change the commutators in (3.7), the proofs of collinear factorisation we have presented become proofs of the complete factorisation of colour-diagonal physics from colour off-diagonal. This is for observables insensitive to the presence of Coulomb exchanges. Since we know that Coulomb exchanges can be factorised onto the initial state [17], this means that there is a complete factorisation of colourdiagonal from colour off-diagonal physics in lepton-lepton, deep-inelastic and Drell-Yan scattering.
Finally, we should remark that it is possible to write down infra-red finite versions of each of the factorised versions of our algorithm, using the procedure in Section 2.5.
-33 -In this section we will first demonstrate how DGLAP evolution emerges. After that, we illustrate the use of the algorithm by calculating thrust, the hemisphere jet mass and gaps-between-jets in e + e − .

DGLAP evolution
We will now show how our algorithm can be used to generate DGLAP evolution, which resums the collinear physics into the running of parton distribution functions. We focus on unpolarised incoming hadrons that collide and produce some high-p T system of interest. We will neglect threshold effects as sub-leading, which is shown carefully in [50][51][52]. The methods employed in this section can readily be extended to other processes, including those dependent on fragmentation functions. DGLAP evolution [34,36] states that where f i (x, µ) is the parton distribution function for partons of type i. P ij (z) are the regularised splitting functions defined at the end of Appendix A. Iterative solutions can be found by expanding the parton distributions: where f (n) i (x, µ) = 0 for all n ≥ 1. This gives which has a separable solution of the form f We can write this in terms of the unregularised splitting functions (e.g. see [36]) where f (n) j (x) = 0 for x > 1 and we have removed factors of n f from P qg . For hadron-hadron collisions, we label the two incoming partons as a and b and their momentum fractions in the hard process as x a and x b . We can take the factorised expression x a corresponding to variant B of our algorithm (3.1) and attach parton distribution functions: (4.6) x am and x bm are the momentum fractions of partons a and b respectively after m collinear emissions generated by Col † m (µ) • Col m (µ); they can be related back to x a and x b by momentum conservation along the collinear cascade. The operator acts to attach parton distributions of the correct flavour/species to partons a and b. There is a technicality relating to parton flavour. That is because DGLAP evolution cares about quark flavour whilst we have defined the splitting operators to sum over quark flavours (in the case g → qq). We could have avoided this technicality by defining the splitting operators per flavour (i.e. set n f = 1 throughout Appendix A). Then we would have to sum over quark flavours throughout the rest of the paper. Instead, we choose to handle quark flavour by keeping track of flavour along the evolution chain, and whenever a g → qq splitting occurs we label the subsequent parton flavour generically, i.e. for two-flavours the relevant set of parton flavours would be {u,ū, d,d, q, g}. Note that since we evolve away from the hard scattering, a g → qq branching from an incoming g actually corresponds to a q → qg (or q →qg) splitting in the usual DGLAP sense. This can be seen in (A.1), where the terms involving δ initial j and s j = ±1 involve the P gq splitting function (the P qg splitting function -35 -appears in the corresponding δ initial terms). With this in mind we have that where α and β label parton type. For n f = 2, we would have α, β ∈ {u,ū, d,d, q, g}, and 2n f f q is the singlet distribution function, i.e. f q = (u +ū + d +d)/(2n f ) for n f = 2. For completeness, we here also attach labels A and B to indicate the type of hadron (we will drop that label elsewhere in this section).
After expanding Col † m (µ) • Col m (µ) in powers of α s , then spin averaging at every vertex, substituting for (4.4) and evaluating the transverse momentum integrals, we find (4.8) Figure 8 illustrates how terms in our algorithm should be grouped in order to generate the iterative relation in (4.5) and so arrive at (4.8). Hence we see that variant B iteratively generates DGLAP evolution up to the hard scale. The derivation of fragmentation function evolution of final-state partons proceeds similarly. For processes where Coulomb exchanges are relevant, DGLAP evolution is generated up to the scale of the last Coulomb exchange. Also note that, in the infra-red finite reformulations of A and B, DGLAP can be found in the B n for n ≥ 1.

Some familiar resummations
In this subsection we show how to resum a number of observables in e + e − → hadrons. The idea is to use well-known results to illustrate the use of the algorithm. The simplicity of the hard process means we can use H(Q) = N −1 c σ H 1 for the hard-scattering matrix. We perform all calculations using the LL recoil from Section 2.4.1.

Thrust
The resummed thrust distribution was initially computed at LL accuracy in [53], then at NLL in [54]. The current state-of-the-art computation is at N 3 LL [55]. Thrust is defined as where the thrust axis n points along the initial hard parton axis at leading-order in the soft and collinear limits; see Section 3.1 in [54]. We will only need to define thrust for 3 partons; of which 2 are hard (p 1 and p 2 ) and one is soft (k). We can work in the dipole zero-momentum frame using p q = E q (1, 0, 0, 1), pq = Eq(1, 0, 0, −1), k = (k ⊥ cosh y, k ⊥ , k ⊥ sinh y), and we can fix 2E q = 2Eq = Q. Thrust is evaluated as (4.10) -36 -When calculated in the hard-collinear limit, we can let T = 1 + O(k 2 ⊥ ) as all partons lie on the thrust axis up to sub-leading contributions. The thrust distribution R T is defined as As thrust is global, the calculation is most readily performed using the manifestly infrared finite version of variant A (see Section 2.5). Using this, all terms with one or more emissions cancel exactly. The measurement function u(k, {∅}) is This is unity for a hard-collinear emission, since |y| → ∞ at LL accuracy. This kills the hard-collinear terms, since they contain a factor (1 − u(k, {∅})) = 0, which is as expected since they contribute no double logarithms. Thus we can immediately write (4.13) After integrating (4.14) Here we used the fact that θ ij (k) restricts the integration so that k 0 < Q. The colour trace can be evaluated to give (4.16)

Hemisphere jet mass
The hemisphere jet mass is subject to non-global logarithms, which greatly increase the challenge of resummation. It was first resummed at LL and LC in [12]. The current stateof-the-art is split between fixed-order computation (α 5 s with leading colour [56] and α 4 s with full colour [57]) and resummation using numerical techniques to introduce full colour dependence with sub-leading logarithms [58]. The measurement function corresponding to the hemisphere jet mass in e + e − → hadrons is where S + 2 and S − 2 are the hemispheres centred on the two primary jets. m ± is the total invariant mass in the S ± 2 hemisphere. The measurement function can be simplified by considering m ± At the order we will perform the calculation we only need to consider one emission, hence where E ± is the energy of the quark/anti-quark defining the S ± 2 hemisphere and Q = 2E ± . Note the similarity between this and the measurement function for thrust, which is expected since it is well known that, at lowest-order, thrust can be expressed as the sum over the two hemisphere jet masses defined by the thrust axis.
Again, we can use the manifestly infra-red finite version of A to find Σ(ρ): From the calculation in the previous section, we can immediately write This gives the global contribution. The non-global contributions are found by evaluating the remaining terms (corresponding to summing over real emissions in (2.61)). This calculation can be found in [25], where the non-global terms are evaluated using the FKS algorithm, which is entirely sufficient in this case. Hence we find

Gaps-between-jets
The LL measurement function in the case of gaps-between-jets is where the 'in' region corresponds to two cones centred on the two leading jets and the 'out' region is the region between these cones. The observable vetoes emissions in the out region that have transverse momentum greater than Q 0 . At order α 5 s this observable is sensitive to super-leading logarithms [13,14]. These will be correctly calculated using variants A, B and their manifestly infra-red finite versions, but not their factorised form unless Coulomb exchanges are interleaved as in Section 3.3). Using the manifestly infra-red finite version of variant A we correctly find that where Y is the rapidity range of the out region and the 1 + O(α 2 s ) factor is the stack of non-global logarithms, which can be computed by considering real gluon emission into the out region, as encoded in (2.61). These were calculated up to O(α 5 s ) in [13,14,59]. We note a kinematic maximum on the rapidity of an emitted gluon, i.e. 2|y| < Y max = ln Q 2Q 0 + Q 2Q 0 − 1 . This means that as Y → Y max all soft radiation goes into the in region. At leading-order in the soft approximation Y max = ln(Q/Q 0 ), i.e. for Y ≥ Y max the observable becomes doubly logarithmic.

Conclusions
Our primary goal in writing this paper is to provide the theoretical basis for the future development of a computer code that is able systematically to resum enhanced logarithms due to soft and/or collinear partons including quantum mechanical interference effects. The algorithm we present, and its variants, are (mostly) Markovian and their recursive nature makes them well suited for the task. First steps towards this goal are under development, using the CVolver code to perform the colour evolution [60][61][62].
The algorithms in this paper correctly account for the leading soft and/or collinear logarithms, though we have been careful to try and present them in such a way as to make -39 -the extension beyond LL possible. For example, we have taken account of the momentum re-mappings that are necessary in order to account for energy-momentum conservation and we have included g → qq transitions which are strictly single logarithmic.

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 lodowska-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 and the Particle Physics Group of the University of Vienna for financial support and kind hospitality. JH thanks the UK Science and Technology Facilities Council for the award of a studentship. SP acknowledges partial support by the COST actions CA16201 PARTICLEFACE and CA16108 VBSCAN. We would like to thank Matthew De Angelis for extensive discussions. We would also like to thank Mrinal Dasgupta, Jack Helliwell and Mike Seymour for helpful discussions and comments on the manuscript. Figures have been prepared using JaxoDraw [63].

A Splitting functions
The splitting operator P ij (see Figure 2), which is explicitly used in variant B of our algorithm, is built from the spin dependent DGLAP splitting functions [38,64]. It is an operator in colour and helicity spaces and is defined using the spinor-helicity formalism [39]. We use the convention v(p, λ) = Cū T (p, λ) where C = iγ 2 γ 0 , which defines our crossing symmetry to have no global minus sign. Using rotational symmetry and parity invariance one generally can write M( where M is a matrix element and {λ i } the set of helicity states on which M depends. Together these define the correct treatment for antiparticles, which should evolve as if they are particles with the opposite helicity. Thus P ij is Here s j is the spin/helicity of parton j and z i is the momentum fraction between parton i and its parent parton, j (as in (2.9)). T g j are the basis-independent colour-charge operators for the emission of a gluon [25,65]. T q j is the colour charge operator for the emission of a qq pair from a gluon. In the colour flow basis it is where τ j exchanges the anti-colour lines associated with the colour line of parton j. For example, let parton j have colour line c 2 and anti-colourc 5 , τ j would exchange anti-colour linesc 2 andc 5 . A full definition of τ j , and other colour flow operators, can be found in [25], where τ j is written s α,β . Note T q j · T q j = T R 1. 10 We have defined S s as the operator that 10 Strictly speaking this is only valid when acting on a physical matrix element.
-42 -adds a parton with helicity s. Just as T g i · T g i = C i 1, it can be shown that S s · S s = 1. We have also defined a 'swap' operator, W ij , which swaps the colour and helicity of particles i and j. Finally, we defined P 1 i as the operator that flips the helicity of parton i and P 2 i that halves the helicity of i. There is some freedom in how we introduce operators to keep track of the evolving helicity state (for example one could have instead made use of (S s ) † , which deletes a parton of helicity s). Examples to illustrate the use of these helicity operators can be found in Appendix B. The (unregularised) collinear splitting functions are It should be understood that P ij always acts as is a generalised spin operator and S (υ j →υ) s i ,s i is a c-number coefficient corresponding to a υ j → υ splitting. For example, when j is a quark The Sudakov factors in variant B can be written in a variety of ways using Note also that there is a subtle factor of two difference between initial state and final state splittings in P ij . The factor arises as partons in the initial state must be convoluted with PDFs which changes the pole structure of splittings and increases the number of diagrams that must be summed over relative to splittings in the final state. In the final state (without fragmentation function dependence), soft poles from real emissions can be found at both z = 1 and z = 0. These poles cancel the poles from loop diagrams. For real emissions in the initial state, the z = 0 poles are absent due to kinematics whilst the z = 1 poles cancel the poles from loops. The factor of 2 ensures the correct pattern of cancellations. Finally, we also note the factors of n f , the number of quark flavours, in P ij . They are present since we sum democratically over flavours whenever there is a g → qq branching. Note that since we always evolve away from the hard process this means that we sum over quark flavours in the case of an initial-state q → gq branching. Care must be taken however, since if the branching cascade terminates with an initial-state quark (or anti-quark) then it is necessary to divide by a factor of n f before convoluting with the corresponding parton distribution function. The same holds in the case where fragmentation functions are needed. In Section 4, we introduced the notation to handle this. Of course, one could set n f = 1 in the above splitting operators, after which it would be necessary to sum over flavours as appropriate.
For variant A, we need the hard-collinear emission operator P ij . This operator is defined at cross-section level through the relation where O is a generalised operator. Note that ...P ij OP † ij ... is not necessarily Casimir in colour. However, as we observed in section 3, ignoring Coulomb contributions, the collinear physics can be factorised and becomes colour-diagonal after taking the trace. Therefore, for processes where Coulomb terms do not contribute (e.g. e + e − and DIS) we could use the emergent colour-diagonal structure to greatly simplify the P ij and P ij operators. For example, we could redefine P ij with a simpler amplitude-level statement. To this end, we can introduce the hard-collinear splitting functions; The newly simplified P ij is equal to the operator found by substituting P → P inside P ij , This new P ij operator is constructed so that when used in the LHS of (A.7) the expression becomes exact after a trace is taken. Additionally, it correctly computes spin correlations after collinear factorisation. Simplifying the collinear emission operators would be very pertinent to an efficient computational implementation of our algorithm. P • υ i υ j and P • υ i υ j are splitting functions used exclusively in our Sudakov factors and they are defined with all colour factors removed: (A.9) In Section 3.1.2 and Section 4 we make use of the plus prescription (see (3.22)). Applying the plus prescription means (1)]. (A.10) The plus prescription is, in our case, is defined by where the structure of the subtraction terms is determined by the corresponding structure of the virtual corrections and this simply means that P V is determined using (A.1) but with parton j always treated as if it is final state. The two splitting functions affected by the plus prescription are are cross-section level decay matrices. D and ρ are calculated from products of amplitude level matrices, D and S respectively. For instance, after n emissions from parton 1: where usual matrix multiplication is implied. The algorithm of Collins and Knowles is able to determine the spin density and decay matrices such that computational time only grows linearly with the number of partons [66][67][68]. Now let us turn to the calculation of splitting functions in our notation. We write M n s 1 ...sn = s 1 ...s n | n , which ignores colour since it is not our focus here, i.e. more correctly we should write M n c 1 ...cn,s 1 ...sn = ( c 1 ...c n | ⊗ s 1 ...s n |) |n . We wish to define an operator P k→ij that adds a new (collinear) particle (j) to |n that is emitted off leg k, i.e. |n + 1 col = k∈{n} P k→ij |n .
As before, we will focus on the q → qg collinear splitting. Note that M n+1 s 1 ...s i ,λ j = s 1 ...s i , λ j | P k→ij |n . Inserting the identity gives This is the definition for S s presented in Appendix A 11 . We will now construct a decay matrix D (j) s j s j , for a final-state hard parton j using the spin operators we have just introduced. Let us first consider the situation where there are 11 Repeating this procedure for the other splitting operators leads us to introduce the operators W ij , P 1 i and P 2 i used in Appendix A.
-47 -no soft interactions and only include emissions from the initial primary leg, j: dσ ∝ n {i} n; j| V † 1,Q P † i 1 j V † 2,1 ...V † n,n−1 P † inj V † 0,n V 0,n P inj V 2,1 ...V 2,1 P i 1 j V 1,Q |n; j , (B.7) where the partons in the set {i} are transverse momentum ordered. V a,b is a Sudakov factor: We can evaluate (B.7) by inserting identity operators and extracting Sudakov factors, which are proportional to identity operators, into a single numerical factor. Hence dσ ∝ n {i} s j s j # s j s j i 1 ...ini 1 ...i n n; j| s j s j | P † i 1 ...P † in P i n ...P i 1 s j s j |n; j , (B.9) where each P † i is a 'pure' colour-helicity operator with no scalar pre-factor. For instance P i = T g j ⊗ S 1 i in the case of a q + 1 2 → q + 1 2 g +1 splitting or P i = T q j ⊗ P 1 j P 2 j S −1 i in the case of a g −1 → q + 1 2 q − 1 2 splitting. # s j s j i 1 ...ini 1 ...i n is a c-number coefficient built from helicity dependent splitting functions and expanded Sudakov factors. We can now make the link with the previous approach, i.e. The expectation value is calculable and equals a product of n Casimir co-coefficients, e.g. if parton j is a gluon and each operator P i corresponds to a g → gg splitting then the expectation value equals C n A (up to a normalisation for the colour evolution). Following this procedure, spin-density and decay matrices can be derived using the algorithm presented in this paper. Let's see this explicitly. Knowles' algorithm calculates spin-density and decay matrices using other intermediate matrices ρ , ρ , D and D [67]. We will calculate ρ using the factorised form of variant B (which we refer to as B-f), with LL recoil. ρ ss describes the distribution of spin states for a give parton after a single collinear emission. It is normalised by the trace of itself so that it maintains a probabilistic interpretation. Knowles begins by defining where ρ is a spin density matrix for a parton in the hard process that is to be inherited by a forwardly evolving shower. In the language of this paper ρ s 1 s 1 ∝ s 1 | H(Q) |s 1 12 . V s 1 s 2 s is the spin-dependent collinear splitting function for the transition s 1 → ss 2 with parton type indices suppressed. Importantly, parton type indices are not summed over in V s 1 s 2 s . When using Knowles' algorithm, it is assumed that the structure of a cascade has already been fully decided; all except the spin that is. Consider a term from B-f corresponding to one collinear emission from a final-state hard parton. Labelling this term P , we have , = s 2 s 2 s, s 2 P 2 1 H(Q)P † 2 1 s , s 2 δ s 2 s 2 Tr(P 2 1 H(Q)P † 2 1 ) .

(B.12)
We have used the LL recoil with variant B and so integrals over the recoil functions were trivial. In the second line, we have labelled the collinear parton as parton 2 and the hard parton as parton 1. We can insert identity operators and evaluate the trace explicitly to find P ss = s 1 ,s 1 ,s 2 ,s 2 s, s 2 | P 2 1 |s 1 s 1 | P † 2 1 |s , s 2 ρ s 1 s 1 δ s 2 s 2 s,s 1 ,s 1 ,s 2 ,s 2 s, s 2 | P 2 1 |s 1 s 1 | P † 2 1 |s, s 2 ρ s 1 s 1 δ s 2 s 2 . (B.13) Now note that ss 2 | P 2 1 |s 1 = V s 1 s 2 s , with the possibilities of parton 2 being a gluon or quark summed over. Hence When comparing with Collins and Knowles it was necessary for us to pick a species for parton 2 as their algorithm is defined for pre-determined decay chains. This is why ρ ss is typically used without a label specifying the species of the partons involved, as their species is always provided by context. We will finish off by calculating ρ ss for a q → qg splitting. ρ s 1 s 1 is hermitian and so can be expressed as ρ s 1 s 1 = 1 + ρ i σ i where σ i are the Pauli matrices. Using (B.13) and normalising correctly gives Similarly, matrices for the other collinear splittings can be found. The most algebraically complex is the g → gg splitting (as usual). In that case P (g→gg) ++ = 1 + where φ is the azimuthal angle to the plane of the splitting. The exact angular dependence depends on the definition of the Weyl spinor products. We have chosen the definition so as to match with the matrices defined in [67], where a factor e i(s 1 −s 2 −s)φ has been pulled out from the definitions of V s 1 s 2 s .