Momentum conservation and unitarity in parton showers and NLL resummation

We present a systematic study of differences between NLL resummation and parton showers. We first construct a Markovian Monte-Carlo algorithm for resummation of additive observables in electron-positron annihilation. Approximations intrinsic to the pure NLL result are then removed, in order to obtain a traditional, momentum and probability conserving parton shower based on the coherent branching formalism. The impact of each approximation is studied, and an overall comparison is made between the parton shower and pure NLL resummation. Differences compared to modern parton-shower algorithms formulated in terms of color dipoles are analyzed.

We present a systematic study of differences between NLL resummation and parton showers. We first construct a Markovian Monte-Carlo algorithm for resummation of additive observables in electron-positron annihilation. Approximations intrinsic to the pure NLL result are then removed, in order to obtain a traditional, momentum and probability conserving parton shower based on the coherent branching formalism. The impact of each approximation is studied, and an overall comparison is made between the parton shower and pure NLL resummation. Differences compared to modern parton-shower algorithms formulated in terms of color dipoles are analyzed.

I. INTRODUCTION
Searches for new physics and measurements of Standard Model parameters at the Large Hadron Collider and possible future colliders require ever increasing precision in the analysis of multi-scale events. Large scale hierarchies in such reactions will generally result in large Bremsstrahlung effects. In order to reliably predict measurable quantities, such as a fiducial cross section, the radiative corrections determined in QCD perturbation theory must be resummed to all orders. Resummation was first performed for energy-energy correlations in e + e − collisions [1][2][3], transverse momentum dependent cross sections in Drell-Yan events [4][5][6] and e + e − hadronic event shapes [7]. Several observables in hadron collisions have also been resummed analytically [8]. Such calculations have been extended to very high precision and used, for example, to extract the strong coupling from experimental data in e + e − annihilation to hadrons [9,10]. Effective field theory methods [11,12] also contribute to rapid progress in this field. General semi-analytic approaches to the problem have been constructed [13][14][15][16][17][18][19] and automated [20] based on direct QCD resummation. They depend only on universal coefficients and are applicable to different processes and a large class of observables. An alternative to analytical calculations is the simulation of events in a Markov-Chain Monte-Carlo known as a parton shower [21][22][23][24]. While the formal precision of this approach is comparable to analytic resummation only in processes with a trivial color structure at the leading order, parton-showers typically give a good description of experimental measurements and are therefore an integral part of the high-energy physics toolkit.
Even in the simplest scenarios the resummation performed by a parton shower is not identical to an analytic computation. This study will investigate the differences in some detail. We first show how a parton shower can be constructed that reproduces the pure next-to-leading logarithmic (NLL) resummed result as obtained by the semianalytic CAESAR formalism [18]. For simplicity we will focus on additive observables in e + e − annihilation to jets. Starting from this algorithm, we successively include effects beyond NLL accuracy that arise from momentum and probability conservation, such that a traditional parton shower in the coherent branching formalism is recovered eventually. To our knowledge this is the first time that a systematic study of this type has been performed. While we focus on a very simple setup, for which parton showers have been shown to achieve NLL accuracy [25], we argue that most differences investigated here will also arise in more complicated scenarios, such as hadron-hadron collisions and processes with a non-trivial color structure at the Born level. They will impact any prediction made for the Large Hadron Collider and possible future colliders, and -while formally sub-leading -they may be numerically large and should be taken into account as a systematic uncertainty. This paper is organized as follows: Section II recalls those parts of the CAESAR formalism and of the parton shower formalism needed in this study. Section III presents the technical details of a modified parton shower reproducing exactly the analytic NLL result. Section IV analyzes the role of NLL approximations in detail by removing them from the previously constructed shower one-by-one. Section V compares the full parton-shower result against a more conventional parton-shower implementation, where soft double counting is removed by partial fractioning of the soft eikonal. Section VI presents our conclusions.

