Colour and logarithmic accuracy in final-state parton showers

Standard dipole parton showers are known to yield incorrect subleading-colour contributions to the leading (double) logarithmic terms for a variety of observables. In this work, concentrating on final-state showers, we present two simple, computationally efficient prescriptions to correct this problem, exploiting a Lund-diagram type classification of emission regions. We study the resulting effective multiple-emission matrix elements generated by the shower, and discuss their impact on subleading colour contributions to leading and next-to-leading logarithms (NLL) for a range of observables. In particular we show that the new schemes give the correct full colour NLL terms for global observables and multiplicities. Subleading colour issues remain at NLL (single logarithms) for non-global observables, though one of our two schemes reproduces the correct full-colour matrix-element for any number of energy-ordered commensurate-angle pairs of emissions. While we carry out our tests within the PanScales shower framework, the schemes are sufficiently simple that it should be straightforward to implement them also in other shower frameworks.


Introduction
Parton showers are ubiquitous tools in high-energy collider physics. In recent years, however, it has become clear that differences between parton showers are among the limiting systematics in many collider physics applications. This has motivated multiple efforts to better understand the consequences of the approximations contained within parton showers, and to exploit that understanding to guide their further development.
The majority of today's most commonly used showers, in particular those of the dipole family [1], make use of the idea of approximating QCD as if it had a large number of colours (N c ). Within this approximation one can view each event as a collection of independent colour dipoles: each gluon in an event functions as the colour-triplet end of one dipole and the colour anti-triplet end of another, while each (anti-)quark is the colour (anti-)triplet end of a single dipole. Those dipoles then radiate independently (and incoherently) from each other. This makes it relatively straightforward, at each stage of the showering, to generate a radiation pattern that is correct across small and large angles in the large-N c limit.
There are several ongoing efforts to include subleading-N c corrections, for example Refs. [2][3][4][5][6][7][8][9][10]. Including subleading colour corrections in full generality turns out to be computationally challenging, because as the parton multiplicity increases one should keep track of a rapidly-growing number of possible colour configurations, with contributions from higher-dimensional colour representations. Here we explore a complementary approach, one which connects the questions of subleading colour and subleading logarithmic contributions.
To help make things concrete, let us recall the PanScales criteria [11,12] for assessing the logarithmic accuracy of a shower 1. We should identify the kinematic configurations for which a shower correctly reproduces tree-level squared matrix elements. Typically it is useful to discuss this as a function of the separation between emissions in a Lund diagram [13]. We return to this in more detail in section 2.
2. We should evaluate the logarithmic accuracy of the shower's predictions for a range of common observables. Suppose we calculate some property P (α s , L) of an event, where α s is the strong coupling at a scale close to the hard scale, Q, of the event, and L, which we take negative throughout this paper, is the logarithm of a ratio of scales. For example this might be the cross section for events whose thrust is larger than 1 − e −|L| , or it might be the number of subjets found when clustering the event with resolution scale Q 2 e −2|L| . There are two ways of classifying logarithmic accuracy.
(1.1) LL stands for leading-logarithmic accuracy, NLL for next-to-leading logarithmic, and so forth. The N k LL functions, α k−1 s g k+1 (α s L), resum terms α n s L n+1−k and may in some cases involve operators rather than numbers. The LL function, α −1 s g 1 (α s L), starts off with a double logarithmic term α s L 2 . Certain observables, such as fragmentation functions and energy flow into a limited angular region, start only from the g 2 function, (in much of the literature, the g 2 function is then called LL; for consistency across the full set of observables, here we still call it NLL).
(b) For other observables, for example subjet multiplicities and certain other jet rates, there is no simple exponentiation of double logarithmic (DL) terms, and one may instead write P (α s , L) = P (α s , 0) h 1 (α s L 2 ) (1.2) where the N k DL function, i.e. α k/2 s h k+1 (α s L 2 ), resums terms α n s L 2n−k . In other work, this classification is often called N k LL (or occasionally N k LL Σ ). We adopt the N k DL nomenclature here to avoid confusion with the N k LL of Eq. (1.1).
For both the matrix element and observable-resummation logarithmic accuracy criteria, one may keep track of powers of the number of colours. Making the number of colours explicit, the leading-colour (LC) part of the LL function, LL-LC, involves terms α n s N n c L n+1 , while the next-to-leading colour (NLC) part, LL-NLC, involves terms α n s N n−2 c L n+1 and so forth. When needed, we will use the abbreviation FC to explicitly denote contributions that include the full colour structure.
Standard dipole showers correctly capture the full set of LL-LC terms (or DL-LC terms, as appropriate for the event property being measured). For exponentiating event properties, it is natural to consider values of the logarithm down to L ∼ −1/α s where NLL terms are of order 1 (cf. Eq. (1.1)). Keeping in mind that, numerically, α s ∼ 1/N 2 c ∼ 0.1, one then concludes that LL-NLC and NLL-LC terms are of comparable importance. 1 For observables that do not exponentiate, one instead considers values of the logarithm down to L ∼ −1/ √ α s , and with the same α s ∼ 1/N 2 c equivalence, one may take DL-NLC terms to be comparable to NNDL-LC terms.
Recently there has been significant progress in designing classes of showers that are NLL and NDL (LC) accurate aside from spin correlations, and in numerically demonstrating that accuracy in practice [12] (Refs. [15,16] instead examine an analytical approach). An approach that bears similarities to one of those shower classes was discussed in Ref. [17]. At this point, to consistently control all first subleading aspects beyond the LL-LC approximation, it becomes essential to identify approaches to construct showers that are correct 2 Angular ordering and Lund diagrams Let us start by elaborating on the first of our two criteria for logarithmic accuracy, i.e. the reproduction of matrix elements in suitably ordered limits. It is convenient to use Lund diagrams as a way of visualising the phase-space (see Ref. [29] for a concrete prescription to construct the Lund diagram from an event's kinematics). At LL accuracy, the tree-level matrix elements should be correct for any number of emissions that are well separated in a Lund diagram in both the logarithm of transverse momentum (k t ) and in rapidity (η = − ln tan θ/2), which one might call double strong ordering; at NLL accuracy, the treelevel matrix elements should be correct for any number of emissions that are well separated in at least one direction in the Lund diagram. Well-separated means that the distance D between points in the Lund diagram corresponding to any given pair of emissions should satisfy e −D 1. The correctness of the matrix element should hold no matter what that direction is, e.g. some may be well separated in rapidity but have similar ln k t , while others may be well separated in ln k t but have similar η values. In this article, the only respect in which we will relax this requirement on (leading-colour) logarithmic accuracy concerns the treatment of azimuthal correlations in collinear splittings, as induced by spin correlations, a topic that we defer to future work.
Throughout this section and the next ones, we will discuss how we attribute the correct colour factor for real emissions. The reader should keep in mind that virtual contributions are also being implicitly corrected at the same time, a consequence of the unitary nature of the showers that we consider in this paper. 3 For subleading colour effects at LL accuracy, one only needs to obtain the correct treelevel matrix element in regions where emissions are all well separated in rapidity. In this limit, for radiation at an angle θ, the question of colour reduces to that of examining the set of partons contained within a cone of aperture θ around the dipole end that is closer in angle. If that cone contains a single net quark (or anti-quark), i.e. |n q − nq| = 1, then the radiation is associated with a C F = (N 2 c −1)/(2N c ) colour factor, while if the cone contains zero net quarks, the radiation is associated with a C A = N c colour factor. 4 In the limit where all emissions are well separated in rapidity, these are the only two possible situations, and there will also never be any partons close to the edge of a given cone (because then the radiation would end up close in rapidity to an existing parton, keeping in mind that g → qq splittings are implicitly collinear up to and including NLL accuracy).
Although it is informative to refer to angular cones, it is more practical, going forwards, to use Lund diagrams. The angular-ordered colour-assignment prescription is illustrated with Lund diagrams in Fig. 1 for two e + e − →qq events, each dressed with additional 3 The discussion of the relation between real and virtual corrections is simple until one has four or more partons at commensurate angles; from that point onwards one should worry about amplitude-level evolution [30], which is beyond the accuracy and scope of this article. Additionally, when considering both initial and final-state emitters, non-trivial iπ terms enter at amplitude level, associated with Coulomb gluons, and these are a source of super-leading logarithms and coherence violation [31,32]. They have been addressed in the case of initial-final showers in Ref. [33], but are not relevant for the final-state showers considered here. 4 Recall that in dipole showers, the CA colour factor will be shared equally between two dipoles. Figure 1: Two example events and associated Lund diagrams. Points in the Lund diagram represent branchings, the triangles represent the two-dimensional logarithmic phase-space, while the coloured lines at the base of the diagram represent the phase-space for individual colour-dipoles (at the corresponding k t value), with dashed and solid segments indicating the parts of the dipole that should be associated respectively with a C F versus a C A /2 colour factor.
radiation. Recall that the main large triangle corresponds to the phase-space associated with the original e + e − →qq event (primary Lund plane), while each "leaf" that comes away from the main plane represents the additional phase-space that becomes available following emissions from that qq system (for example radiation collinear to the gluon in aqqg system). For concreteness we imagine a transverse momentum-ordered shower and consider the state of the event at a value of the ordering variable k t corresponding to the lower edge of the diagram (though our arguments apply to a range of shower ordering choices). The phase-space for emission at that k t is given by the base of the Lund diagram. Let us first consider Fig. 1a. In a normal leading-N c picture, this event consists of three colour dipoles:qg 2 , g 2 g 1 and g 1 q, represented as red, blue and green solid/dashed lines at the base of the Lund diagram. The dashed and solid styles indicate the colour factor based on angular ordering, which we work through in the rest of this paragraph. Along the dashed part of the (red)qg 2 dipole, i.e. the part on the primary Lund plane, any subsequent gluon emission is closer in angle to theq than to g 2 and a cone drawn around theq contains just theq, so the colour factor is C F . Along the solid part of the dipole, i.e. the part on the leaf associated with g 2 , the cone should be drawn around the gluon g 2 , since that is the dipole end that is closer in angle. The only particle contained within the cone around g 2 is the gluon g 2 itself, and one should use a colour factor of C A /2.
Next, we consider the (blue) g 2 g 1 dipole. In the solid blue region, along g 2 's leaf, the cone is to be drawn around g 2 and the only particle that is contained is a gluon, so we have a C A /2 colour factor; the situation is analogous for the solid blue region along g 1 's leaf. For the part of the dipole that is dashed, along the primary Lund plane, we need to separately examine the parts to the left and the right of the vertical dotted line. To the left, g 2 is the end of the dipole that is closer in angle. Writing the angle of the radiation, r, with respect to g 2 as θ rg 2 , and keeping in mind that we are in a region of the Lund plane were θ rg 2 θq g 2 , the cone of angle θ rg 2 that we draw around g 2 automatically containsq as well as g 2 , so the colour factor is C F . The situation is similar to the right of the vertical dotted line, but with a cone containing g 1 and q. Finally the (green) g 1 q dipole can be understood in analogy with the (red)qg 2 dipole.
Based on the above reasoning it is relatively straightforward to see that if one has an arbitrarily large number of gluon emissions (including secondary gluon emissions, but no g → qq splittings), then to determine the colour factor, one should identify whether an emission is on the primary Lund plane and if it is assign a C F colour factor, otherwise assign a C A /2 colour factor. That was the key observation made long ago by Gustafson [24]. A specific scheme (based on boosts) for achieving this was incorporated [25] as a modification of the Ariadne parton shower [13,34].
The next question one may ask is what happens if one allows for g →qq splittings. Whether one should consider this to be part of the LL terms is a matter that one may debate: on one hand g →qq splittings can exist in double-strongly ordered configurations. On the other hand if one integrates over their phase-space (as is relevant for thinking about logarithmically enhanced contributions to common observables), each g →qq splitting contributes only a single logarithm, i.e. an effect that is at most NLL in our counting. Here, we take the view that we should favour the (slightly) more ambitious goal and so aim to obtain the correct colour factors for double-strongly ordered configurations including g →qq splittings. Such a configuration is shown in Fig. 1b and observing the (blue) g 3 q dipole, one sees that it consists of a sequence of solid (C A /2) and dashed (C F ) segments, with the final dashed segment along the q leaf, which is not part of the primary Lund plane. In general the number of segments on a dipole can be anywhere between 1 (e.g. the solid green g 2 g 3 dipole) and infinity, though in practice the average number of segments per dipole turns out to be of order 1 even in high-multiplicity events, a consequence of the smaller number of quarks than gluons produced in the parton shower. This picture differs from the standard approach within dipole showers of breaking every dipole into two parts, each associated with one of the two ends. Note that in this paper we will apply our new segmentation just for the purposes of colour assignment, retaining the default two-part segmentation for the kinematic maps of each shower that we examine. 5

A solution with segments and transition points
Based on the reasoning in section 2, here we propose the first of our concrete schemes for achieving LL full-colour accuracy. We start by introducing the key ideas with the help of a worked example, section 3.1. We then discuss choices we can make that affect aspects beyond LL-FC accuracy, section 3.2, or that arise in occasional special cases, section 3.3 5 One could also imagine constructing a shower with global recoil whose kinematic map follows the pattern of the colour map. Indeed, we wonder whether this might be the most physical approach to constructing a global recoil map, avoiding certain ugly features of existing global maps [12,17,35]. However we leave the investigation of this question to future work. Figure 2: Examples of colour segments and transition points for representative configurations: (a) a q → qg splitting (i.e. the emission of a gluon from a quark segment), (b) a g → gg splitting (i.e. a gluon emission from a gluon segment), and (c) a g → qq splitting.
(some readers may prefer to skip these parts). In section 3.4 we give our full algorithm. Finally in section 3.5 we show how it can be adapted to the Pythia 8 shower [36,37].

A worked example
The key insight from section 2 is that to obtain LL full colour accuracy it is enough to break every dipole into a suitable sequence of C F and C A colour segments. We label those segments by their extremities in rapidity and the colour factor along the segment. For an initial e + e − →qq event we start with a single dipole consisting of one segment stretching from a rapidity of η = −∞ (q end) to η = +∞ (q end), associated with a C F colour factor. We denote this as [−∞, C F , ∞]q q , (3.1) where our notation consists of a sequence of segment boundaries (or, equivalently, transition points) and segment colour factors. When a dipole emits a gluon g 1 , we define the gluon's rapidity within the dipole, η g 1 = ±| ln tan θ/2|, in terms of its angle θ with respect to the dipole end that is closer in the event frame. 6 We assign a positive (negative) rapidity if it is closer to the triplet (anti-triplet) end. Since there is only a C F segment in theqq dipole, the gluon is necessarily emitted with a C F colour factor. When we radiate a gluon g 1 from theqq dipole, we end up with two dipoles,qg 1 and g 1 q, each of which now has one C F and one C A region, cf. Fig. 2a. The segmentation for the two dipoles is as follows where we have highlighted in blue the extra segments that are a consequence of the gluon emission and, for now, we take η L g 1 ≡ η R g 1 ≡ η g 1 . We see here that both new dipoles have a transition point at the rapidity of the emission itself (see Fig. 2a). Note that since gluons are shared across two dipoles, the second segment should really have a colour factor C A /2. For compactness, we suppress the explicit factor 1/2 in our notation.
Next, we consider a second emission, g 2 . It is sufficient to examine the situation where it is emitted from theqg 1 dipole. We use the same procedure to evaluate the rapidity η g 2 of g 2 as given above, but now with respect to theqg 1 dipole, and we examine where η g 2 lies in theqg 1 sequence of segments. Since there are two segments, there are two possible cases: (a) −∞ < η g 2 < η L g 1 , where radiation occurs with a C F colour factor and (b) η L g 1 < η g 2 < ∞ where it occurs with a C A /2 colour factor. As well as differing in the colour factor for emission, the two cases also differ in terms of the resulting set of new segments. Case (a), emission from the C F region, gives the following segmentation of the resulting three dipoles 3) i.e. one splits theqg 1 sequence of segments in Eq. (3.2) into two separate sequences: one is for theqg 2 dipole, which keeps everything to the left of the C F segment where the radiation occurred, plus a closing η L g 2 , C A , ∞] right segment; and one for the g 2 g 1 dipole, which keeps everything to the right of the C F segment, plus a closing [−∞, C A , η R g 2 left segment. This pattern corresponds exactly to what we see in Fig. 1a. Case (b), insertion into the C A region, is shown in Fig. 2b and gives Again we have highlighted the new part in blue. Since we are inserting a gluon into a C A region there are no additional C F −C A transition points. The final case that we need to consider is the splitting of a gluon to q q , cf. Fig. 2c. In a strong ordering (i.e. collinear) limit, as relevant for accuracies up to and including NLL, the |η| associated with the g → q q splitting will always be larger than the last non-infinite η transition point in each dipole's sequence, i.e. the splitting will always be associated with a C A segment that is at the extreme left or right ends of a sequence. We consider aqg 1 q event and the branching of g 1 → q q in theqg 1 dipole, i.e. we have Eq. (3.2) as our starting point. Defining where θ q q is the opening angle of the q q pair in the event frame, we then obtain the following segments i.e. theqg 1 and g 1 q dipoles have respectively becomeqq andq q dipoles and each of those dipoles has additional transition points to a C F colour factor at η L q q = |η q q | and η R q q = −|η q q | respectively. The opposite signs for η L q q and η R q q arise because the transition needs to be at angles close to the q for theqq dipole, i.e. positive rapidity, and at angles close to theq end of theq q dipole, i.e. negative rapidity.
Before formulating a full, concrete algorithm, some details need to be specified, related to O (1) choices for the rapidity boundaries (section 3.2), and the handling of cases where cut θ cut θ cut θ 1 q q g Figure 3: Illustration of the configuration that we consider to examine the equivalence between exact angular ordering and our partitioning of dipoles, as relevant for NDL accuracy in particle and jet multiplicities. We will consider the emission of a further soft gluon, g 2 , with respect to this configuration, at an angle of at least θ cut away from theq, q and g 1 . The dashed/solid coloured lines indicate schematically where we use a C F v. C A /2 colour factor for each of the two dipoles (blue:qg 1 , red g 1 q).
two emissions are close in rapidity, with resulting ambiguities in how to construct the segments (section 3.3).

Specific rapidity definitions and NDL accuracy
There are two aspects to examine concerning the exact choices of rapidity transition points, both relevant for configurations where two branchings occur at similar angles. One involves identifying a choice that can provide full-colour NDL (NDL-FC) accuracy for basic quantities like jet and particle multiplicities, i.e. control of terms α n s L 2n−1 . The control of this class of terms at NDL-FC accuracy has long been one of the strong arguments in favour of angular ordered showers [38]. The other question is more practical: in the segment algorithm we implicitly assume that all transition points are ordered, with a consistent alternating set of C F and C A segments. Those properties are trivially maintained when all branchings are at disparate angles, but that is no longer necessarily the case when two or more branchings are at commensurate angles.
We first consider what is required for NDL-FC accuracy. We start with a systemqg 1 q, with θ g 1 q 1, and examine the emission of a second much softer gluon, g 2 , in a limited range of energy, ∆ln E, and with the following angular constraints: (3.7) The angular limits for the emission of g 2 are represented schematically in Fig. 3 and we work in a limit where θ cut θ g 1 q . The integrated full-colour rate of g 2 emission in this region is well known to be [39,40] where η g 1 = − ln tan θ g 1 q /2, η cut = − ln tan θ cut /2. The parton shower's ability to reproduce this result is a key requirement for NDL full-colour accuracy. In particular, if the evaluation of this integrated emission rate for the shower has an extra O (1) constant in the square bracket, NDL accuracy will not be achieved.
To provide a concrete demonstration of how to satisfy Eq. (3.8), we consider the PanScales showers [12]. The kinematics for branching in these showers is parameterised in terms of the shower ordering variable ln v, a longitudinal variableη (linearly related to the logarithm of a light-cone momentum component) and an azimuthal angle φ. The emission density (working at fixed coupling for illustrative purposes) is where C is the colour factor (C F or C A with C F = C A /2 in the large-N c limit). 7 To integrate the shower over the region of Eq. (3.7), one needs to relateη to the actual rapidity of a soft and collinear emission from a dipole ij with respect to the closer of i and j. The following approximate formula (specific to the PanScales showers) provides the (signed) rapidity with respect to the closer of i and j (when the dipole for which we evaluate η approx is obvious, we omit the (ij) superscript). It is approximate in the following sense: for an emission k from an ij dipole, the true rapidity with respect to the closer of i and j has additional O (1) corrections when the smaller of the θ ik and θ jk is commensurate with θ ij . For a soft emission, those corrections vanish when the ratio of min(θ ik , θ jk ) to θ ij is small. This is the case notably close to the θ cut limits in Eq. (3.7). For the configuration in Fig. 3, with θ g 1 q 1 (θq g 1 π), we havē Note that η approx depends linearly onη except atη = 0, where it is discontinuous for the g 1 q (i.e. small-angle) dipole. With these results one can translate the constraints of 7 Recall, e.g., that the PanLocal family of mappings [12] for emission of p k from a dipole { pi, pj} is where k ⊥ = kt [n ⊥,1 cos φ + n ⊥,2 sin φ], with n 2 ⊥,m = −1, n ⊥,m · p i/j = 0 (m = 1, 2), n ⊥,1 · n ⊥,2 = 0 and Here the parameter β sets the choice of ordering variable, sĩ = 2 pi · pj, sĩ = 2 pi · Q, and Q is the total event momentum. The light-cone components of p k are given by a k = s sĩsĩ kte +η and b k = sĩ sĩs kte −η . The quantity f in Eq. (3.10) determines how transverse recoil is shared between pi and pj, cf. below. The ai, bi, aj, bj are fully specified by the requirements p 2 i/j = 0, (pi + pj + p k ) = ( pi + pj) and pi =pi for kt → 0. These and the shower-dependent choices for f are given in the supplemental material to Ref. [12]. The PanGlobal shower drops the −f k ⊥ and −(1 − f )k ⊥ terms in Eq. (3.10), sets bi = aj = 0, and boosts and rescales the event after each emission so as to conserve momentum.
Eq. (3.7) into constraints onη for each of the two dipoles, For each of the dipoles, the procedure of section 3.1 splits the dipole's rapidity range into a C F piece and a C A /2 piece. The choice that we make for the transition points isη = η g 1 for theqg 1 dipole (using C F forη < η g 1 ), andη = 0 for the g 1 q dipole (using C F forη > 0). Together with the limits in Eq. (3.14), one then finds that the C A /2 pieces from the two dipoles add up to give a totalη range of 2(η cut − η g 1 ), while the C F pieces add up to give 2η cut , giving a total result for the parton shower of , (3.15) in agreement with the full-colour result of Eq. (3.8).
While the discussion in terms ofη is useful for locating an emission with respect to a dipole with a single transition point, to handle more complicated events, it becomes simpler to reason directly in terms of η approx . Thus the transition point for theqg 1 dipole is at η approx = η g 1 . For the g 1 q dipole, since η approx is discontinuous atη = 0, any transition between −η g 1 and +η g 1 is equivalent. One simple convention is to set the transition for the g 1 q dipole to be at η approx = 0, i.e. along the dipole's angular bisector in the event frame.
Together, these observations bring us to the following recipe for soft emission from a configuration with an existing collinear branching. Working in the event centre-of-mass frame, when inserting a gluon emission i into a C F segment, determine its angle θ i to the triplet (anti-triplet) end of the dipole if η approx,i > 0 (< 0) and define η i = ∓ ln tan θ i /2, using a negative (positive) sign when the θ i is determined with respect to the triplet (antitriplet) end. In constructing new transition points, e.g. as in Eqs. (3.2) or (3.3), use When considering the generation of a new emission k, determine its η approx,k from Eq. (3.12) and compare that η approx,k to the dipole's existing transition points to determine the segment to which the emission belongs and the associated C F v. C A /2 colour factor. Similar considerations about the total integral for soft gluon emission are relevant for a g →qq splitting. One can demonstrate that the prescription given in Eqs. (3.5) and (3.6) already ensures the correct integral, an analogue of Eq. (3.8).

Special cases
When inserting a gluon g into a C F segment, and that gluon's angle with respect to the dipole end is close to an existing transition point, it can sometimes happen that the insertion η L g or η R g lies outside the segment range (recall that η approx and η exact can differ by an order 1 amount when the emission angle is commensurate with the dipole opening angle). From a DL logarithmic point of view, e.g. for particle or jet multiplicities, such a situation is beyond NDL. Specifically, to trigger it, one must have two energy-ordered gluons at similar angles. In such a situation, any subleading-N c issue induced by misordering of transition points only affects emission of a subsequent, third, soft gluon that is once again at a similar angle (it also affects related virtual corrections). The requirement of three energy-ordered gluons at commensurate angles corresponds to NNDL terms for quantities such as multiplicities, because one loses a logarithm relative to DL for each of the second and third gluons when requiring its angle to be similar to the first.
Still, we need a definite prescription to handle such cases. We choose to give priority to the requirement that the transition points in the colour-segment sequence should always remain ordered. Consider the configuration [. . . , C A , η i−1 , C F , η i , C A . . .] ab and a situation where we have located a new emission, according to its η approx , as belonging to the central C F segment. We would normally expect to split the sequence as we discard the last two segments of the new ag sequence (η i−1 , C F , η L g and η L g , C A , ∞) and extend the remaining rightmost C A segment to +∞, We use an analogous procedure when η R g ≥ η i , dropping the two first segments of the gb sequence and extending its leftmost C A segment to −∞. Since η L g ≥ η R g , and since the transition points are always ordered prior to insertion, this adaptation is needed at most on one side, never both.
Similar considerations apply to g → qq splittings, though the impact on accuracy is now only expected to be N 3 DL, because g → qq splittings start one logarithm down. Where normally we would have if η L q q ≤ η i , then we remove the last two segments of the new aq sequence and extend its remaining rightmost (C F ) segment to infinity, There is an analogous adaptation when η R q q ≥ η j . In practice, for physical α s values, these adaptations are used in at most a few percent of gluon insertions in C F segments and for up to 10−20% of g → qq splittings.

The algorithm
With the help of the above reasoning we are now ready to formulate a full segment-based algorithm for the PanScales showers. The event should start with Eq. (3.1) for an initial e + e − →qq event and analogues with C A for the two dipoles of an initial H → gg event. Splitting functions and their integrals for Sudakov factors should by default all be evaluated with the replacement C F → C A /2 = N c /2. Dipoles are always labelled by their anti-triplet end followed by their triplet end.
For emission of a gluon g from an ab dipole, use the parent dipole momenta and PanScales kinematic generation variableη to determine η (ab) approx for the gluon, according to Eq. (3.12). Identify the location of η (ab) approx within the ab dipole's sequence of segments. If it is in a C F segment, reject the emission with probability (1 − 2C F /C A ) (or, alternatively, multiply the event weight by (1 − 2C F /C A )). If it is rejected, continue by generating the next value for the shower ordering variable, starting from the value that was been rejected. If the emission is accepted, and η (ab) approx > 0 (< 0) evaluate η g = − ln tan θ gb /2 (ln tan θ ag /2). Determine η L g and η R g from η g according to Eq. (3.16), and replace the ab dipole's sequence of segments with new sequences for the ag and gb dipoles as follows, The left and right ". . ." represent the full sequences to the left and right ends of the original C F factor. They are assigned respectively to the left and right ends of the ag and gb dipoles' sequences. If η L g is such that its insertion would create a disordered sequence of transition rapidities, remove the last two segments in the ag dipole, i.e. use Eq. (3.17b). Proceed analogously for η R g . If the emitted gluon's η approx,g is in a C A segment, do not apply any additional rejection factor (i.e. generate the gluon with its original C A /2 colour factor) and replace the ab dipole sequence of segments as follows.
where the C A on the left corresponds to the segment associated with η approx,g . Finally we consider g → q q splitting, where the gluon belongs to two dipoles ag and gb. We assume the splitting to have been generated with a normal P g→qq splitting function, so there is no rejection factor to apply. Using the opening angle θ q q between the new quark and anti-quark, determine η L q q = −η R q q = | ln tan θ q q /2|, i.e. Eq. (3.5). The g → q q splitting then acts as follows on the segments: In a situation where the insertion of η L q q would lead to a disordered sequence of transition points, remove the last two segments from the aq dipole, as in Eq. (3.18b), and analogously remove the first two segments from theq b dipole if η R q q is disordered with respect to the other transition points.
3.5 Use in other showers, e.g. Pythia 8 The technique outlined above can be applied to almost any dipole or antenna shower. The one adaptation that is required is to identify a suitable expression for η approx , i.e. the analogue of Eq. (3.12), for that shower. For example in the Pythia 8 shower [36,37], for a dipole ab where the emitter has momentum p b and the spectator p a , the expression for η approx is 22) where p ⊥evol is the shower evolution variable (a transverse momentum), the event momentum is Q, and 1 − z is the fraction of p b 's energy carried away by the gluon in the dipole centre-of-mass frame. One can verify that when an emission is soft and collinear to one of p b or p a , this gives the correct signed rapidity with respect to the closer of b or a in the event centre-of-mass frame.

