Semi-inclusive deep-inelastic scattering at NNLO in QCD

Semi-inclusive hadron production processes in deep-inelastic lepton-nucleon scattering are important probes of the quark flavour structure of the nucleon and of the fragmentation dynamics of quarks into hadrons. We compute the full next-to-next-to-leading order (NNLO) QCD corrections to the coefficient functions for semi-inclusive deep-inelastic scattering (SIDIS) in analytical form. The numerical impact of these corrections for precision physics is illustrated by a detailed comparison with data on single inclusive hadron spectra from the CERN COMPASS experiment.


INTRODUCTION
Identified hadron production in hard scattering processes is described in quantum chromodynamics (QCD) through the production of partons (quarks or gluons) which subsequently fragment into hadrons.This partonto-hadron transition is a non-perturbative process that can be parametrized in terms of fragmentation functions (FFs) which describe the probability of a parton fragmenting into a hadron carrying some fraction of its momentum [1,2].These FFs fulfil Altarelli-Parisi evolution equations in their resolution scale [3], which are in complete analogy to the evolution of parton distributions functions (PDFs) in the nucleon.
Production cross sections for a variety of hadron species have been measured in electron-positron, leptonhadron and hadron-hadron collisions.To include these data sets into a global fit of FFs requires the knowledge of the respective parton-level coefficient functions (differential in the momentum of the fragmenting parton) to the desired perturbative order.At present, these coefficient functions are known to next-to-leading order (NLO) for lepton-hadron [4][5][6] and hadron-hadron collisions [7] and to next-to-next-to-leading order (NNLO) for e + e − annihilation [8,9].Consequently, global fits of FFs at NNLO [10][11][12] focus only on data from e + e − experiments, while having to discard any other collider data.
Semi-inclusive hadron production in deep-inelastic lepton-nucleon scattering (SIDIS) has been measured extensively [13][14][15][16][17] for various hadrons.By considering different hadron species, it is possible to single out different flavour combinations of incoming partons, thereby probing the detailed quark and antiquark flavour decomposition of the PDFs.This SIDIS information is largely complementary to inclusive DIS structure function measurements, which allow to determine only a single flavour combination to high accuracy.Moreover, SIDIS measurements play an important role in the determination of spin-dependent PDFs [18] that have to rely on far fewer hadron-collider observables than ordinary (spin-averaged) PDFs.Again, these studies can be performed in a self-consistent manner only up to NLO due to the unavailability of corrections to the SIDIS coefficient functions at higher orders.
It is the purpose of this work to enable precision physics studies with SIDIS observables by deriving the full NNLO QCD corrections to the SIDIS coefficient functions.We provide their analytical expressions for all partonic channels (combinations of initial state and identified final-state partons) and study the impact of the newly derived corrections on a representative data set on SIDIS charged pion production from COMPASS [17].

KINEMATICS OF SIDIS
We consider the observation of a hadron h following the scattering of a lepton on a nucleon.We closely follow the notation of [19], describing semi-inclusive deep-inelastic scattering as ℓ(k) p(P ) → ℓ(k ′ ) h(P h ) X, with X inclusive final-state radiation.The leptons ℓ momenta define the four-momentum q = k−k ′ of the exchanged virtual vector boson and the energy transfer y = (P • q)/(P • k).The usual exclusive variables for describe the momentum fraction of the nucleon carried by the incoming parton (x) and the momentum fraction of the outgoing parton carried by the identified hadron (z).For √ s center-of-mass energy of the lepton-nucleon system we have Q 2 = xys.
As we consider values of Q ≪ M Z only (highly) virtual photons are exchanged, and the triple-differential cross section reads with α denoting the fine structure constant.

arXiv:2401.16281v2 [hep-ph] 29 May 2024
The transverse F h T and longitudinal F h L SIDIS structure functions are given by the sum over all partonic channels of the convolution between the PDF for a parton p (f p ), the FF of a parton p ′ into the hadron h (D h p ′ ), and the coefficient function for the transition p → p ′ (C i p ′ p ): The factorisation theorem that allows the above expression introduces two factorisation scales: µ F for the initial state and µ A for the final state.With µ R we indicate the renormalisation scale.The coefficient functions encode the hard-scattering part of the process and can be computed in perturbative QCD.Their perturbative expansion in the strong coupling constant α s reads ) At LO, only the qq channel (γ * q → q) contributes, with the LO coefficient functions trivially given by where e q is the quark's charge.At NLO instead also the channels qg and gq start to contribute, and results for qq , C i,(1) gq and C i,(1) qg can be found in the literature [19].In this Letter we present results for the NNLO corrections C i, (2) p ′ p to all partonic channels appearing at this order.Following the notation of [20], the 7 partonic channels appearing at O(α 2 s ) are: again for i = T, L. With q ′ (q ′ ) we indicate a quark (antiquark) of flavour different from q, whereas the NS and PS superscripts in the quark-to-quark channel denote the non-singlet and the pure-singlet components respectively.The coefficient functions are computed by applying projectors to extract the longitudinal and transverse components from the respective parton-level subprocess matrix elements with incoming kinematics fixed by Q 2 and x, which are then integrated over the final state phase space.This integration is fully inclusive in the extra radiation X and keeps the final state momentum fraction of the parton p ′ fixed to ẑ.

