Leading-color fully differential two-loop soft corrections to QCD dipole showers

We compute the next-to-leading order corrections to soft-gluon radiation differentially in the one-emission phase space. We show that their contribution to the evolution of color dipoles can be obtained in a modified subtraction scheme, such that both one- and two-emission terms are amenable to Monte-Carlo integration. The two-loop cusp anomalous dimension is recovered naturally upon integration over the full phase space. We present two independent implementations of the new algorithm in the two event generators Pythia and Sherpa, and we compare the resulting fully differential simulation to the CMW scheme.


I. INTRODUCTION
Experiments at high-energy particle colliders have been integral to unraveling the structure of our universe and have confirmed the validity of the Standard Model of particle physics at an unprecedented accuracy. Going beyond the current level of precision and possibly revealing new fundamental particles and forces will require ever more detailed experimental analyses and theoretical calculations. Monte-Carlo simulations by means of event generators play a vital role in this context, as they link experiment and theory through the detailed description of fully exclusive final states [1,2]. They are required to describe the dynamics of a large number of hadrons originating from QCD Bremsstrahlung, which is modeled in the simulation through so-called parton showers. Modern parton showers are typically based on a unified description of QCD radiative effects in a dipole picture, which encompasses both the leading-order spin-averaged collinear radiation pattern, and the leading-order color-averaged soft radiation pattern. The predictions generated by such algorithms accurately describe many experimental measurements. A notable exception to the success of the parton-shower method arises from its limited phase-space coverage. This problem is alleviated by the matching and merging techniques that allow to correct parton showers to any known fixed-order result at limited final-state multiplicity, and that have been in the focus of interest of the theoretical particle physics community in the past decade [3,4].
Currently the most pressing problem in the context of parton-shower simulations is the lack of options to assess the intrinsic uncertainty of the method itself. The precision of fixed-order perturbative QCD calculations is conventionally quantified by varying the renormalization and factorization scales, and the dependence on these scales is reduced at higher orders in the perturbative expansion if the perturbative series converges. No such technique is currently available for parton showers, essentially because parton showers at higher precision do not yet exist or their practical implementation is incomplete. First steps towards the construction of next-to-leading order (NLO) parton showers have been made in [5][6][7][8][9][10][11][12][13][14], but no method has yet been presented that is capable of simulating fully exclusive final states at hadron colliders. At the same time, it should be expected that a difference exists between parton showers and analytical resummation. While it should be reduced at higher perturbative precision, it cannot be completely eliminated due to differences in the treatment of momentum and probability conservation [15].
In this publication we address one of the most important aspects of next-to-leading order parton showers, namely the simulation of the higher-order corrections to soft gluon radiation, and we show how to implement these corrections in a fully differential form in practice. In integrated form, they lead to the well-known two-loop cusp anomalous dimension [16][17][18][19], which is included in improved leading-order parton showers by means of redefining the strong coupling. This is known as the CMW method [20]. At the differential level, the corrections to soft gluon radiation induce spin correlations and sub-leading color corrections that are not included in leading-order parton showers. As part of the extension to the next-to-leading order, we adapt the algorithm in [21] to include these effects. Moreover, the construction of a local modified subtraction procedure as anticipated in [13] mandates the computation of the two-loop cusp anomalous dimension as an endpoint contribution corresponding to the iterated soft times collinear limit. The resulting algorithm will be a key ingredient in the construction of a fully differential technique for matching parton showers to next-to-next-to leading order calculations.
This paper is organized as follows: Section II present an analytic calculation of the local K-factor due to NLO corrections to soft-gluon radiation. Based on this calculation, Sec. III introduces the modified subtraction method and presents our approach to implementing the required changes in the dipole-like parton shower. Section IV presents a numerical validation of the new Monte-Carlo techniques and an assessment of the effect of the fully differential simulation compared to the CMW method. A summary is given in Sec. V.

