Color matrix element corrections for parton showers

We investigate the effects of keeping the full color structure for parton emissions in parton showers for both LEP and LHC. This is done within the Herwig 7 dipole shower, and includes gluon emission, gluon splitting, initial state branching processes, as well as hadronization. The subleading Nc terms are included as color matrix element corrections to the splitting kernels by evolving an amplitude-level density operator and correcting the radiation pattern for each parton multiplicity, up to a fixed number of full color emissions, after which a standard leading color shower takes over. Our results are compared to data for a wide range of LEP and LHC observables and show that the subleading Nc corrections tend to be small for most observables probing hard, perturbative dynamics, for both LEP and LHC. However, for some of these observables they exceed 10%. On soft physics we find signs of significantly larger effects.


Introduction
At present time event generators [1][2][3][4] have developed into indispensable tools for understanding collider phenomenology. At the same time, the high energy available at the LHC has significantly opened up the perturbative phase space available for radiation. This increases the demand of resummation, performed either analytically or, as in the case of parton showers, numerically e.g. via the Sudakov veto algorithm [5,6].
From a QCD perspective, the high number of colored partons due to the large perturbative phase space, as well as due to the fact that the initial state partons carry color, calls for a better understanding of subleading N c effects. This is also in demand to ease the cancellation of infrared singularities in the matching and merging of parton showers with NLO calculations [7][8][9], which often fully include JHEP11(2018)009 subleading N c effects. It should also be stressed that we expect such corrections to account for subleading N c , leading-logarithmic problems which can arise in dipole-type parton shower algorithms due to an ambiguous definition of an emitter's color charge [10].
To achieve greater accuracy and better understanding of parton shower uncertainties, it is therefore time to include a better description of subleading color contributions in parton showers, similar to how parton showers recently have been improved by matrix element merging at leading order [11][12][13][14][15][16], and matching at next-to-leading order [7,[17][18][19][20][21].
The first steps in this direction has already been performed by some of the authors in the case of an e + e − -collider [22]. Others have pursued another road, keeping only a subset of the color suppressed terms [23,24], and detailed studies have been carried out towards systematically expanding virtual and real effects in shower-type evolution algorithms [25,26].
In the present paper we extend the color matrix element corrections, first implemented in [22], by including initial state hadrons as well as g → qq-splittings, subsequent leading color showering and hadronization. This is done within the Herwig 7.1 [27] implementation of the dipole shower algorithm [28], giving us a full-fledged general purpose event generator which can be used for studying color matrix element corrections to any process occurring at the LHC and other colliders, in practice up to a limited number of colored partons, restricted by the fast growing complexity in color space, however still reaching down to relatively soft emissions.
Our method is based on dipole factorization [29,30], which we outline in section 2. The complication brought about by the color matrix element corrections is discussed in section 3, and implementation details, involving evolution of the color structure treatment, the weighted Sudakov algorithm and the density operator, are discussed in section 4, 5 and 6 respectively. Results for various processes, including initial and final state radiation, standard QCD observables, heavy quark production and Z plus jets are discussed in section 7, and concluding remarks are made in section 8.

Essence of dipole factorization and dipole shower evolution
This paper is based on dipole factorization, stating that whenever the next gluon to be emitted from an n-parton configuration becomes either soft or collinear to one of the existing partons, the squared amplitude for the n+1-parton case can be approximated with |M n+1 (. . . , p i , . . . , p j , . . . , p k , . . .)| 2 ≈ k =i,j 1 2p i · p j M n (. . . , pĩ j , . . . , pk, . . .)|V ij,k (p i , p j , p k )|M n (. . . , pĩ j , . . . , pk, . . .) , (2.1) in terms of the old amplitude |M n . In the above, an emitterĩj → i, j whereas a spectator k →k absorbs the longitudinal recoil, such that all partons, before and after emission, stay JHEP11(2018)009 on-shell. For final state radiation we use the standard Sudakov decomposition, with p 2 ij = p 2 k = 0, a spacelike transverse momentum k ⊥ with k 2 ⊥ = −p 2 ⊥ and k ⊥ · pĩ j = k ⊥ · pk = 0. The cases of initial state emitter or spectator are discussed in [28].
Note that the sum in eq. (2.1) only runs overĩj =k. This is possible since the collinear singularity corresponding to the square of the diagram where partonĩj is the emitter, has been rewritten as a sum of interferences between that diagram and every other diagram using color conservation The splitting kernels are given in terms of the standard dipole splitting kernels as where Tĩ j · T k describes the color space effect of exchanging a gluon between partonĩj and partonk, and thus contains an implicit sum over gluon indices. In the massless case, the final-final dipole splitting kernels are given by Note, however, that dipole factorization is valid also for massive particles in the quasicollinear limit, and our implementation is general, using the massive dipole splitting kernels and kinematics available in Herwig, as detailed in [31]. For the dipole configurations involving initial state partons, we use the corresponding expressions as given in [30,32]. In eq. (2.7) the factors of ij in eq. (2.6); in order to use the standard definition of eq. (2.7), T 2 ij is introduced in the denominator of eq. (2.6). In the large N c limit the color correlator becomes where δĩ j ≡ 1 ifĩj is a gluon, and zero otherwise. In this way, the factor C F for gluon emission off a quark is reproduced in the large N c limit by coherent emission from the quark and its color-connected partner, and the factor C A is reproduced for gluon splitting JHEP11(2018)009 by the sum of the coherent emissions from the gluon and its two color-connected partners, hence the factor 1/2 for gluons in eq. (2.8).
We remark that the leading N c version of the above describes how to get a leading N c correct emission pattern. It does not address the issue of assigning a leading N c color flow from which subsequent emission can be performed. The default strategy is to use a color flow where the radiated gluon is inserted between the emitter and the spectator. However, with a gluon splitting kernel which is symmetric in z ↔ 1 − z, the question of which gluon is the radiator and which is radiated arises. In the limit where z is small, it is rather the parton with momentum fraction z which is soft and thus should be seen as radiated and therefore be inserted between the emitter (with the large momentum fraction 1 − z) and the spectator. For this reason, to mimic the swap of color, we swap momenta of the emitter and the emitted parton with probability 1 − z, such that the soft gluon always tends to be inserted between the harder parton and the spectator in the color structure, guaranteeing that we get the correct soft limit. We remark, however that the probability 1 − z is a choice, and that and any function, having the same limits when z → 1 and z → 0, would have been a valid choice. Alternatively, the splitting kernels could have been redefined to only contain the 1 − z singularity (see [33,34] for comparison).
As we want to describe LHC collisions, the color matrix element corrections have to be applied also to initial state hadrons, meaning that we have to deal with initial-initial emissions as well as initial-final and final-initial emissions. The initial state emission cases are treated in a standard backward evolution scheme, meaning that if the emitter is an initial state parton, the backward evolution is done by folding in the parton distribution functions (PDFs), using appropriate splitting kernels, and colors are updated as if the resulting (low energy) parton was emitted. To be more precise, denoting the emitter participating in the hard process by i, the radiated parton with j, and the (initial) parton going into the PDF byĩj, the used splitting kernel is Pĩ j→ij , and the emission probability is evolved using PDF ratios [35], as described in section 6.5 in [2]. The color structure of the full color shower, on the other hand, is, as discussed in section 6, treated as if i →ĩj, j.
In the shower algorithm presented here, we also extend our analysis from [22] by including g → qq splitting. The description of how this splitting fits into the dipole formalism is given in section 6.

Color matrix element corrections
We like to stress that this paper deals with color matrix element corrections to parton showers. We thus correct each emission with the full color correlations, keeping all soft and collinear contributions to the emission, i.e., we use the right "antenna pattern". In this sense, we do more than the standard leading N c showers -where only the leading N c -terms, and the color suppressed term in C F = T R (N 2 c − 1)/N c , are kept -but less than a full matrix element correction.

JHEP11(2018)009
Squaring the amplitude in order to calculate a cross section gives with S n being the scalar product matrix, (S n ) αβ = α n |β n , for color basis vectors |α n and |β n . Upon emission from theĩj,k-pair, the relevant color factor changes to in terms of matrix representations, Tĩ j , Tk ∈ C d n+1 ,dn , of Tĩ j , Tk. While this improves the radiation pattern to fixed order, we should remark that it is not the full story of color suppressed terms in the context of parton showers. To fully include all subleading N c terms in the soft and collinear limits, virtual color rearranging terms associated with the same singularity structure should also be kept. To accomplish this, a full resummation of virtual exchanges is needed. Unfortunately, within the current event generator structure these contributions cannot be included, and we postpone their inclusion for future work. Instead we present a fully functional subleading N c dipole shower, building on the algorithm presented in [22], but including g → qq splitting, hadronization, full mass dependence and initial state hadrons, meaning final-final, initial-final, final-initial and initial-initial dipoles.
The rewriting of the collinear singularities in terms of dipole splitting kernels eqs. (2.6)-(2.7), using color conservation, eq. (2.5), means that the color flow will be replaced by every other color flow, except the flow associated with the color structure of the collinear singularity. In order to be able to continue the full N c matrix element corrected parton shower with a leading N c parton shower (as well as with subsequent hadronization), we do, however, keep the color structure associated with emission from the emitter, as in the standard event record, and we use that for the subsequent leading N c shower and hadronization. Not having one color structure to start from when it comes to hadronization would imply considering a new approach to a hadronization model, which is much beyond the scope of the current paper. On the other hand, it is still interesting to add hadronization in a standard way to enable comparison to data.
For our simulations we use trace bases and the ColorFull implementation [44]. In the trace bases, the color structure is expressed in terms of open and closed quark-lines, of for all possible quark and gluon permutations and all possible number of traces [42,44]. These bases have several advantages. The color structure can trivially be translated into the leading N c color flow, which is needed for subsequent leading N c emission and standard hadronization. The processes of gluon emission, gluon splitting, and gluon exchange can all be very easily described, giving respectively at most two, two and four new basis vectors [39,41,42]. For example, considering gluon emission off an open quark-line we have the two trace basis vectors where we have also illustrated the different color flows in the N c → ∞ limit. 1 If we instead consider gluon splitting, we have One of the disadvantages of trace basis lies in their overcompleteness, meaning that hence -strictly speaking -they are actually not bases, but rather spanning sets.
For a small set of partons, up to approximately five qq-pairs and gluons, the overcompleteness is fairly moderate, but for more partons it becomes significant [49,51]. The number of basis vectors scales as a factorial in the trace basis case, with roughly (N g + N qq )!/e basis vectors [49], whereas the number of basis vectors for finite N c scale only as an exponential. On the other hand, basis vectors which are not needed for a given process can often easily be identified and crossed out for trace bases, in particular, at tree-level, the number of required basis vectors scales approximately as (N g + N qq − 1)!. For example, in eq. The overcompleteness, and the fact that the scalar product matrices are dense, i.e., the scalar product between most pairs of basis vectors does not vanish, is the Achilles' heel of the trace bases. Instead of vanishing for two different basis vectors, the scalar product is suppressed by one or more powers of 1/N c . Thus, when N c → ∞, the bases are orthogonal, JHEP11(2018)009 and the color structure can be replaced by color flows, as in eq. (4.2) and eq. (4.3). Since we keep the full color structure, it is the calculation of scalar products that limits the number of subleading N c emissions that we can keep. At about N g + N qq = 6, the color structure calculations start to take up significant time, and going beyond N g + N qq = 7 is very time consuming within our current setup.
This unfavorable scaling behavior could be circumvented by using orthogonal multiplet bases [49,51], in which case the calculation of the radiation matrices Tĩ j,n would be more time consuming, however not to the same degree [50]. Alternatively color structure could be sampled over. This is a road which we have attempted to pursue. Within our current framework, using the weighted Sudakov veto algorithm from [55] (modified as described in section 5), it is, however, impractical to go beyond 3 subleading N c emissions at LEP, or two subleading N c emissions at LHC, due to the very poor statistical convergence from large weight fluctuations. Therefore, within our current implementation, sampling has proven disadvantageous.

The weighted Sudakov algorithm
The shower described in this paper treats up to N max emissions with the full color correlations. This corrects the emissions to appear first in the p ⊥ -ordered evolution, down to smaller scales, and then the leading color shower handles the subsequent lower p ⊥ emissions. The radiation pattern we aim at describing in the full color shower is where the factor after the multiplication sign is the color matrix element correction, which we will denote by ω ñ ijk . For the cases when the color matrix element correction is negative, the weighted Sudakov veto algorithm from [55] is used, with one modification: in the competition version of the algorithm from [55], the weight for an emission gets a contribution from each veto and accept step, for all trial emissions of all competing pairs. The algorithm will, however, produce the same distribution if the total weight only receives contributions from the accept step of the winning emission and veto steps at scales larger than the winning scale. Discarding all weight contributions below some scale (in this case the scale of the winning emission) makes the algorithm generate another distribution below that scale, but the radiation pattern of the losing trial emissions below the winning scale cannot affect the final distribution, since these emissions are discarded anyway. We therefore choose to discard the weights below the winning scale. This modification significantly improves the convergence. The proof that this modified algorithm generates the same distribution as the algorithm of [55] can be found in appendix A.
Two choices in the algorithm from [55] are: the acceptance probability (denoted in [55]) and the overestimate proposal distribution (denoted R). These free choices were JHEP11(2018)009 used to improve the convergence of the algorithm, the details of the choices can be found in appendix B.
In total, to reproduce the radiation pattern in eq. (5.1), the shower steps given below are repeated until N max emissions have been corrected or no emission is found above the cut-off scale µ.
1. The starting scale Q ⊥ is given by the hard scale. For the first emission, this is taken to be the Z mass for LEP and the average transverse momentum of hard jets in the final state for LHC. For subsequent emissions, Q ⊥ is given by the scale of the previous emission.
2. All processes for all pairs of partons compete with each other and a winning hardest scale is chosen. For each dipole,ĩj,k, candidate emissions,ĩj,k → i, j, k at scales p ⊥,ĩj,k , are chosen according to the Sudakov form factor where P ij,k , in accordance with eq. (5.1), is 3) and z ± (p 2 ⊥ ) follow from the phase space boundaries at fixed transverse momentum. If the color matrix element correction is positive, the standard Sudakov veto algorithm is used (resulting in that the trial emission always contributes a factor 1 to the event weight) and if it is negative, the modified weighted veto algorithm is used (where the weight is, in general, multiplied by the weight in eq. (B.4)). The winning emission defines the details of the kinematics and the recoil is absorbed by the spectatork of the winning dipole, such that all partons are on-shell after the emission.
3. If no scale above the cut-off µ was found, the shower terminates.
4. If this is the emission N max , the leading N c shower will continue showering the event, otherwise the density operator is updated as will be described in section 6.
If N max emissions have been corrected, the leading N c shower continues with the color structure given by the large N c flow associated with emissions from the selected emitters, as discussed in section 3. The leading N c shower then continues until reaching the cutoff scale µ. Finally, the event may, or may not, be hadronized. If the hadronization is performed it starts from the leading N c color flow.

