N3LO Corrections to Jet Production in Deep Inelastic Scattering using the Projection-to-Born Method

Computations of higher-order QCD corrections for processes with exclusive final states require a subtraction method for real-radiation contributions. We present the first-ever generalisation of a subtraction method for third-order (N3LO) QCD corrections. The Projection-to-Born method is used to combine inclusive N3LO coefficient functions with an exclusive second-order (NNLO) calculation for a final state with an extra jet. The input requirements, advantages, and potential applications of the method are discussed, and validations at lower orders are performed. As a test case, we compute the N3LO corrections to kinematical distributions and production rates for single-jet production in deep inelastic scattering in the laboratory frame, and compare them with data from the ZEUS experiment at HERA. The corrections are small in the central rapidity region, where they stabilize the predictions to sub per-cent level. The corrections increase substantially towards forward rapidity where large logarithmic effects are expected, thereby yielding an improved description of the data in this region.


Introduction
Collider experiments have provided a wealth of precision measurements of basic lowmultiplicity production processes and scattering reactions in particle physics. An equally high level of accuracy in the theoretical predictions is required to turn the experimental data into highly accurate determinations of fundamental parameters (e.g., coupling constants and particle masses) or to use them in indirect searches for new-physics effects. The necessary theoretical precision can be obtained by computing the relevant scattering cross sections to a sufficiently high order in perturbation theory.
In this context, one distinguishes predictions for exclusive (sometimes also called fully differential or fiducial, depending on the specific application) cross sections and inclusive cross sections. Fiducial cross sections take account of the kinematical coverage and finalstate reconstruction procedure of the experimental measurement, and can be compared directly with data. Inclusive cross sections are the result of an extrapolation to full kinematical coverage. This extrapolation is usually performed as part of the experimental analysis; it does however require detailed modelling and theory input, thereby introducing additional sources of uncertainty. Wherever possible, experimental measurements of precision observables at the LHC have started to shift their focus towards fiducial cross sections.
The computation of higher-order corrections differs in important technical details between inclusive and exclusive cross sections. In inclusive cross sections, kinematical information on individual final-state products can be fully integrated over, thereby considerably reducing the number of independent scales in the problem under consideration. Results for higher-order corrections for inclusive cross sections can be cast in the form of coefficient functions, which can often be obtained in closed analytical form. In contrast, predictions for exclusive cross sections need to keep track of the full final-state information in all subprocesses relevant at a given order. This is usually realised in the form of a parton-level event generator, which applies the experimental event reconstruction and kinematical cuts (collectively called measurement function) to all subprocesses in order to reconstruct the fully exclusive differential cross section. Owing to these important technical differences, inclusive cross sections can often be computed to a higher perturbative order than exclusive cross sections.
QCD corrections for inclusive cross sections are available to the fourth order (N 4 LO) for e + e − → hadrons [1], and to the third order (N 3 LO) for deep inelastic scattering [2] and Higgs-boson production (integrated over rapidity) in gluon fusion at hadron colliders [3]. Building upon the results for deep inelastic scattering (DIS), corrections to this order have also been inferred for Higgs boson production in vector boson fusion [4]. In processes with hadrons in the initial state, these results for the partonic coefficient functions in principle require parton distributions accurate to N 3 LO, which will be enabled by the ongoing progress in the calculation of the Altarelli-Parisi splitting functions to this order [5].
Infrared singular contributions appear in two forms: either as explicit poles from virtual loop corrections or in the form of implicit poles from real-radiation corrections, which turn into explicit poles only after integration over the phase space associated to the real radiation. An infrared subtraction method extracts the infrared poles from singular real radiation contributions, optimally in a process-independent manner, and allows them to be cancelled against the virtual corrections. The currently available methods require a varying level of preparation and partly allow the re-use of results from lower-order calculations in order to obtain NNLO corrections to exclusive cross sections.
In this respect, the Projection-to-Born (P2B) method [31] is very efficient in repurposing already available calculations: to compute the NNLO corrections to fully differential exclusive cross sections related to a final state X, it combines the next-to-leading order (NLO) calculation for differential cross sections for X + j final states with the NNLO corrections to the inclusive cross section for final state X. In its practical implementation, it requires an extension of the NLO X + j calculation, in which the kinematics of each X + j final state is projected to an equivalent Born kinematics to construct a subtraction term that is subsequently used to regulate the singular behaviour present when the jet becomes unresolved. This method is applicable to all processes where NNLO corrections to the fully inclusive cross section are known differentially in the Born-level kinematics, i.e., the production of a vector boson [39] or a Higgs boson, the latter both in gluon fusion [40] and in vector boson fusion [41], at hadron colliders as well as deep inelastic lepton-hadron scattering [42]. Up to now, it has been applied in the calculation of NNLO corrections to Higgs-plus-two-jet production in vector boson fusion [31].
There has been substantial recent progress in calculating N 3 LO corrections to inclusive cross sections [3,4], and NNLO calculations are now becoming available for fully differential cross sections involving jet final states. By combining these two recent developments, we are now able to extend the P2B method to compute N 3 LO corrections to fully differential cross sections. This is the main purpose of this paper which is organized as follows. Section 2 summarises the different parton-level contributions to cross sections up to N 3 LO in QCD, and describes how the P2B method is applied to cancel their infrared singularities. As a proof-of-principle application of the method, we compute the N 3 LO corrections to single-jet production in the laboratory frame of deep inelastic scattering in Section 3 and compare our predictions to data from the HERA electron-proton collider in Section 4. Finally, Section 5 contains our conclusions.

