Numerical resummation at subleading color in the strongly ordered soft gluon limit

We present a Monte Carlo approach to soft-gluon resummation at subleading color which can be used to improve existing parton shower algorithms. At the single-emission level, soft-collinear enhancements of the splitting functions are explicitly linked to quadratic Casimir operators, while wide angle single-soft enhancements are connected to nontrivial color correlators. We focus on a numerically stable implementation of color matrix element corrections to all orders and approximate the virtual corrections by requiring unitarity at the single-emission level. We provide a proof-of-concept implementation to compute nonglobal event shapes at lepton colliders.

We present a Monte Carlo approach to soft-gluon resummation at subleading color which can be used to improve existing parton shower algorithms. At the single-emission level, soft-collinear enhancements of the splitting functions are explicitly linked to quadratic Casimir operators, while wide angle single-soft enhancements are connected to nontrivial color correlators. We focus on a numerically stable implementation of color matrix element corrections to all orders and approximate the virtual corrections by requiring unitarity at the single-emission level. We provide a proof-ofconcept implementation to compute nonglobal event shapes at lepton colliders.

I. INTRODUCTION
Soft-gluon resummation is one of the most important tools in perturbative QCD, as it allows one to systematically and fairly straightforwardly compute radiative corrections to all orders for a large class of observables [1]. The effect of gluon radiation is typically computed for single or multiple emissions, and recoil effects are approximated at the same level. 1 If the observable is simple, all-order corrections can be obtained by exponentiating these results, and the remaining obstacle is posed by color coherence, which may lead to a soft anomalous dimension in matrix form. A general framework for resumming event shapes based on this concept was developed at next-to-leading logarithmic (NLL) accuracy [3][4][5][6] and at next-to-next-to-leading logarithmic accuracy [7,8].
nonglobal observables require a more sophisticated treatment, which was first discussed in the context of e + e − and deep inelastic scattering (DIS) event shape resummation [9,10] using a Monte Carlo (MC) approach at leading color accuracy [9]. A corresponding evolution equation was derived [11], which enabled the inclusion of subleading color effects [12,13]. Numerical results have subsequently been computed for example for interjet energy flows [14] and for the hemisphere mass distribution in e + e − →hadrons [15]. nonglobal observables have also been investigated using methods of effective field theories [16][17][18].
nonglobal logarithms are particularly important in the context of interjet radiation and in the presence of a jet veto [19,20]. In the context of Large Hadron Collider phenomenology they have therefore received considerable attention [21][22][23][24][25]. Several approaches have been suggested for their numerical resummation, ranging from subleading color parton showers [26][27][28][29][30][31] to evolution at the amplitude level [32][33][34][35][36]. While color-corrected parton showers can exhibit good numerical convergence, the evaluation of the color matrix elements becomes prohibitively expensive at high parton multiplicity, and therefore the approach cannot be used beyond a very limited number of emissions. Amplitudelevel evolution on the other hand will typically suffer from a slow rate of convergence in the Monte Carlo simulation. In order to address the problem of slow convergence, we propose a novel algorithm. Using color conservation, the squared soft-gluon current is rearranged into a soft-collinear contribution proportional to the quadratic Casimir operator, and a collinearly suppressed correction term proportional to the color correlators. Based on the independence of color and kinematics operators, the color matrix elements are integrated with Monte Carlo methods at each step of the evolution. This allows one to reach good precision on the color coefficients, while limiting the run-time of computer simulations at high multiplicity. Appendix A shows how color coherence emerges in our approach.
This manuscript is organized as follows: Section II discusses the resummation formalism, and Sec. III introduces the phase-space mapping needed to implement it away from the exact soft limit. Section IV presents our Monte Carlo technique to compute the color matrix elements. The difference between subleading color and improved leading color evolution used in standard parton showers is analyzed in Sec. V by studying the light jet mass and narrow jet broadening distributions in e + e − →hadrons. Section VI contains an outlook.