A solution with nested ordered double-soft (NODS) corrections
The approach of section 3 had two elements: one was the determination of an effective colour factor for each new emission; the other was to establish whether a new emission should be attributed to a C F segment from the point of view of the colour-factor identification for subsequent emissions.
In this section we consider an approach that retains the second of these elements, but replaces the binary colour-factor choice (C F v. C A /2) with a local matrix-element correction that reproduces the full-N c radiation pattern for configurations involving a pair of energy-ordered soft gluons that are close in rapidity in the Lund diagram, but with all other emissions well separated in rapidity from them and from each other. The procedure will actually give the correct full-N c matrix element even when there are multiple such pairs around, as long as each is well separated in rapidity from all others (in the same sense as in our discussion at the beginning of section 2).
As in the previous section, we will start by explaining our general approach in the case of a simple subset of event structures in section 4.1. Next, in section 4.2, we will consider issues that arise for more general event structures. We will then give our complete algorithm in section 4.3.

Angular-ordered primary-only events
To understand the NODS procedure, consider a situation with a singleqq pair and gluons g 1 . . . g n , each of which is primary in the Lund-diagram sense and well separated in rapidity from all the other gluons, with an ordering in Lund-diagram primary rapidity of η 1 η 2 . . . η n . This configuration will dominantly be associated with a leading-colour dipole structureq1, 12, 23, . . ., nq (up to corrections suppressed by powers of e −|η i −η j | ), which we can represent asq (4.1) using n = 4 for concreteness. The corresponding leading-colour squared matrix element for emission of a soft gluon with momentum k is 8 where we have introduced the shorthand and assigned a C A /2 colour factor to each dipole. In the specific limit that we are considering, i.e. primary emissions that are all well separated in rapidity, the full-colour matrix element reduces to One simple approach to reproducing Eq. (4.4) would be to introduce an acceptance probability for each emission of We will introduce a general shorthand for such expressions .
One approach along these lines was proposed in Ref. [41]. A downside of any approach that uses the full set of momenta in the acceptance is that for an N -particle event it would involve evaluating O (N ) dot products for each emission, leading to a contribution to the total showering time that scales as N 2 . This is not necessarily prohibitive; for example the implementation of the Pythia 8 showering algorithm scales as N 2 , and showers with global recoil such as PanGlobal also scale as N 2 in their current implementation. However it turns out to be possible to formulate an expression for p accept that maintains the same accuracy but can be evaluated in O (1) time.
To understand how this can be done, we consider a dipole in the middle of the chain, say i, i + 1, where both i and i + 1 are gluons. Specifically for this dipole we are free to use rather than the full p accept (q, 1, 2, . . . , n, q). In |M LC | 2 we are justified in dropping all dipoles q1, . . . , (i − 2, i − 1), and (i + 2, i + 3), . . . , nq, because of the requirement that all emissions are well separated in rapidity. To see this, imagine that i is at negative rapidity (θq i 1), i + 1 is at positive rapidity (θ i+1,q 1). Consider an emitted gluon k such that θ ik θq i . Then, because of our requirement that all existing emissions should be well separated in rapidity, we have θ i−2,i−1 θ i−2,k θ i−1,k . Examining Eq. (4.3) one can then see that the (i − 2, i − 1) term is negligible compared to the terms included in |M LC,i,i+1 | 2 , which is simply a consequence of the fact that a small-angle dipole does not substantially emit at large angles. Similarly for the (i + 2, i + 3) term, as well as all the other terms that we have neglected in Eq. (4.7a). Now we turn to Eq. (4.7b): here we have replaced aqq with i − 1, i + 2, which is justified since θ i−1,k − θq ,k θq k and θ k,i+2 − θ k,q θ kq for any momentum k that is likely to be radiated by the (i, i + 1) dipole.
Note that with the truncation adopted in Eq. (4.7a), it would not be sensible to use the (qq) emission factor in the (C F − C A /2) term of Eq. (4.7b), because such a term would have a negative divergence for k exactly collinear toq or q that would not be compensated for by any of the terms in the |M LC,i,i+1 | 2 contribution. In contrast, one can show that the acceptance probability as written Eq. (4.7) is bounded to be in the range i.e. for the physical value of N c it is always positive definite and so straightforward to use in event generation. 9 We have given the justification for writing Eq. (4.7) in a specific frame, one in which it is manifest that one can replaceq with i − 1 and q with i + 2. However the underlying expressions are Lorentz invariant, and the validity of those replacements is ultimately ensured by the condition that all existing emissions are on the primary Lund plane and separated by large differences in rapidity.
We also need to consider dipoles at the end of the chain, for example theq1 dipole in Eq. (4.1). In such a case, we simply write 9) and analogously at the other extremity.
In what follows, when we write it is useful to introduce the following terminology: i − 1 is the auxiliary momentum at the colour-anti-triplet (i) end of the dipole, and i+2 is the auxiliary at the colour-triplet (i+1) end of the dipole. We have emphasised that i−1 and i+2 are auxiliaries by separating their 9 Without loss of generality, for massless particles, one can always rotate and boost an event such that i − 1 and i + 2 are along the negative and positive z axes and k is along the −x direction. The configuration that gives the minimum is then one where i and i + 1 have their momenta along the (sin θ, 0, ± cos θ) directions with tan θ 2 = 1 2 .
indices from the others with semi-colons. One can then represent the event of Eq. (4.1) as where the dotted line beneath each solid black dipole represents the −1/(2N c ) dipole that we use in the matrix-element correction. It shows, for example, that the 23 dipole has auxiliaries 1 and 4, leading to an acceptance factor p accept (1; 2, 3; 4).