II. ANALYTIC COMPUTATION OF DOUBLE-SOFT CORRECTIONS
We employ the formalism for the construction of parton showers at next-to-leading order accuracy originally proposed in [13]. This technique is based on a modified subtraction method combined with a new algorithm for mapping n-particle on-shell momentum configurations to n + 2-particle on-shell momentum configurations and the replacement of explicit symmetry factors by appropriate light-cone momentum fractions that can be identified as "tags" for evolving partons. The extension of this method to soft evolution at next-to-leading order requires the removal of overlap between the explicitly included higher-order corrections in the CMW scheme [20] and the potentially included triplecollinear splitting functions [13]. In this section we will first derive analytic results for the double-soft corrections at next-to-leading order. We define the kinematics in Sec. II A, present the individual corrections in Sec. II B and collect the results in Sec. II C. Based on this calculation, section III introduces a modified subtraction technique and addresses the overlap removal. The leading-order contributions to the soft function, which described the interaction between two hard jets of collinear particles through soft gluon exchange [22][23][24] are shown in Fig. 1. The double solid lines represent the hard legs, and the dashed line indicates the cut. The virtual correction is given by a scaleless integral and vanishes in dimensional regularization [25]. The diagram in Fig. 1(b) and its mirror conjugate generate the eikonal factor Here and in the following we will label the Wilson lines by i and j, while the soft momenta will be denoted by 1 and possibly 2. We also refer to the combined soft momentum as q, where q = p 1 and q = p 1 + p 2 in one-and two-emission configurations, respectively. We restrict our analysis to the improved leading-color approximation typically used in parton-shower simulations. In processes with n possibly color-connected partons, the eikonal term, Eq. (1), is first partial-fractioned [26], and subsequently the color-insertion operator T i T j is approximated by assuming independence of the kinematics. This leads to the replacement As the partial fraction D i,j (q) can be matched to the collinear limit unambiguously, the corresponding color Casimir operator, C i , should indeed be associated with the emission in the soft-collinear limit. This approximation proves to be very accurate in practice. We therefore postpone the exact treatment of the color insertion operators to future work and perform our analysis based on S (0) ij (q). We also point out that including the full next-to-leading order corrections to Eq. (2) requires that the first sub-leading color correction be implemented in the parton shower if the two-loop cusp anomalous dimension is to be recovered in the fully differential calculation. These terms are related to color factors of the form C F − C A /2, where the first contribution is absorbed into the exponentiated leading-order soft result, and the second term becomes part of the genuine two-loop result [27,28]. This will be discussed in detail in Sec. III, and related numerical comparisons will be made in Sec. IV.
The virtual corrections to the single emission have been computed in [29,30]. They are given by (3) Next-to-leading order real-emission contributions to dipole-shower evolution in the soft limit. The double solid lines represent hard (identified) partons i.e. Wilson lines.
The diagrams contributing to the gluonic real-emission corrections are schematically displayed in Fig. 2(a)-(e), while the quark contribution is shown in Fig. 2(f). The vacuum polarization diagrams with gluons have corresponding ghost diagrams, and all terms also occur in the mirror symmetric configuration. Their sum is given by the soft insertion operators computed in [31] S (qq) In the limit of strongly ordered soft emissions, S The full real-emission corrections are obtained by adding the cut vacuum polarization diagrams displayed in Fig. 2(g) and (h), as well as the corresponding terms with the gluons attached to the other Wilson line. They are given by [25] C (qq) To simplify the integration, we define the soft remainder as well as two collinear coefficients 4 The precise meaning of z and φ will be discussed in Sec. III C. In terms of the above functions we can write (1, 2) .

A. Kinematics
We perform the calculation in a scheme that is applicable to both initial-and final-state evolution. We parametrize the final-state momenta using two light-like momenta l and n as The component along l is denoted as p + and the component along n as p − . The reference momenta for the Sudakov decomposition are defined in terms of rescaled hard momenta, where q = p 1 in configurations with one, and q = p 1 + p 2 in configurations with two soft gluons, and where Q 2 = (p i + p j + q) 2 . This implies in particular that 2ln = Q 2 , irrespective of the number of gluons in the final state, and that 0 < α, β < 1 for any of the final-state momenta. We parametrize the integrations over the soft momenta p 1 and p 2 as follows [25] The transverse momentum integrals can be written as and where Ω(n) = 2π n/2 /Γ(n/2) and where we have used the relation p 2 = Q 2 α p β p − p 2 T to perform the integrals over the magnitudes of the transverse momenta. The remaining angular integral has to be carried out differently for different powers of the invariant s 12 = Q 2 (α 1 β 2 + α 2 β 1 − 2 √ α 1 β 1 α 2 β 2 cos φ 12 ) that appears in the expressions of the soft current.
To parametrize the measurement as well as the mapping from four-to three-particle topologies, we introduce the observables B. Contributions at leading and next-to-leading order The leading order momentum space soft function is given by the integral of Eq. (1) To simplify the notation we have defined κ 2 = Q 2 αβ. Next we replace the bare coupling, α 0 s , by the renormalized one in the MS scheme, Thus the leading-order soft function in the dipole shower scheme reads Similarly, the contribution from the virtual corrections, Eq. (3), to the next-to-leading order soft dipole evolution is given by The calculation of the real-emission contributions is straightforward but tedious. We discuss the details in App. A. The contribution from the strong ordering approximation, Eq. (5), reads The contributions from the soft remainder and the collinear terms, Eq. (7), are given by In Sec. III we will devise a modified subtraction method that allows to compute the coefficients of the above functions in four dimensions. The results obtained here are used as a cross-check on the new technique.