Method
Infrared singularities from real-radiation contributions at higher orders in perturbation theory arise from phase-space regions where one or more of the final-state particles become soft and/or collinear. These singularities appear only upon integration over the final-state phase space, and can be extracted from the real-radiation contributions by using an infrared subtraction method. The currently available methods at NNLO can be broadly classified in two types, depending on whether the divergent phase-space integral is regulated by applying cuts to prevent it to encompass the singular regions (cut-based or slicing: [37,38]), or whether it is regulated by introducing counter-terms that render the integrand finite in the singular regions (counter-term-based or subtraction: [31][32][33][34][35][36]). In the latter methods, the construction of counter-terms typically exploits the known factorisation properties of the phase space and the QCD matrix elements in all singular regions.
The P2B method [31] is the simplest possible incarnation of a counter-term-based method. The counter-term is given by the full matrix element itself, which is evaluated at its original phase-space point. The only difference between the real-radiation contribution and its counter-term is in the measurement function: the real-radiation contribution uses the actual measurement function that defines the exclusive cross section, while the measurement function of the counter-term is unity everywhere in phase space (i.e., projected to the Born-level kinematics of the leading-order process), corresponding to a fully inclusive cross section. This method can be applied provided that the complete kinematics of the inclusive cross section for a final state X can be inferred from the momenta of non-QCD particles (i.e., not requiring a recombination or clustering of momenta); this restricts its application in principle to the production of one or more colourless particles at hadron colliders and to DIS processes.
For processes at hadron-hadron colliders, it should be kept in mind that the kinematics of the inclusive process is described by two variables: the mass of the final state X and its rapidity, and that the inclusive cross section is thus not to be confused with the total cross section (which is integrated over rapidity). Conceptually, the P2B method can be extended for these processes to any order in perturbation theory, provided the ingredients (the inclusive cross section at the desired order, and the fully differential cross section with an extra jet at one order below) are available: Here dO abbreviates the kinematical definition (usually multiply differential) of an infraredsafe observable, defined using the actual event-kinematics, while dO B is the limiting value of dO if evaluated for the Born-level production process of X. The essence of the P2B method is to define a kinematical mapping that uniquely assigns It should however be noted that the P2B method only accounts for the infrared subtraction of the most singular parts, i.e. for those contributions that turn from implicit to explicit poles when integrating out the last remaining parton. All other infrared cancellations have taken place already within the construction of the exclusive cross section with an extra jet at the previous order, dσ N k−1 LO X+j /dO. In this part of the calculation, a different infrared subtraction method must be applied, since the kinematics of the process with an extra jet typically does not allow one to define a fully inclusive cross section. In the application of the P2B method to the NNLO corrections to Higgs-plus-two-jet production in vector boson fusion [31], the NLO corrections to Higgs-plus-three-jet production were taken from an existing calculation [43] based on the dipole subtraction method [44].
Fully differential NNLO corrections are now becoming available for a substantial number of jet production processes. These calculations can be used as ingredients to P2B calculations at N 3 LO accuracy, provided the corresponding inclusive cross sections are known to this order. To understand how the P2B method is implemented at this order, we note that the N 3 LO cross section for the production of a final state X (which consists of n particles at Born level) is assembled as follows: where dσ ABC X denotes the parton-level contributions to the cross section from tree-level triple-real radiation (RRR), from double-real radiation at one loop (RRV ), from singlereal radiation at two loops (RV V ) and from the three-loop virtual corrections to the Born process (V V V ). The virtual loop contributions are ultraviolet-renormalised. Processes with incoming partons also contain mass factorisation counter-terms; these are implicitly contained in the virtual loop corrections in the above formula. The mass factorisation counter-terms up to N 3 LO are constructed from the three-loop splitting functions [45]. The function J(O i ) defines the observable under consideration for an i-particle final state, with the Born-level definition J(O B ) ≡ J(O n ), and Φ i denotes the i-particle phase space integration. Each term in the above sum is separately infrared divergent, with infrared singularities arising from the virtual loop integrations and from the phase-space integrations over unresolved particles. The inclusive cross section depends only on the n-particle Bornlevel kinematics, and its N 3 LO contribution can be assembled as The parton-level processes dσ RRR, RRV, RV V X also contribute to the NNLO corrections to the (fully differential) X + j production. Provided that a jet j is resolved in the final state, the infrared cancellations among these three contributions can be accomplished by an NNLO subtraction method. Using the antenna subtraction method [33,34,46], the NNLO cross section schematically reads [34] dσ NNLO where dσ S, T, U X+j denote the subtraction terms, constructed from antenna functions and reduced matrix elements for the three parton-level contributions. Using the factorisation properties of the multi-particle phase space into a phase space of lower multiplicity and an antenna phase space [46], By construction, the integrations do not depend on the definition J of the observable, and can thus be carried out for all antenna functions [33,47,48] in a process-independent and observable-independent manner.
In the N 3 LO contribution to the exclusive cross section for final state X, the final state jet j in Eq. (2.5) is no longer guaranteed to be resolved, but can become soft and/or collinear, thereby resulting in an infrared divergence upon phase-space integration. At this point, the P2B subtraction sets in: For completeness, we also summarise the structure of the NLO and NNLO corrections [31] for the final state X in the P2B method, where the NNLO corrections take the NLO calculation of X + j production using antenna subtraction as an ingredient: The contributions dσ A X denote the corrections to the Born-level process from single-real radiation (A = R), double-real radiation (A = RR) and real-virtual corrections (A = RV ), and mass-factorisation contributions are again implicitly included with the virtual corrections.