METHOD
At NNLO in QCD, three types of parton-level contributions must be taken into account, relative to the underlying Born-level process: two-loop virtual corrections (double-virtual, VV), one-loop corrections to single real radiation processes (real-virtual, RV) and tree-level double real radiation processes (RR), with example diagrams shown in Figure 1.These are accompanied by contributions from QCD renormalization and mass factorization of the PDFs and FFs.
NNLO QCD corrections to processes with identified particles at hadron colliders have recently been derived for identified photons [21,22] and for the production of bottom hadrons in top quark decays [23].These calculations are performed in a fully exclusive manner in the form of parton-level event generators which provide the full kinematical information on all final state partons (and on a single identified hadron or photon), which can then be subjected to the precise final state definition that is used in the experiment.To enable these computations, a method to identify and extract infrared singular real radiation up to NNLO had to be employed.The calculation for identified photons used the antenna subtraction method [24][25][26], while the bottom hadron production relies on a sector-improved residue subtraction [27].
The analytic ingredients to the antenna subtraction method are so-called antenna functions, which encapsulate all infrared singular radiation that emerges between two hard radiator partons.These antenna functions are used to construct the real radiation subtraction terms, and they are integrated analytically over the respective antenna phase spaces to make the infrared pole structure explicit.In the case of fragmentation processes [28], one of the hard radiators is the fragmenting parton, while the other can be in the initial or in the final state.The kinematical situation of fragmentation antenna functions with one radiator in the initial state corresponds exactly to the kinematics of SIDIS [29].Consequently, the RV and RR contributions to SIDIS can be obtained employing exactly the same methods as were used to derive the integrated fragmentation antenna functions with one initial state radiator.These methods were described in detail in [21] and we only provide a brief summary here.
The one-loop squared matrix elements for the RV contributions can be expressed in terms of one-loop bubble and box integrals, which are known in exact form in ϵ, with ϵ the dimensional regularisation parameter for N = 4 − 2ϵ space-time dimensions.For fixed x and ẑ, the associated phase space integral is fully constrained, such that only expansions in the end-point distributions in x = 1 and ẑ = 1 are required to obtain the final result for this contribution.To avoid ambiguities associated with the analytic continuation of the one-loop master integrals, the parameter space of the RV contribution is segmented into four sectors: (x ≤ 0.5, x ≤ ẑ ≤ 1 − x), (ẑ ≤ 0.5, ẑ < x ≤ 1 − ẑ), (x > 0.5, 1 − x < ẑ < x) and (ẑ > 0.5, 1 − ẑ < x ≤ ẑ), where manifestly real-valued expressions for the contributions are obtained [21].The expressions are continuous across the boundaries of the regions.
The RR contributions correspond to integrations over a three-particle phase space, with the momentum fraction of one of the particles fixed by ẑ.They can be expressed as cuts of two-loop integrals in forward kinematics, with ẑ expressed as a linear cut propagator.These integrals are reduced to master integrals using integration-by-parts (IBP) identities [30,31], as implemented in Reduze2 [32].The RR contributions to the SIDIS coefficient functions are expressed in terms of 13 integral families, which contain a total of 21 master integrals.These master integrals are determined by solving their differential equations [33] in x and ẑ, using PolyLogTools [34] and HPL [35].The boundary terms for these differential equations are obtained by integrating the generic solutions over ẑ and comparing to the master integrals relevant to inclusive integrated antenna functions [36] with one initial-state and one final-state radiator.Of the 21 master integrals, 9 were already computed in the context of photon fragmentation at NNLO [21], and derivation of the remaining 12 integrals will be described in detail elsewhere [37].
The VV contributions correspond to the well-known two-loop quark form factor [38] in space-like kinematics.
All contributions are computed using FORM [39] and are assembled to yield the bare SIDIS coefficient functions, which still contain ultraviolet and collinear pole terms.By adding the renormalization and mass factorization counterterms (including convolutions of lowerorder terms using MT [40]), the finite physical SIDIS coefficient functions are obtained.