II. RESUMMATION FORMALISM
Soft-gluon resummation is typically performed for a given, fixed number of hard partons, generated at scales that are widely separated from the scale of additional soft radiation. These partons are assumed to be unchanged after the emission of a soft gluon, leading to the notion of Wilson lines and eventually the exponentiation of the soft anomalous dimension matrix. We will adopt a different approach, based on the physical picture in the strongly ordered soft limit. By the very definition of strong ordering, each radiated soft gluon must be treated as a new Wilson line for subsequent gluon emissions.
We denote the Born matrix element for n partons by |M n and the color insertion operator for parton i as T i [1]. The approximate n + 1-parton squared matrix element in the soft limit then reads where we have defined the squared n-parton soft current The invariants s ij are defined in terms of the parton momenta, p i , as s ij = 2p i p j . Since |m n+1 is an n + 1-parton matrix element, Γ is defined in a higher dimensional color space than Γ n . Equation (1) generalizes to k + 1 emissions as Note that the complexity of Γ n increases rapidly with the number of partons, such that the evaluation of color factors encoded in Eq. (3) becomes increasingly cumbersome. The evolution of the parton ensemble is governed by the differential branching probability where dΦ +1 = d 4 q δ(q 2 )/(2π) 3 is the four-dimensional differential phase-space element for the emission of the gluon with momentum q. We parametrize this phase space as where κ 2 ij is the evolution variable of the parton shower,z i is the splitting variable, φ ij is an azimuthal angle, and J(κ 2 ij ,z i , φ ij,m ) is a Jacobian factor. For details, see Sec. III. Equation (4) describes resolved real-emission corrections. A standard choice for parton-shower algorithms is to define a no-branching probability, Π(κ 2 ), such that virtual and unresolved real-emission corrections are defined in terms of the resolved real-emission corrections by means of unitarity. This approach should be improved by accounting for the different color structure in virtual corrections [30,32,33], as well as for Coulomb phases [32,33,37]. It was pointed out in [38] that one can therefore not yet claim a fully color correct evolution. We will postpone these problems to a future publication. Instead we focus our attention on a suitable rearrangement of color and kinematics factors in the real components, in order to improve the numerical convergence of the simulation. While this is not sufficient for arbitrarily complicated observables, it constitutes an important step toward a complete full-color resummation algorithm.
We define the no-branching probability such as to restore unitarity This equation has the solution where Sketch of the kinematics mapping described in Sec. III. The emitting parton is i, and the reference momentum for definition of the azimuthal angle is j. The emitted gluon carries momentum q. The forward and backward light-cone momenta are given by l and n.
The squared n-parton soft current, Eq.
(2) has a form which is not particularly suitable for implementation in numerical simulations. We use the partial fractioning approach of [39] to rearrange it as where we have defined the splitting operator Note that in general T i ΓT j will not equal T j ΓT i , hence we cannot combine the two terms on the right-hand side of Eq. (9). We use color conservation to rewrite them as Combining the second and the last terms in parentheses, we obtain 2 where we have defined the splitting operator [40] Note thatP tends to zero in the iq-collinear limit. Equation (12) should therefore be viewed as a rearrangement of Eq. (9), where the collinearly enhanced terms are made explicit, and the remainder is singly soft enhanced only. While additional rearrangements would allow one to achieve a further kinematical suppression by combining multiple operators asP i jk +P j il +P k li +P l kj , such rearrangements will produce additional terms proportional to T i Γ T i . We find Eq. (12) to be the most suitable form for a Monte Carlo implementation. Examples for its relation to analytically known soft insertion operators are given in Appendix B.

III. KINEMATICS MAPPING
In order to implement Eq. (12) in a numerical simulation, the operatorP i jk must be well defined. When evaluating the difference between P i k and P i j , we assume that either the underlying Born configurations are identical in both terms, or that their difference gives rise to subleading power corrections. Since the latter may be difficult to prove in the general case, we use a kinematics mapping, which ensures that the underlying Born state is the same for identical i and q. Such a mapping is defined, for example in [41][42][43][44][45][46], and is schematically depicted in Fig. 1. Here and in the following, a tilde denotes momenta before the emission of the soft gluon.
We define the variables in Eq. (5) as follows where the lightlike reference vector n µ is defined in Eq. (15). The lightlike auxiliary vector m µ is given by m µ = K − K 2 /(2 Kn) n and defines the reference axis for φ ij in the transverse plane. Note that the inverse mapping leads to the same n + k-parton momentum configuration for any choice of p µ j and m µ , as long as p µ i and q are identical. This is an important feature needed for the rearrangement of color insertion operators in Sec. II.
The longitudinal recoil generated in the emission of the soft gluon q is absorbed by all partons except the emitter. Due to our choice of evolution and splitting variables, the mapping depends nontrivially on the azimuthal angle. In order to construct the momenta, we first define the lightlike vectors The rescaled invariant γ is given by We can now parametrize the momenta as and Note that x and K T are invariant under the mapping, as the momenta p j and K are completely determined by a boost ofp j andK along the direction of n. The variables x and K 2 T can therefore be computed using the Born kinematics. Solving Eqs. (14) for z and k 2 T then yields The phase-space boundaries are given by The Jacobian introduced in Eq. (5) is given by The fact that the inverse of this mapping yields the same underlying Born kinematics for all configurations where the emitting particle is parton i allows one to rearrange the soft anomalous dimension matrix Γ into Eq. (12) without the need for an additional reweighting away from the exact soft limit.  , 0). In case (a), the newly generated color c differs from ci, and the gluon is assigned the color state (ci, c). In case (b), the newly generated color c differs from ci and the gluon is assigned the color state (c, c). In case (c), the newly generated color is the same as the original color, ci.
The prefactors indicate the weight of each diagram, which originates in the Fierz identity T a ik T a mn = 1/2(δinδ mk − δ ik δmn/Nc) and appears squared in the sampling algorithm. See the main text for details.

