Inclusive-jet and Di-jet Production in Polarized Deep Inelastic Scattering

We present the calculation for single-inclusive jet production in (longitudinally) polarized deep-inelastic lepton-nucleon scattering at next-to-next-to leading order (NNLO) accuracy, based on the Projection-to-Born method. As a necessary ingredient to achieve the NNLO results, we also introduce the next-to-leading-order (NLO) calculation for the production of di-jets in polarized DIS. Our di-jet calculation is based on an extension of the dipole subtraction method to account for polarized initial-state partons. We analyze the phenomenological consequences of higher order QCD corrections for the Electron-Ion Collider kinematics.

cesses, improved the description of the gluon spin distribution, showing that its contribution to the proton spin is not negligible [3], although providing constraints only for a reduced range of proton momentum fractions. Furthermore, the amount of spin carried by the sea quarks is also still an open question [4,5]. In that sense, the future US-based Electron-Ion-Collider (EIC), allowing a much wider kinematical range, and reaching an unprecedented precision for polarized measurements [6], is expected to provide new insights on the spin decomposition of the proton in terms of its fundamental building blocks [7][8][9].
In addition to high-precision measurements for a wider range of momentum fractions, the improvement of our picture of the proton spin will require a consistent increase in the accuracy of the theoretical description of the observables to be measured. It is known that leading order (LO) perturbative calculations O(α 0 S ) in QCD only provide qualitative descriptions, since higher order corrections in the strong coupling constant are sizable. Although a remarkable effort to compute higher order corrections for unpolarized processes has taken place during the last 30 years, setting next-to-next-to-leading order (NNLO) as the standard for Large-Hadron-Collider (LHC) calculations and even reaching the following order for some processes, the picture for polarized calculations is not as developed. Polarized calculations in dimensional regularization necessarily involve dealing with extensions of the γ 5 matrix and Levi-Civita tensor to an arbitrary number of dimension, making the computation much more intricate than its unpolarized counterpart. Until recently, NNLO corrections for polarized processes were only obtained for completely inclusive Drell-Yan [10] and DIS [11], in addition to the helicity splitting functions [12][13][14]. More exclusive observables provide results that can be directly compared to experimental data, and could, in principle, be used to disentangle individual contributions associated to different partons. In particular, jet production in DIS is an extremely useful tool to probe the partonic densities, since it can give a stronger grip on the gluon distribution, while avoiding non-perturbative corrections associated to final-state hadronization. Developments in techniques for flavour and charge tagging in jet production could further improve the potential of jet measurements to disentangle individual flavour contributions in global analysis [15,16].
Higher order corrections are not only necessary to improve the accuracy of the theoretical description. It is also important to check the stability of the perturbative series, that is, how these corrections affect the resulting cross sections and spin asymmetries, since only processes perturbatively well behaved can be used as good probes for parton distributions, and be utilized for its extraction. Furthermore, for the specific case of jet production, it is only at higher orders in QCD that jet structure is fully developed, allowing to realistically match the theoretical description to the experimental data and the cuts imposed in the jet reconstruction.
In this paper we present the NLO calculation for di-jet production in polarized and unpolarized lepton-nucleon DIS, based on an extension of the Catani-Seymour dipole subtraction method [17] to account for polarized initial particles. We analyze the structure of higher order corrections in the Electron-Ion-Collider kinematics, its perturbative stability and phenomenological implications.
Through a detailed study of the polarized cross sections and asymmetries we also identify the most important partonic contributions for different kinematical regions. Additionally, we expand on our previous results [18] for single-exclusive jet production in DIS at NNLO, achieved by combining the aforementioned di-jet result with the inclusive polarized NNLO DIS structure functions [11] through the application of the Projection-to-Born (P2B) method [19]. We analyze the perturbative stability of the higher order corrections to the cross section and asymmetries, as well as the contributions from the different partons to the NNLO corrections. Both the NLO single-and di-jet, as well as the NNLO single-jet calculations are implemented in our code POLDIS [49].
The remaining of the paper is organized as follow: in section II we begin by defining the kinematics of both single-and di-jet production in DIS. In section III we detail both our extension of the dipole subtraction method for polarized QCD processes, and its use in the P2B method in order to achieve polarized jet production at NNLO. In section IV we present the phenomenological results for inclusive NLO di-jet production at the EIC in the Breit-frame, and in section V we do the same for inclusive NNLO single-jet production in the laboratory frame. Finally, in section VI we summarize our work and present our conclusions.

