Implementing NLO DGLAP evolution in Parton Showers

We present a parton shower which implements the DGLAP evolution of parton densities and fragmentation functions at next-to-leading order precision up to effects stemming from local four-momentum conservation. The Monte-Carlo simulation is based on including next-to-leading order collinear splitting functions in an existing parton shower and combining their soft enhanced contributions with the corresponding terms at leading order. Soft double counting is avoided by matching to the soft eikonal. Example results from two independent realizations of the algorithm, implemented in the two event generation frameworks Pythia and Sherpa, illustrate the improved precision of the new formalism.


I. INTRODUCTION
Parton shower algorithms describing QCD and QED multiple radiation have been a central ingredient of simulation programs for particle physics experiments at the energy frontier [1]. After their inception about three decades ago, where the focus was on QCD radiation in the final state [2], efficient algorithms for initial state radiation were developed [3], which amount to evolution back in "time" from the hard scattering to the incoming beam hadrons. The study of quantum interference effects in successive emissions led to the notion of QCD coherence in parton evolution [4,5], and angular ordering was identified as a convenient scheme that incorporates such effects [6][7][8]. As an alternative scheme, the color dipole model [9,10] includes QCD coherence in a natural way. Matrix-element corrections have been investigated as a source of coherence [11,12].
After about a decade of work on matching [13][14][15][16][17][18] and merging [19][20][21][22][23][24][25][26][27][28][29][30] algorithms, the necessity of increased control over the parton shower for a more seamless combination with fixed-order calculations at higher orders triggered a resurgence of interest in improving parton-shower algorithms themselves. As a consequence, new parton showers [31][32][33][34][35][36][37][38], have been constructed that are based on ordering subsequent emissions in transverse momenta, and there were also new constructions with improved and generalized angular ordering parameters [39]. The possibility of including next-to-leading order corrections into parton showers was explored over three decades ago [40][41][42][43], and it was revisited recently in a different framework [44,45]. Next-to-leading order corrections to a single final-state gluon emission off a qq dipole have been presented in [44] as a first higher-order extension of the antenna shower formalism. How this approach maps onto NLO DGLAP evolution was briefly addressed in [45], which furthermore introduced final-state double-gluon radiation into this formalism. In addition to this, NLO splitting functions have been recomputed using a novel regularization scheme [46,47], with the aim to improve parton-shower simulations. The dependence of NLO matching terms on the parton-shower evolution variable has also been investigated [48].
This publication is dedicated to the construction of a parton shower that implements the next-to-leading order (NLO) DGLAP equations up to momentum conserving effects. We employ the non-flavor changing NLO splitting functions in the MS scheme in their integrated form [49][50][51][52][53][54], and we include the flavor-changing NLO splitting kernels fully differentially using the method presented in [55]. We identify the contribution to the NLO splitting functions which is already included in the leading-order (LO) realization of the parton-shower, and correspondingly subtract it from the NLO splitting function. This term is given by the two-loop cusp anomalous dimension, which is usually included at LO using the CMW scheme [56]. After its subtraction, the remaining splitting function is purely collinear, and no double-counting arises upon implementing it as a higher-order correction to the existing splitting kernels of the parton shower. However, the NLO parts of the splitting functions are negative in large parts of the phase space which presents a technical challenge. We overcome this problem through the weighting algorithm first proposed in [57,58]. Our approach can be considered as a first step towards a fully next-to-leading order accurate parton shower and acts as a baseline for further development. Future projects will need to address the leading-color approximation and the simulation of soft emissions beyond the leading order. The Monte-Carlo techniques developed here are expected to become useful in this context as well. A clear phenomenological benefit of the present implementation is that consistency between the parton shower and NLO PDF evolution is achieved for the very first time.
The outline of this paper is as follows. Section II introduces the parton-shower formalism at leading order and establishes the connection to the DGLAP equation in order to identify the correct treatment of the final-state Sudakov factor. Section III outlines the specific implementation in the DIRE parton showers [38]. First results and applications are presented in Sec. IV. Section V contains our conclusions.

