Hadronic Final States in DIS at NNLO QCD with Parton Showers

We present a parton-shower matched NNLO QCD calculation for hadronic final state production in Deep Inelastic Scattering. The computation is based on the UNLOPS method and is implemented in the publicly available event generation framework SHERPA. Results are compared to measurements performed by the H1 collaboration.


I. INTRODUCTION
Deep-inelastic lepton-nucleon scattering (DIS) has presented a formidable challenge to the theoretical high-energy physics community for some time. On the one hand, the reaction provides an extremely precise probe of the nucleon structure for nearly five decades [1]. QCD corrections to the structure functions have been computed through third order in perturbation theory [2][3][4][5][6]. On the other hand, the description of quantities like inclusive jet or di-jet differential cross sections remained difficult, even with computations performed at NLO QCD [7][8][9]. Tremendous progress has recently been made with the fully differential calculation of jet production in DIS at NNLO QCD accuracy [10][11][12], and with the fully differential calculation of inclusive DIS at N 3 LO precision [13]. Several of these results have been used in experimental analyses and provide a much improved description of the measurements [14,15].
In this publication we provide a Monte-Carlo simulation for the estimation of parton shower and hadronization effects that will allow for a yet more precise comparison between theory and experiment. We achieve this by combining the known NNLO QCD perturbative results [2][3][4] with a parton shower [16] based on dipole factorization [17]. We implement this simulation in the general purpose event generator SHERPA [18,19]. The leading-order QCD configuration corresponding to the DIS process is shown in Fig. 1(a). The kinematics are typically parametrized in terms of the exchanged boson's virtuality Q 2 = −q 2 = (k − k ) 2 and the Bjørken variable x = Q 2 /(2 qP ). The jet reconstruction is performed in the Breit frame, which is defined by the condition q + 2x P = 0. At the leading order, a single final-state parton emerges, which carries zero transverse momentum. The kinematical variables x and Q 2 can in principle be inferred from the incoming and outgoing lepton momenta k and k alone, and the only relevant scale in the problem is Q 2 . A measurement of jet production in DIS, however, introduces additional scales of the order of the jet transverse momenta and leads to significant challenges in the description of the corresponding final states. The problems are related to possibly inverted scale hierarchies between Q 2 and p 2 T , the (squared) transverse energy of the jet(s) in the Breit frame. While the proton is probed at Q 2 when considering inclusive DIS, the relevant hardness scale for the QCD Compton process leading to the production of a hard jet is of the order of Q 2 + p 2 T , see the sketch in Fig. 1(b). In the extreme case of two resolved hard jets and small Q 2 , one must even picture the core reaction as a 2 → 2 pure QCD process, followed by the initial-state branching of a virtual photon into quarks, with a relevant hardness scale of H 2 T , the total transverse hadronic energy in the Breit frame. The large available phase space and the increased scales lead to significantly enhanced multi-jet cross sections in the region p 2 T /Q 2 1. Similar effects were also observed in other reactions [20]. The problem was addressed in [21], using multi-jet merging. However, the non-unitary merging technique employed in the simulation made it difficult to simultaneously predict exclusive quantities like jet cross sections and inclusive quantities like the structure functions. In the approach presented here this problem is solved by means of a unitary merging method. Higher-order radiative effects are taken into account to a good approximation by including the complete fixed-order NNLO QCD corrections to inclusive DIS and by choosing an appropriate scale [12]. This paper is organized as follows: Section II presents an introduction to the UNLOPS matching technique and discusses the projection-to-Born method used in our study. Section III presents the numerical validation, assessment of theoretical uncertainties and comparison to experimental data. An outlook is given in Sec. IV.