C. Complete next-to-leading order corrections
The complete Born-local next-to-leading order corrections to the soft function in the dipole approach are given by the sum of Eqs. (18), (19) and (20). The coupling renormalization, Eq. (16), contributes an additional counterterm We finally obtain the fully differential two-loop momentum space soft function Note that Eq. (22) only depends on α and β through κ 2 , which is a consequence of rescaling invariance in the soft limit [32,33]. The constant Γ (2) cusp is the well known two-loop cusp anomalous dimension [16][17][18][19] and the constant Γ (2) soft is the two-loop soft anomalous dimension computed in [34,35], Using Eq. (33) to expand Eqs. (17) and (22) about the poles in the light-cone momenta q + and q − , defined according to Eq. (9), we obtain and In this context we have defined the functions Note that only the two terms proportional to L 1,n in Eq. (26) contribute to the differential radiation pattern as κ > 0. They correspond to a next-to-leading order K-factor modifying the soft eikonal, such that the soft-gluon emission probability becomes In the CMW scale scheme [20] the Γ cusp contribution is absorbed into the definition of the strong coupling as Upon setting µ R = κ we can further eliminate the explicit β 0 term in Eq. (28) [36]. In this scheme, which is commonly used in parton showers and dipole showers [2], the Monte-Carlo simulation correctly accounts for the effects of nextto-leading order soft QCD corrections at the inclusive level, i.e. integrated over all real-emission configurations. This approximation is valid in principle only for finite κ, whereas in the double-soft limit additional corrections arise from the L 0,0 and L 0,n terms in Eq. (26). However, we will detail in the following that the net effect of implementing two-loop soft corrections fully differentially in the parton shower phase space indeed reduces to generating Eq. (28) at the inclusive level, thereby confirming the findings of [20]. The connection to analytic soft-gluon resummation is established in App. B.

III. IMPLEMENTATION OF THE CALCULATION IN FOUR DIMENSIONS
A general scheme to implement higher-order corrections in parton showers in form of a modified local subtraction method was suggested in [13]. Here we proceed to work out the details of the method in the double soft limit. Regarding the divergence structure of the full double real corrections, this is one of the most demanding regions due to the overlap between various singular configurations, and it can be viewed as a part of the complete solution which will include the simulation of higher-order corrections also in all triple-collinear limits.

A. Modified subtraction method
Our technique is based on the modified subtraction method discussed in [37]. We identify the parton-shower splitting kernels with generalized factorization terms in the MS scheme. These terms can be computed by expanding the differential cross section for a particular final state of interest in terms of plus distributions corresponding to lightcone singularities along the directions of the fast partons. Schematically, for a process with no infrared divergences at the leading order, we can use the next-to-leading order factorization formula [26] for real-emission corrections where In this context, |M ij,k n | 2 are the color-correlated Born matrix elements for the n-particle final state, and dΦ (n) is the corresponding differential phase-space element. TheV ij,k are the dipole insertion operators defined in [26]. They reduce to −T ij T k S (0) ik (j) in the soft limit, cf. Eq. (2). The corresponding one-emission differential phase-space element is given by dΦ ij,k analytically to extract the poles in ε. We will instead perform these integrals in a Monte-Carlo fashion. We first parametrize the emission phase space in the collinear limit s ij → 0 in terms of the virtuality t = s ij and the light-cone momentum fraction z = s ik /(s ik + s jk ) for final state radiation and z = 1 − s jk /s ik for initial-state radiation Note the sign of the exponent of the z ∓ε term, which is negative for emissions from final-state particles and positive for initial-state radiation. The integrand in Eq. (31) can now be expanded in in powers of the dimensional regularization parameter, ε, using the relation which is applied to both the t and the z integral. The 1/ε poles generated in this manner will cancel against the virtual corrections and renormalization terms. This produces a non-locality of the finite remainder which is corrected by the resummation, as the first-order expansion of the parton-shower generates the complementary distribution of the real-emission corrections in phase space [37]. In order to compute the finite remainder, we simply need to compute the O(ε 0 ) terms of Eq. (33) applied to Eq. (31). This can be done fully differentially in the remaining phase-space variables, however we need to take into account that the underlying n-particle phase space and matrix element have an ε dependence that contributes finite terms when combined with the poles from real and virtual corrections. This technique was used in [13] to obtain the matching coefficients for the flavor-changing splitting functions. In the following we will describe how it is implemented in the context of the two-loop soft corrections.