II. EXTENSION OF THE PARTON-SHOWER FORMALISM
In this section we will highlight the correspondence between the parton shower formalism and the analytic structure of the DGLAP evolution equations [59], on which the parton shower is based. We will thereby focus on the refinements needed in order to realize NLO accurate parton evolution. This includes the implementation of the complete set of splitting kernels at O(α 2 s ), a subset of which are the flavor-changing kernels discussed in [55]. Another significant change concerns the implementation of symmetry factors. In the computation of the NLO splitting functions [49][50][51][52][53][54], it is assumed that a certain final-state parton is identified, while the parton shower treats all particles democratically. If the full set of splitting functions -both at leading and at next-to-leading order -is implemented naively, the emission probability will thus be overestimated. At leading-order the problem can be solved by adding simple symmetry factors. At next-to-leading order the solution will be different, as the splitting functions include contributions from three-particle final states that have been integrated out. We will show in the following how this problem can be approached [55].

A. Unconstrained evolution with identified final-states
The DGLAP equations are schematically identical for initial and final state. However, their implementation in parton-shower programs differs between the two, owing to the fact that Monte-Carlo simulations are typically performed for inclusive final states. The inclusive evolution equations for the fragmentation functions D a (x, Q 2 ) for parton of type a to fragment into a hadron read where the P ab are the unregularized DGLAP evolution kernels, and the plus prescription is defined to enforce the momentum and flavor sum rules: For finite ε, the endpoint subtraction in Eq. (2) can be interpreted as the approximate virtual plus unresolved real corrections, which are included in the parton shower because the Monte-Carlo algorithm naturally implements a unitarity constraint [60]. The precise value of ε in this case depends on the infrared cutoff on the evolution variable, and is determined by local four-momentum conservation in the parton branching process. For 0 < ε 1, Eq. (1) changes to Using the Sudakov form factor the generating function for splittings of parton a is defined as Equation (3) can now be written in the simple form The generalization to an n-parton state, a = {a 1 , . . . , a n }, with jets and incoming hadrons resolved at scale t can be made in terms of parton distribution functions (PDFs) f , and fragmenting jet functions, G [61,62]. 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. (6), where each term in the sum corresponds to a resolved jet in the final state or a hadron in the initial state. This equation can be solved using Markovian Monte-Carlo techniques in the form of a parton shower [1]. In most cases, however, parton showers implement final-state branchings in unconstrained evolution, which means that final-state hadrons are not resolved. We can use Eq. (3) (which also applies to G [61,62]), to write the corresponding differential decay probability for such an evolution as As highlighted in [60], it is necessary to use the Sudakov factor, Eq. (4), in final-state parton showers beyond the leading order. At the leading order, the factor ζ in Eq. (4) 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. However, at the next-to-leading order it becomes an identifier for the parton that undergoes evolution, which is essential in order to obtain the correct anomalous dimensions upon integration of the NLO DGLAP evolution kernels. We will thus define the final-state Sudakov factor in our implementation according to Eq. (4).

B. Splitting functions
The crucial ingredient of NLO DGLAP evolution are the O(α 2 s ) corrections to the evolution kernels. These corrections depend on the scheme in which PDFs and fragmentation functions are renormalized. We will work in the MS scheme, which allows us to use the results of [49]. Technical challenges in the implementation of the splitting functions in the parton shower include the overlap with the CMW scheme for setting the renormalization scale commonly used in leading-order parton showers [56] as well as the fact that the evolution kernels are negative in large parts of the accessible phase space. We will discuss these problems in a general context in the following and give more details on the implementation in the Dire parton showers in Sec. III.
At O(α s ), the unregularized DGLAP splitting functions are At O(α 2 s ), the quark splitting functions are typically written in terms of singlet (S) and non-singlet (V) components as In the timelike case, their components are The flavor-changing splitting kernels, P qq and P (1) qq first appear at O(α 2 s ). They are new channels contributing to the real-emission corrections to P (0) qq . In order to account for their more involved flavor structure, we must simulate them fully differentially in the 1 → 3 phase space. To this end, we use the method presented in [55]. All other splitting functions have an analogous 1 → 2 topology and are implemented using this topology.
Several new structures appear in the next-to-leading order splitting functions, which require a modification of the branching algorithm used at the leading-order. Firstly, the NLO splitting functions may exhibit new types of apparent singularities, like the term −20/9 C F T F /z contributing to p S(1) qq . Such terms are regulated by the symmetry factor in Eq. (4), which highlights again that without the correct definition of the Sudakov factor one cannot construct a meaningful Monte-Carlo implementation, as the resulting integrals would have unphysical divergences.
In addition, p V (1) qq and P (1) gg include the two-loop cusp anomalous dimension, given by [56] This term is routinely included in standard parton-shower Monte Carlo simulation, typically through a redefinition of the scale at which the strong coupling is evaluated [56]. It must therefore be subtracted from p V (1) qq and P (1) gg before these splitting functions can be included. After the subtraction, no soft enhanced terms remain, and the result is a purely collinear splitting function. This is important to avoid double counting of singular limits in the parton shower [8]. Furthermore, p V (1) qq and P (1) gg also contain a term originating from the renormalization of the strong coupling constant, which is given by the leading-order splitting function times β 0 log z, where The leading contribution from this term upon integration over z is generated in combination with the soft factor 2/(1−z) of the leading-order splitting function, and gives a contribution −β 0 π 2 /3 to the collinear anomalous dimension.

