Building a consistent parton shower

Modern parton showers are built using one of two models: dipole showers or angular ordered showers. Both have distinct strengths and weaknesses. Dipole showers correctly account for wide-angle, soft gluon emissions and track the leading flows in QCD colour charge but they are known to mishandle partonic recoil. Angular ordered showers keep better track of partonic recoil and correctly include large amounts of wide-angle, soft physics but azimuthal averaging means they are known to mishandle some correlations. In this paper, we derive both approaches from the same starting point; linking our understanding of the two showers. This insight allows us to construct a new dipole shower that has all the strengths of a standard dipole shower together with those of an angular ordered shower. We show that this new approach corrects the next-to-leading-log errors previously observed in parton showers and improves their sub-leading-colour accuracy.


Introduction
Parton showers simulate the particle content of scattering events at collider experiments and provide the backbone to modern experimental analyses [1][2][3][4][5][6][7]. Yet questions over their accuracy and on how best to improve them remain. In this paper we present a unified analysis of the two main approaches to formulating parton showers: dipole showers [2][3][4]8] and angular ordered showers [5,6,9]. As a result, we are able to construct a new dipole shower that does not suffer from the next-to-leading logarithm (NLL) and next-to-leading colour (NLC) problems suffered by existing parton showers [10].
In our previous papers [11,12] we introduced an algorithm for amplitude-level parton branching (the PB algorithm). The PB algorithm was designed to capture both the soft and collinear logarithms associated with the leading infra-red singularities of scattering amplitudes without making any approximations on the spin and colour. In [12] we showed how the PB algorithm can be used to derive the resummation of observables at leading-logarithmic accuracy (it has the capacity to be extended to include next-to-leadinglogarithms) and we showed that it gives rise to the collinear factorisation of parton density and fragmentation functions. In [11] we showed that the colour evolution is equivalent to that of other approaches [13][14][15][16]. The PB algorithm is the starting point for the analysis presented here.
In the next section, we present a brief overview of the algorithm before going on to use it to derive both dipole and angular ordered showers. In these derivations we keep close track of the approximations made, with the goal of gaining a solid understanding of the sources for errors in these showers. We focus on deriving showers in e + e − , though much of the machinery necessary to derive showers for hadron-hadron processes is also present in this paper. The full discussion of our derivations is technical and largely handled in Appendix A.
More specifically, in Section 2.2, we derive an angular ordered shower starting from the PB algorithm. In doing so we are able to constrain the recoil functions in the original PB algorithm, since angular ordered showers provide clear constraints on how momentum longitudinal to a jet must be conserved in order to get NLL physics correct. In Section 2.3 we then derive a dipole shower from the PB algorithm, taking particular care over the constraints observed from our angular ordered derivation. The result is a dipole shower that reduces the NLC errors noted in [10] (complete removal of NLC errors requires amplitudelevel evolution). Having pinned down longitudinal recoil, in Section 3 we present a scheme (inspired by [17]) for the transverse recoil. This completes the specification of our shower. We then go on to recreate the fixed order analysis of [10] and show that our shower corrects the NLL errors from incorrect transverse recoil previously observed in dipole showers. In Appendix D we go further and show that our new shower is sufficient for the correct NLL resummations of thrust and the generating functions for jet multiplicity.

Amplitude evolution overview
The PB algorithm defines a sequence of transitions in a Markov chain of amplitude density matrices: A 0 (q 0 ⊥ ; {p} 0 ) → A 1 (q 1 ⊥ ; {p} 1 ) → ... → A n (q n ⊥ ; {p} n ). The sequence is illustrated in Figure 1. We use n to index the number of partons dressing the hard process. Each amplitude is defined at a given scale (parametrised by an ordering variable), this is its first argument. The second argument, after a semi-colon, specifies its full dependence on the relevant parton momenta (which we often choose to omit). The Markov chain uses the initial condition A 0 (Q; {p} 0 ) = H(Q; P 1 , ..., P n H ), where H ≡ |M M| is the hard process density matrix for a process of hard scale Q and with n H hard partons. The hard partons' momenta form the set {P 1 , ..., P n H } ≡ {p} 0 . The Markov chain terminates on the amplitudes A n (µ; {p} n ); µ is an infra-red cut-off and {p} n = {P 1 , ..., P n H , q 1 , ..., q n } where q 1 , ..., q n are the momenta of the n partons that dress the hard process. Steps in the Markov chain are constructed from the action of two operators, D n and Γ n . The D n operators are emission operators; they act as maps from a state A n−1 (q ⊥ ; {p} n−1 ) to a state A n (q ⊥ ; {p} n ), and they describe the emission of the nth parton. Operators Γ n provide a map from a state A n (q ⊥ ; {p} n ) onto some otherÃ n (q ⊥ ; {p} n ). Physically, they dress the density operator with (iterated) virtual corrections. The path-ordered exponent of Γ n is an amplitude level Sudakov factor/operator which we call V a,b : V a,b evolves a state A n (b; {p} n ) to a state at a lower scaleÃ n (a; {p} n ); for a complete discussion of V a,b see [12]. In [12] we presented the PB algorithm in the following form: A n (q ⊥ ; {p} n ) = dR n V q ⊥ ,q n ⊥ D n A n−1 (q n ⊥ ; {p} n−1 )D † n V † q ⊥ ,q n ⊥ Θ(q ⊥ ≤ q n ⊥ ). (2. 2) The algorithm maps the set of partonic momenta prior to the nth emission ({p n−1 }) onto a new set ({p n }), by adding a parton (q n ). In order to conserve energy-momentum, the set of momenta prior to the emission are adjusted after each emission, i.e. {p n−1 } → {p n−1 } and {p n } = {p n−1 ∪ q n }. We achieve this by integrating over delta functions relating the two sets of momenta. This is all hidden inside dR n , which we describe in Appendix A. 1 and give examples of in Section 3. We also provide definitions of each operator involved in the evolution in Appendix A.1. Figure 1. A general term in the Markov chain of amplitude density matrices, A n , constructed by the PB algorithm. H ≡ |M M| is the initial hard process; in this case it has two hard coloured legs, a and b. D n dresses an amplitude with the nth emission that is either soft or collinear. Collinear emissions are emitted symmetrically from the amplitude and conjugate amplitude, such as gluon 1. Soft emissions appear as interference terms, such as gluon 2. Γ n dresses the amplitude after n soft or collinear emissions with a loop.
In this paper, it better suits our purposes to work with the PB algorithm expressed as an evolution equation, i.e. working differentially in the ordering variable, q ⊥ . Broadly speaking, q n ⊥ is the transverse momentum of the nth parton and it is a function of the n-parton phase-space. The precise definition of q n ⊥ is context dependent and is given in Appendix A.1. The evolution equation is + dR n D n (q n ⊥ ) A n−1 (q n ⊥ ; {p} n−1 ) D † n (q n ⊥ ) q ⊥ δ(q ⊥ − q n ⊥ ). (2.3) It is from this equation that we will derive generalised dipole and angular ordered showers. The phase-space measure for the nth parton emitted in the cascade is variously written as We typically parametrise the evolution so that real emissions use the phase-space measure dΠ n and loops d ln q n ⊥ dS (qn) 2 . From each A n we can compute the differential n H +n parton cross section: dσ n (µ) = n i=1 dΠ i Tr A n (µ), (2.5) where µ is either an infra-red regulator that should be taken to zero or the shower cut-off scale. We will focus on e + e − hard matrix elements, in which case observables are computed using Σ(µ; {p} 0 , {v}) = n dσ n (µ) u({p} n , {v}), (2.6) where u({p} n , {v}) is a measurement function for an observable defined by the set of parameters {v}. 1 The formula for processes involving incoming hadrons is given in Appendix A.1.1.