Evolution of the density operator
The parton shower starts from the hard matrix element |M n . However, after emission, the resulting "dipole" color structure from eq. (3.3) cannot, within our framework, be cast

JHEP11(2018)009
into the form of some new amplitude |M n+1 . Instead we see from eq. (3.3) that the relevant n + 1-parton quantity, corresponding to M n ≡ M n M † n for emission fromĩj,k is ∼ Tk ,n M n M † n T † ij,n . For gluon emission (final or initial), keeping all contributions to the emission probability, and using the dipole factorization eq. (2.1) with the splitting kernels eq. (2.6), we see that if we define the matrix element square for n + 1 particles can be written analogously to eq. (3.2) as Thus, when a phase space point has been selected for gluon emission, M n could be updated according to eq. (6.1). Note, however, that the overall normalization of M n+1 is irrelevant, since when used in the n + 2-version of eq. (5.3) to calculate the emission of n + 2 partons, M n+1 enters in both the numerator and denominator. Thus we could ignore any constant factor. In fact, for technical reasons, we only keep the eikonal parts, ∼ p i · p k /(p i · p j p k · p j ), of eq. (6.1). Clearly the dropped hard collinear pieces should not alter the subsequent emission of soft wide-angle radiation.
For the case of g → qq, there is no interference between various possible emitters, and the amplitude is symmetric in all final state gluons, meaning that M n can be updated using only one term where tĩ j,n represents the color space map corresponding to t g qq , i.e., the matrix where element αβ is the transition from basis vector β in the initial (smaller) basis, to basis vector α in the final (larger) basis, where the difference between the color structures is that the gluonĩj has been contracted and replaced by the qq-pair i, j, giving one or two new basis vectors. The possible recoil partners used to set up the gluon splitting into quarks are picked using eq. (5.1) with the g → qq splitting kernel, but with color matrix element corrections as for gluon emission. While this choice is ad-hoc in this case, it has the advantage of nicely fitting into the dipole picture. In the collinear limit, where the splitting becomes relevant, we can use color conservation eq. (2.5), to observe that the sum over the kernels using different spectators is indeed collapsing to the expected collinear splitting function.
Note that the same update of M n+1 , eq. (6.3), and the same recoil strategy, is applied irrespectively of if the splitting gluon is final or initial. The only difference is thus, as for the gluon emission case, the standard convolution with the PDFs for initial parton splitting rates. If we have an initial state quark (antiquark) which evolves backwards into a gluon going into the parton distribution function and an antiquark (quark) which is radiated, the shower has no interference with other diagrams, and the density matrix can be updated according to M n+1 = Tĩ j,n M n T † ij,n , (6.4)