II. NLL RESUMMATION AND THE PARTON SHOWER FORMALISM
We first review the methods used for analytic resummation in CAESAR [18] as well as the parton shower algorithm [21][22][23][24]. They are cast into a common language in order to allow an easy comparison between the two. We focus on the simplest case of resummation of a 2-jet observable in e + e − → jets, i.e. resummation of soft gluons emitted from a pair of two hard quark lines.

A. Prerequisites and notation
Following the CAESAR formalism, we denote the momenta of the hard partons as p 1 , . . . , p n . Additional soft emissions are denoted by k, and the observable we wish to compute by v. In general, the observable will be a function of both the hard and the soft momenta, v = V ({p}, {k}), while in the soft approximation it reduces to a function of the soft momenta alone, v = V ({k}). In the rest frame of two hard legs, i and j, considered to be the radiating color dipole, we can parametrize the momentum of a single emission as We define the rapidity of the emission in this frame as η ij = 1/2 ln(z i,j /z j,i ). The observable, computed as a function of k when radiated collinear to the hard parton, l, can then be written as 1 .
where, in the collinear limit, we have k T,l = k T,lj and η l = η lj for any j ∦ l. We restrict our analysis to the case of additive observables, which can be calculated in the presence of multiple soft gluons as a simple sum, V (k 1 , . . . , k n ) = n i V (k i ). Such observables are of great interest phenomenologically and, while relatively easy to compute, already exhibit most complications associated with the effects of NLL approximations.
The parton shower used in our study will be based on DGLAP evolution [26][27][28][29]. At NLL, for recursive infrared and collinear safe observables, gluon splitting only contributes at the inclusive level and is therefore taken into account effectively by working in the CMW scheme [25]. In analogy to the NLL CAESAR formalism, our parton shower will therefore only implement gluon radiation off the hard partons, and soft double counting will be removed by sectorization of the soft-emission phase space. Technical details are given in Sec. III, and a comparison to more conventional parton showers, which include gluon splitting, is performed in Sec. V. The basis for DGLAP evolution are the collinear factorization properties of QCD matrix elements. With |M n (1, . . . , n)| 2 being the squared n-parton matrix element, the factorization formula in the limit that partons i and j become collinear reads In this context, dΦ n is the n-particle phase space element, and P ij i (z) is the Altarelli-Parisi splitting kernel associated with the branching of an intermediate parton ij into partons i and j. Except for the analysis in Sec. V, the only relevant splitting kernel in our study is the quark-to-quark transition The treatment of gluon radiators is discussed in App. A. We denote the unregularized splitting probability between two scales, t and t , as Following standard practice to improve the logarithmic accuracy of the resummation, the strong coupling is evaluated at the transverse momentum of the gluon [30], and the soft enhanced term of the splitting functions is rescaled by 1 + α s /(2π)K, where K = (67/18 − π 2 /6) C A − 10/9 T R n f [25]. The latter method is known as the CMW scheme. The integration boundaries for z depend on the evolution variable and are given by the constraint that the momentum in the anti-collinear direction must be preserved. For the case of evolution in collinear transverse momentum, The probability for no splitting between two scales can be inferred from a unitarity constraint, i.e. the condition that the parton shower be probability conserving. For final-state evolution the no-branching probability is given by Note that this particular form of the no-branching probability is equivalent to the Sudakov form factor only at leading order, cf. App. A. Since we neglect gluon splitting, the functional form of R is unchanged until the shower terminates, which greatly simplifies the calculation 2 . The parton shower algorithm solves for the scale t , based on a starting scale t and the total branching probability (differential in ln t), It terminates when a cutoff scale t c is reached. Typically, t c is defined such as to mark the transition to the nonperturbative regime, i.e. the region where α s /(2π) ≈ 1.