II. JET PRODUCTION KINEMATICS
We start considering the case of inclusive-single jet production in DIS. Specifically, we study the process e(k) + P (p) → e(k ) + jet(p T , η) + X, where k and p are the momenta of the incoming electron and proton, respectively, and k is the momentum of the outgoing electron. We work in the laboratory frame (L), where single-jet production receives non-vanishing contributions already at O(α 0 S ). We only consider, for the time being, neutral-current processes mediated by the exchange of a virtual photon, with its momentum q = k − k and virtuality Q 2 = −q 2 fully determined by the electron kinematics. The inelasticity y and Bjorken variable x are then defined as usual by In addition to the variables commonly used for fully inclusive DIS, more insight on the underlying partonic kinematics can be obtained through the study of the final-state jet, which can be characterized in terms of its transverse momentum p T with respect to the beam, and its pseudorapidity η.
At higher orders in α S , the production of multiple final-state jets becomes available. Di-jet production can be better studied in the Breit frame (B), where there is no contribution of O(α 0 S ) to the production of transverse jets. Formally, the Breit frame is defined as the one that satisfies 2xp + q = 0. Note that for the O(α 0 S ) process, this implies that the virtual photon and incoming quark collide head-on, completely reversing the momentum of the quark (hence the commonly used nickname brick-wall frame), as its represented schematically in Fig. 1. The first non-vanishing contribution is then obtained at O(α S ), with two final-state partons with opposed transverse momentum.
For di-jet production, we specify the process The availability of a second jet allows for a more in-depth study of the partonic kinematics. As in the H1 [20,21] and ZEUS [22] experiments, and in addition to the jets transverse momentum and pseudorapidities, the di-jet production cross section can be studied in term of the di-jet variables such as the invariant mass M 12 , the di-jet momentum fraction ξ 2 , as well as the average momentum p T 2 and pseudorapidity difference η * in the Breit-Frame, which are defined by It is worth noticing that, at the LO of di-jet production, ξ 2 is the momentum fraction carried by the incoming parton.