B. Separation of iterated double-collinear endpoints
First, we must account for the fact that there is no equivalent of L n,m in the parton shower. The factorized plus distributions are instead replaced by a double-plus distribution and two related endpoint terms. We define the double plus distribution by its action on a test function . 3. Illustration of the kinematical configurations corresponding to the endpoint contributions L0,n in Eq. (35) in oneparticle (left) and two-particle (right) emissions. Note in particular that all partons in the two-particle configuration are forced to be collinear to the same Wilson line, cf. the explanation in Sec. III B.
Using this relation, we can write Eqs. (25) and (26) as Note that because L 0,n is located at either q − = 0 or q + = 0, the corresponding terms are not included in a standard parton shower. In order to add these contributions we will need to implement endpoint terms in the (iterated) double collinear limit. The relevant kinematical configurations are depicted in Fig. 3. They can be explained as follows: Suppose that a single soft momentum q is emitted off the two light-like partons i and j as in Fig 3(a). As q − = 0 or q + = 0, Eqs. (9) and (10) imply q || p i or q || p j , hence all terms L 0,n are related to single soft gluon radiation in the collinear limit. In the case of double-soft radiation of momenta p 1 and p 2 , depicted in Fig. 3(b), the situation is similar, but slightly more involved. Because the radiated partons are both in the final state, s i1 and s i2 must have the same sign. The limit q − = 0 then implies that s i1 + s i2 = 0, which can only be fulfilled if s i1 = s i2 = 0. Therefore, p i || p 1 and p i || p 2 , which leads to p 1 || p 2 , such that s 12 = 0 and s i12 = 0. The conclusion is that all L 0,n terms correspond to the regions where soft emissions are collinear to one of the Wilson lines. If there are two emissions, they must be collinear to the same Wilson line. The change in color flow generated by the soft radiation then reduces the phase space available to subsequent gluon radiation by a factor proportional to the light-cone momentum fraction of the gluon that is color-adjacent to the anti-collinear Wilson line. Let us assume that the corresponding dipole is spanned by p j and p 1 , then the phase space for subsequent gluon radiation is α 1 Q 2 . As we are interested in the soft gluon limit, α 1 → 0, the remaining phase space is typically close to zero. Any further QCD radiation from the j1-dipole will be suppressed by α 1 , and radiation from the remaining dipoles cannot occur because s i1 = s i2 = s 12 = 0. It follows that the effect of the collinear configurations corresponding to L 0,n is to generate a radiator of reduced invariant mass, oriented along the light-cone directions of the original Wilson lines. The phenomenologically relevant branching probability for such configurations cannot be determined in the double soft limit alone, but requires in addition the computation of endpoint contributions in the triple-collinear limit. We therefore postpone the discussion of these terms to a forthcoming publication.