Considerations for general events
Most events do not consist just of angular-ordered primary emissions discussed so far: primary gluons can themselves emit secondary gluons, and any gluon may split to qq.
The first case that we consider is emission from a region that would be associated with a C F colour factor in the segment approach, e.g. the emission of 5 from the 23 dipole, (4.12) When 5 is well-separated in rapidity from both 2 and 3, the acceptance factor for emission of 5 will reduce to 2C F /C A . The new 25 dipole retains the left-hand auxiliary (1) of the parent dipole and acquires a new right-hand auxiliary (3), corresponding to the triplet end of the parent dipole. Similarly, the new 53 dipole acquires a new left-hand auxiliary (2), corresponding to the anti-triplet end of the parent dipole, and retains the parent's right-hand auxiliary (4).
Next, we consider a collinear g → gg splitting. Imagine that we have the configuration in Eq. (4.11) and that dipole 23 branches, with its 2 end splitting to gluons 5 and 6 in a collinear configuration θ 56 θ 15 θ 16 , illustrated as follows (4.13) where the new 56 dipole is highlighted in red. The 15 dipole, which is the successor of the 12 dipole retains the 12 dipole's auxiliaries. Similarly the 63 dipole retains the 23 dipole's auxiliaries. The new 56 dipole, since it is produced far in Lund-plane rapidity from any region with a 2C F /C A correction, does not need any auxiliaries, in the same way that it has no transitions points (and so no C F segments) in the approach of section 3.
For configurations that are not strongly ordered, e.g. emission of 5 from the 23 dipole at an angle that is commensurate with θ 21 , the acceptance factor is unambiguous (and correct) in our approach. However there is an ambiguity in the assignment of auxiliaries to the child dipoles, which translates to an ambiguity in the acceptance for a subsequent third emission at similar angles, a configuration where we do not aim for full-N c accuracy. We resolve that ambiguity as follows: as well as retaining information on auxiliaries, we also retain information on segments and their transitions, as in the approach of section 3. If an emission occurs in a C F segment then we assign new auxiliaries as in Eq. (4.12). Instead, if it occurs in a C A /2 segment, we do not introduce any new auxiliaries, as in Eq. (4.13).
Finally, we consider a g → q q splitting (for which we always assign a normal P g→qq splitting function and use a colour acceptance factor of 1). Consider 2 → q q , (4.14) where we have highlighted the new q q pair in red. To justify the choice for the new auxiliary variables, we consider the radiation of a soft gluon after the 2 → q q splitting. In the segmented approach of section 3 we would have said that that theq 3 dipole now has two C F segments. In our NODS approach, we introduce an acceptance factor for each C F -like segment in the dipole and determine the overall acceptance for the dipole as the product of the individual acceptance factors. E.g. for theq 3 dipole we write where the first C F -like segment starting from theq is associated with a single auxiliary (q ), at its right-hand end, 10 while the second C F segment retains the auxiliaries of the parent 23 dipole. For strongly angular-ordered configurations, i.e. θ q q θ q 1 , only one of the two p accept (. . .) factors in Eq. (4.15) will ever differ substantially from 1. For configurations where strong angular ordering does not hold, for which we do not aim to achieve full-N c accuracy, there will be phase-space regions for a subsequent emission k where both factors will be below 1.