JHEP11(2018)009
without an irrelevant, overall factor. In the standard dipole large N c shower, the momentum recoil is absorbed by the color-connected partner, implying that the emission can be accounted for within the dipole shower formalism. A comment on the update of the density matrix is in order here. In this paper we keep the full color structure, of all possible emitter-spectator pairs which could have contributed to an emission. An alternative strategy is to pairwise sample over emitters and spectators, corresponding to keeping only one term in the double sum in eq. (6.1). While this samples from all color structures, we like to remark that it actually corresponds to a slightly different approximation. In our current implementation, after defining the kinematics of the new emission from the winning pair, all pairs contribute to the color structure of the next emission, and the weight of their color structure, is given by eq. (6.1), implicitly multiplying the same Sudakov factor, the Sudakov factor of the winning pair. In a sampling procedure, the terms in the eq. (6.1) would end up in different events, and each term would be associated with the Sudakov factor of the pair emitting in that winning phase space point. The difference in a sufficiently large Monte Carlo sample, thus lies within the Sudakov factor, coming with different scales for different pairs. We note however, that the two approaches should agree in the soft limit, and we have also checked that the numerical difference is small.

Outline of the simulation
In this section we outline the simulation, and consider various shower distributions with the aim of understanding and validating the shower evolution.
We use the Herwig 7.1 dipole shower, with settings according to the 7.1.3 release [27], with the modified weighted Sudakov veto algorithm outlined in section 5, and with color matrix element corrections as described in section 3, starting from lowest order 2 → 2 processes. The LHC generation p ⊥ -cut is by default put to 30 GeV, and the default jet analysis veto is p ⊥ = 50 GeV. The energy is 13 TeV for LHC and 91.2 GeV for LEP unless stated otherwise. If jet clustering is required, our default choice is the anti-k T algorithm, as provided by the fastjet package [56], with R = 0.4 at LHC. At LEP, no generation cuts are applied, and both at LEP and LHC we use the original Rivet [57] analysis published along with the data in comparison to data.
To ensure statistical convergence, we start with investigating the weight distribution. In figure 1 we show the weight distribution for up to five/three subleading N c emissions at LEP/LHC respectively corresponding to up to six/seven qq-pairs plus gluons in the color bases. We note that although the weight distributions get broader with the number of emissions, it stays sufficiently narrow to ensure convergence for the considered number of subleading emissions. The number of colored partons, which we can practically include, is therefore limited by the evaluation of scalar products in color space, as described in section 4. (We will find, however, that we are able to keep sufficiently many full color emissions for standard hard observables to converge.)  Somewhat against intuition, we see a broader weight distribution for LEP events than for LHC events, despite the fact that we tend to have more colored partons at the LHC. This can be attributed to the fact that the corrections often tend to be negative at LEP (starting from e + e − → qq), due to the negative contribution from coherent emission from the qq-pair. In line with this, we also note that if we separately study qq → qq, qg → qg and gg → gg, we find the largest weight variations for qq → qq, another case where we can expect large negative corrections from qq-pairs.