B. Casting analytic resummation into the parton shower language
To enable a comparison with the semi-analytic resummation framework of CAESAR, we consider the cumulative cross section in an arbitrary observable, v, defined as The calculation is simplified by choosing a parton shower evolution variable, ξ, that (up to a power) corresponds to V (k) This implies that splittings giving the largest contribution to the observable are produced first. Note that here and in the following we use k T = k T,l and η = η l . If the effects of multiple emissions could be ignored, the cumulative cross section in Eq. (8) would be given by the square of the survival probability, Eq. (6), corresponding to the fact that radiation of a single gluon can originate from either of the two hard legs in the two-quark leading-order final state. It would then be sufficient to compute the probability R(v) = R(v, 1) for emissions resulting in observable values larger than v. Already at the level of a single emission this would lead to double counting [22]. The problem can be circumvented by sectorizing the phase space using the requirement η > 0. Note that this constraint is not strictly necessary for the collinear part of the splitting function if the parton shower implementation is capable of handling negative weights. However, this is not the case for most traditional shower algorithms, which prompts us to apply the condition to the entire splitting function. The combined probability for a single emission from any of the two hard legs at ξ > Q 2 v 2/(a+b) can then be written as This should be compared to Eq. (2.17) of Ref. [18], which can be rewritten in our parametrization as A brief summary of semi-analytic resummation based on [18] and using Eq. (11) can be found in App. B. The noemission probability based on R NLL (v) can also be computed in a Markovian Monte-Carlo simulation, by starting from the parton-shower expression, Eq. (10), and performing the following manipulations: • The z-integration in the soft term runs from 0 to 1 − (ξ/Q 2 ) (a+b)/2a , where the upper bound stems from the requirement that η > 0 (the Θ-function in Eq. (10)), eliminating the double counting of soft-gluon radiation.
• The collinear term proportional to (1 + z) is integrated from 0 to 1 in order to produce the collinear anomalous dimension, B q . At the same time, α s is evaluated at ξ.
Note in particular that the z-integration is extended beyond the values z min (and z max in the collinear case) allowed by local four-momentum conservation. This will be one of the effects investigated in Sec. IV. The complete parton-shower prediction of the cumulative cross section, Σ(v), including effects from arbitrarily many emissions, and using the approximation We compare Eq. (12) to the main result of [17], which reads The exponential corresponds to the pure survival probability in terms of Eq. (11). The function F (v) accounts for the effect of multiple emissions. For the simple observables considered here it can be written as [18] 3 : Following the notation of Ref. [18], R NLL (v) is the derivative of R with respect to L = − ln v, excluding all terms formally not relevant at NLL accuracy. Note that Eq. (14) is a pure NLL contribution to Σ NLL (v), as R (v) by itself is sub-leading. If we intend to generate Eq. (14) using a parton shower, the branching probability, Eq. (10), must be modified such as to reflect the differentiation w.r.t. the lower integration limit in Eq. (11), which leads to ξ = Q 2 v 2/(a+b) , as well as the condition that higher logarithmic terms are dropped in R (v). We can satisfy these constraints using the following modifications of the plain parton shower: • The z-integration in the soft term runs from 0 to 1 − v 1/a .
• The strong coupling runs at one loop and is evaluated at v 2/(a+b) (1 − z) 2b/(a+b) .
• The collinear term is dropped.
We can now rewrite Eq. (13) in a form that is similar to Eq. (12) with R given by The choices of α s , z max and µ 2 corresponding to NLL resummation in the CAESAR formalism and in a DGLAP-based parton shower are given in Tab. I. The physical limits on the z-integral in Eq. (10), which are a consequence of local four-momentum conservation, are not easily formulated in terms of ξ and will be investigated separately in Sec. IV. It is interesting to note that F (v) by itself can be extracted from the same formalism by starting the shower evolution at   (12) (parton shower). The effects of switching between the two parametrizations are investigated in the figure referred to in the last column. More details can be found in Sec. IV.