Full algorithm
We are now ready to specify our full matrix-element-based NODS colour-handling algorithm. It retains the core framework of the segmented approach, specified in section 3.4, with the following augmentations: 10 The correctness of the p accept (q , 3; q ) factor can be seen using arguments similar to those of section 4.1.
In the event centre-of-mass frame, angular ordering implies θ q q 1, i.e. we only need to consider collinear g → q q splitting. We imagine boosting the event to a frame in which θ q q is of order 1. In the angularordered limit, all other particles will then be constricted along a single collinear direction, which is well separated from the q andq directions. Next, we temporarily imagine that a CA/2 dipole stretches between theq and q, so as to avoid complications with their actual (quark-like) colour structure. Then, the full acceptance factor correction associated with the q q CF colour structure would be p accept (q , 3, 4, q,q, 1, q ). For emissions along theq end of theq 3 dipole, all terms (34), (4q), . . . , (q1), are irrelevant. Furthermore, one can replace (1q ) with (3q ), since 1 and 3 are collinear and the (1q ) term is relevant only when the emission is emission is far in angle from 1. Thus the acceptance factor can be reduced to p accept (q , 3; q ).

1.
Where the extremity of a segment has finite rapidity, it is associated with an auxiliary, which we denoteā towards the anti-triplet end of the segment and a towards the triplet end.
2. Emission acceptance: (a) The acceptance for gluon emission from a dipole ij is the product of individual acceptance factors for each of the dipole's C F segments. The individual acceptance factor for a given C F segment is (b) For a g → q q splitting, the colour acceptance factor is set to 1.

Auxiliaries update:
(a) If a gluon emission k from an ij dipole occurs in a C F segment with auxiliariesā and a, then the corresponding C F segment in the child ik dipole has auxiliaries a and j, while that in the child kj dipole has auxiliaries i and a. (Where an auxiliary is absent in the parent dipole because the segment stretches to infinite rapidity, it remains absent in the child dipole).
(b) If a gluon emission from an ij dipole occurs in a C A segment, the child dipoles' C F segments retain the auxiliaries of the parent dipole.
(c) For g → q q splitting, the new C F segment in the anti-triplet-q dipole acquires theq as its anti-triplet (left) end auxiliary, while the new segment in theqtriplet dipole acquires q as its triplet (right-hand) end auxiliary.
(d) In those special cases of the algorithm of section 3.4 where two segments are removed from the extremity of a sequence, the corresponding auxiliaries are also to be removed.
When storing the auxiliary for a segment, there is some freedom: for example one can choose the auxiliary's momentum at the time that it is associated with the segment, or its (possibly different) momentum at a later stage in the event when one is evaluating Eq. (4.16). In practice we make a third choice: when gluon emission causes a left-hand auxiliary to be first added to a segment-sequence of an ij dipole, we store the difference in direction between the auxiliary and the anti-triplet end of the dipole, δ iā = dā − d i , where d i is the 3-vector direction of i. For a right-hand auxiliary point, we store the difference in direction between the auxiliary and the triplet end of the dipole δ ja = d a − d j . When we later come to need an auxiliary momentum to evaluate Eq. (4.16) in an m dipole that descends from the original ij, we reconstruct auxiliary directions using differences with respect to the and m directions at that stage of the event, For a g → q q splitting, in the C F segment on the anti-triplet-q dipole we store its lefthand auxiliary (q ) direction difference relative to the q (i.e. triplet) direction, δ q q so that when that segment is eventually evaluated in the context of emission from a descendent m dipole we use and analogously for theq -triplet dipole. We are not aware of a reason for preferring one of these approaches over any other, however that based on direction differences has the advantage of being easiest to use for certain of our tests in section 7 and so it is our main choice.

Colour-factor from emitter (CFFE) algorithms
When examining the numerical behaviour of our algorithms in the next sections, it will be useful to have an illustration of the behaviour of existing algorithms for assigning colour, in order to gauge the practical impact of our proposals.
In standard dipole showers (e.g. the Pythia 8 shower [36] or the Dire v1 shower [42]), each dipole is split into two elements, according to which end of the dipole is deemed the "emitter". Roughly speaking, each element accounts for the rapidity phase-space that extends from zero rapidity in the dipole centre-of-mass frame to the extremity of rapidity phase-space at the emitter end. For gluon emission, each element acquires the full-colour splitting function associated with the flavour of the emitter, i.e. yielding P gq (z) = 2C F /z (P gg (z) = C A /z) for emission of a soft gluon with momentum fraction z 1 if the emitter is a quark (gluon). Accordingly, we refer to this approach as the "colour-factor from emitter" (CFFE) approach. Ref. [11] showed that the CFFE approach yields wrong DL-NLC terms starting from second order, α 2 s L 4 , for some observables, e.g. the thrust. One can also examine the CFFE approach in PanScales-like showers, where the phasespace is again divided between the two elements of the dipole, at an angle that is roughly equidistant between the emitter and spectator ends of the dipole, in the event centre-ofmass frame. 11 Using the approach of Ref. [11], one can show that the method still yields incorrect DL-NLC terms starting from second order, α 2 s L 4 , for some observables. One could fix the second-order issue by dividing the dipole in a way that mimics the transition points that we have discussed in section 3. However, at third order, it becomes impossible to obtain the correct answer for general double-logarithmic observables with a single transition point for the dipole. One can see this by examining Fig. 1a, where the (blue) g 2 g 1 dipole requires two transition points, one from C A /2 to C F and another back to C A /2.