JHEP11(2018)009
We next turn to the convergence of standard observables with respect to the number of subleading N c corrected emissions. As an example, we consider the rapidities of the hardest LHC jets in figure 2. As can be seen, the curves converge as more subleading emissions are added, and in this respect we see the same pattern for all standard hard LHC observables; they all converge when up to three subleading emissions are added, i.e., starting with a JHEP11(2018)009  2 → 2 topology and adding three subleading emissions (followed by leading N c showering) gives results very similar to adding just two subleading N c emissions (followed by leading N c showering). LEP observables show a similar convergence pattern. Only when explicitly considering very many jets, does the convergence fail. This strongly suggests that for standard hard observables, subleading N c corrections can be well approximated by color correcting the first few emissions.
The convergence can also be underpinned by studying the evolution scale at which the N c = 3 parton shower terminates, and further evolution only is given by the leading N c shower. This is investigated in figure 3 where the transverse momentum, as measured in the frame of the emitting dipole, is shown in orange for the last full color corrected emission, i.e., while keeping up to three subleading emissions at the LHC and up to five subleading emissions at LEP. For comparison we also show the distribution of the first emission not to be corrected.
We find that the typical scale of the last color corrected emission is about 1 GeV at LEP and about 5-10 GeV at LHC, using a 50 GeV cut. Increasing the cut to 1 TeV gives a much harder last subleading N c corrected emission, as expected.
The convergence of observables is also in line with results from a very recent paper on subleading N c corrections at LEP [58], where the authors claim to observe good convergence while keeping subleading corrections down to a variable cut-off at around 3 GeV.
In [58] a Monte Carlo sampling over color is advocated, and the point is made that this avoids the factorial scaling in color space. In view of figure 3, we note, however, that we typically go further down in p ⊥ despite keeping the full color structure, although clearly there is no scale down to which we guarantee that we always keep the N c = 3 color structure. In fact, even if we limit ourselves to four subleading N c emissions, for which the time penalty due to color structure treatment is negligible, we also tend to go well below 3 GeV. Therefore, in the case of hard observables at LEP, color structure sampling seems to be of no benefit for the convergence of observables.
We have also performed a number of standard shower variation checks, including shower scale variation and infrared cutoff variation. While varying the shower scale, we find effects very similar to the leading N c case. Increasing the infrared cutoff from 1 to  Figure 4. Parton level plots for fraction of LEP-events containing n jets with E > 5 GeV (left), thrust (middle) and aplanarity (right), the jets have been clustered with the anti-k ⊥ -type generalized e + e − clustering algorithm with R = 0.7. For hadronized events the effect on thrust vanishes, whereas the effects on number of jets and aplanarity are somewhat reduced. Note that the case of one full N c emission should agree with the leading N c shower, as is seen.
2 GeV likewise results in differences well comparable to the leading N c shower, and so does adding multiple interactions.
Finally, we have checked what happens if the shower is turned off completely after one to three subleading emissions. Here, we find large differences for the LHC, even if we keep up to three subleading N c corrected emissions. We therefore conclude that it is essential to keep showering beyond three color corrected emissions, but it is -for standard hard QCD observables, or likely any observable which is mostly sensitive to the hardest jetsnot important to keep color correcting the subsequent, softer and softer, emissions. For observables depending on soft physics, the situation may be different, as indicated below.