III. IMPLEMENTATION IN THE DIRE PARTON SHOWER
Our numerical simulations are based on the DIRE parton shower, presented in [38]. This section first presents a brief overview of the model as implemented at leading order, before moving to the modifications needed for an implementation of the next-to-leading order contributions.
A. Parton-shower model at leading order The evolution and splitting parameters κ and z j used in DIRE for splittings of a combined parton ij to partons i and j in the presence of a spectator k are given by In this context, Q 2 plays the role of the maximally attainable momentum squared, which is defined as Q 2 FF = 2(p i + p j )p k + 2p i p j for final-state splittings with final-state spectator, Q 2 FI = Q 2 IF = 2(p i + p j )p k for final(initial)-state splittings with initial(final)-state spectator, and Q 2 II = 2p i p k for initial-state splittings with initial-state spectator. The splitting functions for initial-state branchings are given by the modified DGLAP splitting functions [38] where z = 1 − z j . It is interesting to note that the dimensionless quantity κ 2 plays the role of the IR regulator in the very same fashion as the principal value regulator δ 2 introduced in Eq. (3.13) of [49]. In our algorithm, κ has a physical interpretation, as the scaled transverse momentum in the soft limit. As such, it also sets the renormalization and factorization scale through µ 2 R/F = κ 2 Q 2 . For final-state branchings, the matching to the differential cross section in the soft limit requires the replacement where the j-soft part of the splitting function is given by In a similar fashion we have The two terms in Eq. (15) correspond to different color flows in the parton shower. For the first term partons i and k are considered radiators and j is the soft gluon insertion, while for the second term partons j and k are the radiators and i is the soft gluon insertion. Therefore, in the first term gluon j is color-connected to the spectator parton, while in the second term gluon i is color-connected to the spectator. The two contributions are evolved using the two different variables κ 2 j,ik and κ 2 i,jk . Following standard practice to improve the logarithmic accuracy of the parton shower, the soft enhanced term of the splitting functions, Eqs. (14), is rescaled by 1 + α s (t)/(2π) Γ (2) [56]. We do not absorb this rescaling into a redefinition of the strong coupling, as this would generate higher-logarithmic contributions stemming from the interaction with the purely collinear parts of the splitting functions.