II. COMPUTATIONAL SETUP
The starting point of our simulation is a fully differential calculation of the inclusive DIS process at NNLO QCD using the projection-to-Born method [22]. This technique relies on a map F b that uniquely assigns any given flavor and momentum configuration in the one emission phase space Φ 1 and in the double-real emission phase space Φ 2 to a point in the Born phase space Φ 0 . Given such a map, we define the Born-differential NNLO cross section as In this context, B n are the Born differential cross sections for DIS plus n partons, and V n and VV 0 are the corresponding UV renormalized virtual and double virtual corrections, including the appropriate collinear mass factorization counterterms. This cross section is free of divergences if F B is an infrared safe observable, i.e. if it maps a Born phase space point supplemented by infinitely soft gluons and/or collinear parton branchings to the same Born phase space point. In DIS, the construction of F B is straightforward: We require that the mapping preserve the lepton momenta. This choice uniquely determines the kinematics of the corresponding Born configuration, where the momentum of the incoming QCD parton is then set to p = xP and the momentum of the outgoing QCD parton is given by p + q.
In terms of the Born-differential NNLO cross section, any infrared-safe observable O can now be calculated as follows: While the first line of Eq. (2) generates the observable dependence correctly in the Born phase space, the second and third line correct for it's dependence in the single and double emission phase space. They are generated by events in the single and double emission phase space with the appropriate flavor and momentum configurations (corresponding to the O(Φ 1 ) and O(Φ 2 ) terms). For each event, a duplicate event with inverted weight and Born-projected flavorkinematics structure is added (corresponding to the O(F B (Φ 1 )) and O(F B (Φ 2 )) terms). Since the events in the single and double emission phase space correspond to a regular NLO calculation of the jet-associated Born process, any of the well established techniques for the computation of virtual corrections and infrared-subtraction at NLO can be used. In our work, we employ the BLACKHAT library [23][24][25] for the computation of the virtual corrections and Catani-Seymour dipole subtraction [17] as implemented in the matrix element generator AMEGIC [26,27]. Our implementation of the generic NNLO corrections in Eq. (1) is based on the two-loop DIS structure functions available in the literature [2][3][4].
We match the fixed-order computation in the projection-to-Born method to a parton shower using the UNLOPS algorithm [28,29]. The effect of additional emissions generated in the parton shower approach can be described using a generating functional, which is recursively defined for an n-parton final state and the observable O as wheret = t(Φ 1 ), and where the parton-shower no-branching probability is given by Here, dΦ 1 is the differential one-emission phase space, which is parametrized in terms of the evolution and splitting variables t and z as dΦ 1 = dt dz dφ/(2π)J(t, z, φ), with J a possible Jacobian factor. The starting scale of the evolution is given by t n , while t c denotes the cutoff scale. K n is the evolution kernel for the n-parton state. In the case of DGLAP evolution of the DIS process, it can be written as [30,31] where the first term corresponds to initial-state radiation and the second term to final-state radiation. The matching to a fixed-order higher-order calculation is achieved by exploiting the factorization of tree-level matrix elements, which implies schematically that B n+1 → B n K n in the soft or collinear limit. Note that in processes with a more complicated color structure or external gluons, spin and color correlations between the underlying Born configuration, B n , and the splitting kernels in K n must be taken into account. In our simulation of the DIS di-jet topologies, these correlations are included by means of an NLO matching in the S-MC@NLO method [32,33]. The NNLO matching in the UNLOPS method proceeds in two steps. In order to reproduce the logarithmic coefficients of the parton shower resummation, the real emission terms in the fixed-order calculation are reweighted. The nominal accuracy of the fixed-order NNLO calculation is then restored by subtracting the fixed-order expansion of the reweighted result to the second order in the strong coupling such as to remove the overlap with the exact NNLO result. In Eq. (9) we quote only the final formula. The details of the matching procedure are given in [28,29]. The only modifications of the original method that lead to Eq. (9) are due to our NNLO fixed-order input being computed in the projection-to-Born method, rather than the q T -cutoff technique.
We start with the one-jet differential NLO cross sections for standard and hard events as defined in the S-MC@NLO method [32,34], The generating functional for matching the one-jet process at NLO is given in terms of dipole subtraction terms, S 1 , and the underlying Born differential cross sections, B 1 , as The no-emission probability is defined as in Eq. (4) with K 1 → S 1 /B 1 . Events are generated above the parton-shower cutoff scale, t c , below which the DIS process is considered to be inclusive. We introduce a regular and an exceptional part of the S-MC@NLO hard remainder function Exceptional contributions appear in regions of phase space for which no ordered parton shower history can be identified.
The prime example are configurations where the transverse momenta of jets in the Breit frame are much larger than Q 2 , cf. the sketch in Fig. 1(b). The final UNLOPS matching formula at NNLO accuracy reads We have definedB R 1 =B 1 − B 1 and introduced the UNLOPS matching weight, w 1 , which is given by [28] f a (x a ) and f a (x a ) denote the PDFs associated with the external and intermediate parton, respectively, and the constant K = (67/18 − π 2 /6) C A − 10/9 T R n f is the 2-loop cusp anomalous dimension, which restores the physical coupling in soft-gluon emissions [35,36]. The subtraction terms for the no-branching probability of the parton shower, and for the weight w 1 , are given by where b aa = exp{−δ aa K/β 0 }. Finally, the resummation scale in Eq. (9) is given by µ 2 Q = max(t 1 , Q 2 ). The importance of scale choices in higher-order perturbative QCD calculations has been in the focus of interest recently [37,38]. In order to reflect the dynamics of the DIS di-jet and tri-jet processes in the high transverse momentum region, we select a scale similar to the one proposed in [12]. Instead of the jet transverse momenta we employ the total transverse hadronic energy in the Breit frame, H T . The central renormalization and factorization scale in our fixed-order calculations is then given by Equation (12) smoothly interpolates between the regions of normal scale hierarchies, where Q is larger than the transverse momenta of jets in the Breit frame, and the regions of inverted scale hierarchies, where the transverse momenta are much larger than Q. This corresponds to the two situations sketched in Fig. 1. Our scale choice Eq. (12) effectively selects largest scale in the reaction in both cases.