Parton level analyses
We first recapitulate the exercise from [22] and run an e + e − simulation at √ s = 91.2 GeV, without hadronization, but this time including subsequent leading N c showering beyond the (up to) five subleading N c emissions, as well as g → qq splittings. 2 Again we find that the corrections to most LEP observables are small. As examples of observables which show some effect, we show the fraction of events containing n jets with E > 5 GeV, the thrust distribution, and the aplanarity in figure 4, in all cases showing corrections below 10%. On the other hand, effects can be significant in tailored situations. For example, considering the average transverse momentum with respect to a thrust axis defined by the three hardest partons, we find corrections above 10%. We caution, however, that without modification this is not an observable. Nch MC/Data Figure 5. Out-of-plane p ⊥ w.r.t. the thrust and thrust major axes (left), light hemisphere mass (middle) and fraction of events containing N ch charged particles, using to data and Rivet analyses from [59,60].

Hadron level
We have studied hadronized LEP events for a large class of observables from [59,60]. For planarity, sphericity, oblateness, and in-and out-of-plane p ⊥ w.r.t. sphericity axis, we find small differences of a few percent or less, whereas C-parameter shows some effects of 5-10% in the low C-parameter region. As examples, in figure 5, we show the total out-of-plane p ⊥ (w.r.t. the plane defined by the thrust and thrust major axes) and the light hemisphere mass in comparison to data from [59].
In general the deviation of simulated results compared to data from [59][60][61] is clearly dominated by other factors, and the overall description of data does not change visibly. Turning, on the other hand, to observables sensitive to soft physics, we find larger differences. As an example we show the charged multiplicity distribution in figure 5 (right), compared to data from [60], but we likewise see large effects for Durham jet resolution variables in the soft region, on the color singlet cluster masses for the Herwig hadronization model [2] and on individual hadron multiplicities.
While it is tempting to interpret the charged particle multiplicity plot in figure 5, as improved data description, we remark that the hadronization model in use, is the standard Herwig cluster hadronization model [2], with a leading N c color flow, as described in section 3. We therefore caution that the differences should only be seen as an indication of effects on soft physics, from the altered particle kinematics entering the hadronization. While this can be interpreted as a need to retune the full color parton shower, retuning shower parameters is beyond the scope of the present paper.