B. Extension to the next-to-leading order
We now describe the extensions of the DIRE parton shower that are necessary to construct a simulation which describes the DGLAP evolution of parton distributions and fragmentation functions at next-to leading order precision. As an important construction paradigm, we consider contributions at different orders in the strong coupling as separate evolution kernels, and we restrict ourselves to the inclusive radiation pattern where possible. The latter implies that in general we do not attempt to simulate the emission of an unordered pair of partons according to the triple collinear splitting functions. The notable exception to this is the treatment of flavor-changing splitting functions, where the implementation of a 1 → 2 rather than a 1 → 3 transition is not possible due to local flavor conservation. The generation of these contributions is described in detail in [55] 1 . The main remaining complication in the implementation of the integrated NLO splitting functions arises from the fact that they assume negative values in large regions of phase space, hence we potentially need to generate branchings based on negative "probabilities". To this end we use the method developed in [14,57,58].
We start by formally replacing the leading-order splitting functions of Eq. (14) with the combined leading-order plus next-to-leading order evolution kernels. 2 As described in Sec. II, the soft enhanced part of P (1) qq and P (1) gg matches the term α s /(2π)Γ (2) 2C a /(1 − z), at leading order, which is included in the implementation of the leading-order parton shower by rescaling the soft enhanced part of the splitting functions. We therefore subtract this contribution from the NLO splitting kernel and define In addition, we include in the soft enhanced part of the leading-order splitting function the three-loop coefficient Γ For final-state gluon evolution this requires the independent modification of both terms in Eq. (15). Scale variations can be performed in the DIRE showers by using a method similar to [64]. When varying the argument of the strong coupling, i.e. replacing α s (t) → α s (c t) f (c, t), with c a constant, the appropriate counterterm at O(α 2 s ), which multiplies the leading-order splitting functions, P We use the multiplicative threshold matching described in [64], as the additive matching generates artificially large deviations in the case of two-loop and three-loop running of the coupling. The product in Eq. (21) runs over the number n th of parton mass thresholds in the interval (t, c · t) with t 0 = t, t n th +1 = c · t and t i are the encompassed parton mass thresholds. If c < 1, the ordering is reversed to recover the correct sign. β 0 (t) is the QCD beta function coefficient, which depends on the scalet through the number of active parton flavors. 1 The contribution from triple collinear splitting functions of type q → q and q →q to the overall NLO corrections is numerically small. A more detailed discussion can be found in [55]. 2 For a complete list of the NLO splitting functions see App. A. Note that we do not require the knowledge of p V (1) qq in our approach, because flavor-changing splittings are generated fully differentially in the 1 → 3 phase space. 3 The normalization differs by factor four between our notation and that of [63]. 4 Note that the lowest-order DGLAP kernels, P (0) ab , are defined at O(αs), and we use a strict order counting. The scale variations in our approach are therefore more conservative than the ones presented in [64].

IV. DIRE PREDICTIONS
We have implemented our new algorithms into the DIRE parton showers, which implies two entirely independent realizations within the general purpose event generation frameworks PYTHIA [3,65] and SHERPA [66,67]. This section presents a first application of our new algorithm to the simulation of the reactions e + e − →hadrons, pp → e + ν e and pp → h. We compare the magnitude of the next-to-leading order corrections and the size of their uncertainties to the respective leading-order predictions. Note that we only quote the renomalization scale uncertainties, which are the ones that can be expected to decrease when moving from leading to next-to-leading order evolution. There are of course other uncertainties, for example those related to the kinematics mapping and the choice of the evolution variable in the parton shower. However, these effects arise identically both at leading and at next-to-leading order, and they are therefore not included in the uncertainty bands. In addition, nonperturbative effects will contribute their own uncertainty, which is somewhat harder to quantify. However, it is expected that a reduced perturbative uncertainty will lead to a more realistic extraction of nonperturbative model parameters, and that the uncertianties on those parameters can therefore be reduced as well. Figure 1 shows predictions from our new implementation compared to leading-order results from the DIRE parton shower for differential jet rates in the Durham scheme compared to experimental results from the JADE and OPAL collaborations [68]. Results have been obtained with DIRE+SHERPA using the default settings of SHERPA version 2.2.3. The perturbative region is to the right of the plots, and y ∼ 2.8 · 10 −3 corresponds to the b-quark mass. The simulation of nonperturbative effects dominates the predictions below ∼ 10 −4 . In the perturbative region, the results are in excellent agreement with the experimental measurements. The shapes of distributions receive only minor changes compared to the leading-order result, however, the uncertainties are greatly reduced. Figure 2 shows a comparison for event shapes measured by the ALEPH collaboration [69]. The perturbative region is to the right of the plot, except for the thrust distribution, where it is to the left. We notice some deviation in the predictions for jet broadening and for the C-parameter, which are largely unchanged compared to the leading-order prediction. These deviations are mostly within the 2σ uncertainty of the experimental measurements, and they occur close to the nonperturbative region, which indicates that they may be related to hadronization effects.
In Figure 3, we illustrate the effect of NLO kernels on differential jet resolutions in Drell-Yan lepton-pair production as well as on Higgs-boson production in gluon fusion. In both cases, the impact of varying the renormalization scale in the parton shower is greatly reduced upon inclusion of NLO corrections, and shape-changes of O(10%) can be observed. It is interesting to note that these shape changes have the opposite effect in Drell-Yan lepton pair and Higgs boson production. This effect could not have been obtained by changing the argument of α s at leading order only, as in Eq. (21). Figure 4 confronts DIRE with Drell-Yan transverse momentum spectra measured by ATLAS [70]. We limit the comparison to the soft and semi-hard region of transverse momenta, p T < 30 GeV. Parton shower predictions are insufficient in the hard region, and the shower is usually supplemented with fixed-order calculations through matching or merging in order to improve upon this deficiency. Note that no tuning of DIRE+PYTHIA has been performed, neither in the default version nor for the present publication. Our results have been obtained with PYTHIA 8.226, using the NNPDF 3.0 (NLO) PDF set [71], α s (M Z ) = 0.118 throughout the simulation. The ISR/FSR shower cut-off has been set to 3 GeV 2 , and a primordial transverse momentum of k ⊥ = 2 GeV was used. All other parameters are given by the default tune of PYTHIA 8 [72]. NLO corrections improve the agreement with data particularly in the region where resummation has a large impact.