III. RESULTS
In this section we present numerical results of our calculation, quantitative assessments of the theoretical uncertainties, and comparisons to experimental data. We use the publicly available event generation framework SHERPA [18,19], modified to include the changes described in Sec. II. Parton showers are generated using the DIRE model [16]. Comparing the radiation pattern of DIRE in higher-multiplicity events to fixed-order predictions we find that the simulation can be improved by redefining the evolution variables for dipoles with initial-state emitter and final state spectator or vice versa as t = Q 2 u (1 − z)/z, as opposed to the definition t = Q 2 u (1 − z) given in [16]. We use the CT14nnlo PDF set [39] and choose the strong coupling accordingly. Analyses of the simulated events are performed with the help of RIVET [40,41].
In order to validate our implementation of the Born differential NNLO cross section, we compare fixed-order predictions for the reduced cross section to results obtained with the publicly available APFEL library [42]. The reduced cross section is defined as where Y + = (1 + (1 − y) 2 ) and y = Q 2 /(sx) with s the center-of-mass energy of the collider. As shown in Figure 2, where we show a comparison differentially in Q 2 and x, we observe good agreement within the statistical uncertainties. Note that for this cross-check we set the renormalization and factorization scale to µ 2 R/F = Q 2 . Figures 3 and 4 show the inclusive jet cross section, the di-jet cross section and the tri-jet cross section in the high and low Q 2 region as predicted by our simulations in comparison to experimental data from [14,15]. The shaded uncertainty bands display the estimated theoretical uncertainties at fixed order and are obtained from a correlated variation of the renormalization and factorization scale by factors of two up and down. The hatched uncertainty band is a combination of the fixed-order uncertainties and the estimated parton-shower uncertainty, which is obtained by varying the scales at which the strong coupling and the PDF are evaluated in the parton shower by factors of √ 2 up and down, using the technique in [43]. The parton-shower starting scale µ Q is set to Q 2 . Variations of this scale have a minor effect on our predictions, because the dominant hierarchy in the measurement phase space is such that p T > Q. The predictions shown in green are obtained from a fixed-order NLO calculation for inclusive jet production. They have the same fixed-order accuracy as our parton shower matched simulations because the Born configurations of inclusive DIS do not contribute to the observables at fixed order. In the fixed-order predictions, we account for hadronization effects using the correction factors tabulated in [14,15]. They are obtained in the usual fashion, i.e. by comparing leading-order parton shower Monte Carlo simulations before and after hadronization. The corresponding ratios are typically applied as multiplicative corrections to fixed-order calculations at the level of the observables. The predictions shown in red are obtained using our UNLOPS matched calculation supplemented with the string hadronization model [44,45]. We observe good agreement between the matched calculation and the fixed-order prediction, which indicates that the QCD evolution and hadronization effects are well under control. Our results test and confirm, for the first time at this level of theoretical precision, the validity of the approach outlined above, where hadronization corrections are extracted from parton shower Monte Carlo simulations and then multiplicatively applied to the fixed-order calculation.
We quantify the size and uncertainty of the hadronization corrections in our simulations in Fig. 5. We plot the ratio between the particle-level predictions, with hadronization and hadron decays applied, to the prediction obtained after parton showering. We show two different results, one obtained using the Lund string fragmentation model [44,45] Sherpa MC e + p → e + X x B = 3.2 · 10 −4 (×1.1) x B = 5 · 10 −4 (×1.04)  FIG. 2. Comparison of NNLO cross sections, differential in x and Q 2 , between our dedicated implementation in SHERPA and the publicly available APFEL library [42].
as implemented in PYTHIA 6.4 [46] and one obtained using the cluster fragmentation model [47,48] as implemented in SHERPA [49]. The perturbative input to the two simulations is identical, hence we quote their difference as the estimated hadronization uncertainty. We note that the hadronization corrections in our approach agree well with the results computed by H1, which are based on computations using DJANGOH [50] and RAPGAP [51]. These two generators both include the exact expressions for the QCD Compton process at leading order, but DJANGOH uses the Lund dipole cascade model of ARIADNE [52] to simulate higher-order radiative corrections, while RAPGAP is based on collinear parton evolution. Both make use of the Lund string fragmentation model. It is encouraging that in our approach, using both a higher-order perturbative input and a different parton-shower model, the size of the hadronization corrections is very similar. In addition, the hadronization uncertainties are small, except for the very low transverse momentum region in the tri-jet cross section.