III. MARKOV-CHAIN MONTE CARLO IMPLEMENTATION
As described in Sec. II, the NLL resummation is nearly equivalent to a parton shower at the single-emission level. The differences lie in the treatment of the collinear term and of the lower integration boundary on z. These differences also introduce a change in the scale of the running coupling in Eq. (15). The choice of integration boundaries in the analytic resummation implies that the splitting function turns negative in parts of the phase space. To deal with this situation in the Monte Carlo simulation, we use the methods discussed in [37,38]. Splittings are generated according to an overestimate of the strong coupling and the splitting kernel with an in principle arbitrary constant γ. For practical calculations we choose γ = 2. Note that the values of z max soft and z max coll are overestimated by a common value in P max , which we have made explicit by writing z max . Splittings are vetoed with a constant probability 1/C and are associated with a weight This correction accounts in particular for the negative sign of the integrand, Eq. (19), in the region z > z soft max . In addition, it is possible to veto emissions violating the condition i V (k i ) < v, which would contribute with zero weight, to improve numerical accuracy [38]. The value of C determines how many emissions are proposed, and thus potentially vetoed. It can again in principle be an arbitrary constant larger than 1, but is relevant for the speed of convergence. We choose C = 2 in our implementation. The kernel eventually used for NLL resummation is given by with z max and µ 2 chosen according to Table I. For multiple emissions P res explicitly depends on v. We therefore first choose a value for v and then run the parton shower, implementing the z integration bounds and the scale of the strong coupling as defined in Tab. I. This is a highly inefficient procedure to compute the cumulative cross section. If probability was conserved, the same distribution could be obtained by running the parton shower, computing v, filling the histogram in each bin with lower edge larger than v, and filling the histogram in the bin containing v with weight (v max − v)/∆v, where v max is the upper bin edge and ∆v is the bin width. This will be the method used to compute the predictions in Fig. 6 and Sec. V. While at the level of accuracy we are interested in, it is sufficient to set the cutoff scale of the parton shower to some numerically small value in ξ, exact agreement with the analytic calculation is expected only if the calculation is performed for a finite , and the parton-shower cutoff is set to ξ c = v. We can verify that in this situation we reproduce the analytic result for finite in Eq. (14) and investigate the convergence towards the analytic result for → 0. Figure 1 presents the corresponding comparison for different values of in the case of the thrust (a) [39], a BKS observable (b) [40,41] and a fractional energy correlation (c) [18]. The definitions of the observables and related resummation coefficients are listed in App. C.