IV. COLOR ALGEBRA
The color insertion operators T i . . . T j in Eq. (12) are computed using Monte Carlo summation in the color-flow basis. In the following, we denote the emission of a gluon off parton i by the color branching (c i ,c i ) → (c i ,c i )(c g ,c g ), and the absorption on parton j by the color recombination (c j ,c j )(c g ,c g ) → (c j ,c j ). We choose to sample the color configuration in the emission according to the quadratic Casimir operator, T 2 i . In the case of gluon emission off a quark, we have T 2 i = T a ik T a kj = C F δ ij . In order to allow for a fully differential sampling of the gluon color index, we insert a partition of unity in the form δ ab = 2Tr(T a T b ), leading to 2(T a ik T a mn )(T b nm T b kj ). The first parentheses are associated with the emission of the gluon, and the second with its absorption. With the help of the Fierz identity, T a ik T a mn = 1/2(δ in δ mk − δ ik δ mn /N c ), we identify three distinct color topologies, which are depicted in Fig. 2. Diagram (a) squared corresponds to contributions from δ in δ mk only, and carries a weight of 1/2. For each given quark color (c i , 0), there are N c − 1 such configurations, leading to a relative weight of (N c − 1)/(2C F ). Diagram (b) squared corresponds to contributions from δ ik δ mn only, and carries a weight of 1/(2N 2 c ). For each given quark color (c i , 0), there are N c − 1 such configurations, leading to a relative weight of (N c − 1)(1/N 2 c )/(2C F ). Diagrams (c) squared correspond to contributions from both δ in δ mk and δ ik δ mn , and carry a weight of (1 − 1/N c ) 2 /2. For each given quark color (c i , 0), there is exactly one such configuration, leading to a relative weight of (1 − 1/N c ) 2 /(2C F ).
Based on these considerations, and a similar decomposition of the quadratic Casimir operator for gluons, we are eventually led to the following color sampling algorithm 1. If the emitter is a quark or an antiquark, assign a weight C F to the emission The complete operator T i . . . T j is restored by sampling over all possible recombinations of the intermediate gluon upon insertion of T j . The recombination algorithm proceeds as follows 1. If the absorber is a quark or antiquark (a) If j is a quark and c j =c g orc g = c g , set the merged color to (c g , 0), else assign weight zero (b) If j is an antiquark andc j = c g or c g =c g , set the merged color to (0,c g ), else assign weight zero (a) If c g =c j andc g = c j , assign weight 2 and set the merged color randomly to either (c j ,c g ) or (c g ,c j ) (b) Else if c g =c j /c g = c j , set the merged color to (c j ,c g ) / (c g ,c j ), else assign weight zero

If the absorber is a gluon
Note that arbitrarily many insertions may happen before the gluon emitted by T i is annihilated via T j , as required by Eq. (12). The correctness of the above algorithm follows directly from the decomposition of the generators and the structure constants of SU (N c ) in the color-flow basis [47]. In the context of numerical resummation it is important to note that the color matrix elements in Eq. (8) can be evaluated as a Monte Carlo integral with more than one point per event. This can be used in practice to improve the convergence of the overall simulation. We validate the above algorithm numerically by computing the color coefficients for gluon webs within a quark loop. They can be systematically reduced to maximally non-abelian coefficients, which are related to the quadratic, quartic and higher-point Casimir operators. This leaves a small number of nontrivial intermediate gluon web configurations,   [48] in e + e − →hadrons at √ s = 91.2 GeV. We compare predictions from two independent implementations of our algorithm, labeled "Code 1" and "Code 2". The infrared cutoff is set to which need to be evaluated. Table I lists some of these configurations up to four gluon insertions and compares the analytic results to Monte Carlo predictions from our algorithm at high statistical accuracy. In this context we define F a bc = if abc , where f abc are the SU(3) structure constants. Note in particular, that the fourth coefficient in the table is related to the quartic gluon Casimir operator, leading to an additional contribution of 2/N 2 c which arises from double singlet gluon exchange between two gluons.