III. CALCULATION OF HIGHER ORDER CORRECTIONS
Calculations beyond the leading order in QCD necessarily involves cancellations between the individually divergent pieces coming from infrared real emission and virtual diagrams, in addition to the factorization contributions. In the dimensional regularization scheme the number of dimensions is set to D = 4−2 , and those divergences then appear as poles in . The cancellation between those poles can only be achieved after the integration of each of the divergent parts over its appropriate phase space, thus impeding a direct numerical calculation.
Several methods to numerically compute higher order corrections were developed over the last three decades. The two main approaches are based on either limiting the phase space integration (phase space slicing) in order to avoid the divergent regions, or generating appropriate counterterms (subtraction) to cancel the singularities in each of the pieces of the calculation. For the latter, the proposed counter-term should have the same divergent behaviour as the real and collinear parts, while being simple enough to be integrated analytically in order to cancel the poles coming from the virtual diagrams.
Many general methods for constructing NLO counterterms were proposed. Among them, the dipole subtraction method developed by Catani and Seymour, and based on the dipole factorization formula, allows to calculate any jet production cross section at NLO accuracy. The landscape for the following order is complicated due to the appearance of many more singular configurations, but several methods of varying generality are also available for the computation of NNLO calculations [19,[23][24][25][26][27][28][29][30]. In particular, for processes where the Born kinematic can be inferred from external non-QCD particles, the P2B method allows to obtain N k LO differential calculations for a jet observable O, given that the N k LO inclusive cross section and the differential N k−1 LO for O + jet are known. Consequently, given a NLO DIS di-jet calculation and the NNLO structure functions, we can then compute the NNLO exclusive single jet cross section.
This exclusive NNLO single-jet calculation is implemented in our code POLDIS for both polarized and unpolarized DIS. It allows to compute any infrared safe observable related to single-jet production at NNLO accuracy in the laboratory, as well as to single-and di-jet production in the Breit-frame with NLO precision. The code is partially based on DISENT, which implements the Catani-Seymour dipole subtraction method to obtain the NLO single-and di-jet cross sections in unpolarized DIS. Mayor modifications were made in order to include the polarized di-jet computation, using an extended version of the dipole subtraction to account for initial-state polarized particles, as well as the implementation of the P2B subtraction in order to obtain NNLO results.
Note that the previously reported bug in DISENT in the gluon channel [31][32][33][34] was fixed along with the modifications (see Appendix A).
Both the extension of the dipole subtraction as well as the P2B method will be discussed in more detail in sections III A and III B.
A. The dipole subtraction method for polarized processes For processes involving (polarized) unpolarized initial-state hadrons, QCD calculations necessarily involve convolutions between partonic cross sections and (helicity) parton distribution functions, (p)PDFs, codifying the (spin) momentum distribution of partons inside that hadron. In the case of DIS scattering, the (polarized) unpolarized hadronic cross sections can be written perturbatively as: where ... denotes higher order corrections. The helicity pPDF for a parton a carrying a fraction z of the proton's momentum p is defined as ∆f a (z, µ 2 ) denoting the density of partons of type a and momentum fraction z, with their helicities aligned (anti-aligned) with that of the proton. On the other hand, the polarized partonic cross section ∆σ ≡ 1 2 [σ ++ −σ +− ] is defined in terms of the difference between the cross sections with the incoming lepton and hadron polarized parallel and antiparallel. Up to NLO, the (polarized) unpolarized m-parton cross section is given by where d(∆)σ B is the (polarized) unpolarized partonic Born cross section, and d(∆)σ R and d(∆)σ V stand for the NLO partonic real-emission and virtual matrix elements, respectively. The last term in Eq. (5) is associated to the collinear factorization that must be introduced in the case of cross sections involving initial hadrons, to account for the divergences arising from initial-state radiation.
It is worth noticing that we are working in D = 4 − 2 dimensions, and that each of the integrals in Eq. (5) is separately divergent in the limit → 0. The calculation of polarized cross sections in dimensional regularization is more involved than its unpolarized counterpart, since the extension of the γ 5 matrix and the Levi-Civita tensor µνσρ appearing in the helicity projection operators in D dimensions is far from trivial. One way to consistently treat γ 5 and µνσρ is in the HVBM scheme [35,36], in which the D-dimensional space is separated in the standard four-dimensional subspace, and a (D − 4)-dimensional subspace. In this scheme, µνσρ is treated as a genuinely four-dimensional tensor, while γ 5 is such that {γ 5 , γ µ } = 0 for µ = 0, 1, 2, 3, and [γ 5 , γ µ ] = 0 otherwise.
An alternative to numerically compute the partonic cross section in Eq. (5) is the so-called dipole subtraction method, introduced by Catani and Seymour [17] as a general framework for the calculation of NLO jet cross sections. This is the method used to compute the NLO corrections of jet observables in both DISENT and POLDIS. As in other subtraction-based approaches, the idea behind the procedure is to cancel the infrared singularities that appear in the real, virtual and collinear-factorization pieces of the (polarized) unpolarized cross section, which are integrated in different phase spaces (m particles for the virtual diagrams and m + 1 for the real-emission diagrams), already at the integral level. That cancellation of divergences is achieved through the introduction of a counterterm d(∆)σ A that has the same infrared behaviour (in D dimensions) as d(∆)σ R . By adding and subtracting this term, the NLO calculation can be rewritten as In Eq. (6) the first integral can be numerically performed in four-dimensions since d(∆)σ A acts as a local counter-term of d(∆)σ R . In the second term the cancellation of poles requires the integrals to be performed analytically.
Clearly, the key of the subtraction method lies in the construction of d(∆)σ A , which in addition to reproduce the divergent behaviour of d(∆)σ R should be simple enough to be analytically integrated. In this case the term is constructed by the use of the dipole factorization formula in the collinear and soft limits, where ⊗ stands for the appropriate phase space convolution and sums over color and spin indices. The (∆)V dipole are the universal dipole factors that match the infrared singular behaviour of d(∆)σ R . Note that these terms need to be analytically integrable if D-dimensions over the single-parton spaces related to soft and collinear divergences in order to make use Eq. (6). The construction of these dipole factors for the unpolarized case was already outlined in detail in Catani and Seymour's paper. We now discuss the extension to the particular case of cross sections with one initial-state polarized parton, required for the calculation of the polarized DIS process.
Following the same notation introduced by Catani and Seymour, the complete polarized local counterterm d∆σ A a is: where the terms ∆D ij,k , ∆D a ij and ∆D ai k represent the dipole subtraction terms for final-state singularities with a final-state spectator, final-state singularities with an initial-state spectator, and initial-state singularities, respectively. The sum is performed over all the possible m + 1 finalstate partons configurations, with dφ m+1 denoting the corresponding phase space. Additionally, the 1/n c accounts for the average over the initial-state colors, Φ(p a ) is flux factor, and S {m} is the Bose symmetry factor for identical particles in the final-state. In N in the rest of the QCD-independent factors are included.
It is important to note that to create local counter-terms for the polarized DIS NLO cross section, only the polarization of the initial-state parton needs to be considered. In this case, instead of taking the average of its polarizations, the difference between them is used. Spin states of final-state parton are summed over and therefore they are treated as unpolarized. Thus, the dipole subtraction terms ∆D ij,k and ∆D a ij , associated to final-state singularities, are constructed as in ref [17] (using the corresponding polarized Born cross section). New expressions for the dipole formulas are therefore only needed in the case of initial-state singularities with one initial-state parton, represented by ∆D ai k . As in the case of the unpolarized cross sections, the terms ∆D ai k can be obtained from the dipole factorization formula. In the limit p a · p i → 0, where p a is the moment of the initial-state parton and p i a final-state one, the dipole factorization formula for the polarized (m + 1)-parton matrix element can be expressed as where |1, ..., m; ∆a m,a represents an m-particle state in the color and helicity space, with ∆a denoting that the difference between the incoming parton polarizations is considered. The k ∆D ai k stands for the sum of the polarized dipole contributions, in which the partons a and i act as a single initial-state parton ai, the 'emitter', and the final-state parton k acts as the 'spectator'k.
The ... stands for the other non-singular terms in the p a · p i → 0 limit. Each dipole contribution is given by The T are the color charge operators corresponding to each parton. The emitter and spectator momenta are given respectively by p µ ai = x ik,a p µ a andp µ The splitting functions ∆V ai k are the only new blocks needed for the extension of the dipole subtraction formalism to the polarized case. They are constructed so that they give the correct eikonal factors in the soft limits, and the correct D-dimensional polarized Altarelli-Parisi splitting functions ∆P ij in the corresponding collinear limits. Similarly to ∆P ij , ∆V ai k are matrices in the helicity space of the emitter parton ai, and their expression in given by: where Notice, however, that these expressions of the splitting functions ∆V ai k as matrices in the helicity states of ai are not really needed in the polarized case since the spin structure is trivial for both quarks and gluons. This is due to the fact that the spin correlation terms cancel out due to parity conservation in polarized processes (See Appendix B). Therefore, only the difference between the possible spin states of the emitter parton ai are required to perform the subtraction. In the case of a gluon emitter, this accounts for the contraction with the tensor i ρσµν p ρ a n σ /(2 p a · n), where n is any light-like vector that satisfies n · p a = 0, while for a quark emitter the tensor δ s s /2 is used.
The resulting kernels are: In order to integrate the dipole subtraction term, m+1 dσ A , the D-dimensional integrals of the ∆V ai k terms over the dipole phase space dp i (p k ; p a , x) are needed. The procedure to obtain them is the same one outlined by Catani and Seymour. The resulting expressions ∆V a,ai are: where x is the phase space convolution variable and the ∆P ab (x, 0) are the aforementioned polarized four-dimensional Altarelli-Parisi kernels. In the HVBM scheme they are given by [37]: A final remark must be made about the polarized factorization counterterms d∆σ C a in Eq. (5). These counterterms are explicitly written as: where the value of ∆K ab F.S determines the factorization scheme. We work in the conventional polarized MS factorization scheme in which one needs to compensate for the difference between the polarized and unpolarized quark splitting functions (∆P qq (x, ) and P qq (x, ), respectively) in D−dimensions. Since the difference between the two kernels is given by B. The Projection-to-Born method As it was mentioned, the P2B method allows to obtain the N k LO calculation for a differential observable, provided that its inclusive cross section at that order, as well as the differential cross section for the observable plus a jet are known at N k−1 LO. The idea behind the method is to cancel the most divergent parts by simply using the full matrix element at each phase space point as a counterterm, but binning it in an equivalent Born-projected kinematics of the leading order process (hence the name "Projection-to-Born"). That is, for each event with weight w, a counterterm with weight −w is generated, but with the measurement function evaluated in the kinematics of an equivalent leading order process. Note that this requires the Born kinematics to be fully determined by external non-QCD particles.
The differential cross section for an observable O at N k LO accuracy can be written as: where in dσ N k−1 LO O+jet infrared cancellation at the N k−1 LO level has already taken place (numerical implementations beyond leading-order thus require the use of an additional subtraction method).
It should be noted that as the final-state partons approach the most singular regions, they are So the mapping to the Born kinematics is given by using these parton momenta to evaluate the measurement function for the born-projected counterterms. Note that this mapping only works in jet production in the laboratory frame, since in the Breit-frame the first non-vanishing contributions starts at order O(α s ), with two final-state partons (and hence no mapping is possible in terms of P , x, and q).
In the particular case of single jet production in unpolarized (polarized) DIS at NNLO, the corresponding counterterms are generated from the double-real and one-loop real radiation matrix elements. The combination of those counterterms with the two-loop matrix elements is then equal to the unpolarized (polarized) DIS inclusive cross section at NNLO [11,38,39]. As mentioned, a numerical implementation of the calculation has yet to deal with the sub-leading divergences coming from the single-real radiation and one-loop diagrams contributing to the unpolarized (polarized) di-jet cross section at NLO. Those missing blocks can then be calculated with the implementation of the Catani-Seymour dipole formalism, whose extension to the polarized case was discussed in III A. We can then re-write Eq. (29) for the production of jets in unpolarized (polarized) DIS at NNLO in terms of the counterterms of Eq. (6) as: where we have used that the inclusive part can be expressed in terms of the P2B counterterms and the double-virtual matrix element for the observable O as: In addition, the complete expression for the counterterm d(∆)σ A O+jet is that given by Eq. (8).

IV. RESULTS OF POLARIZED NLO DI-JET PRODUCTION
The first step to reach NNLO accuracy for jet production in DIS lies in the calculation of the NLO di-jet cross section. Precisely, in this section we present our results for polarized inclusive di-jet production at NLO in the Breit frame (B). We consider the Electron-Ion Collider kinematics, with beam energies of E e = 18 GeV and E p = 275 GeV, and reconstruct the jets with the anti-k T algorithm and E-scheme recombination (R = 1). Furthermore, for di-jet production we fix the normalization and factorization scales central values as with the η cut imposed in the laboratory frame, while the lepton kinematics is restricted by 0.2 < y < 0.6, 25 GeV 2 < Q 2 < 2500 GeV 2 .
The parton distributions sets used were the NLOPDF4LHC15 [40] and DSSV [3,41] for the unpolarized and polarized case, respectively. momentum, the availability of a second jet allows to define more kinematical observables to analyze the underlying partonic kinematics in detail. In that sense, it is instructive to study the unpolarized cross section as a function of the usual di-jet kinematical observables p B T 2 , M 12 , η * and ξ 2 , defined in section II, as presented in Fig. 3.
As it was noted for the kinematics of HERA [42], higher order corrections are sizable for all the variables under consideration. The scale variations of the NLO calculation are as large as the LO ones, or even larger, in the lower bins of the M 12 , p B T 2 and ξ 2 distributions, as the infrared limit is approached. As mentioned, this behaviour is mainly due to the asymmetrical cuts in p T imposed to the two jets. In the Breit frame, LO kinematics implies that the two outgoing partons generating the jets have opposite transverse momentum, and therefore the region with p B T 2 < 5 GeV is not accessible at that order. A similar argument can be used to show that new regions of M 12 < 10 GeV and low ξ 2 become accessible only at NLO. This discrepancy in the available phase space at different orders is known to cause instabilities in the perturbative expansion [43]. Actually, for that forbidden phase space region the calculation is effectively a LO one. Note, however, that the use of symmetric cuts in p T leads to even worse perturbative problems, due to the enhancement of large logarithmic contributions related to the back-to-back configuration that can completely spoil In Fig. 4 we show the same distributions of Fig. 3 but for the polarized cross section. Compared to the unpolarized case, for low M 12 , p B T 2 , η * and ξ 2 it can be seen that while the NLO corrections follow the same pattern, they are generally milder, with lower K-factors. There is also a difference in the behaviour of the second order corrections for higher values of M 12 , η * and ξ 2 , resulting in stronger suppressions than the ones observed in the unpolarized case. The ξ 2 distribution is particularly shifted towards higher momentum fractions. The same considerations regarding theoretical uncertainties apply to the polarized case, leading to the strong NLO scale-dependence.
The somewhat big NLO corrections, and the differences between the unpolarized and polarized cases, can be better understood by analysing the previous distribution at different values of Q 2 .
As an example, in Figs. 5 and 6 we present the unpolarized and polarized double-differential distribution, i.e, in bins of Q 2 and log 10 (ξ 2 ), respectively. Regarding the unpolarized distributions of Fig. 5 it can be noted that, as expected, lower Q 2 values are correlated to smaller momentum fractions, from which the cross section receives its most important contributions. Di-jet production measurements at the EIC are therefore expected to explore the mid-x region, 10 −2 < x < 10 −1 .
The NLO cross sections for the high Q 2 bins are in good agreement with the LO calculations and show small scale dependence, indicating good convergence of the perturbative series. In addition to the complementary constraints on the quarks polarized and unpolarized distribution functions, restrictions coming from this region on the gluon helicity distribution, which is mainly probed down to x ∼ 5 × 10 −2 by RHIC data, will be specially important. On the other hand, in Fig. 5 it can be seen that both the K-factors and theoretical uncertainties increase as lower Q 2 values are considered. This is consistent with the aforementioned population of the new phase space region at low ξ 2 becoming available at NLO.
Compared to the unpolarized case, the polarized distributions of Fig. 6 present two striking features: they decrease at lower Q 2 , and they also display significant differences in shape between LO and NLO results in that region. Both features can be explained by the analysis of the contributions from the quark and gluon channels to the polarized cross section. In Fig. 7 we present, precisely, the di-jet double-differential polarized distribution as a function of Q 2 and log 10 (ξ 2 ), distinguishing the contributions initiated by the quark and gluon channels. In this case, the lower insets in the plot depict the ratio between the gluon-and quark-initiated differential cross sections.
The peculiar behaviour of the polarized cross sections as a function of Q 2 can be traced back to the negative sign of the gluon contribution below Q 2 = 600 GeV, which becomes more significant for lower values of Q 2 , as shown in the ratio between the gluon and quark contributions. The cross section leads to an enhancement in the asymmetry at lower momentum fractions, albeit the very small values of the asymmetry in that region.
Once again, the behaviour of the asymmetries can be better understood by studying the doubledifferential Q 2 dependence of the distributions. Fig. 9 depicts the double spin asymmetry as a function of both Q 2 and ξ 2 . The reduction of the polarized cross sections for low values of Q 2 due to the negative gluonic contribution leads to a sizable suppression of the asymmetry in those bins for ξ 2 10 −1 . It is worth mentioning that, for the first to bins of Q 2 , the significant shift in the NLO quark contribution towards lower momentum fractions shown in Fig. 7

V. POLARIZED NNLO INCLUSIVE-JET PRODUCTION
Having discussed our NLO di-jet production calculation, we can now turn to the NNLO corrections for single jet production, obtained through the application of the P2B method. In this section, we present our results for polarized single-inclusive jet production at NNLO in the laboratory frame (L), for the Electron-Ion-Collider kinematics. Similarly to [18], the default distributions are obtained reconstructing the jets with the anti-k T algorithm and E T -scheme recombination, using a jet radius R = 0.8, and fixing the normalization and factorization scales central values as µ 2 F = µ 2 R = Q 2 ≡ µ 0 . As in the previous section, α s is evaluated at NLO accuracy with α s (M z ) = 0.118. The reconstructed jet in the laboratory frame is then required to satisfy: while on the leptonic side we impose the additional cuts: 0.04 < y < 0.95, The lower cut in Q 2 was chosen to avoid differences in the phase space available at different orders. Note that at LO the transverse momentum of the jet in the laboratory frame is given by , and thus the region Q 2 25 GeV 2 is kinematically forbidden for the specified cuts in p L T . Since there is no NNLO global fit of polarized PDFs available, the parton distributions sets used were, once again, the NLO extractions NLOPDF4LHC15 [40] and DSSV [3,41] for the unpolarized and polarized case, respectively.
In Fig. 10 we present the cross section for single-inclusive jet production in polarized DIS, as a function of the jet transverse momentum p L T , its pseudorapidity η L , and in terms of Q 2 and x, calculated at LO, NLO and NNLO accuracy. The lower insets in Fig. 10 show the K-factors, defined as the ratios to the previous order, that is, K NNLO = σ NNLO /σ NLO and K NLO = σ NLO /σ LO . As in the case of di-jet production, the theoretical uncertainty bands were obtained performing a sevenpoint independent variation of the renormalization and factorization scales as µ R , µ F = [ 1 2 , 2]µ 0 . The uncertainty associated to the polarized parton distributions was estimated using the DSSV set of PDFs replicas from [41]. Note that due to the unavailability of proper polarized NNLO PDF, this bands should be taken only as a first attempt to quantify the non-perturbative errors in the NNLO cross section. The same NLO PDFs were used at all orders so as to quantify only the variations arising from the perturbative calculation.
As it can be seen in Fig. 10, the main effect of higher order corrections is to shift these distributions toward higher values of pseudorapidity and lower values of transverse momentum, since more jets originating from the emission of additional partons become available in those regions. In the case of the pseudorapidity distribution, this is translated into high values of the NLO K-factor in the forward region (η L > 1), while a strong suppression in the backward region (η L < −1) is observed. NNLO corrections have the same behaviour, albeit with lower values of K-factor. Similar comments can be made regarding the transverse momentum distribution, which is enhanced for lower values of p L T . For the p T distribution, the NNLO corrections are typically of order 10%, while for the η distributions they are of order 5%. It should be noted that while there is good agreement between the NLO and NNLO calculations, with overlapping bands throughout the kinematical range, anticipating convergence of the perturbative series, the scale bands for the NNLO distributions are still somewhat large in certain bins compared with those of the NLO. This effect is associated with the Even though the growth of the uncertainty bands at NNLO in the p L T and η distributions originates from the difference in the available phase space at each order, the sizes of the bands in this region are further enhanced in the polarized case compared to the unpolarized one. This results in bigger NNLO bands in p L T and η L distributions, as observed in [18]. This enhancement is related to the fact that in the polarized case there are cancellations between processes initiated by different partons. To highlight this point, in Fig. 11 we present the contributions of the most  relevant parton channels to the polarized cross section. As in the unpolarized case, for most of the explored Q 2 and x values, the cross section is dominated by initial u quark contributions. However, as lower values of both Q 2 and x are reached, there are significant cancellations between the u quark channel and the negative contribution of the d quark and gluon channels, which accounts for higher relative uncertainties once the sum over of all the initial parton contributions is taken (the s quark also has a negative contribution, but it is negligible). Since low Q 2 and x correlate with low p L T and η L 0, those same cancellations are translated into the sizable NNLO scale bands in Fig. 10 in those ranges.
It is worth noticing that, even though it is expected to have a greater gluon contribution at low p L T , since that region correlates with low Q 2 , the first bin of the p L T distribution is very small and slightly positive (as opposed to the u and d quarks contributions). This is related to the fact that the gluon contribution to the structure function is positive below x ∼ 2 × 10 −2 . Since the structure function is obtained by the integration over all the p L T range, as lower values of p L T are reached the p L T distribution must become positive at some point. Regarding the uncertainty associated to the PDFs, it is typically of order 5%−10% for the region of {p L T , η L } studied. Though this uncertainty is comparable to the NNLO corrections for most of the kinematical range, it should be noted that for the low p L T region, it becomes smaller than the NNLO corrections, highlighting the relevance that NNLO extractions will have in order to match the accuracy of the perturbative side. As in the case of the scale-variations bands, the PDF uncertainty becomes larger as lower values of x and Q 2 are approached, since the cancellation between the different partonic channels for those bins is sensitive to changes in the partonic distributions.
Another feature associated to cancellation between partonic channels in the polarized cross section is the reduced dependence on the parameters of the jet-reconstruction algorithm, compared to the unpolarized case. To emphasize this point, in Fig. 12, we present the NNLO cross sections as a distribution of both p L T and η L , for different values of the jet radius R = 0.5, 0.8, 1 used in the anti-k T algorithm. In both cases, higher values of jet radius correspond to larger cross sections in the whole kinematical range due to the inclusion of more jets that satisfy the imposed cuts.
However, the polarized case shows a reduced dependence in R at low p L T and the intermediate  η L values, precisely where the stronger cancellations between channels take place. This results in an overall reduction of the dependence of the polarized cross section on the jet parameter. It is worth noticing that while the total cross section is affected by these strong cancellations between channels, with the use of jet tagging techniques [15,16] it could be possible to noticeably modify the shape of the distributions, enhancing the contributions from different partons.
The difference of sensitivity to changes in the jet radius will in turn modify the behaviour of the double spin asymmetries. In Fig. 13 we present the NNLO double spin asymmetries in the p L T and η L distributions for the R values used before. As expected, a larger dependence on R is obtained in those regions where the cancellation between channels for the polarized cross sections are more important. For those regions, the increase in R leads to a relative increase of the unpolarized cross section, and consequently to a reduction in the spin asymmetry. Conversely, lower values of R produce an increment of the asymmetry in the same regions. Fig. 13 also shows the LO and NLO asymmetries for R = 0.8. Regarding higher order corrections, it is worth mentioning that the relative higher NNLO contributions to unpolarized cross section lead to an important suppression of the asymmetry in the high pseudorapidity region, with milder corrections for intermediate η L .
However, note that for η L 1 and p T 10, the variations with the jet radius are greater than those coming from the perturbative series. The jet parameters are therefore expected to have sizable impact in the double spin asymmetries in regions where cancellation between partonic contributions take place in the polarized cross section.

VI. CONCLUSIONS
In this paper we have presented the NLO calculation for the production of di-jets in polarized and unpolarized lepton-nucleon DIS in the Breit frame, for the EIC kinematics. Our calculation is based in a generalization of the dipole subtraction method to handle the polarization of initial-state particles, which is discussed in detail. The cross sections were studied as functions of the leading jets transverse momenta p B T,1 and p B T,2 , invariant mass of the jets M 12 , the mean transverse momentum p B T 2 , the difference in pseudorapidities η * and the di-jet momentum fraction ξ 2 . Additionally, the double-differential distributions in Q 2 and ξ 2 were analyzed. Both for the polarized and unpolarized cross sections, the differential distributions show important NLO corrections, particularly for low values of M 12 and ξ 2 , and higher values of η * , associated to differences in the phase space available at each order. While the NLO corrections obtained show good agreement with the LO calculations and reduced dependence on the choice for the factorization and renormalization scales, for values of Q 2 above 250 GeV, anticipating convergence of the perturbative expansion, the distributions for lower values of Q 2 present sizable corrections as well as a strong dependence on the scale choice. We noted that this effect is further enhanced in the polarized cross sections, due to the non-negligible negative contribution of the gluon-initiated channel, producing noticeable differences between the polarized results and their unpolarized counterparts. This difference in behaviour is translated to the double spin asymmetries, with significant suppression in M 12 , η * and p B T 2 . Once again, the corrections are more significant as lower values values of Q 2 are approached.
The di-jet calculation was in turn used to obtain the polarized NNLO single-inclusive jet production cross-section in the laboratory frame via the P2B method [18], which combines the exclusive NLO di-jet cross section along with the inclusive NNLO polarized structure function. We expanded on our previous results to include a better estimate of the theoretical uncertainty, as well as the dependence on the jet radius. Good agreement was found between the NLO and NNLO results for the range studied in p L T and η L . The somewhat large size of some of the NNLO uncertainty bands was linked to a combination of the effects due to the difference in phase space available at LO at low Q 2 and x, also present in the unpolarized case, as well as the cancellation between partonic channels in the polarized cross section. This channel cancellation also leads to a reduced dependence of the polarized cross section in the jet radius R, which in turn produces a more noticeable dependence of the double spin asymmetries in R in the regions of low p L T and intermediate values of η L . This hints towards a sizable dependence of the polarized cross section and asymmetries with the jet parameters in those regions, as well as important sensibility to the recently proposed jet-tagging techniques.
The results presented on this paper highlight the relevance that higher order QCD corrections will have in the precise description of the jet observables to be obtained in the future EIC, as well as the potential of those measurements to further improve our understanding of the spin structure of the proton and, particularly, in the precise extraction of polarized parton distributions.

Appendix A: Dipole bug in DISENT
The presence of a bug in the gluon channel in DISENT was reported long ago in [31][32][33][34], particularly while studying the event shape distributions in DIS. After a careful analysis, along with an extensive comparison with DISASTER [45] (a code which showed good agreement with resummed event shape calculations), and also by writing independent codes, we found that the Born matrix element used in one of the dipole subtraction terms in the gluon channel had the momentum of two final-state partons interchanged, leading to the reported discrepancies. Due to the nature of the bug, it turns out to produce noticeable differences only in certain extreme regions of the phase space, and remains within the typical statistical uncertainties of the calculations in many others. We have checked that the fixed counterterm actually corrects the reported disagreement between DISENT and DISASTER in the event shapes, as well as the differences between DISENT and the analytical calculation for logarithmically enhanced terms. As an example, we present in Fig. 14 the difference between the O(α 2 S ) coefficient for the fixed-order Monte Carlo calculation and the expansion of the resummed calculation for the gluonic contribution to τ z E , using DISASTER, the v0.1 version of DISENT and its fixed version. The event shape presented was calculated (at x bj = 0.0039, Q 2 = 7.5 GeV and y = 0.001) with the programs Dispatch and DISresum, written by Salam et al [31,32,46]. Similar results are obtained in the case of τ z Q . We also found agreement between DISASTER and the modified version of DISENT for the quark channel, and for other event shapes. It should be noted that, even after fixing all the external particles helicities, the factorization at the amplitude level involves the summation over the helicity states λ e of the intermediate parton.
The case in which p e is a quark is trivial, since helicity conservation at the vertex implies that one of the terms in the sum over λ e is zero. The case with an intermediate gluon is, however, more involved. The exact factorization is lost at the squared-amplitude level, through the appearance of interference terms between the different helicities in the propagator. Those interference terms give rise to, precisely, the spin correlations noted in the dipole factorization formula. The exact form of the correlation terms can be easily obtained by squaring Eq. (B1): where |N n (λ 1 , λ 2 , λ a , {λ X f })| 2 = g 2 S C |S λ 1 λ 2 + 1g (z)| 2 |M n−1 (+, λ a , {λ X f })| 2 and the interference term is given by In Eqs. (B3) and (B4) we introduced the short-hand notation for the color factor with 1/N c 1 denoting the average over the initial parton colors. For the relevant cases, and using the normalization from [48], C can take the values 2C A and C F , for an initial gluon and quark, respectively. Notice that the interference term depends on the initial parton helicity λ 1 only through the spin-dependent kernels S λ 1 λ 2 λe 1g . In the calculation of the unpolarized (polarized) cross section, we can then write: where the helicity factor (λ a ) should only be considered in the polarized case. Using Eq. (B2), the unpolarized (polarized) cross section can in turn be expressed as (z)| 2 |M n−1 (+, λ a , {λ X f })| 2 where we have used that the polarized Altarelli-Parisi kernels for z < 1, ∆P < 1j (z), can be obtained from the helicity-dependent kernels as