Single-inclusive jet production in electron-nucleon collisions through next-to-next-to-leading order in perturbative QCD

We compute the ${\cal O}(\alpha^2\alpha_s^2)$ perturbative corrections to inclusive jet production in electron-nucleon collisions. This process is of particular interest to the physics program of a future Electron Ion Collider (EIC). We include all relevant partonic processes, including deep-inelastic scattering contributions, photon-initiated corrections, and parton-parton scattering terms that first appear at this order. Upon integration over the final-state hadronic phase space we validate our results for the deep-inelastic corrections against the known next-to-next-to-leading order (NNLO) structure functions. Our calculation uses the $N$-jettiness subtraction scheme for performing higher-order computations, and allows for a completely differential description of the deep-inelastic scattering process. We describe the application of this method to inclusive jet production in detail, and present phenomenological results for the proposed EIC. The NNLO corrections have a non-trivial dependence on the jet kinematics and arise from an intricate interplay between all contributing partonic channels.


INTRODUCTION
The production of hadrons and jets at a future Electron Ion Collider (EIC) will play a central role in understanding the structure of the protons and nuclei which comprise the visible matter in the universe. Measurements of inclusive jet and hadron production with transversely polarized protons probe novel phenomena within the proton such as the Sivers function [1], and address fundamental questions concerning the validity of QCD factorization. Event shapes in jet production can give insight into the nuclear medium and its effect on particle propagation [2]. The precision study of these processes at a future EIC will provide a much sharper image of proton and nucleus structure than is currently available. Progress is needed on both the experimental and theoretical fronts to achieve this goal. Currently, much of our knowledge of proton spin phenomena, such as the global fit to helicitydependent structure functions [3], comes from comparison to predictions at the next-to-leading order (NLO) in the strong coupling constant. Theoretical predictions at the NLO level for jet and hadron production in DIS suffer from large theoretical uncertainties from uncalculated higher-order QCD corrections [4] that will eventually hinder the precision determination of proton structure. In some cases even NLO is unknown, and an LO analysis fails to describe the available data [5]. Given the high luminosity and expected precision possible with an EIC, it is desirable to extend the theoretical precision beyond what is currently available. For many observables, a prediction to next-to-next-to-leading order (NNLO) in the perturbative QCD expansion will ultimately be needed.
An important step toward improving the achievable precision for jet production in electron-nucleon collisions was taken in Ref. [4], where the full NLO O(α 2 α s ) correc-tions to unpolarized lN → jX and lN → hX scattering were obtained. Focusing on single-inclusive jet production for this discussion, it was pointed out that two distinct processes contribute: the deep-inelastic scattering (DIS) process lN → ljX, where the final-state lepton is resolved, and γN → jX, where the initial photon is almost on-shell and the final-state lepton is emitted collinear to the initial-state beam direction. Both processes were found to contribute for expected EIC parameters, and the shift of the leading-order prediction was found to be both large and dependent on the final-state jet kinematics.
Our goal in this manuscript is to present the full O(α 2 α 2 s ) NNLO contributions to single-inclusive jet production in electron-nucleon collisions, including all the relevant partonic processes discussed above. Achieving NNLO precision for jet and hadron production is a formidable task. The relevant Feynman diagrams which give rise to the NNLO corrections consist of two-loop virtual corrections, one-loop real-emission diagrams, and double-real emission contributions. Since these three pieces are separately infrared divergent, some way of regularizing and canceling these divergences must be found. However, theoretical techniques for achieving this cancellation in the presence of final-state jets have seen great recent progress. The introduction of the Njettiness subtraction scheme for higher order QCD calculations [6,7] has lead to the first complete NNLO descriptions of jet production processes in hadronic collisions. During the past year several NNLO predictions for processes with final-state jets have become available due to this theoretical advance [6,[8][9][10][11][12]. In some cases the NNLO corrections were critical in explaining the observed data [12]. We discuss here the application of the N -jettiness subtraction scheme to inclusive jet produc-tion in electron-proton collisions. Our result includes both the DIS and photon-initiated contributions, and allows arbitrary selection cuts to be imposed on the final state. Upon integration of the DIS terms over the finalstate hadronic phase space we compare our result against the known NNLO prediction for the inclusive structure function, and we find complete agreement. We present phenomenological results for proposed EIC parameters. We find that all partonic channels, including new ones that first appear at this order, contribute in a non-trivial way to give the complete NNLO correction. We note that the NNLO corrections to similar processes, massive charm-quark production in deeply inelastic scattering and dijet production, were recently obtained [13,14].

LOWER-ORDER RESULTS
We begin by discussing our notation for the hadronic and partonic cross sections, and outlining the expressions for the LO and NLO cross sections. We will express the hadronic cross section in the following notation: where the ellipsis denotes neglected higher-order terms. The LO subscript refers to the O(α 2 ) term, the NLO subscript denotes the O(α 2 α s ) correction, while the NNLO subscript indicates the O(α 2 α 2 s ) contribution. For the partonic cross sections, we introduce superscripts that denote the powers of both α and α s that appear. For example, the leading quark-lepton scattering process is expanded as Here, the dσ indicates the O(α 2 α s ) term. The leading-order hadronic cross section can be written as a convolution of parton distribution functions with a partonic cross section, Here, f q/H (ξ 1 ) is the usual parton distribution function that describes the distributions of a quark q in the hadron H carrying a fraction ξ 1 of the hadron momentum. f l/l (ξ 2 ) is the distribution for finding a lepton with momentum fraction ξ 2 inside the original lepton. At leading order this is just f l/l (ξ 2 ) = δ(1 − ξ 2 ), but it is modified at higher orders in the electromagnetic coupling by photon emission. dσ (2,0) ql is the differential partonic cross section. At leading order only the partonic channel q(p 1 ) + l(p 2 ) → q(p 3 ) + l(p 4 ) and the same process with anti-quarks instead contribute. The relevant Feynman diagram is shown in Fig. 1. It is straightforward to obtain these terms. At the next-to-leading order level several new contributions first occur. The quark-lepton scattering channel that appears at LO receives both virtual and realemission corrections that are separately infrared divergent. We use the antennae subtraction method [15] to regularize and cancel these divergences. Initial-state collinear divergences are handled as usual by absorbing them into the PDFs via mass factorization. A gluonlepton scattering channel also contributes at this order. The collinear divergences that appear in these contributions are removed by mass factorization. Example Feynman diagrams for these processes are shown in Fig. 2. Figure 2: Representative Feynman diagrams contributing to the following perturbative QCD corrections at NLO: virtual corrections to the q(p1) + l(p2) → q(p3) + l(p4) process (left); real emission correction q(p1) + l(p2) → q(p3) + l(p4) + g(p5) (middle); the process g(p1) + l(p2) → q ( p3) + l(p4) +q(p5) (right). We have colored the photon line red, the lepton lines green, the gluon lines blue and the quark lines black.
The processes discussed above exhaust the possible NLO contributions when the final-state lepton is observed. However, for single-inclusive jet production a kinematic configuration is allowed where the t-channel photon is nearly on-shell, and the final-state lepton travels down the beam pipe. The transverse momentum of the leading jet is balanced by the additional jet present in these diagrams. This kinematic configuration leads to a QED collinear divergence for vanishing lepton mass, since the photon can become exactly on-shell in this limit. While it is regulated by the lepton mass, it is more convenient to obtain these corrections by introducing a photon distribution function in analogy with the usual parton distribution function. The collinear divergences that ap-pear in the matrix elements computed with vanishing lepton mass can be absorbed into this distribution function, which can be calculated in perturbation theory. Representative diagrams for the photon-initiated processes are shown in Fig. 3. Figure 3: Representative Feynman diagrams contributing to the q(p1) + γ(p2) → q(p3) + g(p4) (left) and g(p1) + γ(p2) → q(p3) +q(p4) scattering processes.
The full expression for the NLO hadronic cross section then takes the form where we have abbreviated f k i/j = f i/j (ξ k ). The contributions dσ (2,1) gl , dσ (2,1) ql and dσ (2,1) ql denote the usual DIS partonic channels computed to NLO in QCD with zero lepton mass. The terms dσq γ,LO and dσ gγ,LO denote the new contributions arising when Q 2 ≈ 0 and the virtual photon is nearly on-shell. The photon distribution can be expressed as (5) where the splitting function is given by This is the well-known Weizsäcker-Williams (WW) distribution for the photon inside of a lepton [16]. The appearance of the renormalization scale µ indicates that an MS subtraction of the QED collinear divergence is used in the calculation of the gl and ql scattering channels, and consequently in the derivation of the photon distribution function.

CALCULATION OF THE NNLO RESULT
The calculation of the full O(α 2 α 2 s ) corrections involves several distinct contributions. The quark-lepton and gluon-lepton scattering channels receive two-loop double virtual corrections, one-loop corrections to single real-emission diagrams, and double-real emission corrections. These contributions necessitate the use of a fullfledged NNLO subtraction scheme. We use the recentlydeveloped N -jettiness subtraction scheme [6,7]. Its application to this process is discussed here in detail. In addition, the photon-initiated scattering channels receive virtual and single real-emission corrections. The calculation of these terms follows the standard application of the antennae subtraction scheme at NLO.
There is in addition a new effect that appears at the NNLO level. The initial-state lepton can emit a photon which splits into a qq pair, all of which are collinear to the initial lepton direction. In the limit of vanishing fermion masses there is a collinear singularity associated with this contribution. This divergence appears in the quark-lepton, gluon-lepton, and photon-initiated scattering channels. It can be absorbed into a distribution function that describes the quark distribution inside a lepton. Treating the collinear singularity in this way leads to new scattering channels that first appear at NNLO: qq → qq, qq → q q , qq → gg, qq → qq , and qg → qg. For our numerical predictions for these channels we need the quark distribution in a lepton. We obtain this by solving the DGLAP equation, which to the order we are working can be written as where the two needed splitting kernels are This expression for the NLO splitting kernel can be obtained from upon replacement of the appropriate QCD couplings with electromagnetic ones. To derive the full quark-in-lepton distribution we use as an initial condition f q/l (ξ, m 2 l ) = 0. Solving Eq. (7) with this initial condition gives With this distribution function it is straightforward to obtain numerical predictions for these partonic channels. Representative Feynman diagrams for these processes are shown in Fig. 4. We can now write down the full result for the O(α 2 α 2 s ) correction to the cross section: The most difficult contribution is the quark-lepton scattering channel. It receives contributions from twoloop virtual corrections (double-virtual), one-loop corrections to single real emission terms (real-virtual), and double-real emission corrections. Sample Feynman diagrams for these corrections are shown in Fig. 5. These are separately infrared divergent, and require a full NNLO subtraction scheme to combine. We apply the N -jettiness subtraction scheme [6,7]. The starting point of this method is the N -jettiness event shape variable [17], defined in the one-jet case of current interest as with Here, p B and p J are light-like four-vectors along the initial-state hadronic beam and final-state jet directions, respectively. 1 The q i denote the four-momenta of all final-state partons. Values of T 1 near zero indicate a final state containing a single narrow energy deposition, while larger values denote a final state containing two or more well-separated energy depositions. A measurement of T 1 is itself of phenomenological interest. It has been proposed as a probe of nuclear properties in electron-ion collisions [2,18], and has also been suggested as a way to precisely determine the strong coupling constant [19]. We will use T 1 to establish the complete O(α 2 α 2 s ) calculation of the quark-lepton scattering channel. Our ability to do so relies on two key observations, as first discussed in Ref. [6] for a general N -jet process.
• Restricting T 1 > 0 removes all doubly-unresolved limits of the quark-lepton matrix elements, for example when the two additional partons that appear in the double-real emission corrections are soft or collinear to the beam or the final-state jet. This can be seen from Eq. (11); if T 1 > 0 then at least one q i must be resolved. Since all doublyunresolved limits are removed, the O(α 2 α 2 s ) correction in this phase space region can be obtained from an NLO calculation of two-jet production in electron-nucleon collisions.
• When T 1 is smaller than any other hard scale in the problem, it can be resummed to all orders in perturbation theory [20,21]. Expansion of this resummation formula to O(α 2 α 2 s ) gives the NNLO correction to the quark-lepton scattering channel for small T 1 .
The path to a full NNLO calculation is now clear. We partition the phase space for the real-virtual and doublereal corrections into regions above and below a cutoff on this definition is dimensionless, unlike the choice in previous applications of N -jettiness subtraction [6].
T 1 , which we label T cut 1 : We have abbreviated θ < 1 = θ(T cut 1 − T 1 ) and θ > 1 = θ(T 1 − T cut 1 ). The first three terms in this expression all have T 1 < T cut 1 , and have been collectively denoted as dσ The remaining two terms have T 1 > T cut 1 , and have been collectively denoted as dσ . The double-virtual corrections necessarily have T 1 = 0. We obtain dσ (2,2) ql (T 1 > T cut 1 ) from a NLO calculation of two-jet production. This is possible since no genuine double-unresolved limit occurs in this phase-space region. We discuss the calculation of dσ (2,2) ql (T 1 < T cut 1 ) using the all-orders resummation of this process in the following sub-section. We note that only the quark-lepton and gluon-lepton scattering channels have support for T 1 = 0. For the other processes there are two final-state jets with non-zero transverse momentum. Such configurations necessarily have T 1 > 0, and therefore only receive contributions from the abovethe-cut phase-space region. We only need these contributions to at most NLO in QCD perturbation theory, which can be obtained via standard techniques. We use antennae subtraction. The only non-standard aspect of this NLO calculation is the appearance of triple-collinear QED limits associated with the emission of a photon and a qq pair which require the use of integrated antennae found in Refs. [22][23][24]. A powerful aspect of the N -jettiness subtraction method is its ability to upgrade existing NLO calculations to NNLO precision. Previous applications of N -jettiness subtraction [6,8,9] have used the NLO dipole subtraction technique [25] to facilitate the calculation of the above-cut phase-space region. This work demonstrates that it can also be used in conjunction with the NLO antennae subtraction scheme.
Below T cut 1 An all-orders resummation of the T 1 event-shape variable in the DIS process for the limit T 1 1 was con-structed in Refs. [18,19]: We have allowed the index q to run over both quarks and anti-quarks. x denotes the usual Bjorken scaling variable for DIS: where P is the initial-state nucleon four-momentum. Φ 2 denotes the Born phase space, which consists of a quark and a lepton. The derivation of this result relies heavily on the machinery of soft-collinear effective theory (SCET) [26]. A summary of the SCET functions that appear in this expression and what they describe is given below.
• H is the hard function that encodes the effect of hard virtual corrections. At leading order in its α s expansion it reduces to the leading-order partonic cross section. At higher orders it also contains the finite contributions of the pure virtual corrections, renormalized using the MS scheme. It depends only on the Born-level kinematics and on the scale choice.
• J q , the quark jet function, describes the effect of radiation collinear to the final-state jet (which for this process is initiated by a quark at LO). It depends on t J , the contribution of final-state collinear radiation to T 1 . It possesses a perturbative expansion in α s .
• S is the soft function that encodes the contributions of soft radiation. It depends on k S , the contribution of soft radiation to T 1 , and has a perturbative expansion in α s .
• B is the beam function that contains the effects of initial-state collinear radiation. It depends on t B , the contribution of initial-state collinear radiation to T 1 . The beam function is non-perturbative; however, up to corrections suppressed by Λ 2 QCD /t B , it can be written as a convolution of perturbative matching coefficients and the usual PDFs: where we have suppressed the scale dependence of the PDF, and i runs over all partons.
The delta function appearing in Eq. (13) combines the contribution of each type of radiation to produce the measured value of T 1 . The ellipsis denotes power corrections that are small as long as we restrict ourselves to the phase-space region T 1 1. The hard, jet, and soft functions as well as the beamfunction matching coefficients all have perturbative expansions in α s that can be obtained from the literature [27][28][29][30]. Upon expansion to O(α 2 s ) and integration over the region T 1 < T cut 1 , Eq. (13) will give exactly the cross section dσ (2,2) ql (T 1 < T cut 1 ) that we require. We can match the beam function onto the PDFs and rewrite the cross section below the cut as We have introduced the schematic notation J q ⊗S ⊗H q ⊗ I qi for the integrations of the SCET functions over T 1 , t J , t B , and k S ; dΦ i Born represents all other terms for the given index i: the parton distribution functions, the integral over the Born phase space, and any measurement function acting on the Born variables. We denote the expansion of these functions in α s as where the superscript denotes the power of α s appearing in each term. With this notation, we need the following contributions to obtain the O(α 2 s ) correction to the cross section below T cut 1 : To simplify this expression, we first note that the hard function has no dependence on the hadronic variables T 1 , t J , t B , and k S . It depends only on the Born phase space and is a multiplicative factor for the hadronic integrations. Next, we note that the leading-order expressions for the SCET functions are proportional to delta functions in their respective hadronic variable: for X = J q , S, or I qi . This simplifies the integrals involving an X (2) . Using the J (2) q term in Eq. (18) as an example, we have where the δ qi comes from the I (0) qi term. Using the fact that the jet function can be written in the form where the a i denote coefficients that can be found in the literature [28], we can immediately derive We have used the standard definition of the plus distribution of a function: This analysis shows how to analytically calculate any term containing one of the NNLO SCET functions. It remains only to calculate contributions containing two NLO SCET functions. We focus on J qi as an example. Using the LO expression for the beam function we can immediately derive the equation The NLO results for the jet and soft functions can be written as Using these expressions it is straightforward to derive the result The vertical bar in the last term indicates that we should take the α m β n coefficient of the series expansion of the bracketed term. We have introduced the abbreviations Using these results it is straightforward to analytically compute all of the necessary hadronic integrals in Eq. (18). The remaining integrals are over the Born phase space and parton distribution functions, and are simple to perform numerically. This completes the calculation of the T 1 < T cut 1 phase space region. We note that the cross section below T cut 1 will contain terms of the form ln n (T cut 1 ), where n ranges from 0 to 4 at NNLO. An important check of this framework is the cancellation of these terms against the identical logarithms that appear for T 1 > T cut 1 . We must also choose T cut 1 small enough to avoid the power corrections in Eq. (13) that go as T cut 1 /Q. Both of these issues will be addressed in our later section on numerical results.
Above T cut 1 We now briefly outline the computation of the quarklepton channel in the phase space region T 1 > T cut 1 . As discussed above this can be obtained from a NLO calculation of the DIS process with an additional jet. In addition to the virtual corrections to the q l → q l g partonic process, there are numerous radiation processes that also contribute: q l → q l g g, q l → q l q q and q l → q l qq. In addition to the usual ultraviolet renormalization of the strong coupling constant, the real and virtual corrections are separately infrared divergent. Remaining divergences after introducing an NLO subtraction scheme are associated with initial-state collinear singularities, and are handled via mass factorization.

VALIDATION AND NUMERICAL RESULTS
We have implemented the NNLO cross section of Eq. (10), as well as the LO and NLO results of Eqs. (3) and (4) in a numerical code DISTRESS that allows for arbitrary cuts to be imposed on the final-state lepton and jets. We describe below the checks we have performed on our calculation.
The antennae subtraction method provides an analytic cancellation of the 1/ poles that appear in an NLO calculation. We are therefore able to check this cancellation of poles for all components computed in this way. This includes the entire σ NLO , as well as the following contributions to the NNLO hadronic cross section: dσ (1,2) gγ and dσ (1,2) qγ . The various contributions dσ (0,2) ij that occur in σ NNLO are finite, and simple to obtain. Our NLO results for the transverse momentum distribution and pseudorapidity distribution of the leading jet are compared against the plots of Ref. [4]. We find good agreement with these results.
This leaves only the validation of the dσ contributions, which both utilize the N -jettiness subtraction technique. There are two primary checks that these pieces must satisfy. First, they must be independent of the parameter T cut N . This checks the implementation of the beam, jet and soft functions, which have logarithmic dependence on this parameter. It also determines the range of T cut N for which the power corrections denoted by the ellipsis in Eq. (13) are negligibly small. Second, upon integration over the final-state hadronic phase space we must reproduce the NNLO structure functions first determined in Ref. [27]. This is an extremely powerful check on our calculation, which essentially cannot be passed if any contribution is implemented incorrectly.
We show in Fig. 6 the results of these checks for the ql and gl scattering channels. We have set the total centerof-mass energy of the lepton-proton collision to 100 GeV. For the purpose of this validation check only we have imposed the phase-space cut Q 2 > 100 GeV 2 , and have integrated inclusively over the hadronic phase space. We have equated the renormalization and factorization scales to the common value µ = Q, and have used the CT14 NNLO parton distribution functions [31]. T cut 1 has been varied from 5 × 10 −6 to 1 × 10 −4 , and the ratio of the NNLO correction to the LO result for the cross section is shown. The solid lines show the prediction of the inclusive structure function. We first note that correction is extremely small, less than 1% of the leading-order result. Nevertheless we have excellent numerical control over the NNLO coefficient, as indicated by the vertical error bars. Our numerical error on the NNLO coefficient is at the percent-level, sufficient for 0.01% precision on the total cross section. The N -jettiness prediction for the ql scattering channel is independent of T cut 1 over the studied range, while the gl scattering channel is independent of T cut 1 for T cut 1 < 10 −5 . Both channels are in excellent agreement with the structure-function result. We have also checked bin-by-bin that the transverse momentum and pseudorapidity distributions of the jet have no dependence on T cut 1 . Having established the validity of our calculation we present phenomenological results for proposed EIC run parameters. We set the collider energy to √ s = 100 GeV and study the inclusive-jet transverse momentum and pseudorapidity distributions in the range p jet T > 5 GeV and |η jet | < 2. We use the CT14 parton distribution set [31] extracted to NNLO in QCD perturbation theory. We reconstruct jets using the anti-k t algorithm [32] with radius parameter R = 0.5. Our central scale choice for both the renormalization and factorization scales is µ = p jet T . To estimate the theoretical errors from missing higher-order corrections we vary the scale around its central value by a factor of two. The transverse momentum distributions at LO, NLO and NNLO are shown in Fig. 7. The K-factors, defined as the ratios of higher order over lower order cross sections, are shown in the lower panel of this figure. The NNLO corrections are small in the region p jet T > 10 GeV , changing the NLO result by no more than 10% over the studied p jet T region. The shift of the NLO cross section is slightly positive in the low transverse momentum region, and become less than unity at high-p jet T . Both the NNLO corrections and the scale dependence grow large at low-p jet T . The large scale dependence arises primarily from the partonic channels qq and gq. These channels are effectively treated at leading order in our calculation, since they first appear at O(α 2 α 2 s ), and they are evaluated at the low scale µ = p jet T /2 in our estimate of the theoretical uncertainty. It is therefore not surprising that their uncertainty dominates at low-p jet T These channels do not contribute at NLO, and consequently the NLO scale uncertainty is smaller. This is an example of the potential pitfalls in using the scale uncertainty as an estimate of the theoretical uncertainty. Only an explicit calculation can reveal qualitatively new effects that occur at higher orders in perturbation theory. We next show the pseudorapidity distribution in Fig. 8, with the restriction p jet T > 10 GeV. There are a few surprising aspects present in the NNLO corrections. First, the scale dependence at NNLO in the region η jet < 0 is larger than the corresponding NLO scale variation. Although the corrections are near unity over most of the studied pseudorapidity range, they become sizable near η jet ≈ 2, reducing the NLO rate by nearly 50%. To determine the origin of these effects we show in Fig. 9 the breakdown of the NNLO correction into its separate partonic channels. This reveals that the total NNLO correction comes from an intricate interplay between all contributing channels, with different ones dominating in different η jet regions. Only the gluon-lepton partonic process is negligible over all of phase space. For negative η jet , the dominant contribution is given by the quarkquark process. As discussed before, this appears first at O(α 2 α 2 s ). It is therefore effectively treated at leadingorder in our calculation, and consequently has a large scale dependence. We note that the quark-in-lepton distribution from Eq. (9) is larger at high-x than the corresponding photon-in-lepton one, leading to this channel being larger in the negative η jet region. At high η jet , the distribution receives sizable contributions from the gluonphoton process. No single partonic channel furnishes a good approximation to the shape of the full NNLO correction.  Figure 9: Breakdown of the NNLO correction to the ηjet distribution into its constituent partonic channels, as a ratio to the full NLO cross section in the bin under consideration. Also shown is the total result obtained by summing all channels. The bands indicate the scale variation.

CONCLUSIONS
We have presented in this paper the full calculation of the O(α 2 α 2 s ) perturbative corrections to jet production in electron-nucleus collisions. To obtain this result we have utilized the N -jettiness subtraction scheme introduced to allow NNLO calculations in hadronic collisions. We have described the application of this method to inclusive jet production in detail, and have shown that upon integration over the final-state hadronic phase that we reproduce the known NNLO result for the inclusive structure functions. Our results have been implemented in a numerical program DISTRESS that we plan to make publicly available for future phenomenological studies.
We have shown numerical results for jet production at a proposed future EIC. Several new partonic channels appear at the O(α 2 α 2 s ) level, which have an important effect on the kinematic distributions of the jet. No single partonic channel furnishes a good approximation to the full NNLO result. The magnitude of the corrections we find indicate that higher-order predictions will be an important part of achieving the precision understanding of proton structure desired at the EIC, and we expect that the methods described here will be an integral part of achieving this goal.