V. NUMERICAL RESULTS
In this section we apply our new algorithm to e + e − → jets at the Z pole, √ s = 91.2 GeV. We use a two-loop running coupling defined by α s (Q 2 ) = 0.118, and the quark mass thresholds m c = 1.3 GeV and m b = 4.75 GeV. We crosscheck our predictions using two entirely independent Monte Carlo implementations based on [49], which was validated independently against [50] at high precision. Our simulations do not include collinear contributions to the splitting functions, and they are carried out at the parton level. They can therefore not be compared directly to experimental data. However, they serve as a first proof-of-concept that color matrix element corrections at arbitrary multiplicity can be performed in a numerically stable fashion that enables their application to relevant physics problems in current or past collider experiments. We follow the approach in [51] and terminate the subleading color evolution at a scale t c,F C that is insignificantly larger than the typical parton-shower infrared cutoff of √ t c ∼ 1 GeV. All distributions presented here are generated with √ t c,F C = 3 GeV. We claim that this is not a problem for practical applications, since hadronization effects typically influence numerical predictions up to a scale of the order of the b-quark mass, and the details of the fragmentation model have a much larger impact on measurable distributions in this range than the details of the parton shower. In order to provide a smooth transition to improved leading-color evolution below t c,F C , we choose a leading color state according to the probability for a leading-color matrix element to have produced the partonic final state at scale t c,F C . This is similar to how leading-color configurations are chosen in matching and merging techniques [52,53]. Figure 3 shows a comparison of predictions for the Durham k T -jet rates [48] in e + e − →hadrons at √ s = 91.2 GeV.
Our two independent numerical implementations of the resummation are statistically compatible and show good convergence, even in regions of large k T and for higher jet multiplicity. Figure 4 displays numerical predictions at √ t c = 1 GeV for the Durham 2 → 3 and 3 → 4 jet scales, and for two nonglobal shape observables, the narrow jet broadening, B N , and the light jet mass, M L . We find that the impact of subleading color evolution on all these observables is less than 10%, which agrees with the intuitive notion that corrections to improved leading color evolution should be of order 1/N 2 c . This can be taken as a strong indication that the typically excellent agreement of modern parton-shower predictions with measured nonglobal shape observables is not entirely accidental. A variant of the Durham n → (n + 1) jet scales has been resummed recently and matched to next-to-leading order (NLO), to achieve NLO + NLL accuracy in [54], including a quantification of prospective subleading color contributions. Our results are compatible with the smallness of the effects observed there, in particular when noting that the results in [54] are matched to a fixed order NLO calculation relative to the Born process while we present pure parton shower results here.

VI. CONCLUSIONS
We have presented a novel Monte Carlo method for soft-gluon resummation that allows one to generate parton-level events and can be incorporated into existing parton showers in order to improve their formal precision. Along with this manuscript, we provide a proof-of-concept implementation that can be used for numerical studies in e + e − →hadrons. We find that the impact of subleading color evolution on the Durham k T -jet scales, on narrow jet broadening, and on the light jet mass agrees with the naive expectation that corrections to existing parton-shower approaches should be suppressed by O(1/N 2 c ). It will be interesting to investigate the impact on dedicated observables, which probe nontrivial color correlations, as for example in [55]. In this work we have neglected the matrix structure of soft virtual corrections. A future study will address the feasibility of including these contributions as well.

ACKNOWLEDGMENTS
We thank Joshua Isaacson for numerous fruitful discussions on the stochastic sampling of color configurations. We are grateful to Simone Marzani for his comments on the manuscript, and for many stimulating discussions on soft-gluon resummation. DR thanks Steffen Schumann for his support.
We are particularly grateful to Dave Soper and Zoltan Nagy, as well as Jeff Forshaw, Simon Plätzer and Jack Holguin for pointing out that our original claim to achieve resummation at full color accuracy does not hold due to the unitarity argument in Eq. (7)  A general framework for resumming a large class of observables has been developed in the context of the program Caesar [3][4][5][6]. This formalism, in particular including the computation of color structures for in principle arbitrary multiplicities, was automated in [56] and recently applied to resummed calculations in e + e − → jets [54,57]. In this formalism, color coherence is an integral part in treating multiple real emissions. We therefore find it relevant to test how this is reproduced in our algorithm.

Four radiators
We choose T 2 1 , . . . , T 2 4 , T 1 T 4 and T 1 T 3 to be the independent elements of the color algebra. The remaining insertion operators can be expressed in terms of these operators as Based on Eq. (12), the complete soft insertion operator for the four parton case with 1,2 and 3,4 being the same type of parton (either quark or gluon), then reads = C 1 w 12 + C 3 w 34 + T 1 T 3 w 12 + w 34 − w 13 − w 24 + T 1 T 4 w 12 + w 34 − w 14 − w 23 . (B4)