Matrix element tests
In order to demonstrate how the methods presented above effectively incorporate subleadingcolour effects, we examine parton shower results at fixed order for 2 → 3 + g and 2 → 4 + g final states (e + e − → X + g), where X ∈ {qg 1 q,qq q q,qg 1 g 2 q} and g is an additional soft gluon that is radiated from system X. We compare the parton-shower (PS) result dσ PS /dηdψ (σ PS for brevity) to the exact ratio of full-colour squared matrix elements dσ FC /dηdψ = dΦ g |M X+g | 2 /|M X | 2 . Here, η and ψ are the rapidity and azimuthal variables as obtained from a Lund declustering procedure (see below). We do this for each of the segment, NODS and CFFE parton-shower colour schemes.
In sections 6.1 and 6.2, we show differential results for soft-collinear configurations, where all partons in X are strongly ordered in energy and rapidity, and the additional soft gluon g can be at angles commensurate with any other particle (while still strongly ordered in energy compared to them). To test NDL-FC accuracy, in section 6.3 we compare the integrated soft-gluon emission rate I PS , determined numerically from dσ PS /dηdψ, to the expected value I FC , see Eq. (3.8). We do this for three different configurations in which the gluon g 1 is either soft and collinear, hard and collinear, or soft and large-angle with respect to the initial qq dipole.
Note that the approach that we develop below for matrix-element tests could also be adapted to provide interesting measurements within Lund-diagram type jet substructure analyses.
6.1 Differential matrix element: 2 → 3 + g configuration All differential plots shown below are produced with the PanGlobal shower, with β = 0 (see [12] and Eq. (3.11)). To set up the event (e + e − →qg 1 q, see Figs. 2a and 3), we let the parton shower split the initialqq dipole, which emits a gluon g 1 with predefined kinematics, corresponding in our soft-collinear case to z g 1 = 10 −8 , η g 1 = 5, ψ g 1 = π , (6.1) in the event frame. The emission of the additional gluon, g, is then also performed by the parton shower: we fix the value of the evolution variable ln v (≡ ln k t for the PanGlobal β = 0 shower), while sampling over the two remaining shower degrees of freedomη and φ. We cluster the event with the Cambridge/Aachen algorithm [43,44], and work backwards through the clustering history to determine the effective Lund diagram [29] (we call this Lund declustering). We then log the event weight in a histogram in variables (η, ψ) defined with respect to the emission's parent Lund leaf. We fill separate histograms for emission on the primary and secondary Lund leaves. Starting from the same initialqg 1 q configuration, we also sample the exact full-colour analytic tree-level matrix element, see Eq. (4.4), neglecting any recoil from the last emission (since the impact of the emitted soft gluon g is negligible). Finally, we take the ratio of the parton shower and exact full-colour histograms.
The results are collected in Fig. 4 for the CFFE (section 5), segment (section 3) and NODS (section 4) colour schemes. The labels on each panel, q and g 1 , refer to the par-  Figure 4: Tree-level fixed-order expansion of the parton shower density for emission of a soft gluon from aqg 1 q system (as represented in Fig. 2a). The results are shown for the CFFE (top), segment (middle) and NODS (bottom) colour schemes, as implemented for the PanGlobal shower with β = 0. The first row in each pair of plots depicts the ratio σ PS /σ LC (green-yellow colour bar). The second row shows the relative deviation (σ PS − σ FC )/σ FC between the parton shower and the analytic matrix element (blue-red colour bar). The left-hand panels corresponds to the primary-emission phase-space region, while the right-hand panels corresponds to emission from g 1 , as emphasised in the diagrams at the top.
ton with which the gluon g is clustered by the C/A algorithm (i.e. the emitter according to the Lund declustering sequence). At the top of each column, a diagram indicates the corresponding part of the event that the Lund rapidity η refers to (primary Lund plane, along the quark direction, or secondary plane along the gluon direction). Note that in the following, we do not show differential results at theq-end since, in this collinear configuration, the emission correctly gets an overall factor of C F for all colour schemes that we consider (to within corrections suppressed by powers of e −|η−ηg 1 | ). To help see the effective colour factor applied by the parton shower, for each method the upper row in Fig. 4 depicts ratios, σ PS /σ LC , of the various parton-shower colour schemes to the leading-colour result σ LC . For the latter, we use a constant colour factor C F = C A /2 = 3/2 irrespective of the position of the emission. The ratio σ PS /σ LC is 1 (yellow in Fig. 4) when the effective colour factor is C A /2, while it is 8/9 = 2C F /C A (green) when the effective colour factor is C F . In the lower row, we plot the relative deviation (σ PS − σ FC )/σ FC between the parton-shower weight and the analytic full-colour matrix element. Regions where the shower agrees with the full-colour matrix element come out in white (to within statistical fluctuations).
The two collinear limits of interest, g q and g g 1 , are always mapped to the large-η region of the corresponding panels, i.e. the q-leaf (left column) and g 1 -leaf (right column). Rapidities cover the range η p < η < ∞, where η p corresponds to the opening angle between the parent parton and the next particle in the declustering sequence (with η p = 0 for radiation on the primary plane). Holes correspond to the point where the emission moves across leaves, from a primary to a secondary plane (from the q-leaf to the g 1 -leaf, in our example).
Let us examine the three colour methods in turn: • The CFFE method incorrectly assigns a colour factor C A /2 to the region 0 < η < η g 1 = 5. The fact that it extends from zero up to η g 1 means that the wrong colour factor is being used for the emission of the soft gluon g in a double-logarithmically enhanced region (for fixed kinematics of the soft gluon g, the double log arises from the integration over the transverse momentum and rapidity of g 1 ). This ultimately will lead to incorrect subleading N c contributions to LL and DL terms.
• Instead, with the segment method, colour factors are clearly separated into a C F region for emissions belonging to the q-leaf, and a C A /2 region for those belonging to the g 1 -leaf, as expected from Fig. 2a. A residual deviation from the exact matrix element is present in a region of rapidity localised around η g 1 , and reaches the level of +15% for ψ = 0, as shown in the lower row of plots (blue-red colour scale). If one integrates over azimuthal angle and rapidity, the blue and red regions will compensate each other. In section 6.3 below, we shall verify that that compensation is exact for large η g 1 , so that the method reproduces the correct total rate of soft emission, Eq. (3.8), as needed for NDL-FC accuracy in observables such as multiplicities.
• The NODS procedure reproduces the full squared tree-level matrix element (up to statistical fluctuations associated with our Monte Carlo sampling), as it should, since in this kinematic region the method is effectively using that full tree-level result to correct the leading-colour shower matrix element.

Differential matrix element: 2 → 4 + g configurations
Next, we consider tests for 2 → 4+g configurations, first for e + e − →qq q q +g and then for e + e − →qg 1 g 2 q + g. A first gluon g 1 is emitted off theqq dipole with the same kinematics  Colour factor assignment and relative deviation to the squared tree-level matrix element, as in Fig. 4, for the 2 → 4 + g configurations (a)qq q q + g and (b)qg 1 g 2 q + g, corresponding respectively to the Lund diagrams in Figs. 2c and 1a (with g 2 there moved to the right of g 1 ). The results have been obtained with the β = 0 PanGlobal shower algorithm.
as in Eq. (6.1). The second splitting is performed with qq q q configuration : zq = 1/4, η q q = 10, ψ = 0, (6.2a) These configurations are such that the second splitting happens at a much smaller angle than the first gluon emission. For the first configuration (g 1 →q q ), we choose a z fraction reflecting the absence of soft enhancements. For the second configuration (emission of g 2 from the quark, well separated in rapidity from g 1 ) we focus on a case where g 2 is much softer than g 1 , though the conclusions are unchanged if we take g 1 and g 2 to have commensurate transverse momenta. Results are displayed in Fig. 5. They have features similar to those of Fig. 4, albeit with a richer structure.
Theqq q q case is shown in Fig. 5a. The Lund diagram for the corresponding phase space, and expected colour assignments, was shown in Fig. 2c. There are several possible ways for the additional gluon to cluster with the rest of the event. It is useful to organise them from a Lund-plane viewpoint. The three different possibilities correspond to the left, centre and right panels in each plot of Fig. 5a. First, primary emissions (left-hand panel, labelled q) include cases where g 2 clusters with eitherq, q, or the full g 1 (→ q q )q system, corresponding respectively to η < 0 (not shown), η 5 and 0 < η 5. Next, secondary emissions (middle panel, labelled q ) include clusterings with either q (the harder of q andq ) or the g 1 → q q system. These correspond to the regions η 10 and 5 η 10 respectively. Finally, tertiary emissions (right-hand panel, labelledq ) correspond to clusterings of the soft gluon g withq . The hole observed in the primary (secondary) plane corresponds to the region where emissions are clustered in the secondary (tertiary) plane. The correct colour assignment involves a factor C F everywhere except in the region 5 η 10 of the q -labelled plot, where the emission of g is sensitive to the net colour-octet charge of the whole g 1 → q q system. The CFFE method tracks neither the intermediate particles nor segments, and therefore blindly applies a factor C F across phase-space. Our segment and NODS schemes, in contrast, display the correct behaviour along the intermediate gluon segment. In the case of the segment method, azimuthal deviations from the exact matrix element are observed to be of similar size as shown above in Fig. 4. Note also the discontinuity in the segment-methodq panel at η = 10. This is a consequence of our choice to make discontinuous transitions in the segment approach. Similar features are present in the other plots, but are less immediately visible because they coincide with the phase-space boundaries between primary and secondary Lund planes.
Next, we turn to theqg 1 g 2 q configuration (Fig. 5b), where the gluons g 1 and g 2 are primary emissions, and the holes both appear on the (primary) q-leaf. The corresponding Lund diagram would be Fig. 1a with g 2 there moved to the right of g 1 . In this case, the three panels correspond to primary emissions, emissions from the secondary g 1 Lund leaf and emissions from the secondary g 2 leaf. In the CFFE scheme, the region on the q-leaf extending from η = 0 to the gluon's position η = η g 1 = 10, is assigned a wrong colour factor of C A /2. For the segment and ME methods, the same region gets a correct factor of C F . In all cases a factor of C A /2 is correctly applied to emissions on the g 1 -and g 2 -leaves.
Note that the discrepancies in the CFFE approach are log-enhanced, and therefore as we discussed in section 6.1, can affect subleading-colour contributions to DL terms. For the segment method, the discrepancies are localised around the rapidities of the emission in the parent event and they are designed to integrate to zero. We will explicitly test (and confirm) this below in section 6.3.
The results in this section demonstrate that even though our segment and NODS corrections only ever consider the structure of double-emission matrix elements, their iteration reproduces the full-colour soft matrix elements for higher multiplicities, in the appropriate angular-ordered limits.
For a study of the NODS method for parent events that are not angular-ordered, see Appendix A, which shows that discrepancies do then arise relative to the correct full-colour matrix element, as is to be expected. For some configurations they reach 15%. The results are shown for several parton showers and colour schemes for emission from aqg 1 q system. The kinematics of g 1 are indicated in each column, and the integral bounds are set to η cut = 30. The results are colour-coded green and red, respectively, according to whether the result is consistent with zero, or not.

Integrated results
As noted in section 3.2, for observables such as the multiplicity, a shower that gives the correct full-colour rapidity and azimuth integral for the rate of soft emissions is expected to be correct at NDL-FC accuracy. Here we evaluate such integrals explicitly. Keeping in mind Eqs. (3.8) and (3.15), we compute the parton shower's integrated rate of emission, I PS , forqg 1 q + g, with the collinear cut-off set to η cut = 30 (i.e. we integrate over the whole rapidity and azimuth phase-space, but exclude cones of half-angle corresponding to η cut around each parton). Fig. 6 summarises the integral results for the CFFE, segment, and NODS colour schemes, for each of three different kinematic regimes for g 1 : a soft-collinear regime, a hard-collinear regime and a soft large-angle regime (cf. values of z g 1 and η g 1 as labelled on the plots). The known full-colour result is: which extends Eq. (3.8) beyond the small-angle limit for g 1 (i.e. Eq. (6.3) is valid for any η g 1 ). The plot shows the difference between the shower result and Eq. (6.3), multiplied by This normalisation is chosen so that N (I PS − I FC ) is equivalent to the effective net extent in rapidity over which there is a C A /2 versus C F discrepancy.
For soft-and hard-collinear configurations (left and middle columns), one observes a deviation in the integral computed with the CFFE scheme due to the mis-assignment of a factor C A /2 on the q-leaf, see Fig. 4. The difference between the PanScales shower family and the Pythia 8 parton shower is explained by the different separation of dipole elements mentioned in section 5. In comparison, the segment and NODS colour schemes reproduce the expression expected from angular ordering for all the parton showers we considered.
The last configuration that we examine is when g 1 is soft and at large angles (rightmost column of Fig. 6). In the case of the CFFE scheme, there remains a disagreement for all parton showers, although it is again smaller for the PanScales family than for Pythia 8.
Interestingly, the segment method also shows a residual deviation from the analytic expression. This expected deviation stems from our use of the small-angle limit to determine transition points between segments also in the large-angle region. It can be calculated explicitly for PanScales showers following the same procedure as in section 3.2 for arbitrary values of θ g 1 q : where η t = 1 2 ln 1+ √ 5 2 0.241 solves η t = 1 2 ln(1 + e −2ηt ). The fact that Eq. (6.5) is of order α s /N c for η g 1 close to zero, together with its exponential decrease as η g 1 increases, ensures that its impact on observables such as the multiplicity is NNDL-NLC.
As expected, the NODS scheme reproduces the full-colour analytic integral even when g 1 is at large angles.