Parton level analyses
We now turn to the LHC and start with reconsidering figure 2. Considering the leading two jets, with our standard p ⊥ > 50 GeV analysis cut, we see that they tend to have slightly different rapidity distributions, with the second jet being more central. For most other standard observables we find small differences, of a few percent or less. Nevertheless it is illustrative to separately consider scattering involving different partons. Doing so we find that for qq → qq and gg → gg, the subleading N c corrections are small for all studied

JHEP11(2018)009
Leading Nc 3 Nc = 3 emissions observables, whereas for qg → qg, they can be more sizable. In figure 6 we therefore revisit the rapidity distributions of the two leading jets, and find corrections going in opposite directions for the two jets. Since qg-induced scattering contributes with a large fraction of the cross section for the applied cuts, and since LHC data contains qg → qg and gq → gq (along with all other processes), it can be concluded that significant cancellation of subleading N c corrections is present at LHC. Using the cuts of [62], 400 GeV < M 12 < 600 GeV , 3.8 < |y 1 + y 2 | < 5.2 , (7.1) where y 1 and y 2 are the rapidities of the hardest and second hardest jet and M 12 is the invariant mass of the two hardest jets, events dominated by the hard process qg → qg can be statistically enhanced. These cuts select events with one of the two hardest jets being forward (in either direction) and one central. With these cuts we find differences of 5 − 10% for the rapidity distributions of the hardest three jets, as illustrated in figure 7 for the hardest jet. From figure 7, we see that in the N c = 3 shower, the hardest jet tends to be central less often as compared to the leading N c shower. The rapidity distribution of the second hardest jet shows that it is forward less often. There are also 5 −10% differences in ∆φ ij = φ i − φ j , ∆η ij = η i − η j and ∆R ij = ∆φ 2 ij + ∆η 2 ij , for i = 1, 2 , j = 3. As an example ∆φ 13 is also shown in figure 7. In general, with these cuts subleading N c effects show sizable corrections for many standard QCD observables.

Hadron level analyses
We now turn our attention to hadronized events and to comparisons with LHC data. We have compared the subleading N c corrected parton shower to experimental data for a wide range of QCD observables, using data from [63][64][65][66][67][68][69]. First, in figure 8, we consider the event

JHEP11(2018)009
Leading Nc 3 Nc = 3 emissions where p ⊥,i is the transverse momentum of the central jet i, having pseudorapidity η < 1.3, andn T is the direction perpendicular to the beam axis, which maximizes the sum. We also consider central thrust minor, defined again in terms of central jets with η < 1.3, and compare to data from [63]. As can be seen we find relatively small corrections, at the 5% level or below. A comparison to the jet shapes and jet masses for high p ⊥ jets, from [64], shows yet smaller subleading N c corrections, typically below a few percent.
We have also compared data to the so-called color coherence effects from [67]. In figure 9 we show the distribution of the β-angle, defined in terms of the pseudorapidities  Figure 10. Fraction of events having no additional jet with p ⊥ above Q 0 within a rapidity interval |y| < 0.8 (left) and fraction of events where the scalar sum of transverse momenta within |y| < 0.8 does not exceed Q sum (right) for tt events at √ s = 7 TeV. Data and Rivet analysis are taken from [65]. η 2 and η 3 and the azimuthal angles φ 2 and φ 3 , as Experimental data is first compared to shower predictions using a 2 → 2 hard matrix element, and then using a 2 → 3. We clearly see that the use of a 2 → 3 hard matrix element significantly improves the description of data, whereas adding subleading N c showering to the 2 → 2 process changes the distribution compared to the leading N c shower very marginally. This casts doubt upon the description of this observable as probing color coherence, and rather illustrates its dependence on the hard matrix element (or alternatively on other details of the shower algorithm).