Application to jet production in deep inelastic scattering
The basic parton-level process in deep inelastic lepton-proton scattering (DIS) is the elastic scattering at large momentum transfer of the lepton and a quark, mediated by a virtual photon, W-or Z-boson. The outgoing quark forms a jet, of which the momentum can be inferred using momentum conservation using the initial state momenta and the outgoing electron momentum. Viewed in the laboratory frame of the lepton-proton system (which can be either of fixed-target or collider type), this jet always has a non-vanishing transverse momentum and a finite rapidity. The inclusive single-jet cross section (at fixed kinematics of the jet) at Born level is therefore identical to the inclusive structure function (at fixed kinematics of the lepton).
Jet production in DIS has been measured [49] in the laboratory frame [50][51][52][53][54] and in the Breit frame [55,56]. In the latter, the virtual photon and the proton collide head-on and the basic lepton-quark scattering process always yields a quark at zero transverse momentum. Consequently, the first non-trivial jet production process in the Breit frame is two-jet production. Owing to the higher final state multiplicity, this process has a higher sensitivity to QCD dynamics; the vast majority of jet production studies at HERA were therefore performed in the Breit frame [55,56].
Next-to-leading order QCD corrections to jet production in DIS have been available for a long time for single-jet production in the laboratory frame [57] as well as for two-jet [58] and three-jet production [59] in arbitrary frames. Very recently, the calculation of NNLO QCD corrections to two-jet production in DIS [29] was completed. This calculation uses the antenna subtraction method [33,34] and is implemented in a flexible parton-level event generator (NNLOjet), which allows to compute any infrared-safe observable derived from the process under consideration. While [29] and its subsequent application to a precision determination of the strong coupling constant [60] discuss only jet production observables in the Breit frame, the very same NNLOjet implementation can be applied to compute two-jet production in DIS in the laboratory frame.
The Born level kinematics of single-jet production in DIS in this frame can be reconstructed entirely from the lepton kinematics: denoting the incoming and outgoing lepton momenta by p a and p 1 (with q = p a − p 1 ), and the incoming proton momentum by P , the incoming and outgoing quark momenta are determined by p b = xP and p 2 = xP − q with x = −q 2 /(2q · P ). The Born-level value of any observable O B ≡ O 2 derived from single-jet production in DIS is therefore expressed in terms of the momenta (p a , p 1 , P ). As a result, we can define a mapping which, starting from any momentum set of final-state multiplicity m > 2, maps onto to this Born kinematics: {p i } m → {p i,B } 2 . It is then given by With this mapping at hand and the identification X = 1j and X +j = 2j, Eqs. (2.9)-(2.11) provide the N 3 LO corrections to single-jet production in DIS in the laboratory frame, based on the existing NNLO calculation of two-jet production processes [29] and the inclusive N 3 LO DIS structure function [2]. It should be noted that this reasoning exploits the specific constraints on the Born-level kinematics in the laboratory frame, and that the P2B method can not be applied in the Breit frame.
In order to validate our implementation of the P2B method up to NNLO, we have performed an independent calculation of single jet production in the laboratory frame based entirely on the antenna subtraction method. In Fig. 1, we compare results from both methods for the inclusive jet pseudorapidity η jet , the inclusive jet transverse energy E jet T , the momentum transfer Q 2 = −q 2 , and the Bjorken scaling variable x. The setup of the calculation follows the ZEUS measurement [54] described in Sect. 4.1 below. The bottom panels in each plot display the ratios between the P2B and antenna subtraction method for the NLO and NNLO coefficients and the central choice µ R = µ F = Q of the renormalization and mass-factorization scales. The error bars correspond to the numerical integration error of the P2B prediction, while the filled area shows the corresponding error of the antenna-subtraction results. The dashed lines are for illustrative purposes and show ±1% for NLO-only and ±5% for NNLO-only.
We observe agreement between the two methods for the O(α s ) NLO coefficient, which is well below 1% across almost the entire kinematic range and fully consistent within the respective Monte Carlo errors. For the O(α 2 s ) NNLO coefficient, the statistical uncertainties are typically below 1% in regions which contribute the bulk of the cross section; they can increase to a few percent in the tails of the distributions. Again, we observe excellent agreement between the two independent calculations. With respect to the full NNLO prediction, the agreement between the two methods was validated at the sub-permille level. A similar level of agreement is also observed for the other scale settings of the seven-point scale variation described below.