Numerical tests for observables
In this section, we follow the approach proposed in Ref. [12] for testing the logarithmic structure of showers with specific observables in e + e − collisions. Such tests implicitly probe virtual as well as real corrections and so provide a powerful verification of the correct overall assembly of the different shower elements. As the reader will recall from the introduction, the tests can be separated into two broad classes, according to whether one is probing a double-logarithmic structure (DL, NDL, etc.) or an exponentiated (LL, NLL, etc.) structure. We start by outlining the structure of the tests, keeping in mind that we will adjust the details for each specific observable.
For example to test N k DL accuracy we will study a quantity such as 12 where V N k DL is the known N k DL prediction from resummation and V PS is the result from the parton shower. For a parton shower that is correct to N k DL accuracy, δV N k DL should be zero. Values of ξ for different momentum ranges are shown in table 1. In practice we will often use ξ = α s L 2 = 5, which is towards the upper end of the phenomenologically relevant combinations of α s and L accessible at the LHC. We perform such studies for multiplicities (section 7.1) and event shapes (section 7.2.1). For observables whose logarithmic prediction exponentiates, Eq. (1.1), we can study ln V (α s , L), taking the limit of α s → 0 with fixed λ = α s L . (7.3) To test N k LL accuracy we can examine where V N k LL is the known N k LL prediction from resummation. For a shower that is correct to N k LL accuracy, δ ln V N k LL should be zero. For a shower that is incorrect at LL accuracy, the unnormalised difference between ln V PS and ln V LL diverges for α s → 0, hence Eq. (7.4) multiplies that difference by α s to give a finite, non-zero result. 13 In practice we will often use λ = α s L = −0.5. This corresponds to a slightly narrower range of logarithm than our choice for ξ, in part to help mitigate some of the technical difficulties of the α s → 0 limit. We perform such studies for event shapes (section 7.2.2) and non-global logarithms (section 7.3).
Generation with very small α s and fixed ξ or λ is often difficult. Many of the techniques that we use were outlined in the supplemental material to Ref. [12]. For the work presented here we added three main new advances: 1. We implemented a weighted generation technique that is equivalent to evolving multiple replicas of an event, discarding a replica when it emits into a region of phase-space that we wish to veto, and then adjusting the number of replicas and their weights so as to continue generating with the original effective number of replicas (cf. section 3 of Ref. [49]). For the combinations of α s , shower and event-shape that were most challenging in Ref. [12], this enabled us to save about an order of magnitude in computing time, associated with accessing regions with very strong Sudakov suppression. It also enabled us to reach small α s values that were simply not feasible in Ref. [12], facilitating the extrapolation to α s = 0.
2. We adjusted the shower implementation so that it can track differences in directions between neighbouring particles in the dipole chain. This works around issues that arise in normal shower implementations where it becomes difficult to determine angles between particles (and dot products, etc.) when those angles go below machine precision . This, together with the next point, was especially useful in allowing for smaller α s and larger values of the (absolute) logarithm in double-logarithmic tests, though it also facilitated cutoff dependence tests in the NLL event-shape studies. It has a small ∼ 30% speed penalty, and some implementation overhead, but avoids the need for the double-double and quad-double types [50] used, with substantially larger speed penalties, in Ref. [12].
3. To allow for momenta across such disparate scales that the logarithm of the ratio of scales is truly large (|L| 1 2 log µ 354 where µ 1.8 × 10 308 is the maximum number that can be represented in double precision), we implemented a new floating type that supplements a normal double-precision number with a 64-bit integer to store the exponent, replacing the usual 11 exponent bits in a double precision number. This came with a speed penalty relative to double precision numbers of about ×4, but was substantially faster than solutions we investigated based on the MPFR [51] or Boost multiprecision libraries. It was particularly valuable for double-logarithmic tests and useful also for tests of non-global observables.
Throughout we will run with the physical colour factors, N c = 3, C F = (N 2 c − 1)/(2N c ) = 4/3, C A = N c and n f = 5 (except for leading-N c comparison results, where we use C F = 1 2 C A = 3 2 ).

Particle multiplicity
One of the most powerful tests of a parton shower is its prediction for the particle or subjet multiplicity. To reproduce the NDL multiplicity requires correct modelling of a variety of aspects of a parton shower, including nested g → gg and g → qq splittings. Ref. [12] demonstrated agreement between a range of dipole showers and NDL-LC predictions. We expect the segment and NODS schemes of sections 3 and 4 to bring NDL-FC agreement. For example, in the segment method, the general pattern of C F dipole segments associated with the Bornqq pair should guarantee DL-FC terms, while the specific choices of transition points for the segments, together with the C F segments for g → q q splittings, should ensure NDL-FC accuracy. Strictly speaking, the particle multiplicity is infrared and collinear (IRC) unsafe. A closely related, but IRC-safe, quantity is the subjet multiplicity [52] in the k t jet algorithm [53]. Ref. [12] directly computed that subjet multiplicity comparing to the results of Ref. [52]. Here we take a slightly different approach that avoids the need for jet clustering, but can still be directly compared to the results of Ref. [52]. We run the shower as normal, but set the running coupling to zero below some threshold transverse momentum k t,cut . At NDL accuracy, the multiplicity with such a procedure should be identical to that with a clustering threshold for the subjets of y cut = k 2 t,cut /Q 2 . Only starting from NNDL do we expect differences to arise between the multiplicity for a shower definition of a k t cutoff and the k t algorithm definition. 14 We have explicitly verified that this is the case at leading colour.
The analytic expression for the multiplicity in e + e − → qq events, up to and including NDL terms, can be straightforwardly extracted from Ref. [52], (7.5) in the approximation that |L| Fig. 7 shows the multiplicity N for ξ = 5 at a value of α s = 5 × 10 −12 , which is sufficiently tiny that one can neglect subleading corrections to within an accuracy of √ α s 2.24 × 10 −6 . The figure compares several shower algorithms (the main PanScales showers and our implementation of the Pythia 8 shower), and three colour schemes, with the analytic DL 14 Those differences can for example arise because of details of the specific definition of transverse momentum in the collinear limit, and because of effects related to the jet algorithm's clustering properties. result. We see that the CFFE scheme, the default in many showers including Pythia 8, differs by up to ∼ 3% from the analytic result. A 3% difference is not a huge effect, but it remains a subleading-N c DL difference and gives a measure of the practical impact of such contributions, i.e. ∼ 1/(3N 2 c ). For comparison, the differential matrix element plots, Figs. 4 and 5, showed ∼ 10% differences for the CFFE approach in extended phase-space regions. Our new segment and NODS approaches coincide well with the DL analytic multiplicity result for all showers, including when applied to the Pythia 8 shower.
To test the NDL multiplicity terms, we need to examine δN NDL , as defined in Eq. (7.2). For this comparison, we need α s to be large enough for √ α s terms in the multiplicity to be visible after summing over a finite number of Monte Carlo events, while also small enough that we can safely extrapolate away O (α s ) terms. To this end, Fig. 8a shows a quantity conceptually similar to δN NDL in Eq.

as a function of
√ α s . It shows results for the NODS method for several showers and one sees that most of the results go to zero as α s → 0, with the exception of the red curve, to which we return below (the segment method, not shown, gives almost identical results). The result of the actual α s → 0 extrapolation (using a cubic polynomial extrapolation based on α s = {0.000005, 0.00032, 0.00128, 0.00512}) is shown in Fig. 8b for each shower and for each of the segment and NODS methods. All results are consistent with the (full-colour) NDL result, to within the small statistical errors. α s → 0, as required to achieve full-colour NDL accuracy. The curves depict a fit that is a polynomial in powers of √ α s . (b) Explicit extrapolation of N PS −N NDL N NDL −N DL to α s = 0, for a range of showers, with the segment and NODS colour schemes. In both plots we also illustrate the NDL-level discrepancy that arises in a scheme where quarks from g → qq splittings erroneously continue to emit with a C A /2 colour factor. Fig. 8 also includes results (red points) obtained using a deliberately erroneous prescription that omits the insertion of new C F segments following g → qq branchings, resulting in those quarks emitting with a C A /2 colour factor. While this prescription gives the correct DL-FC result, one should expect it to fail to reproduce the NDL-FC results, because g → qq splittings start to contribute to the multiplicity from NDL terms onwards, as can be verified by inspecting the n f terms in Eq. (7.5). 15 We see in Fig. 8a that with this incorrect treatment of g → qq splittings, the limit √ α s → 0 fails to converge to the NDL expectation. Instead the extrapolated NDL coefficient is ∼ 3% larger than the analytic expectation.
We have also carried out similar tests for the multiplicity in H → gg events for all showers shown in Fig. 8 and found a similar level of agreement with the full-colour NDL predictions for both the segment and NODS colour prescriptions.

Event shapes
The main event shapes that we consider here are obtained by considering all primary Lund declusterings [29], and for each declustering evaluating, where k ti is the transverse momentum of the declustered subjet i with respect to its partner direction, and η i = ln tan θ i /2, with θ i the angle of the declustered subjet with respect to the partner. The parameter β obs determines the relative weighting of different rapidities and we will consider three values, β obs = 0, 1 2 , 1. We will study two combinations of the u i , For each event shape, we define Σ(α s , L) to be the fraction of events for which that event shape has a value smaller than a threshold defined to be e −|L| . Up to NLL accuracy, Σ(α s , L) for M β obs =0 coincides with that for the Cambridge √ y 23 resolution scale [43], and S β obs =1 with that for one minus the thrust (below, in section 7.2.2, we will consider those explicitly as well).
Ref. [11] showed that the CFFE procedure led to spurious subleading-N c terms starting from order α 2 s L 4 in standard dipole showers. With similar reasoning it is straightforward to show that the issue is present also for the PanScales showers with the CFFE approach. Those issues are caused by mis-attribution of a C A /2 colour factor to emissions that should be seen as coming from the primaryqq pair. Once that issue is fixed, one should obtain NLL-FC accurate results for all PanScales showers that already give NLL-LC accuracy.
Given that the CFFE issue arises at order α 2 s L 4 , one expects to be able to observe it numerically in both DL and LL-style tests. As we shall see below, numerical LL and NLL tests require us to prune the shower branchings, so as to keep multiplicities under control in the limit α s → 0 for fixed λ = α s L. That pruning can be delicate in situations with DL discrepancies. Accordingly we first carry out DL tests, for which we can run complete, unpruned showers.
When it helps to limit notational ambiguity with respect to β obs , in some cases below we will write β PS instead of β for the parameter that determines the shower ordering variable (cf. Eq. (3.11)).

Double-logarithmic study
To enable us to study the subleading colour double-logarithmic contributions, Fig. 9 shows Σ PS /Σ DL as a function of ξ = α s L 2 , in the limit of α s → 0, for M β obs observables with three different values of β obs . It includes results for two PanScales showers and for Pythia 8, for all three colour schemes. To understand the plots, it is useful to recall that the PanScales showers are ordered in a variable ∼ k t e −β PS |η| (we show PanGlobal with β PS = 0 and PanLocal with β PS = 1/2). We expect (and see) that double logarithmic discrepancies are present in the CFFE scheme for any β obs > β PS . 16 For Pythia 8, which is k t ordered, we expect DL discrepancies to be present in the CFFE scheme for any β obs = 0, though our tests only consider β obs ≥ 0. Discrepancies at DL accuracy are clearly visible as deviations of the curves from 1. In the physically accessible region, ξ 7, those discrepancies are fairly modest, no more than 2−3 percent. The results for the segment and NODS colour scheme are all consistent with 1. While shown only for a subset of showers and for the M β obs observables, the results for other PanScales showers are the same for any given β PS and for the S β obs observables.
In the case of the PanScales showers, one can also calculate the analytic form of the CFFE discrepancies, cf. Appendix B. The analytic results agree well with the observed shower discrepancies. It is interesting to consider the large-ξ limit of the discrepancy, which will be relevant also in understanding the results of the next subsection. For example, for 16 For β obs ≤ βPS, a given contour of fixed M β obs is traversed by the shower in order of increasing absolute rapidity, which means that relevant emissions occur either on the primaryqq dipole, or on aqg or gq dipole. Taking the example of a gluon emitted at some rapidity yg > 0 and forming a gq dipole, emissions that are closer to the quark than to the gluon, i.e. along the quark with rapidity y > yg, are emitted with the correct CF colour factor. Since the only emissions that can increase the event shape value are those that occur at rapidities larger than previous gluon emissions, one obtains the correct full-colour double logarithmic result.
the thrust determined with any of the PanScales showers with β PS = 0, we find that However, the approach to this limit is very slow: for ξ = 5 the ratio is about 1.009, to be compared to an asymptotic value of 1 + 1/17 1.0588.