C. Differential subtraction terms
We will now derive the modified subtraction terms needed for implementing the two-loop soft corrections in the parton-shower. Due to the non-abelian exponentiation theorem [27,28] it is sufficient to consider the gluon splitting functions and include the complete soft eikonal instead of the partial-fractioned term, Eq. (2). However, we must include sub-leading color configurations, corresponding to double soft-gluon radiation off the hard Wilson lines in order to account for coherence effects. The subtraction terms related to the functionsL n,m in Eq. (35) can be defined as S (qq)ct ij where S are given by Eqs. (5) and (7), respectively. In the collinear limit, Eqs. (37) reduce to where z 1 is the light-cone momentum fraction of p 1 in the direction of p 1 + p 2 − s 12 /(s n1 + s n2 )n, with n an auxiliary light-like vector satisfying (p 1 + p 2 )n = 0. The spin-dependent DGLAP splitting functions, P µν (z), are given by The soft gluon current J µ ij is given by the standard expression in the eikonal limit.
Note the minus sign in this expression, which arises from color conservation along the hard Wilson line, i.e. T i = −T j .
In processes with a non-trivial color structure this condition holds only at leading color. We rewrite Eq. (40) such that the transversality of the current becomes manifest: The transverse momentum in Eq. (39) can be parametrized as k µ ⊥ = j µ 12,⊥ (n). We can now prove that φ ij 12 defined in Eq. (7) is indeed an azimuthal angle, as cos φ ij 12 = k ⊥ j ij,⊥ (p 1 + p 2 ), and we can replace z 1 z 2 → s n1 s n2 /(s n1 + s n2 ) 2 . In order to obtain the correct differential radiation pattern in the leading-order simulation, we implement 2 cos 2 φ ij 12 as a correction factor applied to the purely collinear parts of the g → qq and g → gg splitting functions, see Sec. III D for details.
The pure soft terms of Eq. (37) can be rewritten as where The first contribution in Eq. (43) can be interpreted as the eikonal expression for emission of the combined soft-gluon cluster 12 from the hard Wilson lines i and j, and the subsequent radiation of gluon 2 off the leading-color dipoles spanned by i1, j1 or the sub-leading color dipole spanned by ij. The second term describes the same situation with the two gluons interchanged. The last term is a negative contribution arising from the dipole spanned by i and j. This contribution is sub-leading in the global 1/N c expansion, but it contributes at leading color in the double-soft limit and must therefore be included in the parton-shower simulation as the first correction to leading-color evolution. Partial fractioning Eq. (43) following the approach in [26], we find where Equation (45) can be interpreted as the soft enhanced part of the dipole shower splitting function in the limit where partons i, 1 and 2 become triple-collinear, with parton j defining the anti-collinear direction. Note that the in the i1-collinear limit, Eq. (45) develops an integrable singularity that vanishes upon azimuthal integration. This problem will be discussed in Sec. III D. The only remaining two-particle singularity is approached as partons 1 and 2 become collinear. The integrals of Eq. (43) have been computed in Eqs. (19) and (A19). They combine to give Upon defining approximate virtual corrections as we would readily obtain the desired result, Eq. (22), at O(ε). We have verified that the corresponding subtracted real-emission contribution could be computed directly in four dimensions and cross-checked the finite term against the difference between Eqs. (22) and (46). Nevertheless, S ij,A is not a suitable local subtraction term for Monte-Carlo simulation, because the difference to the full real-emission correction contains integrable singularities. In the following section, we will therefore devise a technique to simulate the complete soft subtraction term, Eq. (42), by reweighting the leading-order parton shower.