IV. EFFECTS OF APPROXIMATIONS
This section is dedicated to the detailed investigation of the effects of local four-momentum conservation and approximations made in the NLL calculation compared to the parton shower. In order to cover different choices of the parameters a and b, we again present results for the thrust, a BKS observable (x = 1/2) and a fractional energy correlation (x = 1). All distributions are shown for Q = 91.2 GeV, and for a strong coupling defined by α s (Q 2 ) = 0.118 and a fixed number of flavors, n f = 5. We have cross-checked all of our predictions using two independent Monte-Carlo implementations based on [42].
We first investigate constraints arising from momentum conservation in the anti-collinear direction at single emission level, which reads This induces both a lower and an upper bound on z given by z min/max = (1 ∓ 1 − 4k 2 T /Q 2 )/2. Figure 2 shows a comparison between the pure NLL predictions and those where this constraint has been implemented. The effect on the cumulative distributions is moderate, about 5% in the medium and low-v region. In addition, we investigate the effect of choosing the scale in the collinear term to be k 2 T . This alters the slope of the thrust and BKS 1/2 distributions in the small-v region, due to additional sub-leading logarithmic terms in R(v).
The upper bound z max is generally weaker than the constraint arising from the condition η > 0, listed in Tab. I. Figure 2 displays the additional effect on the NLL prediction when this constraint is applied in form of z max ≶v,coll as used in typical parton showers (cf. Eq. (10)). The effects are about 10% on all observables, and they lower the prediction for Σ(v) due to an increased branching probability. Again, we also investigate the effect of choosing the scale in the collinear term to be k 2 T , which generates the same slope differences at small v observed before. Next we investigate the effect of lifting the restriction on the z integration in the calculation of F(v), i.e. removing the constraint z < 1 − v 1/a if ξ < Q 2 v 2/(a+b) and replacing it by the constraint η > 0. In this case R (v) must be computed down to very small scales in Eq. (14), (except for FC 1 ) and it becomes mandatory to introduce an additional cutoff, as one would otherwise need to evaluate α s at values where perturbation theory is no longer valid. We choose to implement this by adding the requirement k min T = 0.5 GeV. The difference to the pure NLL result is shown in Fig. 3. Independent of the observable, this change is one of the largest differences observed in this study. The large relative difference between the pure NLL result and the modified prediction at small v shows that sub-leading logarithmic effects become important.
We also study the effect originating in the evaluation of the running coupling at Q 2 v 2/(a+b) (1 − z) 2b/(a+b) if ξ < Q 2 v 2/(a+b) . Again we implement the constraint k min T = 0.5 GeV. Figure 3 shows that the predictions for all observables single emission   exhibit relatively large changes. They also show convergence issues at small values, which arise from the lower cutoff in k T , leading to an insufficient sampling of F at higher number of emissions. This effect is most pronounced in FC 1 , where it starts to appear around 10 −1.5 . Note, however, that practical measurements of FC 1 would be impacted by non-perturbative corrections in this regime. The problem is therefore purely academic in nature, hence we do not attempt to solve it here. Formally the CMW scheme is the key to achieving NLL accuracy in a parton-shower computation of the observables considered here [25]. The numerical impact on R(v) is investigated in Fig. 4. Figure 5 displays the effect of replacing 1-loop by 2-loop running couplings and of using the CMW scheme in sub-leading terms of the NLL calculation (cf. Tab. I). The red line is computed by making the replacements only in the soft-enhanced part of the splitting function for ξ < Q 2 v 2/(a+b) , and the red dotted line corresponds to not using the CMW scheme if ξ < Q 2 v 2/(a+b) . It is evident that the effects are sizable over most of the observable range, and most pronounced at small v. The use of the CMW scheme has the biggest impact. Note in particular that not using the CMW scheme in the computation of F(v) has nearly the same impact as not using the CMW scheme in the computation of R(v). Figure 6 shows the cumulative effect of all changes discussed so far. In addition we present results from a simulation where the observable is computed using its definition in terms of four-momenta rather than using the soft approximation in Eq. (2) (see App. C for details). In this context it becomes important to take into account that emissions away from the strict soft limit inevitably change the momenta of the hard partons. Subsequent emissions are then computed based on the momenta of the quark lines with recoil effects taken into account. This can have a significant impact on the result, depending on the precise definition of the transverse momentum and momentum fraction. The magenta line in Fig. 6 corresponds to the conventions of [43], while the green line corresponds to the conventions  of [44]. In the latter case the transverse momentum coincides with Eq. (2.5) of [18]. 4 Note that the phase-space sectorization constraint, η > 0, generates a different restriction on z once recoil is taken into account, and that this condition depends on the choice of evolution and splitting variable.

V. COMPARISON WITH A DIPOLE-LIKE PARTON SHOWER
This section presents a comparison of our previous results with predictions from a dipole-like parton shower. In such parton showers the soft enhanced part of the collinear splitting function is typically replaced by a partial fraction of the soft eikonal matched to the collinear limit [43,45]. At the same time, the phase-space sectorization is removed, i.e. the restriction η > 0 is lifted. A complete description of the parton-shower algorithm employed here can be found in [44]. Figure 7 shows a comparison between results from the dipole-like parton shower in its default configuration (including gluon splitting) and from a modified version, tailored to match the settings of the parton shower used in Sec. IV, Fig. 6. It is interesting to observe that the dipole-shower prediction lies between the parton-shower result and the analytic result for all observables, and in the case of thrust agrees very well with the analytic prediction. In the measurable range at LEP energies, the predictions for FC 1 also agree fairly well between the dipole-shower and the analytic result. Figure 8 displays a cross-check on the logarithmic terms implemented by the dipole shower as compared to the parton shower and the analytic result. We extract R(k T /Q) for a fixed value of the strong coupling, α s = 0.118, using the technique described in [38]. The slope of the distribution corresponds to the leading logarithm, while the offset of the analytic result corresponds to the next-to-leading logarithm. Any parton-or dipole-shower prediction must approach the analytic result as k T → 0, which is verified by the convergence of the predictions at small k T .