LL and NLL studies
Our next set of tests is to carry out LL and NLL studies for the same set of event shapes discussed above, supplemented with the total and wide-jet broadenings (B T , B W ) [54], the Cambridge y 23 jet resolution parameter [43], thrust [55] and fractional energy-energy correlation moments (FC x ) [56]. We consider the limit α s → 0 for fixed λ = α s L.
A fundamental difficulty with these LL and NLL studies is that the logarithm of the multiplicity, ln N , scales as √ α s L 2 , cf. Eq. (7.5), and when we fix λ = α s L and take the limit α s → 0, the logarithm of the multiplicity blows up as 1/α s . Since generation time and memory consumption scale at least in proportion to N , event generation becomes prohibitively expensive with too small a value of α s .
To address this problem we adopt a strategy where, at each stage of the shower, branchings that are guaranteed to be irrelevant to the observable are vetoed. Specifically, for a given β obs , we track the maximum value of O approx = k t,approx e −β obs |ηapprox| that has occurred so far in the showering of the event (using Eq. (3.12) for η approx and k t,approx equivalent to the k t in Eq. (3.11)). We accept a given new branching only if it has O approx > e −∆ O approx,max , with ∆ a parameter to be chosen appropriately. The logarithm of the multiplicity is then expected to scale roughly as √ α s ∆ 2 +ln |λ∆|. The value of ∆ should be small enough that the multiplicity remains under control, and large enough that recursively IRC safe observables [56] are not affected by it, which is ensured by the requirement e −∆ 1. For shower-observable combinations where there is a double-logarithmic subleading-N c effect, we also need to be aware that the mechanism generating that effect may only be fully operational if α s ∆ 2 1, which is in tension with the constraint on the multiplicity. In practice it is difficult to ensure that the full DL discrepancy is captured in LL tests, especially when considering the interplay between the α s → 0 and ∆ → ∞ limits. Our approach will instead be to ensure that the combinations of α s and ∆ are such that the presence of any DL issues is correctly diagnosed in our LL tests, even if we do not reproduce the exact value of any discrepancy. This is illustrated in Fig. 10, which shows LL tests for the CFFE colour approach, using ∆ = 18 and a quadratic polynomial extrapolation to α s = 0 based on runs at α s = {0.0025, 0.005, 0.01}. There, one clearly sees that the set of shower-observable combinations that fails the LL test is consistent with the expectations from the DL tests in Fig. 9 and from Appendix B. However, taking the example of the β obs = 1 observables, for the PanGlobal shower with β PS = 0, one sees a discrepancy of 0.056, to be compared to the expectation of 0.068 in table 2 of Appendix B, as obtained when one first takes the limit ∆ → ∞, and then α s → 0 (this 0.068 is the running coupling analogue of Eq. (7.9)). The LL tests for the segment and NODS colour schemes (not shown in Fig. 10) are consistent with the analytic LL results for all observables and showers, including Pythia 8. Accordingly in Fig. 11, where we show the full-colour NLL tests, i.e. examining δ ln Σ NLL , Eq. (7.4), we include results just for those two schemes. 17 For all showers that were in agreement with NLL predictions at NLL-LC in Ref. [12], we now see that these two new colour schemes ensure agreement with the NLL-FC predictions. This is to be expected, because, beyond the recoil issues that were relevant for NLL-LC accuracy, NLL-FC requires the correct treatment of the colour factor only for emissions that are widely separated in rapidity, which the segment and NODS schemes both accomplish by design.
Note that we have only tested event-shape observables that vanish in the 2-jet limit. One could also envisage testing event shapes such as the thrust minor and D-parameter, which vanish in the limit of 3 narrow jets, and whose NLL-FC resummations have long been known [57,58] for planar events. Our expectation is that the NODS scheme (but not the segment scheme) will yield the correct NLL-FC results also for these observables, as long as the appropriate matching is included for the 3-jet matrix element, and the segment variables and auxiliary momenta are properly initialised.

Energy flow in a rapidity slice
The final observable that we consider is the probability, Σ(α s , L), for the energy in a given central slice of rapidity to be less than e −|L| Q. Such an observable is of interest because its resummation involves non-global logarithms, single-logarithmic terms α n s L n that involve configurations with an arbitrary number of (soft, large-angle) gluons in the neighbourhood 17 One can also carry out NLL tests for observable-shower combinations that have the correct LL-FC result in the CFFE approach. For the PanScales showers, where LL-FC results were correct for β obs ≥ βPS, one finds that most observables are correct at NLL-FC only for β obs > βPS, with the exception of those in the max-type class, which remain correct for β obs ≥ βPS. NLL accuracy tests NODS method Figure 11: NLL global event-shape tests of the segment and NODS colour schemes, showing NLL agreement for β = 1/2 PanScales showers and for the β = 0 PanGlobal shower. In contrast to the NLL-LC tests of Ref. [12], the Pythia 8 β obs > 0 results here are coloured green rather than amber, because our colour code does not incorporate the information about failure of exponentiation in fixed-order tests, tests that we have not explicitly repeated for this paper.
of the slice [22,59] (see also Ref. [60]). The full-colour resummation for such observables is sensitive to arbitrarily complex colour correlators, both in the real emissions and the virtual corrections, which need to be evaluated at amplitude level. The resulting subleadingcolour single-logarithmic corrections go far beyond the scope of the colour schemes that we introduced in sections 3 and 4. In particular, we expect the segment scheme to be correct at full colour only up to order α s L, and the NODS scheme to be correct at full colour up to order α 2 s L 2 . Recall that leading-colour all-order single-logarithmic accuracy for PanScales showers was demonstrated in Ref. [12].
The question we ask in this section is to what extent our schemes differ from a fullcolour computation. Such a calculation was performed for the energy flow in a slice for the Z →qq process, by Hatta and Ueda [28]. Their approach was based on a refinement of a proposal by Weigert [61], whereby the problem was reduced to the simulation of associated Langevin dynamics in the space of Wilson lines, and solved on a two-dimensional angular grid. It was also applied, with Hagiwara, to the calculation of a hemisphere observable [62] and, very recently, to H → gg decays and a number of 2 → 2 scattering processes [63]. The approach was formulated in such a way that it included only the single logarithms, α n s L n . Accordingly, in using it as a reference to which to compare our α s → 0 limit, we can be sure that any difference is exclusively associated with subleading colour effects.
Two other calculations are based on parton showers at finite α s and a truncation of the 1/N 2 c series: Nagy and Soper examined rapidity gaps between dijets at hadron colliders [6,33], while De Angelis, Forshaw and Plätzer examined the energy in a slice for the Z →qq and H → gg processes [8]. Only the processes examined in that latter paper are within the scope of our work here, but their simulation at finite α s precludes a meaningful direct comparison, because it would be impossible to know whether any differences between their results and ours are associated with subleading-colour effects or instead subleading logarithmic (α n s L n−1 ) effects. The specific definition of the energy in a rapidity slice is as follows: we cluster each event with the e + e − Cambridge algorithm [43] with the y cut parameter set to 1 and then undo one step of the resulting clustering sequence to obtain two back-to-back jets. Defining rapidity with respect to those two jet axes, we examine the total energy contained in the rapidity region |y| < − ln tan θ cut /2, with θ cut setting the boundary of the slice. We will make the choice θ cut = π/3, which coincides with that of Hatta and Ueda [28]. As in section 7.2.2, if we carry out a full run of the shower in the α s → 0 limit for fixed α s L, parton multiplicities will grow too large to handle. Here, we solve this problem by vetoing emissions whose rapidity with respect to the closer end of the emitting dipole is larger than η max , which by default we take to be 10 (we will verify that reducing it to 8 does not modify the results).
It is commonplace in studies of non-global logarithms to show the results as a function of and where the result on the right-hand side is given for a one-loop running, which is the only contribution that survives in the limit α s → 0 at fixed λ = α s L. Values of τ for various momentum ranges are shown in table 1 and, at the LHC, the largest accessible value is τ 0.4. Fig. 12 (left) shows the fraction of events, Σ(τ (α s , L)), whose energy flow in the slice is less than e −|L| Q. It includes several results: a leading-N c result with C F = 1 2 C A = 3/2, the Hatta-Ueda full-N c prediction, 18 and our predictions with the segment and NODS 18 We are very grateful to Hatta and Ueda for rerunning their code with higher statistics than in their The agreement is potentially surprising given that our schemes do not achieve NLL-FC (α n s L n ) accuracy for non-global observables. The thin band for our results represents the statistical uncertainty added in quadrature to estimates of systematics obtained using the difference between our default runs (η max = 10 and α s = 0.7 × 10 −8 ) and runs with η max = 8 and α s = 1.4 × 10 −8 . Our results for other showers with the same colour schemes are very similar, as is to be expected.
methods. Recall that those methods are not expected to work beyond order λ and λ 2 respectively. However in Fig. 12 (left) they are indistinguishable from the full-N c Hatta-Ueda result. To further probe this observation, the right hand plot shows ratios to a reference, which we take to be the PanLocal-antenna β = 1/2 NODS (the specific choice is largely immaterial, since our aim is to compare different predictions on this ratio plot). One sees that the difference between the full-N c Hatta-Ueda result and our leading-N c result is about 23% at τ = 0.4. Remarkably, both our segment and NODS methods seem to be in good agreement with the Hatta-Ueda result across the full range of τ : the whole range is within two standard deviations of the Hatta-Ueda result, and in much of the range the agreement is within one standard deviation. Some caution is needed in interpreting these results: firstly, they correspond to one specific choice of slice size. Secondly, when using a finite-resolution angular grid (as in the Hatta-Ueda approach), there are inevitably original paper and providing us with the corresponding numerical results. some residual systematic effects associated with that finite resolution, and in this case we cannot exclude the possibility that they are comparable to the statistical error. 19 A further observation is that the segment and NODS results come out largely identical, even though the former (latter) is FC-accurate only to α s L (α 2 s L 2 ). We attribute this to a partial cancellation that arises after azimuthal integrations. In particular, in Appendix C we demonstrate that if one breaks the azimuthal symmetry by considering a patch in azimuth and rapidity instead of a slice (which covers all azimuths), a difference between the segment and NODS scheme does appear at the statistically significant level of a couple of percent.
A final comment concerns the same slice observable for H → gg events. For a given value of τ , both our segment and NODS methods are identical to the leading-N c (2C F = C A = 3) result in the α s → 0, fixed τ limit, because g → qq branchings only affect NNLL terms, α n s L n−1 . Subleading-N c effects start for the slice only from order α 3 s L 3 onwards, because soft-gluon emission from a ggg system (i.e. α 2 s L 2 ) is given exactly by the sum of emission from three C A /2 dipoles. Therefore we have the same formal subleading-N c accuracy as the NODS method for Z → qq: exact up to α 2 s L 2 , and only leading-N c for terms α 3 s L 3 onwards. Very recently, Hatta and Ueda have shown that an FC calculation for the slice observable in H → gg is in good agreement with the large-N c result, to within statistical errors [63]. The findings of Ref. [8] appear to be consistent with this conclusion. 20 Accordingly, the NODS and segment methods agree with the FC results for H → gg, to within the precision of the predictions for the latter.