V. SUMMARY
In this paper we have presented an extension of the parton shower formalism to include, for the first time, the DGLAP evolution at next-to-leading order precision for both initial and final state radiation. The new terms are of order α 2 S , and they fall into three categories: Soft terms ∝ 1/(1 − z) which are multiplied by the two-loop cusp anomalous dimension, and which are routinely included in parton shower simulations through a suitable rescaling of the argument of the strong coupling. In addition there are genuine, non-trivial higher-order terms which modify already existing leading order terms. Although they are negative over a wide region of phase space we can include them as separate terms through existing reweighting techniques. Finally, there are new structures which correspond to flavor-changing transitions of the type q → q or q →q and which originate in genuine 1 → 3 transitions. The algorithm for their simulation is detailed in a separate publication. Including all these terms corresponds to adding the process-independent collinear enhanced NLO corrections present in standard DGLAP evolution into the parton shower. The overall effect of this increased precision is twofold. While for e − e + annihilations to hadrons the central values of distributions experience only marginal shifts, the situation is different for distributions at hadron colliders. This is exemplified by the transverse momentum distribution of Z-bosons produced at tree-level in qq annihilation   Results for leading and next-to-leading order DGLAP evolution in comparison to LEP data from [68].

MC/Data
Dire PS

MC/Data
Dire PS  Results for leading and next-to-leading order DGLAP evolution in comparison to LEP data from [69].
dσ/d log 10 (d01/GeV)  Predictions for leading and next-to-leading order DGLAP evolution for the differential kT -jet resolution parameters in pp → e + e − + X (LHC √ s = 7 TeV) and pp → h + X (LHC √ s = 8 TeV). Results for leading and next-to-leading order DGLAP evolution in comparison to ATLAS data from [70]. and the differential jet rates in Z and Higgs boson production through gluon fusion, respectively, which experience some shifts of up to about 10% relative to corresponding leading order distributions. In both cases, the uncertainty from variations of the renormalization scale by factors of two is significantly reduced when going from leading to next-to-leading order precision. For the first time, we are able to quote a realistic renormalization scale uncertainty as we only add renormalization counterterms which appear at the perturbative order to which we control the expansion of the splitting functions. While the work presented here represents a significant improvement over existing parton showers, it includes only parts of the higher-order corrections. We did not improve upon the leading color approximation typically used in the parton shower. Ways to include such corrections have been discussed in [73]. Furthermore, we did not include the effect of higher-order soft terms, i.e. the effect of multiple unordered soft emissions. We expect these terms to have only limited impact on inclusive observables such as standard event shapes or the transverse momentum of singlet particles produced at hadron colliders. They will mostly contribute to the further stabilization of perturbative predictions for these observables. However, we appreciate that they will certainly impact on non-global observables such as out-of-cone radiation which in turn renders their inclusion an important task for the future.