D. Monte Carlo implementation details
We employ the techniques described in [13,21] to generate the final-state momenta, and we evaluate the splitting functions directly in terms of the kinematic invariants s nm with n, m ∈ {1, 2, i, j}. The kinematics mapping in 2 → 3 branchings is based on [26,38] and is summarized in App. A of [21]. The kinematics mapping and (D-dimensional) phase-space factorization in 2 → 4 splittings was derived in [13], App. A. Note that in both cases we use the Lorentz invariant and numerically stable technique of [13] to construct the transverse components of the momenta.
In order to simulate Eq. (45) in the parton shower, we must correct the leading-order soft radiation pattern. First we need to account for the fact that the eikonal generated by the leading-order parton shower is not identical to s ij /((s i1 + s i2 )(s j1 + s j2 )) if the soft-gluon emission is followed by a subsequent branching of any of the emerging momenta. In the transition ( ı12,) → (ĩ, 12, j) followed by ( 12,ĩ) → (1, 2, i), we obtain instead the following probability for the emission of the final soft cluster 12 The similarity of the kinematics mapping in final-state splittings with final-and initial-state spectator [21,26] implies that Eq. (48) holds for both final-and initial-state Wilson lines, i (corresponding to the ± sign). Note that the term proportional to s 12 in the denominator cannot be neglected in the double-soft limit. We can correct the mismatch between Eq. (48) and the target distribution s ij /((s i1 + s i2 )(s j1 + s j2 )) in Eq. (43) by applying a reweighting factor in the branching of the soft gluon ( 12,ĩ) → (1, 2, i) In the transition ( ı12,) → ( ı1,2, j) followed by ( ı1,2) → (i, 1, 2), with i in the final state, we obtain the following probability for the emission of the final soft cluster 12 If the radiator ı1 is in the initial state, we obtain instead Eq. (48) with ± → −. The weight factor arising from Eq. (50) generates pseudo-singularities in the parton shower phase space, which is undesirable in a Monte-Carlo simulation. We will therefore choose to implement a different strategy in the leading order parton shower. The kinematics in the soft enhanced part of the i1-collinear emission will be chosen according to the identified particle prescription of [26]. Note that due to our definition of the evolution and splitting variable in final-state evolution with final-state spectator [21], the Jacobian factor related to this modification is unity. Eventually, all kinematical correction factors are then given by Eq. (49). Upon including the phase-space correction factors, the collinear terms in the gluon splitting functions implementing the spin correlations present in Eq. (38) read The remaining phase-space effects leading to S (sct) i,j,A are taken into account by multiplying the 1-soft parts of the i1and 12-collinear splitting functions by w 12 ij . Finally, we need to account for the additional strongly ordered term in Eq. (42). This is achieved by means of the identities Using Eq. (48) we can then write where we have defined the weight factor Note that the negative contribution in Eq. (53) does not have a parton-shower correspondence. At the same time, we have so far omitted the squared leading-order contribution arising from Eq. (1). We can correct for both mismatches by adding a subleading color contribution to the i1-collinear terms of the splitting function of the Wilson lines. This term reads The weight factor of theC ij term in Eq. (55) was derived by considering its diagrammatic representation, which arises from the Abelian parts of Figs. 2(a) and (b) [25]. We may also consider the result of the integration in Sec. II and its Fourier transform in impact parameter space, cf. App. B. In fact, Eqs. (A10) and (A18) generate the exact same result up to O(1) as the square of the leading-order term, Eq. (B4), hence proving that Eq. (55) is a valid form in the double-soft region that will allow us to reproduce the squared leading-order term at the integrated level. In our numerical implementation we include Eq. (55) in the i1-collinear sector with spectator 2. This means that we mis-identify in principle the related evolution variable, which should be κ 2 in the notation of [21], and hence proportional to s i1 s 1j instead of s i1 s 12 . We correct for this effect by reweighting with a ratio of strong couplings, taken at the current vs. the correct evolution variable, and by setting Eq. (55) to zero as the evolution variable falls below the parton-shower cutoff. A second sub-leading color contribution is given by the difference It accounts for the fact that the second soft emission off the Wilson lines occurs with the color charge C A due to the interference with a color octet. We can now define the combined sub-leading color contribution to the parton-shower evolution as Note that P (slc) ij vanishes in the i1-collinear limit, such that the correct color factor is recovered in collinear evolution. The remaining parts of the improved leading-order, fully differential splitting functions related to the 12-collinear, 1-soft final-state singularities are given by the leading-color expressions Note that we omitted the collinear parts of the splitting functions related to the Wilson lines i and j, as these are unchanged by the double-soft corrections. The notation is such that the subscripts indicate the partons which are color-adjacent to the splitting products, while the superscript indicates one that is color-adjacent to the adjacent parton. In particular, the color connection in (P qq ) k 2 (1, i) would be 2 ↔ 1 ↔ i ↔ k. Note that the ordering of the arguments and lower indices is important. This is apparent in the 12-collinear limit, where the gluon-to-gluon kernel receives a second contribution, (P gg ) l ji (2, 1), which is related to a different evolution variable in the dipole shower approach because it corresponds to the 12-collinear, 2-soft singularity [21]. The leading-order g → gg splitting function in the collinear limit is recovered only upon adding (P gg ) k ij (1, 2) and (P gg ) l ji (2, 1). The pure collinear term in (P qq ) k i (1, 2) could in principle be modified by the weight, Eq. (49), but there is no indication, based on the double-soft limit, as to whether this would constitute an improvement of the parton shower or not. We postpone the analysis of this term to a future publication.
We emphasize that, despite the reweighting of the leading-order parton shower to the full real-emission pattern of the double-soft limit, a hard correction remains to be computed using the techniques of [13]. This correction arises because the leading-order parton shower does not fill the complete two-emission phase space, see for example [39]. The correction is given bỹ where t c is the infrared cutoff of the parton shower. The two terms in the Θ-function correspond to the ordering variables in the first and second emission, respectively. To simplify the notation we have defined S i,12 , which is given as s i12 in the case of final-state Wilson lines and 2p i (p 1 + p 2 ) in the case of initial-state Wilson lines. Correspondingly, S ij,q is given by Q 2 for two final-state Wilson lines, 2p i p j for two initial-state Wilson lines, and 2p i (p j + q) if i is in the initial, and j is in the final state [21]. Following the discussion in Sec. III A, the O(1) remainder in Eqs. (35) is implemented as an endpoint contribution at s 12 = 0 and s i1 = 0, s i2 = 0, s j1 = 0, s j2 = 0 for all κ > 0. This allows us to simulate the radiation pattern fully differentially at the next-to-leading order. The related endpoint terms are given bỹ gg,ij (1, 2) = δ(s 12 ) wl,ij (1, 2) = − δ(s i1 )  [pb]  [pb] Pythia vs Sherpa   FIG. 4. The effect of the phase-space weights w 12 ij , andw 12 ij , defined in Eqs. (49) and (54), on the leading-order parton-shower evolution, limited to two emissions. We show the differential jet rate y34 in the Durham algorithm [48] as a proxy for the rate change, and the angle defined in Eq. (7) as a proxy for the impact on differential distributions. The process considered is e + e − →hadrons at LEP I energies.
The factor 1/2 inS (cusp) gq,ij removes the double counting of soft-collinear regions when swapping the role of i and j. It would in principle be desirable to work with partial fractions of the eikonals s ij /(s i1 s j1 ) and s ij /(s i2 s j2 ). However, these partial fractions cannot be defined unambiguously in the exact limits s i1 → 0, s j1 → 0, s i2 → 0, s j2 → 0. One possible solution would be to introduce an additional rapidity regulator, similar to [40] or [41,42]. We leave the investigation of this possibility to future work. We implement the contributions proportional to the beta function as a double endpoint which contributes an additional term to the soft enhanced parts of the leading order splitting functions. 1