Angular ordered shower
In this section we give an overview of the derivation of an angular ordered shower, starting from Eq. (2.3). The unabridged derivation is given in Appendix A.2. Angular ordering is derived after averaging over the azimuth of each emitted parton, as measured relative to their parent parton (and neglecting all subsequent azimuthal correlations). After performing this averaging in Eq. (2.3), the colour structures can be greatly simplified (a manifestation of QCD coherence). We exploit this to re-write the evolution in terms of squared matrix elements, |M n | 2 . What follows is a little more detail of the key steps.
1. The D n operators in Eq. (2.3) describe the emission of soft gluons from dipoles (via eikonal currents) and the emission of hard-collinear partons. The probability for the emission of a soft gluon is partitioned as n in · n jn n in · n n jn · n = P injn + P jnin , where 2P injn = n in · n jn − n in · n n in · n n jn · n + 1 n in · n , n in = p in /E in and n = q n /E qn , and E is an energy in the event zero-momentum frame. Note that P injn only has a pole when the emission is parallel to i n . When integrated, this term gives rise to a theta function that enforces angular ordering.
2. We average over the emitted parton's azimuth, ... 1,...,n , such that (for some quantity f ) The relevant angles are defined in Figure 2. We use this operation on both sides of Eq. (2.3) and spin-average, see Appendices A.2 and B for details. It is at this point we see that P injn n ∝ Θ(θ jn,in − θ n,in ).
3. We perform a change of variables, q n ⊥ → ζ n,jn = 1 − cos θ n,jn , so as to make the angular ordering explicit. We merge the soft and hard-collinear emission kernels; expressing them in terms of collinear splitting functions. We also must sort out recoil so that the longitudinal component of the total momentum in a 1 → 2 splitting is conserved. Finally, using kinematic variables defined in the event zero-momentum frame 2 allows us to saturate the Θ(θ jn,in − θ n,in ) angular ordering constraint for emissions originating from the primary hard partons (which are anti-parallel to each other). For all other emissions, it is necessary to approximate Θ(θ jn,in − θ n,in ) ≈ 1. This approximation (which corresponds to strong ordering in angles) is equivalent to assuming the angle of the current emission is smaller than the opening angle 2 i.e. for e + e − → qq, zn =pi n · n/pi n · n and n is chosen so that n||Pq for all emissions in the quark jet and vice versa for the anti-quark jet.
θ in,jn θ n,jn i n j n φ n Figure 2. The angles used to derive angular ordering by azimuthal averaging. φ n is the azimuth that is averaged over. In some equations two azimuths are present, in these situations we give φ n a second index, e.g. φ n,jn . Angular ordering corresponds to θ in,jn > θ n,jn . of every other dipole, not just the opening angle of its parent dipole. This is the familiar angular ordering used in both resummations [19,20] and parton showers when showering from an e + e − → qq hard process [5]. Strong ordering in angles simplifies the colour structures, so that all colour-charge operators can be reduced to Casimir, i.e. C F for a quark and C A for a gluon. The simplified colour reduces the evolution equation to an evolution of matrix elements, |M n | 2 .
The final result is The angular ordering variable ζ n,jn = 1−cos θ n,jn . P υυ jn (z n ) are the usual collinear splitting functions, e.g. P qq (z n ) = C F 1+z 2 n 1−zn . Here we have used υ jn to label the species of parton j n and υ to label the species j n transitions to; if υ jn = q then υ = q and if υ jn = g then υ = q, g. z n is the momentum fraction of parton n, i.e. if we have a collinear splitting that induces j n−1 → j n n then p jn ≈ z n p j n−1 and q n ≈ (1 − z n )p j n−1 . Θ on shell is a product of theta functions that ensures each parton is integrated over the phase space corresponding to a real particle (see Section A.2.2). In the first term, Θ on shell is a function of ζ, z and the n-parton phase space. In the second term Θ on shell is a function of ζ n,jn , z n and the (n − 1)-parton phase space. |M n (ζ; {P 1 , ..., P n H , (z 1 , ζ 1,j 1 ), ..., (z n , ζ n,jn )})| 2 1,...,n is the azimuthally averaged, squared matrix element for a hard process dressed with n stronglyordered partons with a unique branching topology; each emitted parton is specified by a pair (z m , ζ m,jm ) and parton j m is the corresponding parent. The delta function enforces longitudinal momentum conservation; |M n | 2 depends on the momentum after the emission, p jn , and |M n−1 | 2 depends on the momentum before the emission, p jn .
Observables in e + e − are computed after summing over emission topologies: where µ should be taken to zero (or the shower cutoff) and for hadron-hadron collisions see Appendix A.1.1. 3 There are several noteworthy points involved in this derivation: • In order to reduce the colour structures to being diagonal, we made the approximation Θ(θ jn,in − θ n,in ) ≈ 1 for emissions from partons other than the two primary hard particles. The approximation is generally only good to LL accuracy (though angular ordered showers are able to go beyond this when combined with the CMW running of the coupling [20], e.g. to compute thrust at NLL [19]). Moreover, modern angular ordered showers retain information on the hard-process, leading N c colour flows by working in the dipole frames of initially colour-connected partons. This improves the approximation for hard processes with greater than two hard jets, since it is then only required to assume Θ(θ jn,in − θ n,in ) ≈ 1 for emissions from partons other than the primary hard partons. During the subsequent evolution, traditional angular ordered showers lose the information on QCD colour flows 4 , while dipole showers retain it to all orders at leading N c . We will exploit this in our dipole shower construction. Appendix A.2 and A.3 give more details on this point.
• The shower does not yet fully conserve energy and momentum. Rather it only conserves energy-momentum longitudinal to a jet. Accounting fully for energymomentum conservation is formally sub-leading in many observables. However, it is phenomenologically important and necessary for shower unitarity. Furthermore, if total energy-momentum conservation is handled incorrectly it can spoil the NLL accuracy of a shower for some observables [10]. We will return to this in Section 3.
• We averaged the azimuthal dependence of the matrix elements. However, this ignores possible azimuthal dependence of the observable. Really one should compute |M n | 2 u({p} n , {v}) 1,...,n . It is therefore important to ask whether is a good approximation. In other words, are the azimuthal dependencies of the matrix element and the observable correlated? This is clearly an observable dependent statement. Despite this we can make some progress; we can remove the approximation and find |M n | 2 u({p} n ) 1,...,n = |M n | 2 1,...,n u({p} n ) 1,...,n + n m=1 σ m ( |M n | 2 1,...,n ) σ m ( u({p} n ) 1,...,n ) Cor m ( |M n | 2 1,...,n , u({p} n ) 1,...,n ) + higher order correlations, (2.9) where σ n (x) = x 2 n − x 2 n and Cor n (x, y) = . The first order correlation term (the second line of Eq. (2.9)) acts as a switch. If it is suppressed relative to the uncorrelated term then all higher correlations will be too. If it is not suppressed then higher order correlations may not be. In Appendix 2.2 we show that the higher order correlations are subdominant in the computation of NLL thrust. This is because the observable is two-jet dominated 5 and exponentiates, and so at NLL accuracy σ m ( u({p} n ) 1,...,n ) ≈ 0. However, we also find that the correlation term can 4 Some azimuthal correlations due to colour correlations can be re-instantiated in coherent branching algorithms [21,22]. 5 Observables, such as thrust, for which the leading logarithms quantify small deviations from the two-jet limit or, more generally, the n-jet limit in the case of n-jettiness [23] provide a formally leading contribution to non-global logarithms. In Appendix A.2 we observe that the correlation terms contribute leading logarithms to observables like gaps-between-jets, for which α n s L n logs are leading. The miscalculation of non-global logarithms by angular ordered showers has previously been subject to numerical study in [24,25], where it was observed that leading non-global logarithms are incorrectly computed by angular ordered showers. However, [24,25] also observed the error to be a phenomenologically small effect.