N 3 LO results and comparison to HERA data
The first observation of jet production in DIS was made by the fixed-target E665 experiment [50]. Shortly thereafter, the experiments H1 and ZEUS at the electron-proton collider HERA embarked on a large program of jet production studies in deep inelastic scattering, with measurements both in the laboratory frame [51][52][53][54] and in the Breit frame [55,56]. Early measurements established the jet production process and determined total jet rates as function of the resolution parameter [51][52][53]. Subsequent studies on a larger dataset led to the determination of differential distributions of jets [54][55][56] in the kinematical variables. The resulting HERA legacy dataset provides important constraints on QCD dynamics and the parton distributions of the proton.
In the following, we will compare the N 3 LO results to single-jet measurements in the laboratory frame that were performed by the ZEUS collaboration for differential distributions [54] and jet rates [53]. The theory uncertainties on the predictions are obtained from a seven-point scale variation of renormalisation and factorisation scales around a central value of µ 2 F = µ 2 R = Q 2 , independently varying µ R and µ F up and down by factors [1/2 , 2] and discarding the two combinations with the largest separation of both scales. Parton distribution functions fitted at N 3 LO accuracy are not yet available; the results presented below are obtained using the NNLO NNPDF3.1 set [61] with α s (M Z ) = 0.118. The same set is also used to evaluate the LO and NLO expressions.

Differential distributions
The ZEUS measurement [54] of differential distributions for jet production in the laboratory frame is based on data that were taken with a proton beam energy E p = 820 GeV and an electron beam energy E e = 27.5 GeV. We compare our results to the measurement performed in the 'global' region 1 defined by ZEUS through the fiducial cuts where E e denotes the energy of the outgoing electron. Jets are reconstructed using the k T clustering algorithm [62] in the longitudinally invariant mode (E T -weighted recombination scheme) [63] and are required to satisfy  [54] for the inclusive jet pseudorapidity η jet , the inclusive jet transverse energy E jet T , the momentum transfer Q 2 , and the Bjorken scaling variable x. We observe that, for the first time, the scale-uncertainty bands overlap across the full kinematic range, when going from NNLO to N 3 LO. The inclusion of the N 3 LO corrections further reduces the scale uncertainties, by typically a factor of two or more.
The low-x region as well as the first Q 2 bins are kinematically suppressed at LO, hence the perturbative accuracy is effectively reduced. As a consequence, we observe larger higher-order corrections with residual N 3 LO scale uncertainties at the level of about ±10% in this region. Large perturbative corrections and correspondingly large uncertainties are also observed in the forward region, η jet 1. In this region, the LO process is again suppressed kinematically, and most of the jet activity here arises from final states containing several jets. Forward jet production in DIS has been studied extensively [49] with a view on establishing large logarithmic corrections arising from the high-energy limit [64] that potentially require resummation. We observe that in the forward region these corrections are built up by including consecutive perturbative orders. The overlap of the scale uncertainty bands at NNLO and N 3 LO indicates a stabilization of the expansion at the present order.
Overall, the N 3 LO QCD predictions provide an excellent description of the ZEUS data. An improvement over the NNLO description is observed in particular in those kinematical regions where higher-order corrections are large (low x, low Q 2 , forward region: η jet 1): Here the N 3 LO corrections induce changes to the shape which bring the central predictions in line with the measurements. The shape of the E jet T distribution, especially in the region 8-20 GeV, is also affected by higher-order corrections; however, its experimental accuracy is systematically limited by relatively large jet-energy-scale uncertainties.