IV. NUMERICAL RESULTS
In this section we present numerical cross-checks of our algorithm, and we compare the magnitude of the corrections generated by the double-soft splitting functions to the leading-order parton shower result in the CMW scheme [20]. We restrict the analysis to pure final-state evolution, but we stress that the formulae relevant to initial-state evolution have also been presented in Sec. III D. We have implemented our algorithm into the DIRE parton showers, which implies two entirely independent realizations within the general purpose event generation frameworks PYTHIA [43,44] and SHERPA [45,46] that are cross-checked point by point and in the full simulation at high statistical precision. We ij,A compared to uncorrelated parton-shower evolution, including the phase-space suppression investigated in Fig. 4. All simulations are limited to two emissions. The left and middle panels show the impact on the pure z(1 − z) contribution (top panels) and on the complete g → gg and g → qq splitting function (bottom panels). In both cases the production of the gluon is described by the eikonal part of the q → qg splitting function only. The right panel shows the effect of spin correlations on the complete two-emission pattern. In order for the results to be as similar as possible, the weight w 12 ij from Eq. (49) is included. The process considered is e + e − →hadrons at LEP I energies.
use the strong coupling according to the CT10nlo PDF set [47]. The process under investigation is e + e − →hadrons at LEP I energy (91.2 GeV). We choose to exemplify the effects of the double-soft corrections using the k T jet rates y 23 and y 34 in the Durham algorithm [48] and the angle α 34 between the two softest jets [49]. Figure 4 shows the impact of the phase-space weights, w 12 ij andw 12 ij , defined in Eqs. (49) and (54). These weights generate a strong suppression of the radiation probability. The effect is eventually compensated by other corrections (see Fig 8), such that a fairly good agreement with the leading-order approximation is obtained. The lower panels in Fig. 4 show a comparison between the results from Pythia against those from Sherpa. The two predictions agree up to statistical fluctuations, providing a strong cross-check on the consistency of our implementation. Figure 5 shows the impact of the spin correlations implemented by the cos 2 φ ij 12 dependence of Eq. (7) compared to a spin averaged simulation. While the related effects are striking when investigating the z(1 − z)-dependent parts of the splitting functions in isolation, they are greatly diminished in the complete calculation. Figure 6 displays the impact of the generic sub-leading color corrections in Eq. (57). The effects are generally smaller than expected based on a naive estimate (i.e. O(1/2N c )), because Eq. (57) is suppressed in the collinear region, cf. the discussion in Sec. III D. Figure 7 shows the impact of the subtracted real-emission corrections, Eq. (59), and the endpoint terms, Eq. (60), on the radiation pattern in q → q(gg), q → q(q q ) and q → q(qq) splittings, where the particles in parentheses are the soft emissions. We have verified that exact agreement between our implementations is obtained also in the case of g → g(gg) and g → g(qq). The numerical impact of these corrections is similar to the quark-induced case. Note that the 3 → 4 jet rates receive corrections from the subtracted real emission only, while the 2 → 3 jet rates are impacted by both the subtracted real-emission and the endpoint terms. Figure 8 compares the results from a leading-order simulation according to [21] to our complete next-to-leading order prediction. In the leading-order case we present the calculation with and without the CMW scheme [20]. We observe that the CMW prediction matches the rates of the full next-to-leading order result for the Durham jet rate y 23 fairly well in the intermediate-y region, but there are some discrepancies in the low and high-y region. In addition, there is a considerable rate change in y 34 . The angular observable α 34 shows deviations between the CMW prediction  and the full next-to-leading order result at very small and very large values. In all cases, the scale uncertainty is greatly reduced at the next-to-leading order, and the next-to-leading order predictions lie within the leading order uncertainty bands. Note in particular that our result presents the first genuine estimate of the perturbative uncertainty in a parton-shower simulation. Some earlier attempts, despite generating variations of the same order, treated the problem in an approximate manner [50]. Other techniques [14,51] assumed the scale variations on collinear parts of the splitting functions to be identical to the soft parts, and therefore generate artificially small uncertainty bands, which may not reflect the true perturbative precision at the order to which the computation is performed.