Dipole shower
In the PB algorithm, the mechanism for energy-momentum conservation is unspecified. This is because interference terms make it difficult to see how recoil should be distributed.
There are no such issues in angular ordered showers. In this case, the naive guess for how to conserve momentum longitudinal to a jet is correct and is sufficient for the computation of NLL DGLAP evolution and jet physics [26][27][28][29][30][31]. We can exploit this to constrain the form of the recoil ( dR n ) so that the PB algorithm is consistent with an angular ordered shower. In this section, we will derive a dipole shower with this constraint in place from the outset. The resulting dipole shower is very similar to the dipole showers that are commonplace in event generators [2,3]. However, it has a crucial difference: it does not use Catani-Seymour dipole factorisation [32].
To derive the dipole shower proceed as follows.
1. Expand Eq. (2.3) in powers of the number of colours N c and keep only the leading terms, which go as α n s N n c , see [11,33]. This is necessary as only in the leading colour limit can we write evolution equations for |M n | 2 . For the same purpose, spin average the evolution, see Appendix B for details.
2. The colour expansion reduces the evolution equation so that it only depends on dipoles formed by colour connected partons. We use the form of dR n to partition each dipole into two parts, introducing longitudinal momentum conservation to each part of the dipole in such a way that it is exactly consistent with the angular ordered shower. This is similar to how dipoles are usually partitioned using Catani-Seymour dipole factorisation. This partitioning allows us to exchange the sum over dipoles with a sum over emitting parton colour lines.
3. Use the dipole partitioning to restore the (full-colour) hard-collinear physics that is correctly computed by an angularly ordered shower. This is uniquely determined by how longitudinal recoil is assigned. The result is a dipole shower with NLC and NLL sensitivity to collinear physics.
In Appendix 2.3 the complete proof is presented. The final result, expressed in the colour flow basis, is where σ is a colour flow and σ/n is the same colour flow but with the nth colour line removed. We use i c n to index the (anti-)colour line(s) of parton i in a final state dressed with n soft or collinear partons, i.e. if parton i is a quark it has a single colour line and so i c n = i q n , if parton i is a gluon it will have a colour and an anti-colour line so i c n = i g n , iḡ n respectively. i c n is the (anti-)colour line connected to i c n . Momenta with colour line indices are the momenta of the partons associated to that colour line, i.e. p i c n = p in . The shower is ordered in dipole p T , defined as (2.11) The dipole splitting functions are These splitting functions are related to those in the previous section according to P gg (z) = P gg (z) + P gg (1 − z), and P qq (z) = P qq (z). Note that to simplify Eq. (2.10) we have omitted the possibility of g → qq transitions, which is sub-leading in colour and only contributes a leading logarithm to single-logarithm, collinear-sensitive observables. In Appendix 2.3 we present Eq. (2.10) with this splitting included. Being explicit, we would write the squared matrix element as As for the angular ordered shower, this is the squared matrix element for a hard process dressed with n strongly-ordered partons with a unique branching topology, i.e. each emitted parton is specified by a triplet (z m , q (i c m ,i c m) m ⊥ , φ m ) and is emitted from the parton with colour line i c m . The dipole recoil function is given by where and where , and T = in p in . (2.14) Note, in the limit that q n is collinear to p i c n , Asym i c n i c n (q n ) = 1/2. Thus, in this limit R dipole i c n → R i c n , as required. Our expression for R dipole i c n should be compared to the recoil function one would find using Catani-Seymour dipole factorisation: i c n if we were to make the replacement T → p i c n + pīc n . Observables are computed after summing over emission topologies: (2.16) There are several noteworthy points involved in this derivation: • This shower was built around preserving the beneficial features of an angular ordered shower. In fact, azimuthally averaging the dipole shower reinstates an angular ordering. Angular ordered showers provide a sufficient framework to resum global two-jet dominated observables up to α n s L 2n−1 terms with full colour. This will also be achieved by the dipole shower presented here, fixing much of the sub-leading colour errors noted in [10].
• Traditional angular ordered showers fail to correctly compute α n s L 2n−1 logarithms for n > 2 jet observables. This is because soft, wide-angle physics is miscalculated because of the Θ(θ jn,in − θ n,in ) ≈ 1 approximation, as previously discussed. 6 It is never necessary to make this approximation in the dipole shower since we can use the underlying colour flows to define variables in the relevant dipole frame, for which Θ(θ jn,in − θ n,in ) = 1 is always true. Thus we expect the dipole shower to have the capacity to re-sum α n s L 2n−1 logarithms at leading colour. 7 • In the soft limit the dipole shower generates iterative solutions to the BMS equation [16] (the proof is as in Section 3 of [11]). This demonstrates that the shower computes non-global logarithms at leading colour correctly. 6 Modern implementations of angular ordered showers do use colour flow information from the hard process, allowing them to compute α n s L 2n−1 terms at leading colour for global n > 2 jet dominated observables by deriving appropriate initial conditions from the respective large-N colour flows of the hard process [5]. 7 Eq. (2.10) as it stands only provides a sufficient framework for this resummation. It is not yet sufficient in itself: one would need to enhance the shower with a running coupling and, possibly, higher order splitting functions.
• At this point in our theoretical development, the dipole shower does not completely conserve energy and momentum. Rather it only conserves momentum longitudinal to the emitting parton. Accounting for total energy-momentum conservation is not needed to compute some observables to NLL accuracy, e.g. thrust. Regardless, it is an important effect that if handled incorrectly can spoil the NLL accuracy of the shower [10]. Addressing this is the focus of the next section.

Improving recoil in dipole showers
In this section we will address the problem of energy-momentum conservation in a dipole shower, though our approach is simple to map onto an angular ordered shower. The mechanism for energy-momentum conservation (or recoil scheme) we present lacks a formal derivation. Rather it is inspired by the study of recoil by Bewick et al. [17]. Bewick et al. analysed several approaches to recoil in angular ordered showers, reproducing some of the fixed-order checks of [10] and performing further numerical checks. They observed that among the better performing recoil schemes are globally defined schemes; schemes that redistribute momentum across an entire jet or event. From our perspective, a global Boost to ZMF to conserve energy rescale momentumP J Figure 3. A summary of the dipole shower global recoil scheme (a scheme for energy-momentum conservation). In words: A new particle is emitted which leaves some momentum unbalanced (in the direction of the colour connected parton and in the plane transverse to the dipole); perform a Lorentz boost to the new ZMF, and re-scale the jet momenta in such a way that the rescaling does not change the k ⊥ of the emission. This leaves an n-parton ensemble with the same total energy and total momentum as the n − 1-parton ensemble.
scheme is also preferable, as it is more simply implemented in a dipole shower. Momentum conservation on an emission-by-emission basis is also desirable when it comes to matching to fixed-order and merging of hard processes of different jet multiplicity. In the two-jet limit, our scheme becomes that which is analysed in [17] and implemented in HERWIG's angular ordered shower [5]. For comparison, in Appendix C we summarise the implementation and limitations of a spectator recoil scheme, as implemented in [2,3,34].
We start with an observation that is key to all global recoil schemes: when a parton is emitted from another, the parent parton must have been off-shell. We parametrise the amount by which it is off shell by giving it a virtual mass. A parton shower approximates the sum over the multiplicities of QCD radiation dressing a given hard process. Each term in the sum should have the same total energy and the same zero-momentum frame (ZMF). Naively adding a parton to an n − 1 on-shell parton state changes the total energy and ZMF. We will redistribute parton momenta as simply as possible in order to restore the ZMF and total energy. We will do this using a single global Lorentz boost and a single rescaling that preserves the transverse momentum ordering. This procedure is illustrated in Figure 3. Below we will spell out how to implement this recoil scheme. The simplicity of the scheme can get lost in its mathematical definition and so we encourage the reader to keep Figure 3 in mind.
Let us now make Figure 3 quantitative. We require that energy is conserved, and that momentum is conserved where, in the ZMF, P J is the 3-momentum of Jth jet amongst the n − 1 jets constructed from an n − 1 parton ensemble, i.e. P J = p in for J = i n (recall that i n labels parton i in an n-parton ensemble). We introduce the extra notation because it is the momentum of jets that we particularly focus on conserving.P J is what we wish to find; it is the momentum of the Jth jet now constructed from an n parton ensemble after the necessary redistribution of momenta (all jets contain a single parton except for one which contains two partons; the original parton and the newly added parton). m i is the mass of parton i, and m i = 0 since we consider only massless partons.P 2 J is the virtual mass squared of the Jth jet, it also is zero for all jets other than the jet built of two partons. We can achieve our desired redistribution by a Lorentz boost, Λ µ ν , from the ZMF of the n − 1 parton ensemble to the ZMF of the n parton ensemble. Once in this frame, we re-scale all the jet momenta by a global factor κ i c n (the index will prove necessary later on) so as to preserve the centre-ofmass energy. In all, we wish to findP J µ = κ i c n Λ ν µP j ν whereP j is the four-momentum of the Jth jet constructed from the n parton ensemble before the redistribution of momenta. We place a hat on all intermediary kinematic variables (i.e. those after the emission but before redistribution of momenta). We denote the 3-momentum ofP J asP J = κ i c n ΛP J . Λ µ ν is specified by solving Eq. (3.2) and κ i c n is specified by solving which comes from requiring E before = E after = Q. We will express this recoil scheme in terms of the shower kinematics and solve forP J . We use the following Sudakov decomposition for a 1 → 2 (p i c n →p i c nq n ) parton transition: We label the jet in which the splitting takes place as P J emit , so that P J emit = p i c n . From Eq. (3.4): For all jets other than "J emit"P J = P J andP 2 J = 0. The Lorentz boost, Λ µ ν (i c n , i c n ), can now be found. The boost is in the direction of p i c n and is given by the boost velocity Finally we must solve for κ i c n using Eq. (3.3), Note that in both the soft and collinear limits κ i c n → 1. So now we have everything we need to computeP J = κ i c n ΛP J . We can put this in the dipole shower by introducing a recoil function where δ 4 J (f (p i c n )) is a delta function normalised against its Jacobi factor: where X is the (unique) solution to f (X) = 0. Note that in an implementation of the algorithm there is never any need to invert the argument of the delta function to solve for p i c n sincep i c n is what is needed going forwards. In Eq. (2.10), the delta functions simply kill all of the integrals over p jn . For the sake of being explicit, the emitted parton has momentum Note that both z n and dq are Lorentz and jet scaling invariants. This means that all of the emission kernels remain unchanged and so the implementation of this recoil scheme only enters so as to ensure that the real emissions continue to be integrated over the correct phase-space (and through the corresponding Θ on shell for the virtuals). Thus, the complete dipole shower is defined by Eq. (2.10) 8 , Eq. (2.12) and Eq. (3.7).