RESULTS
The results for the full set of coefficient functions up to NNLO order are too lengthy to be presented here and thus are given as an ancillary file attached to the arXiv submission of this Letter.Our results include the full scale dependence (µ R , µ F , and µ A ) that was crosschecked with the solution of the renormalisation group equation for all channels.In this section we discuss the comparison of our results with the literature and numerical results.
By means of the threshold resummation formalism for SIDIS [19], approximate corrections for the qq channel have been derived at NNLO [41], and even at N 3 LO [42].Such approximate NNLO corrections have been adopted in the context of a global QCD analysis of light fragmentation functions [43,44].Moreover, partial results for the qq NNLO longitudinal coefficient function are also available [20].Most recently the leading colour contribution to the qq non-singlet channel was computed in [45].
Concerning the qq channel, we compare our results against the ones of [45].The longitudinal components are in perfect numerical agreement.Regarding the leading color transverse ones, we find analytical agreement for all terms involving endpoint distributions as well as perfect numerical agreement for the regular part in the region (ẑ ≤ 0.5, ẑ < x ≤ 1 − ẑ).We are also in agreement with the threshold expansion terms of [41], which predict all double distributions in the partonic variables, and have been confirmed by [45] as well.Figure 2 illustrates the numerical impact of the newly computed NNLO corrections and assesses the relevance of different partonic channels.Using selected kinematical bins from the COMPASS measurement [17] of SIDIS pion production in muon-nucleon scattering, which is described in more detail in the following section, we compute the K-factors at NLO and NNLO and decompose the cross sections according to different channels.We use the NNPDF3.1 PDF set [46] and the FF set from [43] at NNLO throughout, with α s (M Z ) = 0.118 and with N F = 5 light quarks.The central scales are fixed at µ R = µ F = µ A = Q, with scale variations determined through variations by a factor 2 around the central scale.We further fix µ F = µ A .
We observe moderate NNLO corrections to the Kfactors, which reinforce the tendency of the NLO corrections of an increase of the K-factor with increasing  2. QCD K-factors up to NNLO and fractional contribution of individual channels for selected kinematical bins studied by the COMPASS experiment [17].The g → g, q → q, q → q ′ and q → q′ channels are not shown in the channel decomposition as they are found to give negligible contributions in the kinematical bins considered.
z.The non-uniformity of the NNLO corrections in x and z clearly highlights the phenomenological relevance of the NNLO contributions.In the smallest x-bin (corresponding to the lowest Q avg ), NNLO corrections are somewhat larger, and the overlap of NLO and NNLO uncertainty bands is only marginal.At larger x, the NNLO corrections are within the NLO uncertainty bands and their inclusion leads to considerably smaller uncertainties from 20% at NLO to well below 10% at NNLO.The predictions are largely dominated by the quark-to-quark channel, the gluon-to-quark and quark-to-gluon channels both yield small negative corrections to the SIDIS cross section, especially at small x, with the gluon-to-quark channel being typically larger due to the larger magnitude of the respective fragmentation function.All new channels appearing at NNLO are found to give negligible contributions.

COMPARISON WITH DATA
The COMPASS experiment performs deep-inelastic scattering measurements on various fixed targets with a high-energy muon beam at CERN.In their SIDIS study [17], momentum spectra for charged pions and for unidentified charged hadrons are measured with a 160 GeV muon beam scattering off an isoscalar target, corresponding to a center of mass energy √ s ≈ 17.35 GeV.Events are accepted if Q 2 > 1 GeV 2 and W > 5 GeV with W = (P + q) 2 the invariant mass of the hadronic system.
The measured hadron multiplicities dM h /dz are given by the ratio of the differential cross section for hadron production and the differential inclusive DIS cross sec-tion.Therefore we compute the ratio For brevity, we focus on the h = π + spectra.In our numerical implementation we compute the denominator of ( 7) using the APFEL++ code [47,48].We apply the same experimental cuts in our numerical implementation and we integrate (7) over x and y, according to the given bin ranges.In Figure 3 we present the ratio of data and theory predictions over the NLO result.The uncertainty on theory predictions is estimated by varying the scales in an uncorrelated way between the numerator and denominator of (7).We observe that inclusion of the NNLO corrections modifies the shape of the predictions, in general improving the description of the experimental data.For the lowest values of Q avg ≤ 2 GeV, no reduction of the scale uncertainty is observed.Moving to higher Q avg , this reduction becomes clearly significant, with NNLO uncertainties usually being half the size of their NLO counterparts.

CONCLUSIONS
Semi-inclusive deep inelastic scattering processes will be among the key observables of the physics program at the BNL Electron-Ion Collider (EIC).To enable precision studies with SIDIS data, higher order perturbative corrections are crucial.To prepare the precision SIDIS program at EIC, we derived the analytical expressions for the NNLO QCD corrections to the SIDIS coefficient functions.The NNLO corrections are non-uniform in  the kinematical variables.They lead to a substantial reduction of the uncertainty on the theory predictions at sufficiently large values of Q, where the perturbative expansion is applicable.In comparison with COMPASS results on π + SIDIS production, we observe an improved description of the experimental data.
Our newly derived results allow precision determinations of the quark flavour decomposition of nucleon PDFs and of hadron FFs in SIDIS at the EIC.A natural extension of our work could be towards the polarized SIDIS coefficient functions, thereby enabling precision SIDIS studies in the EIC spin physics program.

FIG. 1 .
FIG. 1. Example Feynman diagrams contributing to C i gq at the RV level (top), and to C i,NS qq and C i,PS qq at the RR level (bottom left and right).