V. CONCLUSIONS
We have presented a calculation of the next-to-leading order corrections to soft-gluon radiation, differentially in the one-emission phase space. This is a crucial ingredient in the construction of a next-to-leading order parton shower. We have demonstrated, for the first time, that the soft next-to-leading order contribution to the evolution of color dipoles can be obtained in a modified subtraction scheme, such that both one-and two-emission terms are amenable to Monte-Carlo integration in four dimensions. The two-loop cusp anomalous dimension emerges naturally in this method. We observe fair agreement between the results of the fully differential simulation and the approximate treatment using the CMW scheme, where the two-loop cusp anomalous dimension is included in an inclusive manner. The similarity of the results is reassuring, because the individual higher-order contributions have kinematical dependencies that can differ strongly from the iterated leading-order result. Our calculation can be seen as a confirmation that the existing leading-order parton showers developed over the past decades have been amended by the dominant effects arising from the higher-order soft corrections, but it also confirms that the higher-order corrections do have an impact beyond a simple K-factor. We are now in place to compute these effects without the need for approximations, and to include them in phenomenological studies as well as experimental analyses at the particle level. This allows in particular to obtain meaningful estimates of the renormalization scale uncertainty.  59), and endpoint terms, Eq. (60) on the radiation pattern in e + e − →hadrons at LEP I energies. We show the contributions from q → qgg (left) and q → qq q (right) to the differential 2 → 3 (red) and 3 → 4 (right) jet rates in the Durham algorithm.  8. Scale variations in the leading-order and next-to-leading order (soft) parton shower simulation of e + e − →hadrons at LEP I energies at parton level. We compare to both the plain leading-order predictions (green) and the result in the CMW scheme (red).

Strong ordering approximation
The real corrections in the strong ordering approximation, Eq. (5), lead to one non-trivial integral over transverse momenta, which is given by Here we have defined the intermediate variables χ = (1 + cos φ)/2 and K = 2 √ α 1 β 1 α 2 β 2 /(α 1 β 2 + α 2 β 1 ). The last line was introduced in [52] 2 . It is obtained by the following three transformations of the hypergeometric function ij,A (q) = Here we have split the integration range in applying the second transformation in Eq. (A2) in order to guarantee that that the argument of the hypergeometric function is between zero and one in the entire range of integration. Light-cone momentum conservation takes the form α = α 1 + α 2 and β = β 1 + β 2 . This suggest the parametrization α 2 = α s and β 1 = β t, with s and t ranging from zero to one. We remap the integration variables to the unit hypercube using the additional transformations in the first and second integral, respectively. We obtain where we have defined κ 2 = Q 2 αβ. Changing variables to x = uv and y = (1 − v)/(1 − uv) we can write We first perform the y-integration and obtain (A7)