NLC and NLL accuracy of the global recoil
In this section we will discuss the colour accuracy of our new dipole shower and test its logarithmic accuracy.
Firstly, the sub-leading colour contained in the shower is inherited from its link to angular ordered showers. It is reasonably simple to see that in the two-jet limit of an exponentiating, global observable the dipole shower assigns colour factors exactly as an angular ordered shower would (that is to say correctly for NLLs). Outside of this limit, leading colour accuracy is guaranteed but it is observable dependent as to how much is achieved beyond this. This is an improvement on existing dipole showers, which have been noted to incorrectly compute NLC even at LL accuracy [10]. Further improvements on sub-leading colour, for more general observables, require amplitude evolution. We doubt that substantial further improvements in the accuracy of sub-leading-colour effects can be achieved in either the dipole shower or coherent branching frameworks. There is already a body of literature exploring possible resolutions to the NLC errors in dipole showers [35][36][37]. Our approach of using angular ordering to improve dipole evolution is similar to that of [35,36], though there it was largely explored only in the context of hadronisation and the computation of jet multiplicity observables. We also note that, by partitioning dipoles so as to identify a unique parent, we expect the sub-leading logarithms associated with unresolved soft and collinear radiation to be captured using the CMW scheme for the running coupling [18,20].
We will now proceed to evaluate the logarithmic accuracy of the recoil scheme discussed in the previous section. We do so in two ways. Firstly by re-creating the analysis of Section 4.2 in [10]. In this analysis, several event shape observables, defined by functions V ({p}), are considered at fixed order. The analysis tests the sub-leading contributions from the soft region found in the limit that the transverse momentum of the second emission is of similar magnitude to that of the first but both are small relative to the hard scale. This limit is considered because it is the limit where dipole showers have previously been shown to mishandle recoil. Specifically, we calculate the difference between the α 2 s LC, NLL contribution to the observable using the fixed-order amplitude, and the shower contribution: As the observables exponentiate, we are looking for differences of the form α 2 s N 2 c L 2 at fixed coupling since these terms contribute to the NLL exponent.
Our second check of logarithmic accuracy is to compare against two known NLL resummations: Thrust and generating functions for jet multiplicity. This is done in Appendix D.
Let us proceed to compute δΣ(L) in the doubly-soft limit in e + e − → qq. We label the quark as parton a and the anti-quark as parton b. In the same way that we label partons with indices i n , each parton label is given a subscript stating the 'current' multiplicity of radiated partons (since a parton's momentum changes to conserve momentum as more partons are radiated). From Eq. (2.10) we can compute the first two soft emissions and find where σ n H is the hard process cross section. θ injn is the product of theta functions defining the on-shell requirements for emission from dipole i n j n (previously given without indices as Θ on shell ). {p} correct are the momenta used to compute Σ Before considering any specific event shape, we can simplify our expressions further by using the recoil delta functions to perform some of the integrals. These fix the final state momenta: q 1 and q 2 are defined with respect to the rescaled momentap a ,p b and so have appropriately modified limits on their phase space. We employ the 'equally soft' limit (Q q 1 ⊥ , q 2 ⊥ ; q 1 ⊥ q 2 ⊥ ) which reduces the complexity of the phase space limits and removes dependence on longitudinal recoil. In total, we find that In the 'equally soft' limit we are considering The κ dependence in the shower integrals (lines 1 and 2 of Eq. (3.11)) causes potentially incorrect O(q 2 ⊥ /2Q 2 ) terms in the phase space limits. 9 These integrate to give dilogarithms in q 2 ⊥ /2Q 2 which do not contribute α 2 s L 2 terms but rather α 2 s L 0 terms that go to zero in both soft and collinear limits. 10 Thus, with NLL accuracy, Eq. (3.11) reduces to Note that δΣ(L) is only non-zero because {p} 2 = {p} correct . We will now consider several specific observables, still following [10]. Dasgupta et al. first consider the two-jet rate in the Cambridge algorithm. They argue that for this observ- and that the α 2 s N 2 c L 2 terms are correctly computed. Similarly, V ({p} 2 ) is also equal to the correct measurement function (up to neglected poly-logs) for the 'fractal moment of energy-energy correlation' (FC 1 ) which, in the soft-collinear limit, is given by V ({p i } correct ) = i p i ⊥ . In the limit we 9 The algebra to show this is awkward but as κi n is simply a ratio of energies, we can argue that it must be an even polynomial when expanded in small q ⊥ . 10 The recoil terms in these integrals are reducible to a few general forms. One such form is where a parametrises the observable, x ∼ q ⊥ /Q and parametrises the coefficients to the O(q 2 ⊥ /2Q 2 ) effects from our recoil scheme; = 0 gives the leading log result. Note that all terms other than the LL result are not enhanced in the a → 0 limit. See Appendix D.1 for more details.
. In fact, it will be the case that for all observables built from Lorentz invariant and jet rescaling insensitive quantities 11 our recoil scheme is sufficient for the computation of α 2 s N 2 c L 2 terms. This being because the scheme is constructed by a Lorentz boost and a formally sub-leading re-weighting. We expect that for suitably simple observables this accuracy will also extend to higher orders, see the resummations in Appendix D. This discussion should be contrasted with that in Appendix C, where we perform the same tests with a spectator recoil scheme [2][3][4]. In agreement with [10], we find that with such a recoil scheme these . This generates NLL errors.

Conclusions
Starting from a general algorithm designed to capture both the soft and collinear logarithms associated with the leading infra-red singularities of scattering amplitudes, we have derived an angular ordered shower and a dipole shower. Our dipole shower is novel in the way that it partitions each dipole in order to account for longitudinal momentum conservation. This partitioning is constructed so as to ensure that the shower implements longitudinal momentum conservation in precisely the same way as the angular ordered shower does. This new dipole partitioning is similar to, but not the same as, Catani-Seymour partitioning. We complete our dipole shower by specifying the transverse recoil. The result is a new dipole shower that formally represents an increase in accuracy when compared to the traditional parton shower models employed by many current event generators [2-6, 8, 38, 39]. For example it will correctly compute the effects of coherence and jet recoil, and the leadinglogarithm, leading-colour contribution from non-global logarithms. It should also reproduce the correct leading-colour, wide-angle, soft radiation pattern beyond the two, three, and four-jet limits. To our knowledge this is not achieved by other parton shower models.
However, our shower still has substantial limitations. In large part that is because it is based on a cross-section-level, semi-classical picture. Operating at cross-section level necessitates that the shower be defined only at leading colour. Full-colour resummation means a more complicated, amplitude-level, approach [12,40,41]. Certainly it would be of considerable interest to compare a parton shower defined at amplitude level, such as the CVolver shower that is currently under construction [42,43] or the Deductor shower [44], with the improved dipole shower we present here.
UK Science and Technology Facilities Council for the award of a studentship. Thanks also to Jack Helliwell, Mike Seymour and Andy Pilkington for comments and discussions.

Additional comment on PanScales
Whilst this paper was being finalised, a study of NLL accuracy in parton showers was released by the PanScales collaboration [45]. There would seem to be a fair degree of similarity between the dipole shower we derive and the PanGlobal shower with β = 0 presented in [45]. Our recoil schemes in particular are similar, and the dipole partitioning they employ appears to be a limit of the one we derive (they are equal when both partons in the dipole have the same energy in the event ZMF).

A.1 Amplitude evolution detailed definitions
Before we proceed with the technical details of the PB evolution, it is necessary that we properly introduce the notation we will later be relying on. In these appendices we will often find ourselves manipulating expressions relating states of differing parton multiplicities (for instance Eq. (2.3) relates an n H +n−1 state to a an n H +n) state. We must label partons and the multiplicity of state they come from carefully since the state's multiplicity determines both the dimension of the colour-helicity space in which the state resides and the momenta of the constituent partons. To this end, we label partons with indices i n , j n , k n , ... which run as i n , j n , ... ∈ {1 H , 2 H , ..., n H } ∪ {1, 2, ..., n − 1}, where {1 H , 2 H , ..., n H } is the set of hard partons and {1, 2, ..., n − 1} the set of partons emitted during the evolution. We use υ in ∈ {q, g} to label the species of a parton i n . The momentum of the i th parton in a state of multiplicity n H + n − 1 is p in ∈ {p} n−1 = {P 1 , P 2 , · · · P n H , q 1 , · · · q n−1 }. The emission operator, D n , adds a new (nth) parton, of four-momentum q n , to the state. After considering energy-momentum conservation, the parton momentum, q n , is added to the set {p} n−1 , to produce the set {p} n . dR n acts in conjunction with D n to map {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). Following this, {p} n = {p} n−1 ∪ {q n } is the set of n momenta including the last emission, q n .
Many of the objects used in this paper carry complicated dependencies. To simplify some lengthy expressions, we will only provide the full list of arguments in an object's definition. In definitions, we will write every object as some f (x; {y}), where x is the evolution variable on which f depends and the set {y} itemises the complete dependences of f . In all expressions subsequent to the definition we will drop the {y} dependence and only write f (x). We can do this safely as, following the initial definition of an object, each object can always be uniquely determined by the subscripts and superscripts we provide.
In Section 2.1 we gave an overview of the roles of D n , dR n and Γ n . Let us now define these operators more carefully 12 where O is some operator in the colour-helicity space and where we have used a shorthand notation to help save space Delta functions of this form are used to carry the frame dependence of the ordering variable in a compact form. Physically, S in n emits a soft parton from the parton labelled i n . These soft partons take the form of interference terms in the evolution. C in n emits a collinear parton from the parton labelled i n . The following two definitions for transverse momenta are used as ordering variables for soft and collinear emissions respectively, where n is a light-like reference vector. The choice of n is determined by how recoil is handled in the evolution and is often taken to be in the backwards direction relative to p in . Strictly speaking, recoil cannot be entirely factorised from each D n however the way in which it acts in each D n follows a simple pattern. Thus we have used the recoil measure dR n as an abridged notation. It is defined to act by the following rules R soft injn and R coll in contain the necessary delta functions and kinematic pre-factors needed to account for energy-momentum conservation. They are discussed in Section 3 and further examples are given in [12]. 13 Explicit expressions defining S in n and C in n are lengthy and can be found in [12]. Finally, indicates that the enclosed operators should act on any incoming partons as if they were in the final state (see Eq. A.1 in [12], which defines the operators from which D n+1 12 For pedagogical reviews of the colour-helicity operators relevant in the definition of these operators see [11,12,46]. 13 In [12] R soft injn and R coll in are written as R soft * n jn R soft n in and R coll * n in R coll n in respectively. . Θ on shell is our short-hand notation for the inclusion of the theta functions necessary for restricting the range of integration to the phase-space for an on-shell parton. These are also specified fully in [12] (see functions θ ij and θ i in Section 2).δ i n+1 j n+1 = 1 if both partons i, j are incoming or both outgoing andδ ij = 0 otherwise.

A.1.1 Computing observables
In the main text our focus is on dressing e + e − → qq. The formalism is more general and can be used to compute observables in hadron-hadron collisions using where f υ i (x i , µ) are the parton distribution functions (PDFs) with momentum fractions x i and u n ({p} n , {v}) is the (n H + n)-body measurement function for an observable described by parameters v i ∈ {v}. Note that Σ is differential in hard process momenta, and that it should be multiplied by the necessary flux factors as necessary. The star operation is defined in Section 4 of [12] but in essence assigns PDF type to a given partonic leg (gluon or quark). In this paper, every concrete use of our formalism concerns the showering of an e + e − hard process and so we will not expand further on the treatment of DGLAP evolution.

A.2 Derivation of the angular ordered shower
This section derives an angular ordered shower from Eq. (2.3). It is split in three parts. Part one forms the main derivation, however it will state some results without proof (when these results are themselves laborious to prove). The subsection following presents the limitations of this derivation. Finally the last subsection fills in the gaps. We will focus on e + e − → qq as the hard processes, and at the end we will sketch the extension to other hard processes.
We begin with the amplitude evolution equation, Eq. (2.3), and introduce an azimuthal averaging operation 1,...,n which averages the lab frame dipole azimuths of partons 1 to n, i.e.
Implicit in this operation is also spin averaging when acting on spin-dependent operators, as discussed in Appendix B. To keep things simple, we will proceed in this section without discussing any dependence on the observable, which means we are implicitly assuming the observable is not a function of the parton azimuths. We devote the next sub-section to addressing this. After averaging, Eq. (2.3) becomes where s in,jn and c jn are the spin-averaged kinematic factors associated with a soft or collinear emission respectively (they will be manipulated into the form of collinear splitting functions shortly). They are defined through the relations We observe that c jn = c jn n provided R coll in is independent of the emission's azimuth (spin correlations provide the only azimuthal dependence for collinear emissions). In Section A.2.2 we show that 14 δq (injn) n ⊥ (q n ⊥ ) s in,jn n = − in d 4 p in P injn φ n,in Θ on shell φ n,in + P jnin φ n,jn Θ on shell φ n,jn R soft injn + O(1), where P injn φ n,in = Θ(θ jn,in − θ n,in ) 1 − cos θ n,in . (A.10) The angles in Eq. (A.10) are defined in Figure 2. Θ on shell φ n,in contains the necessary theta functions to constrain the phase-space of parton q n so that it is real and on-shell, encoding the phase-space limits for energy conservation. Its lengthy definition can also be found in Section A.2.2. The presence of the functions P injn φ n,in enforces an angular ordering, secondary to the k ⊥ ordering. To bring this ordering to the fore, we now change variables: 14 Under the assumption that R soft injn is independent of the azimuth up to O(1) terms, which is true for the two recoil schemes we discuss in this paper. and define ζ n,in = 1 − cos θ n,in . In these new variables Eq. (A.7) becomes 15 ζ ∂ A n (ζ) 1,...,n ∂ζ ≈ −Γ n (ζ) A n (ζ) 1,...,n − A n (ζ) 1,...,n Γ † n (ζ) jn P jn n Θ on shell φ n,jn R col jn T jn A n−1 (q n ⊥ ) 1,...,n−1 T † jn ζ n,jn δ(ζ − ζ n,jn ). (A.11) Here we have used c jn n ≈ P jn n Θ on shell φ n,jn R col jn , where P jn n is a sum over collinear splitting functions with the soft divergences subtracted away, e.g. when j n is a quark, P jn n (z) = (1 − z)P qq /2 where P qq (z) = −(1 + z). The details can be found in Appendix A of [12]. We will formulate the evolution in terms of the full splitting functions once equations have been reduced enough that it becomes convenient to do so. Using the strongly ordered approximation, ζ 1 ζ 2 ... 16 , Also using strong ordering, the leading part of Θ on shell φ n,in does not depend on j n and R soft injn can be chosen so that its leading part can be factorised from the sum over j n as 1 ζ n,in Θ on shell φ n,in R soft injn ≈ 1 ζ n,in Θ on shell φ n,in R col in .
Using these simplifications we can apply colour conservation and, by re-labelling indices, write ζ ∂ A n (ζ) 1,...,n ∂ζ ≈ − Γ n (ζ) A n (ζ) 1,...,n − A n (ζ) 1,...,n Γ † n (ζ) + in d 4 p in × jn P jn n Θ on shell φ n,jn + 2 1 ζ n,jn Θ on shell φ n,jn R col jn × T jn A n−1 (q n ⊥ ) 1,...,n−1 T † jn ζ n,jn δ(ζ − ζ n,jn ). (A.13) 15 Γn(ζ) is defined as Γn(q ⊥ ) after the change of variables has been made rather than naively swapping out the argument. 16 When working in a frame that ensures i and j are back to back, the theta function is saturated without approximation. In this derivation we are concerned with e + e − → qq. Thus we can saturate the theta function for emissions from the primary hard partons, so that they are handled without approximation. This means we pick the backwards direction (n) (used to define kinematic variables for emissions in a jet) to be in the direction of the other jet. This in turn fixes the definition for the momentum fraction used in later equations: zn =p j ·n p j ·n . When working beyond the two-jet limit, tricks can be played to further saturate the theta function using knowledge of the hard process colour flows.
By recognising the evolution will become entirely colour-diagonal once the trace is taken, we can diagonalise the colour structures. In turn this allows us to group the soft evolution kernels and the collinear ones into splitting functions. We find × Θ on shell φ n,jn R col jn Tr A n−1 (q n ⊥ ) 1,...,n−1 ζ n,jn δ(ζ − ζ n,jn ). (A.14) P υυ jn (z n ) are the usual DGLAP splitting functions, e.g. P qq (z n ) = C F 1+z 2 n 1−zn . Here we have used υ jn to label the species of parton j n and υ to label the state j n transitions to; if υ jn = q then υ = q and if υ jn = g then υ = q, g. z n is the momentum faction of parton n, i.e. if we have a collinear splitting that induces j n−1 → j n n then p jn ≈ z n p j n−1 and q n ≈ (1 − z n )p j n−1 . We specifically require that z n =p jn ·n p jn ·n where n is a light-like vector pointing along the primary axis of the jet from which parton j n does not stem.
We can make connection to squared matrix elements by letting from which we find the evolution equation for a final-state angular ordered shower with a conventional phase-space for a coherent shower in dz. After which, Eq. (A.14) can be written as in Eq. (2.7) after |M n | 2 1,...,n → j 1 ,...,jn |M n | 2 1,...,n .
(A. 16) We can start by considering the effects of only averaging over the nth parton and use the following identity where σ n (x) = x 2 n − x 2 n and Cor n (x(φ n ), y(φ n )) is the correlation function of x and y under the variation of φ n . Both |Cor n (|M n | 2 , u({p} n ))| and σ n (u({p} n ) are smaller than unity 17 . Next we can consider averaging over both the nth and (n − 1)th partons: |M n | 2 u({p} n ) n−1,n = |M n | 2 n u({p} n ) n n−1 + σ n (|M n | 2 ) σ n (u({p} n ) Cor n (|M n | 2 , u({p} n )) n−1 , This can be iterated to give We have been slightly lazy with notation; it is implicit that The important question is whether the correlations can provide a logarithmic enhancement to the observable. This is obviously an observable dependent statement. To progress we will place some assumptions on the observable. If the observable is such that the correlation term's contribution to the cross section is suppressed relative to |M n | 2 m u({p} n ) m , we can approximate |M n | 2 u({p} n ) 1,...,n by only keeping the first order correlations, since second order correlations will necessarily be even further suppressed. The approximation assumed by coherent branching is to neglect correlation terms altogether. Let us look at the n = m = 1 term for thrust. At this order u({p} n ) is not a function of the azimuth and so σ 1 (u({p} 1 )) = 0. As the observable exponentiates [18,19], this is sufficient to guarantee that it can be computed to NLL using the coherent branching formalism (these last two sentences are an abridged form of the argument in [19]). For contrast, let us look at the n = m = 2 term in the computation of gaps-between-jets, with the same hard process. The pertinent measurement functions are where Θ in/out (q m ) is unity when parton m is in/out the rapidity region between the two highest p T jets and zero otherwise. In the following subsection, we compute all the ingredients for σ 2 ( |M 2 | 2 1 ). It is reasonably easy to argue (though less easy to compute) that, unless suppressed by multiplicative factors in σ 2 ( u({p} 2 ) 1 ) and correlation functions, σ 2 ( |M 2 | 2 1 ) terms can contribute fourth-order, infra-red poles and with them leading logarithms. By considering the variation of φ 2 , it is also simple to convince oneself that the correlation function must be finite and positive. So, if angular ordering is to adequately describe this observable, it must be the role of σ 2 ( u({p} 2 ) 1 ) to screen against contaminating logarithms. This means we only need to test to see if σ 2 ( u({p} 2 ) 1 ) is non-zero: Furthermore, not only is this non-zero but it contains non-vanishing terms in Θ in (q 1 )Θ out (q 2 ).
While these terms do screen against fourth order poles and logarithms, they are crucial for the computation of the α 2 s L 2 non-global logarithms. As such, a coherence branching algorithm (that makes usage of azimuthal averaging) cannot compute the leading logarithms to gaps-between-jets, as it certainly gets the numerical coefficient to non-global pieces incorrect. This is a general feature: coherent branching will fail to capture leading, non-global logarithms (though in most cases these logarithms are sub-leading in the computation of the overall cross section). This has been previously observed in [24,25], where the effect of the missing correlations was computed numerically to all-orders. They found that, though the missing correlations are a formally leading effect, phenomenologically their effect is < 10%. As is widely known, we observe that coherent branching is always capable of calculating logarithms up to α n s L 2n−1 in observables for which α n s L 2n is the leading logarithm.

A.2.2 Azimuthal averaging
In this appendix we will fill in the details on the azimuthal averaging of the evolution kernels. The general procedure for azimuthal averaging is well known [20] textbook material [26,47]. However, the procedure is less widely discussed taking into account phase-space limits and momentum maps. In this section we provide a more careful treatment than the textbook one. We begin by looking at the following integral (which corresponds to the integrated soft emission spectrum), where E qn is the energy of parton q and dΩ qn is solid angle in the frame which E qn is measured. We can regroup the dipole kinematics as 2P injn = n in · n jn − n in · n n in · n n jn · n + 1 n in · n , (A. 24) where n in = p in /E in . The two terms in this integral are symmetric under the exchange of i and j and so we shall focus only on the first: To compute this the integral we take n in = (1, 0, 0, 1), n jn = (1, sin θ jn,in , 0, cos θ jn,in ), and n = (1, sin θ n,in cos φ n,in , sin θ n,in sin φ n,in , cos θ n,in ). In this basis The textbook treatment would set Θ on shell = 1 here. For us, (1 − sin θ n,in cos φ n,in sin θ jn,in − cos θ jn,in cos θ n,in ) which bounds the φ n,in integration to the range |φ n,in | ∈ [φ − q,i , φ + q,i ). The solutions for the boundaries, φ ± q,i are given by cos φ ± q,i = ± min |α ± |, 1 for α ± > 0 and cos φ ± q,i = 0 otherwise, Note that the expression under the square root is always positive. The usual approach to azimuthal averaging is to employ the soft limit and set Θ on shell = 1, after which the φ n,in integral can be performed by contour integration. However, in our case this is not viable, due to the boundaries on the φ n,in integral. Instead we will write the integral as where Cor(x, y) is the correlation function between two variables x and y, in context the correlation over variation of the azimuth. Firstly note that the usual result from azimuthal averaging. We can also note that Θ on shell φ n,in ∈ [0, 1] and |Cor(P injn , Θ on shell )| ∈ [0, 1]. By brute-force evaluation and noting Θ on shell is binomially valued, we find whereθ on shell = Θ on shell φ n,in =φ crit , and cos φ crit = sign(f )min (|f | , 1) , f (θ n,in , θ jn,in , E in , E jn , q ⊥ ) = The exact angular ordered result is obtained when Θ on shell φ n,in =θ ij = 1, which is the case in the strongly ordered, q ⊥ /Q → 0, and collinear, θ n,in → 0, limits (here Q stands in for any other harder invariant). The remainder of this section is used to show that the correlation term can be neglected at least at α n s L 2n−1 accuracy (and for NLL thrust). It can be skipped if the reader does not need convincing.
injn φ n,in = dφ n,in 2π P 2 injn = dφ n,in 8π n in · n jn − n in · n n in · n n jn · n + 1 n in · n 2 , = 1 (n in · n) 2 dφ n,in 8π cos θ n,in − cos θ jn,in 1 − sin θ n,in cos φ n,in sin θ jn,in − cos θ jn,in cos θ n,in + 1 2 , (A.32) using the substitution z = exp(iφ n,in ) this integral equals cos θ n,in − cos θ jn,in 2z − sin θ n,in (z 2 + 1) sin θ jn,in − 2z cos θ jn,in cos θ n,in + 1 2z z(cos θ n,in − cos θ jn,in ) sin 2 θ n,in sin 2 θ jn, where z ± = 1 − cos θ jn,in cos θ n,in sin θ n,in sin θ jn,in ± 1 − cos θ jn,in cos θ n,in sin θ n,in sin θ jn,in Only the z = z − and z = 0 poles are in the unit circle: (cos θ n,in −cos θ jn,in ) 2 (1−cos θ n,in ) 2 (cos θ n,in −cos θ jn,in ) + (cos θ n,in −cos θ jn,in ) 2 (1−cos θ n,in ) 2 (cos θ n,in −cos θ jn,in ) + 1 4(1−cos θ n,in ) 2 when θ n,in < θ jn,in , This has a collinear divergence that is suitably screened in Eq. (A.30) by the accompanying phase space factor, Θ on shell φ n,in (1 − Θ on shell φ n,in ), as is the soft divergence from the q ⊥ pre-factor in Eq. (A.25). Cor(P injn , Θ on shell ), is bounded above and below by 1 and −1 so at most further dampens the effect of the σ 2 P injn term. As a result it is a finite non-logarithmic correction at order α s and its contribution is suppressed at higher orders (to be seen explicitly one could repeat the analysis of Appendix D.1). Hence, for α n s L 2n−1 accuracy, we need only take the first term on the right hand-side of Eq. (A.30).

A.3 Derivation of the dipole shower
In this section we will derive from Eq. (2.3) an evolution equation for a dipole shower for final-state coloured radiation in e + e − . The extension to an initial state shower does not add complexity but lengthens equations. To derive the dipole shower we will spin average the evolution and make the leading colour approximation. To approximate the colour, we express amplitude density matrices and colour charge operators in the colour-flow basis. We manipulate the colour-flow basis using the mathematical machinery introduced in [11].
Before we begin the derivation let us look at Eq. (2.3) in more detail and apply some of the knowledge we have gained from deriving an angular ordered shower. Angular ordering is most powerful when applied to the two-jet limit in e + e − , the mono-jet limit of DIS and Drell-Yan. In these cases, angular ordering does not approximate the soft radiation pattern at all. Instead, the soft radiation is colour diagonal. The diagonalisation of soft radiation renders the conservation of momentum longitudinal to a jet unambiguous. Matching to the angular ordered limit is sufficient to completely constrain the leading component of momentum conservation in Eq. (2.3) (it must respect the partitioning defined by P injn as given in Appendix A.2). It is required that where T = in p in is a vector for projecting out the energy of a parton in the event ZMF and where E n is the energy of q n in the ZMF. This can be rearranged to give As previously stated in our discussions on angular ordering, This recoil function is ready to use in Eq. (2.3). Now, let us begin computing the leading colour evolution of A n (q ⊥ ). We intend to compute where A (0) τ σ n is the leading colour amplitude for colour flows τ and σ, see [11,33] for details on this procedure. Term by term in Eq. (2.3) we can apply this operation and find and where The sum over "i n , j n c.c. in σ" standards for performing the sum over partons dipoles i n , j n which are colour connected in the colour state σ. P υ in →υ,υn ≡ P υ,υ in are the hard-collinear splitting functions defined in Appendix A of [12]. They are the usual collinear splitting functions with soft poles subtracted away, i.e. P qq = −C F (1 + z n ). Note that as we are working in the strict leading colour limit C F = N c /2. The constants λ in andλ jn are defined in Table 1 of [11], in the situations we will use them (the LC limit) λ inλjn → 1/2. We can observe that the first term on the RHS of Eq. (A.44) is of the same form as the standard dipole type term. Next we can take the leading colour part of the emission operators. We spin average emission kernels, see Appendix B for details, and place carats on objects to remind us that they are spin-averaged. We find Note thatγ (σ) = γ (σ) as the loops do not depend on spin.
For now we will ignore the single logarithmic, hard-collinear pieces as they are easy to introduce later on (they are uniquely attributed to delta functions of the form δ 4 (p jn − z −1 np jn ) in the recoil). This means that for now our final state will simply be the qq pair plus n gluons. It is also typical in the strict LLA to let R soft injn = 1; this will prove to be exact with the recoil scheme given in Section 3 though only approximately so with the spectator scheme in Appendix C. Thus the evolution equation is This is a modified version of the equation for dipole evolution found in [11] that was shown to reproduce BMS evolution [16]. It has been modified to allow for the possibility of kinematic recoil and to account for the phase-space effects from energy conservation. By taking the leading colour limit, the colour evolution has been made diagonal. We can trivially make the connection with squared spin-averaged matrix elements; for a given colour flow, σ, whereM is a dimensionless, spin-averaged and leading-colour matrix element, up to global factors of 2 and π which have been absorbed into the definition of our phase-space measure 18 . Thus This is a generalised leading-colour dipole shower evolution equation with fixed coupling. Commonly one would introduce a running coupling with q ⊥ as its argument, allowing for the computation of further NLLs. At this point this would be a simple extension. We have omitted the running coupling as it does not effect our discussion. From this point on we drop the carat denoting spin averaging, leaving it implicit that the equations are spin averaged.
To manipulate our new dipole construction into the more usual form we now define a recoil function based on colour flows: The usual dimensionful matrix element is retrieved by multiplying with a factor i n+1 2π −1 q −2 i n+1 ⊥ .
where, just as in Section 2.3, we use i c n to index the (anti-)colour line(s) of parton i in a final state dressed with n soft or collinear partons. Using this we can now return to Eq. (A.50) and manipulate the dipoles so that emissions from each half of a dipole are separated: We can now include the sub-leading logarithms from the hard-collinear limit along with full-colour Casimir invariants. The Casimir invariants and collinear logarithms are each uniquely associated with longitudinal recoil and so a single R dipole i c n . We note that Asym i c n i c nn (q n ) gives no logarithmic enhancement in the hard-collinear region, rendering the inclusion of hard-collinear pieces simple (including the re-inclusion of g → qq transitions). Thus we arrive at Eq. (2.10). 19 We can explicitly include the g → qq transitions by extending Eq. (2.10): where P qg (z n ) = n f T R z 2 n . The inclusion of Casimir factors and collinear physics in this fashion ensures our shower correctly computes everything an angular ordered shower can compute. Outside the region of validity for an angular ordered approach, NLC and NLL errors will emerge. However, the usual LC accuracy of a dipole shower is preserved. Also note that at no point in this derivation did we restrict ourselves to a qq final state for the hard process. In Section 2.3 we made this restriction as it allows Eq. (2.16) to be written more simply. For more complex hard-process topologies one should sum over showers originating from each distinct hard-process colour flow (dipole).
So far we have still not constrained the O(q ⊥ /Q) pieces in the recoil function associated with recoil in the backwards direction. These pieces are important for the computation of NLLs. Specifying them is the purpose of Section 3 and Appenidx C. In these sections we 19 When constructing Eq. (2.10) we chose to multiply each matrix element by a phase-space factor so that n | 2 and separate sums over emission topologies, |M This ensures the standard dipole shower phase space can be used [2][3][4]10]. study their effect on NLLs. For contrast, in Section 2 of [12] we considered various recoil functions that specify the O(q ⊥ /Q) pieces. We ensured each possible recoil prescription would consistently produce all leading physics, however we did not check sub-leading effects. One of the prescriptions we considered was based on the spectator recoil commonly employed in modern dipole showers [2,34]. This approach involves partitioning the dipole using Catani-Seymour dipole factorisation [32] and distributing the longitudinal recoil in accordance with this partitioning. The remaining transverse recoil is then given to a third parton, not in the dipole but colour connected to the emitting parton. In [12] we give the functional form of R soft injn necessary to implement this recoil. Using this recoil function instead of the one we present here gives us an evolution equation similar to that governing Pythia8 [2].
In [10] it was shown that the standard spectator recoil prescriptions used in conjunction with Catani-Seymour dipole type showers are subject to errors computing NLLs and miscalculate next-to-leading colour. The errors in NLC occur because of the misattribution of longitudinal components of recoil and so colour factors. The errors in NLLs occur as unphysical artefacts from the shower construction do not cancel when one properly considers the effects of recoil after multiple emissions. It is for this reason that we have taken so much care to ensure consistency between our dipole shower and angular ordered showers, and why we take great care implementing recoil in Section 3.

B Spin averaging
In the derivation of an angular ordered shower and a dipole shower we had to spin average the evolution from Eq. (2.3). We can introduce spin averaging safe in the knowledge that the spin-correlated evolution can be computed from the spin averaged by re-weighting with the algorithm of Collins, Knowles et al [22,48]. In our previous paper [12] we showed that, given collinear factorisation, the evolution of our algorithm is consistent with that of Collins and Knowles et al. We also showed that complete collinear factorisation can be achieved in the PB algorithm (neglecting Coulomb exchanges, which cancel in the leading colour limit). In this appendix we will summarise the spin averaging procedure. We will do so in the leading colour limit, as this is the limit of interest in the dipole shower case and this limit reduces the number of indices on objects. Real emissions in the leading colour limit without spin averaging give rise to where h L/R i is the helicity of the parton with label i on the left/right hand side of the amplitude. In Eq. (B.2) we again used the abbreviation "i n , j n c.c. in σ" to mean that we sum over pairs i n , j n that are colour connected in σ. Note we have been a little sloppy by omitting sums over trivial spin indices of partons not involved in the splittings induced by where we denoted the spin averaged objects with a carat. We have assumed R soft injn and R col in are chosen such that they are not spin dependent, otherwise they too should be averaged in the same fashion.

C Dipole shower with spectator recoil
It is commonplace to use local 'spectator' recoils in dipole showers rather than the global approach we have opted for [2,3]. In this appendix we introduce one such recoil scheme and show that, despite the other improvements to our dipole shower, it suffers the NLL errors pointed out in [10].
Following the approach of [34], we can treat each transition from an n−1 to an n parton matrix element as being generated by a 2 → 3 parton splitting which locally conserves momentum. The splitting is defined such that the parton with colour line i n under goes a where {p} correct is the set of correct momenta for the 4-body system and where θ correct injn = θ injn ({p} correct ). From this we find The kinematics are encapsulated by {p} 2 , just as in the global scheme given in Section 3. They are in fact exactly the same kinematics as those specified in Section 3.3 of [10] and we have arrived at the same expression as B.5 of [10]. Thus, we can follow their argument from Appendix A and Section 4 and conclude that our local dipole prescription does suffer the same NLL errors as other local dipole prescriptions. For example, we can consider the two-jet rate using the Cambridge algorithm, for which V ({p i }) = max i {p i ⊥ }.
In the limit we have considered, this reduces to V ({p} correct ) = q ) since the recoil scheme does not ensure that q at all points in the phase-space for parton 2's emission. [10] show that this error generates a incorrect NLL (N 2 c α 2 s L 2 ). This was expected, as in our local dipole scheme we have only made modifications to fix the NLC of the usual dipole shower procedure. It would be unexpectedly fortuitous if this also fixed the NLL problems.

D.1 Thrust with NLL accuracy using global recoil
Thrust has a long history. It was first resummed to leading log accuracy in 1980 [49] and then later at next-to-leading in 1993 [19]. More recently, it was resummed to N 3 LL [50]. In this section we will analyse the consistency of the dipole shower and recoil scheme we present in Sections 2.3 and 3 with the NLL computation found in [19]. Crucially, the calculation of NLL thrust was performed using a coherent branching algorithm [20] (or equivalently by analytic computation of an angular ordered shower). The coherent branching algorithm employed in the resummation was not strictly momentum conserving and effectively only conserved the momentum longitudinal to the two back-to-back jets. In [19] they show that neglecting the other components is a valid approximation in the computation of NLLs for thrust (see their expansion of the correct phase-space). However, in [10] it was observed that incorrect handling of complete momentum conservation in dipole showers can induce NLL errors in thrust from O(α 3 s ) onwards. These two papers are not inconsistent with each other, the situation is simply that the incorrect re-inclusion of negligible terms can induce large errors.
Our dipole shower algorithm was built around consistency with an angular ordered shower. It returns an angularly ordered shower with correct longitudinal momentum conservation when azimuthally averaged. However there is one main difference between the coherent branching resummed in [19] and our algorithm; ours conserves momentum completely. Thus the only question is whether this complete momentum conservation breaks the NLL accuracy. We can compute the difference between our algorithm's computation of thrust and [19]. As thrust is dominated by the two-jet limit, we initially focus on emissions from the primary hard legs (which is sufficient for NLL accuracy in the approach of [19] by assuming inclusivity over jets from secondary jets). Afterwards we will briefly consider the effects of secondary emissions, i.e. possible recoil effects from the multi-jet limit. Firstly note that thrust can be defined as where P n (Pn) is the total four-momentum in the hemisphere centred on the forwards (backwards) thrust axis. From this definition, it is clear that thrust is invariant under boosts along the thrust axis and is invariant under global jet energy rescaling. Following the notation of Section 3, the difference in the two-jet limit between our dipole algorithm and the NLL result due to recoil is of the general form − ln(κnQ/q n ⊥ ) dy n ...
where each transverse momentum is defined relative to the thrust axis and C n is a constant coefficient.
It is most beneficial to us if we evaluate the logarithmic order of δΣ(L) by starting more generally and then applying the result to thrust. As previously stated, each κ n = 1 − O(q 2 n⊥ /2Q 2 ). We will parametrise this as κ n = 1 − q 2 n⊥ /2Q 2 where is order unity. Note that when = 0, δΣ(L) = 0. Eq. (D.1) is built from repeated sums over elementary integrals of the following type where a parametrises the observable dependence (for thrust a ∼ 1 − T ), x i ∼ q i ⊥ /Q and Θ(f (a, {x i })) parametrises any residual more complex observable dependence. Note that both terms in the square bracket are monotonically decreasing as x i → 0 and that the second is always of smaller magnitude than the first. Thus I evaluates to having the largest possible magnitude when Θ(f (a, {x i })) = 1, as every point in the domain of the integrand adds constructively to the integral. Therefore we will work assuming Θ(f (a, {x i })) = 1 in order to place an upper limit on the order of logarithms produced. With this assumption applied, I is dominated by the term which is in turn proportional to g 2n−2 (a, ) − g 2n−2 (a, 0) where For large n, g n is difficult to evaluate. However we can navigate this by constructing a generating function for g n , so that g n = (∂ ν ) n GF | ν=0 and GF (a, , ν) = a ν − 1 ν 2 + 2 F 1 1, ν 2 + 1; ν 2 + 2; 2 − a ν+2 2 F 1 1, ν 2 + 1; ν 2 + 2; where L = ln(1 − T ), andÃ i,n andB n are order unity constants. Hence for T ≈ 1, the limit in which we resum, δΣ(L) n α n s Cn n! ln(1 − T ) 2n−2 where C n are also order unity coefficients. Also note that the first logarithmic enhancement from our recoil scheme occurs as ∼ α 4 s L 2 . Finally, we note that this argument applies to recoil distributed along any chain of strongly ordered emissions. Therefore recoil from emissions off secondary legs also contributes terms to δΣ(L) that are much less than n α n s Cn n! ln(1 − T ) 2n−2 . 20 We have shown that the recoil scheme for a dipole shower presented in Section 3 does not introduce incorrect next-to-leading logarithms into the resummation of thrust in e + e − . We did this using a very general approach, leading us to believe that for other exponentiating two-jet dominated observables the same result will also be found. We conjecture that our recoil scheme will never contribute incorrect α n s L 2n−2 or greater terms. As we stressed at the beginning of this section, our dipole shower was built around consistency the angularly ordered limit. Therefore, one would only need to add a running coupling and next-to-leading splitting functions and the shower could be used to compute the NLL resummation of thrust. 20 In fact, following the epsilon expansion arguments of [19], recoil from secondary legs will contribute terms less dominant than n α n s Cn n! ln(1 − T ) 2n−4 .

D.2 Generating functions for jet multiplicity using global recoil
We will now use our algorithm with our new recoil scheme (as presented in Section 3) to compute the integral equation defining the spin-uncorrelated generating function for the multiplicity of subjets in the final state of e + e − → hadrons. The generating function was first computed at NLL accuracy (i.e. including all α n s L 2n−1 terms) in [51]. The methodology has since seen a variety of applications [28,31] (and references therein) and can be found in graduate texts [26,47]. We will compute the generating function at LL accuracy, though taking care to include all α n s L 2n−1 logs from recoil. The generating function is defined by Here F is some flux factor for the hard process and P Σ (n, Q) is the probability of finding n partons/subjets in the final state of a process with centre-of-mass energy (or hard-scale) Q. N is the number of partons in the hard process and n Σ is the mean number of subjets in Σ. For e + e − → qq, i.e. computing φ qq (u, Q), it is a textbook result that at our accuracy generating functions factorise as φ qq (u, Q) = φ q (u, τ )φq(u, τ ) where φ a (u, τ ) is the generating function for subjet multiplicity within the jet from a single parton a. τ = 2E sin(θ/2) is the scale of an individual jet and can be thought of as its maximum transverse momentum, E is the energy of each jet and θ the opening angle of the jet, e.g. φ qq (u, Q) = φ q (u, Q)φq(u, Q) as θ = π and E = Q/2 [28,47].
We will now construct an integral equation for φ a (u, τ ). To do so consider also computing φ e + e − →qq[g] (u, q ⊥ 1 ), where the next hardest jet (if one occurs) is a gluon jet of scale q ⊥ 1 . For the computation of φ e + e − →qq[g] (u, q ⊥ 1 ), the hard process can be approximated as H (e + e − →qq[g]) (q 1 ⊥ ) = A 0 (q 1 ⊥ ) + uA 1 (q 1 ⊥ ). Hence φ e + e − →qq[g] (u, q ⊥ 1 ) = F  Born is the Born phase-space for the qq pair 21 . We can rewrite 21 The Born phase-space on the momenta of partons after momentum conservation has been taken into account and includes the momentum conserving delta function δ 4 (Pq + Pq) as well as a delta function fixing the energy. this as φ e + e − →qq[g] (u, q ⊥ 1 ) =φ q (u, q ⊥ 1 )φq(u, q ⊥ 1 )Tr(V q ⊥ 1 ,Q · V q ⊥ 1 ,Q ) + dΠ Born dΠ 1 × dR 1 d 4 P g dφ q (u, q ⊥ 1 ) d 4 P q dφq(u, q ⊥ 1 ) d 4 Pq dφ g (u, q ⊥ 1 ) d 4 P g × Tr V q ⊥ 1 ,Q · V q ⊥ 1 ,Q D † 1 · D 1 1 δ 4 (P g − q 1 ), (D.12) where we have employed azimuthally averaged result of Appendix A.2 since the equation is independent of the azimuth of the gluon. We have also spin averaged at this step. We also note that Eq. (D.12) is equal to φ qq (u, Q) by necessity, i.e. φ qq (u, Q) = φ e + e − →qq[g] (u, q ⊥ 1 ) as within the strong ordering approximation the next hardest jet of an e + e − → qq process must be a gluon jet. After a little work, φ qq (u, Q) = 1 2 φ q (u, q ⊥ 1 )∆ q (q ⊥ 1 , Q) φq(u, q ⊥ 1 )∆q(q ⊥ 1 , Q) dz P qq (z)φ q (u, q ⊥ )φ g (u, q ⊥ ) ≈ φq(u, q ⊥ ), φ g (u, q ⊥ ) = dφ 1 2π d 4 P g dφ g (u, q ⊥ 1 ) d 4 P g δ 4 (P g − q 1 ) ≈ φ q (u, (1 − z)q ⊥ ), (D.14) and where the recoil functions, using the same definitions as Section 3, are given by R primary q 1 = δ 4 J P q 1 − zκ q Λ(q,q)p q , R secondary q 1 = δ 4 J κ q 1 Λ(q,q)P jn −P jn , i.e. eachφ is simply related to each φ by momentum conservation. At our accuracy, momentum conservation simply maps E q → z 1 E q and E g = (1 − z 1 )E q since κ q 1 and the Lorentz boost are unity at our desired accuracy (noting the argument for neglecting the changes in phase-space due to our recoil scheme given in the previous subsection also holds for this resummation as the measurement function is unity and we are resumming logs up to α n s L 2n−1 accuracy). The limits on the z integrals capture angular ordering at NLL accuracy whilst still using a k ⊥ ordering variable. ∆ c (a, b) is a Sudakov factor ∆ c (a, b) = exp We can factorise Eq. (D.13) as φ qq (u, Q) = φ q (u, q ⊥ 1 )∆ q (q ⊥ 1 , Q) dz P qq (z)φ q (u, q ⊥ )φ g (u, q ⊥ ) × (q ↔q) + O(α 2 s ). (D. 16) keeping only terms first order in α s 22 . From this, we can identify φ q (u, Q) =φ q (u, q ⊥ 1 )∆ q (q ⊥ 1 , Q) dz P qq (z)φ q (u, q ⊥ )φ g (u, q ⊥ ). (D. 17) This expression is correct at LL accuracy with complete colour and only requires the coupling to run as α s (z(1 − z)q ⊥ ) in order to capture the full NLL (α n s L 2n−1 ) result. We also can note that the correct NLL resummation might not have been achieved using the local dipole prescription presented in Appendix C. This is because the recoil could introduce a correction in the n > 3 jet limit of the form φq(u, q ⊥ 1 ) φq(u, |q ⊥ 1 − q ⊥ 2 |) (the wavy arrow implying that it will approximately go to). This correction prevents both the usage of naive azimuthal averaging and the factorisation φ qq ≡ φ q φq (which naturally emerged between Eq. (D.13) and Eq. (D.16)), though it is possible that these features could reemerge once the phase space of each jet has been inclusively integrated over. Due to the other known NLL limitations of this recoil scheme, we did not think it worthwhile further proceeding to evaluate the order of these errors but rather conjecture that NLL errors will also be likely here.