Jet rates
Earlier ZEUS measurements [52,53], based on a smaller dataset taken with E p = 820 GeV and E e = 26.7 GeV, determined the jet production rates, i.e., the fraction of events with a certain jet multiplicity. These studies applied the JADE clustering algorithm [65] with the four-momentum recombination scheme. The jet rates were measured as a function of the JADE clustering parameter y cut in the fiducial region defined by [53]: 160 < Q 2 < 1280 GeV 2 , 0.04 < y < 0.95 , 0.01 < x < 0.1 . In this particular measurement, the two-jet rate is defined as where N 1+1 and N 2+1 are the number of recorded (1 + 1)-and (2 + 1)-jet events, with the extra "+1" denoting the proton remnant forming the beam jet. The normalisation is chosen such that  Figure 3. Comparison of ZEUS data to theoretical predictions up to N 3 LO in QCD for the two-jet rate, normalised according to Eq. (4.5). Data are corrected to give parton level results and were extracted from Fig. 9 of Ref. [53]; the error bars show the statistical uncertainties.
In Fig. 3 we present the results for the 2-jet rate R ZEUS (2+1) up to N 3 LO and compare them to the measurement in Ref. [53]. We refrain from showing the one-jet rate separately, as it is trivially related to R ZEUS (2+1) via Eq. (4.5). It should be noted that the errors on the data [53] correspond only to the statistical uncertainties. Systematic uncertainties from jet acceptance corrections (which amount to up to 20%, and whose uncertainty is not quantified) are not provided on a bin-by-bin basis in [53].
It can be seen that the N 3 LO corrections result in a substantial reduction of scale uncertainties. As the value of y cut is lowered, less of the final-state radiation is clustered, thereby resolving more jet structure. As a result, the fractions of two-jet events increases towards lower values of y cut . The low-y cut region is also where the largest scale dependence is observed. At N 3 LO the scale band starts to overlap in the low-y cut region with that of the previous order for the first time. This is accompanied with a substantial reduction of scale uncertainties, and the y cut -dependence of the data is well-described by the predictions at NNLO and N 3 LO.  Using the same kinematical cuts, we have also calculated the jet rates normalising to the total hadronic cross section, i.e. not applying the unusual normalisation in Eq. (4.5). In this setting, the zero-, three-, and four-jet rates are also well-defined and can be evaluated, as shown in Fig. 4. Similar features as in Fig. 3 are observed concerning the reduction of the scale uncertainty and the y cut -dependence. Moreover, it can be seen that the fractions of one-and two-jet events dominate over events with higher jet multiplicities by at least an order of magnitude.

Conclusions
In this paper, we demonstrated how the P2B method for infrared subtractions can be applied to compute fully differential third-order (N 3 LO) QCD corrections to observables with sufficiently simple Born-level kinematics. The implementation of the method for a given process X requires the knowledge of the fully inclusive N 3 LO cross section for the production of X, and a fully differential NNLO calculation for X+jet final states.
As a first application, we have computed the N 3 LO corrections to jet production in deep inelastic lepton-proton scattering. Predictions at this order lead to very small residual scale uncertainties (often at sub per-cent level in the bulk of the phase space) and account properly for enhancements of distributions from multiple radiation near boundaries of the phase space. The N 3 LO results provide an excellent description of data on distributions and rates in jet production, obtained by the ZEUS experiment [53,54], also in those kinematical regions where the NNLO predictions were insufficient to describe their behaviour.
The P2B method could be used in the near future to evaluate more processes at N 3 LO accuracy in particular in proton-proton collisions.