IV. CONCLUSIONS
We have presented the first parton-shower matched calculation of hadronic final state production in DIS at NNLO QCD precision. The techniques needed to perform the simulation are implemented in the publicly available program SHERPA and can be used for event generation at the particle level. In contrast to earlier calculations using non-unitary multi-jet merging techniques, we are able to predict both jet production rates and inclusive quantities like structure functions within a single calculation. The agreement with recent analyses by the H1 collaboration is good. Since our calculation is based on collinear factorization, it can provide the basis for the reliable extraction of hadronization corrections needed in the comparison to fixed-order calculations. FIG. 3. Inclusive jet, di-jet and tri-jet cross section differential in Q 2 as a function of pT,j compared to experimental data from the H1 collaboration [14]. We show separate error bars for the reported statistical and systematic uncertainties. We compare NLO fixed-order predictions corrected for hadronization effects (green) and parton-shower matched NNLO predictions at the particle level (red). The light green and light red uncertainty bands are obtained from a correlated variation of the renormalization and factorization scales. The hatched blue uncertainty band combines the fixed-order uncertainties and partonshower uncertainties in quadrature. See the main text for details.
Sherpa MC Inclusive jet selection 5.5 < Q 2 < 8 GeV 2 (×1)  dσ/dp T,j [pb/GeV] Sherpa MC Di-jet selection 5.5 < Q 2 < 8 GeV 2 (×1)   FIG. 4. Inclusive jet cross section differential in Q 2 as a function of pT,j compared to experimental data from the H1 collaboration [15]. See Fig. 3  FIG. 5. Hadronization corrections determined using the UNLOPS simulation of DIS. Predictions from the Lund string model (red) are compared to results from the cluster fragmentation as implemented in SHERPA (blue) and to results computed by the H1 collaboration [14,15].