JHEP11(2018)009
We have also investigated the effects from subleading color contributions on top-pair production, and compared to data from [65,66]. Here we find that the jet shapes from [66] are essentially unaltered, whereas the measurement of the additional jet activity in tt events [65] show intriguing effects, displayed in figure 10. In particular we note that the data description improves in the region of a modest ratio between Q 0 /Q sum and the hard scale. For gap observables, like in figure 10, with a large scale hierarchy, we remark that we can expect effects from resummation of virtual gluons, which we do not include in this paper. These subleading N c corrections may be sizable [26,[71][72][73].
As a clean test of initial state radiation, we have also compared the N c = 3 shower to event shapes in leptonic Z decay events, [69]. Here, as expected, having fewer colored particles in the hard process, we find very small corrections, at the percent level or below.

Conclusion and outlook
In this paper we have investigated the effect of keeping the full color structure in parton shower emissions in realistic simulations of LHC and LEP events. This is pursued within the dipole shower of the Herwig 7.1 framework [27] as color matrix element corrections.
The N c = 3 color shower corrects the first few (five for LEP and three for LHC) emissions using the full N c = 3 emission pattern. All hard observables we have studied, with the exception of observables explicitly considering very many jets, have converged with respect to keeping additional subleading N c color corrected emissions. The convergence can also be underpinned by noting that the last emission to be corrected, corresponds to a low evolution scale, in comparison to the studied observables, see figure 3. In particular for LEP, we typically go down to evolution scales around one GeV.
Our results show that, for most QCD variables at the LHC, the subleading N c effects are, similarly to at LEP, of the order of a few percent. At the LHC we can, however, see larger effects, 10 − 20%, for the tails of the rapidity distribution of the second hardest jet, figure 2, as well as for more tailored situations, cf. figure 7.
In a gedanken experiment, where quarks and gluons are collided, larger differences can be found, indicating that cancellation of subleading N c effects are present at the LHC. To capture this situation, we consider LHC events while requiring that of the two most energetic jets, one is central and one is forward. In this case we find differences of 5 − 10% also for standard QCD observables. These differences are not only in the tails of the distributions, but over the whole range of several observables.
Turning to soft observables the situation is different. In many cases, including jet resolution variables at low scales, charged particle multiplicities (figure 5), individual hadron multiplicities and the number of very soft jets at LEP, we find large effects, of several 10%. While we cannot expect to make accurate predictions for any of these cases, due to sensitivity to hadronization, multiple interaction and resummation effects, it should be stressed that subleading N c effects can be expected to play an important role for the final state of the shower, entering the hadronization. An immediate extension of this work is therefore retuning of the N c = 3 shower. Indeed, an improved description to (most) observables cannot be expected until retuning is performed.

JHEP11(2018)009
Another natural next step is to include virtual corrections. Virtual gluon exchanges would allow rearrangement of the color structure without emission. This would be expected to have an effect on gap fraction observables, such as in figure 10, and we therefore caution that our conclusions regarding the magnitude of the subleading N c corrections may not be applicable in these cases.
More precisely, virtual gluon exchanges should be the mechanism underlying (perturbative) color reconnection effects. In the longer perspective, it would be desirable to update the standard hadronization model to encompass a subleading N c shower.
For these reasons, we like to stress that this work should be considered as the start of subleading N c corrections at the LHC, not the end. Indeed much work remains to be done.
Note added. While this work has been finalized, a similar approach has been reported in [58]. In [58] color matrix element corrections, as well as the weighted Sudakov algorithm for final state evolution is used, but a sampling of color structure is advocated.

Acknowledgments
We thank Johannes Bellm and Mike Seymour for discussions on color flows in dipole showers and Torbjörn Sjöstrand for comments on the manuscript. Johan Thorén wishes to thank Universität Wien for their hospitality. This work was supported by the Swedish Research Council (contract numbers 2012-02744 and 2016-05996), and in part by the MC-netITN3 H2020 Marie Curie Initial Training Network, contract number 722104, as well as the European Union's Horizon 2020 research and innovation programme (grant agreement No 668679), and the COST action ("Unraveling new physics at the LHC through the precision frontier") No. CA16201.

A Proof of the modified weighted veto algorithm
In this appendix we will recapitulate the weighted Sudakov veto algorithm of [55], and prove that our modified algorithm generates the same probability distribution for competing processes.
The algorithm of [55] generates splittings at the scale q with d additional splitting variables x from a splitting kernel P (q, x) (which can be negative), according to the distribution where Q is the starting scale, µ is the cutoff scale, x µ is a parameter point associated with the cutoff, θ is the Heaviside step function, where we use the convention θ(0) = 1, and ∆ P (q|Q) is the Sudakov form factor,

JHEP11(2018)009
We note that the distribution S P is normalized to unity, in the sense that Q µ dS P (µ, x µ |q , x |Q) = ∆ P (µ|Q) + 1 − ∆ P (µ|Q) = 1, where the integrations are over q and x . As in [55] we use an integrable splitting kernel R(q, x) for which the primitive function can be inverted. This function is required to be positive everywhere, but actually does not need to be an overestimate of P (q, x). We also define an acceptance probability (q, x), such that Now we add the modification to the algorithm, instead of a single weight that is updated upon each trial emission, we keep track of the scale and weight of every veto and accept step. The algorithm starts at a scale q 0 = Q and a step counter i = 0. The algorithm is as follows, for every competing splitting kernel: 1. i is incremented by one.

2.
A trial splitting is generated, at the scale q i and with splitting variables x, from S R (µ, x µ |q i , x|q i−1 ).
3. If q i = µ, the algorithm terminates without an emission, and returns the scales q i and weights ω i for all i, that may have been acquired in the steps below.
4. The trial splitting is accepted with the probability (q i , x). If accepted, the weight is defined and the scales q i and weights ω i for all i, along with the splitting variables x, are returned.

Otherwise, the weight
is stored and the algorithm repeats, starting at step 1.
In the algorithm of [55] the weight associated with this trial splitting is where n steps is the number of veto and accept steps. Our modification will use the weight where q w is the scale of the winning splitting, i.e. we ignore weights below the winning scale.

JHEP11(2018)009
Now we consider the competition algorithm. The standard proof of the competition algorithm does not depend on how the probability distribution, S P , is generated, nor does it depend on the sign of the splitting kernels, hence the standard competition proof applies for the weighted Sudakov veto algorithm. To prove that the weight eq. (A.8) generates the same distribution as the weight eq. (A.7) when using the competition algorithm, we will consider the probability density for a specific scale q and splitting variables x, from the splitting kernel P α , competing with some other splitting kernels P β , β = α. For a splitting kernel P α , the probability density for the algorithm defined by steps 1 to 5 above, starting at scale q 0 and splitting variables x 0 , to traverse a sequence of scales and splitting variables of length n, before emitting at the scale q with the splitting variables x (using the short-hand notation α, where R α is the "nice" splitting kernel corresponding to the splitting kernel P α . The weight associated with this is Pα,Rα,¯ α (µ; q, x|q n , x n | . . . |q 0 , where the upper element of the square brackets correspond to the weight from [55] and the lower element corresponds to the modified weight, and ω i is a function of q i and x i as seen from eq. (A.5) and eq. (A.6). The total probability density for the splitting variables q, x from the splitting kernel P α , competing with the kernels P β , is ∞ n=0 q 1 ,x 1 ,...,qn,xn × ω (n β ) P β ,R β ,¯ β (µ; q (β) , x (β) |q (β) n β , x (β) n β | . . . |q 0 , x 0 )θ(q − q (β) ), (A.11) where the product is over every competing kernel and the θ-function is to enforce that the

JHEP11(2018)009
where the second equality uses eq. (A. 16). We now see that eq. (A.15) and eq. (A.17) are equal, so the modified weights give the correct Sudakov form factor and the distribution in eq. (A.11) will be the same.
To conclude, the physical interpretation of the argument in this appendix, is that emissions below the winning scale do not affect the physics above it. Hence, we are free to use the weights in eq. (A.8), which offer a better convergence for the algorithm.

B Choices for the modified veto algorithm
This appendix contains the details for the choice of acceptance probability and overestimate kernel for the modified weighted Sudakov veto algorithm. The arguments of the splitting kernels, p 2 ⊥ , z, pĩ j and pk, have been suppressed for clarity in the following equations. We use the modified version of the weighted sudakov veto algorithm from [55], as described in section 5. The splitting kernel we aim at generating, P ij,k , is given in eq. (5.3). Defining P Nc→∞ ij,k to be the leading color limit of P ij,k , assumingĩj andk are color connected (cf. eq. (2.8)), and letting R Nc→∞ ij,k be the standard overestimate used in Herwig, as obtained from the ExSample [6] adapted grid proposal, we can write P ij,k = sign(ω ñ ijk )|ω ñ ijk |(1 + δĩ j )P Nc→∞ i.e. we use the same acceptance probability as for the leading N c shower (ifĩj,k were color connected). The acceptance weight is then sign(ω ñ ijk ) = 1 for ω ñ ijk > 0 −1 for ω ñ ijk < 0  Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.