VI. CONCLUSIONS
We have performed a detailed comparison between pure NLL resummation and parton showers for additive observables in e + e − annihilation to hadrons. We have isolated their differences, which can broadly be classified as related to probability or momentum conservation. While a different treatment of these effects leads to formally subleading corrections on the resummed prediction, it can have a numerically sizable impact (20% or more) in the region where experimental measurements are performed. Similar effects can reasonably be expected to arise in other observables, as well as in processes with hadronic initial states and with a more complicated color structure at the Born level. When comparing analytic resummation to parton showers it should be kept in mind that such differences may exist, in which case they should be taken into account as a systematic uncertainty. We have shown in a simple scenario that the differences can be assessed quantitatively by casting analytic resummation into a Markovian Monte-Carlo simulation and introducing momentum and probability conservation. Conversely, parton showers can be modified to violate momentum and probability conservation to reproduce pure NLL resummation. From the practical point of  10) and (11). The plot shows the average number of emissions per bin as a proxy observable [38].
view this approach is disfavored, as it leads to numerically inefficient Monte-Carlo algorithms. unitarity constraint [47]. For 0 < ε 1, Eq. (A1) changes to Using the Sudakov form factor the generating function for splittings of parton a is defined as Equation (A3) can now be written in the simple form The generalization to an n-parton final state, a = {a 1 , . . . , a n }, resolved at scale t can be made in terms of fragmenting jet functions, G [48,49]. If we define the generating function for this state as F a ( x, t, µ 2 ), we can formulate its evolution equation in terms of a sum of the right hand side of Eq. (A6). For unconstrained evolution, we can use Eq. (A3), to write the differential decay probability as Thus, as highlighted in [47], it is generally necessary to use the Sudakov factor, Eq. (A4), in final-state parton shower evolution. At the leading order, the factor ζ in Eq. (A4) simply replaces the commonly used symmetry factor for g → g splitting and it also accounts for the proper counting of the number of active flavors. 5 However, this reasoning applies only if the boundaries of the ζ-integration are defined by momentum conservation, and are therefore symmetric around ζ = 1/2. In our analysis we attempt to extend the lower integration limit to zero, which would generate a spurious singularity arising from the symmetry of the gluon splitting function. Therefore, the commonly used technique of implementing the symmetrized gluon splitting function without an additional factor ζ cannot be used, and the only correct way to treat the problem is to work with Eq. (A4).

Appendix B: Analytic results at NLL accuracy
This section summarizes the components of the CAESAR formalism [18] that are needed for our analysis. The resummed cumulative cross section at NLL is given in this formalism by Σ NLL (v) = e −RNLL(Q 2 ,v) F(v), cf. Eq. (13). The unregularized branching probability R(v) follows from Eq. (11). It is typically written in terms of λ = α s β 0 L, where L = − ln v. One obtains where r(L) is separated into a leading and a sub-leading logarithmic piece as r(L) = Lr 1 (α s L) + r 2 (α s L).
The beta function coefficients and the two-loop cusp anomalous dimension in the MS scheme are given by β 0 = 1 2π T R n f , T R n f .

(B3)
The sub-leading logarithmic term T (L) is defined as The F-function, Eq. (14), for additive observables is given by .
Since both T (L) and r 2 (L) are sub-leading in L, we have Using Eq. (B6) it can be verified that the combination of z max <v,soft and µ 2 <v,soft listed in Tab. I generates the correct value of r (L), and therefore the correct value of the F-function.