Timing assessment
One of the motivations behind this work is to provide a method to include logarithmically relevant subleading colour effects in parton showers that is simple and computationally efficient. In this section we show that the two methods highlighted in sections 3 and 4 bring only a modest time penalty.
To check this, we have run the PanLocal shower, in its antenna variant, setting β = 1 2 , using a fixed-coupling prescription, α s = 0.1, and varying the ordering variable (ln v) cutoff scale between −24 and −10. At leading-N c , the PanLocal antenna shower is currently one of the fastest of our showers, and since the time penalty of the colour schemes is largely independent of the shower, the choice of PanLocal-antenna provides a worst-case estimate of the relative impact of the colour schemes. Following a similar logic, we use a fixed coupling for these tests, because our running-coupling shower implementation is not 19 We base this statement on a run of our PanScales showers, assigning each particle to a bin on an 80 × 60 grid in cos θ and ψ (for extremal cos θ bins, we choose to map all ψ values to a single bin). We then accept emissions only if they are in a distinct bin of the grid from both the emitter and spectator of the parent dipole. Comparing this to our standard continuum-limit runs, we see effects of the same order of magnitude as the Hatta-Ueda statistical uncertainty. 20 Specifically, in its Fig. 7a, the subleading colour curve ("all, d = 2") is in good agreement with the leading-colour curve ("LCV+R"), after rescaling the latter by 8/9 to account for N 2 c − 1 rather than N 2 c in the Born H → gg colour sum, notably in the region ρ ≡ e −|L| 0.1, where statistical fluctuations seem to be under control at the few percent level. yet fully optimised, and so would produce a misleadingly high timing baseline to which to compare the timing penalties of the new colour schemes. Fig. 13a shows the average event showering time divided by the average number of emissions, for each value of the ln v cutoff, plotted as a function of the average multiplicity N for that ln v. It includes results for four colour schemes. One sees that the timing per emission for the C F = C A /2 = 3/2 and CFFE schemes is about 1µs, almost independently of N . The penalty for the segment and NODS schemes is at most 0.2µs and 0.5µs respectively. 21 It decreases with increasing N , which is probably a consequence of the fact that the number of segments per dipole decreases with increasing N , cf. Fig. 13b.

Conclusions
In this work on subleading-colour effects in final-state dipole and antenna showers, we have paid particular attention to the interplay between colour and logarithmic accuracy. Insofar as the accuracy of current parton showers is at best NLL, we argued that an essential requirement is subleading-colour accuracy for LL terms, with potentially missing subleading-N c NLL contributions being on the same footing as leading-colour NNLL terms, which are currently not available in parton showers.
In sections 3 and 4, we outlined two schemes that are simple and efficient to implement in a range of parton showers (including, e.g. the Pythia 8 shower) and that provide fullcolour LL accuracy (LL-FC). One of them, the segment method, varies the colour factor along a dipole, using a Lund-diagram type classification to identify regions that should have either a C F or a C A /2 colour factor, according to whether an angular-ordered picture implies emission from a quark or a gluon. 22 The other method, dubbed NODS, nests fullcolour energy-ordered double-soft matrix-element corrections for emissions from any dipole that, according to the segment approach, contains at least one C F segment.
In practice it was possible to engineer our approaches so as to provide full colour accuracy beyond the LL (or DL) approximation for a range of observables: all global event-shape variables, as well as observables that examine particle or jet multiplicities. For these, both of our schemes achieve NLL-FC or NDL-FC accuracy, as appropriate for the specific observable (aside from spin correlations, which do not affect the observables studied here at NLL-FC accuracy, and whose study we postpone to future work). This involved care with g → qq splittings and with the specific locations of C F to C A /2 transitions in the segment method. Numerically demonstrating the resulting accuracy relied on a range of techniques for taking the α s → 0 limit of our showers while maintaining fixed α s L or α s L 2 , techniques that we have improved relative to earlier work [12].
The one context in which we do not achieve NLL-FC accuracy is for non-global observables (in our definition, LL for these is zero, and they start at NLL, i.e. α n s L n ). There we achieve full-colour accuracy only for a finite set of n: for our segment method, n ≤ 1, and for our NODS method, n ≤ 2. Remarkably, despite this limitation, our schemes still give good numerical agreement with the full-colour NLL results for the energy in a slice, as obtained by Hatta and Ueda [28], to within the 1−3% statistical error of the latter. In future work, it would be interesting to gain a better understanding of why this is the case.
One of the important points related to the computational efficiency of our colour schemes, is that it makes it straightforward to obtain high statistical accuracies. The time penalty for the more sophisticated scheme, NODS, was below a microsecond per emission. That has enabled us to draw robust conclusions about non-trivial sub-leading colour effects, which were typically at the level of a few percent. This provides a target for future work on subleading colour schemes, which should straightforwardly be able to provide the high statistical precision needed to conclusively determine the size of any subleading colour effects that exist beyond those accessible in simple schemes such as those developed here.  Figure 14: Results from the NODS scheme for the PanGlobal (β = 0) shower, in an initial configurationqqg 1 g 2 + g where the two gluons g 1 and g 2 are not strongly angular-ordered. We show three configurations with the gluon g 2 brought increasingly closer to g 1 : for a fixed η g 1 = 5, ψ g 1 = π, the second gluon g 2 is placed at η g 2 = 7 (top), η g 2 = 6 (middle) and η g 2 = 5 (bottom) at an azimuthal angle ψ g 2 = 0. Black dots convey the largest negative deviation between the NODS scheme and the exact tree-level matrix element.
rate of emission, σ PS , and the full-colour tree-level matrix element, σ FC , is observed to be relatively small for ∆η g 1 g 2 = 2, where the largest negative deviation is identified by a black dot and reaches −4% at ψ g = π/2 on the q-leaf. As the second gluon g 2 moves to an angle close to that of the first gluon g 1 , the discrepancy becomes larger, reaching a value of −9% for ∆η g 1 g 2 = 1, and −14% for ∆η g 1 g 2 = 0. Note that the regions where the matrix element is not reproduced are, as expected, localised around kinematic configurations for which the three emissions are at commensurate angles. These configurations are beyond the scope of the accuracy we are aiming for.
The initial 4-partonqg 1 g 2 q configurations studied here differ from the hadron-collider 2 → 2 processes such as those studied in Ref. [63], in that g 1 and g 2 are both soft. Still, our observation that non-trivial subleading-colour effects decrease as g 1 and g 2 become more separated (which is to be expected based on coherence arguments), may be relevant also more widely. Specifically, they suggest that if one seeks to maximise subleading-colour effects in hadron-collider processes, it may be of interest to consider 2 → 2 configurations where the outgoing jets are both at the same rapidity.

B Analytical results for CFFE subleading-N c effects in event shapes
Here, we develop a (semi-)analytical approach to determining the expected DL and LL deviations in event shapes for the CFFE colour scheme as applied to the PanScales family of showers. The intent is to help provide an independent validation and understanding of some of the numerical results shown in section 7.2. The analysis is common to all of the PanLocal-dipole, PanLocal-antenna and PanGlobal showers. In contrast traditional dipole showers such as Pythia 8 are, we believe, more complex to analyse and, as can be seen from Fig. 9, they clearly yield different CFFE DL results.
Let us consider an observable with β obs > β PS and compute Σ PS (L). Whenever we have an emission at a scale x = ln v and rapidity η, say with η > 0, s.t. x − β obs η < L, subsequent emissions (real or virtual) at rapidities between 0 and η will have a factor C A /2 instead of C F . This gives the following expression where the r in and r out Sudakov exponents correspond to the light blue and pink regions in Fig. 15, respectively. When writing this expression, the second term describes the case where the first emission below the observable boundary is at a scale ln v = x. Compared to the expected result, one gets a veto for emissions below the boundary with ln v > x -the r out contribution -as well as a change C F → C A 2 in the region corresponding to r in . (The r out factor is minus the derivative of r out with respect to x.) The first term in the square bracket describes the situation where there is no emission below the observable boundary that affects the observable Sudakov. Finally, the overall square accounts for both hemispheres.
The exponents r in and r out can be easily computed including one-loop running coupling effects. Here we just give their fixed-coupling expressions:   Fig. 16, illustrating the perfect agreement between the two.
The next limit we want to study is the one where we keep λ = α s L fixed and take α s → 0, corresponding to the studies presented in section 7.2.2. Here, running-coupling effects can no longer be neglected. We first change the integration variable in (B.1) to ν = α s x, which should be integrated over a finite range going between λ and 1+β PS 1+β obs λ. For fixed λ, the Sudakov exponents are proportional to 1/α s and so become exponentially suppressed  Table 2: Value of the LL discrepancy for α s → 0 with fixed λ = α s L = −0.5 for various shower-observable combinations in the CFFE colour scheme.
when α s is taken to 0. One can thus evaluate (B.1) in the saddle-point approximation which becomes exact in the limit α s → 0. If one includes running-coupling effects, the saddle-point equation can only be solved analytically for β PS = 0 and we have used a numerical evaluation for β PS > 0. For simplicity, let us quote the analytic result obtained in the fixed-coupling approximation, where we find We note that this result is equivalent to taking the limit ξ → ∞ in (B.3). In particular, it reproduces Eq. (7.9). For completeness, we list in table 2 the values of ln Σ PS ln Σ FC − 1 with fixed and running coupling.

C Segment versus NODS results for patches
When comparing our subleading-colour results in section 7.3 for the energy in a rapidity slice to the full-colour ones obtained by Hatta and Ueda [28], the apparent agreement between the segment and NODS methods might come as a surprise, given that the NODS method reproduces the full-colour double real energy-ordered soft matrix element, while the segment method does not. The difference between the methods is clearly visible in our matrix-element tests, see e.g. Fig. 4. For the energy in a rapidity slice, this means that the segment method should deviate from the full-N c result from O(α 2 s L 2 ) onwards, one order earlier than the NODS method for which deviations are expected from O(α 3 s L 3 ). One possible explanation for the similarity between the segment and NODS methods for the rapidity slice observable is that the O(α 2 s L 2 ) deviation in the segment method largely disappears when integrating over the azimuthal angle of the two emissions. To test this hypothesis, we study the energy deposited in a rectangular patch of both finite rapidity and azimuthal extent, |η| < η cut = − ln tan θcut 2 and |ψ| < ψ cut , thus breaking the azimuthal symmetry of the slice. Fig. 17a shows the ratio between the segment and NODS methods for a patch of the same rapidity width as used in section 7.3 and different azimuthal extent. We see that while the difference between the two methods is smaller than 0.5% for the rapidity slice (0.2% at τ = 0.4), it can reach a few percent for a rectangular patch with finite ψ extent. Among the patch extents we studied, |ψ| < π 2 showed the largest effect. Finally, Fig. 17b shows the ratio between the segment and NODS schemes measured at τ = 0.4 for different rapidity and azimuthal extents. We notice that the deviation from 1 increases as the rapidity width of the patch increases, reaching ∼ 5% for a patch with |η| < 2 and |ψ| < π 2 . The results here suggest that for future comparisons of non-global logarithms in methods implementing subleading-N c effects, it could be desirable to study not just slices but also patches.