Approximate NNLO QCD corrections to semi-inclusive DIS

We determine approximate next-to-next-to-leading order (NNLO) corrections to unpolarized and polarized semi-inclusive DIS. They are derived using the threshold resummation formalism, which we fully develop to next-to-next-to-leading logarithmic (NNLL) accuracy, including the two-loop hard factor. The approximate NNLO terms are obtained by expansion of the resummed expression. They include all terms in Mellin space that are logarithmically enhanced at threshold, or that are constant. In terms of the customary SIDIS variables $x$ and $z$ they include all double distributions (that is,"plus"distributions and $\delta$-functions) in the partonic variables. We also investigate corrections that are suppressed at threshold and we determine the dominant terms among these. Our numerical estimates suggest much significance of the approximate NNLO terms, along with a reduction in scale dependence.


Introduction
Data taken in the semi-inclusive deep-inelastic scattering (SIDIS) process p → hX offer powerful insights into QCD and hadronic structure. Among their main uses are extractions of fragmentation functions [1][2][3][4][5], (polarized) parton distributions [6,7], or even combinations thereof [8,9]. Today, modern "global" analyses of parton distributions are customarily carried out at nextto-next-to-leading order (NNLO) accuracy of QCD perturbation theory. Although SIDIS might in principle offer important complementary information on, for example, the flavor structure of the sea quarks, the analyses usually do not include information from SIDIS. One reason for this is the fact that the NNLO partonic hard-scattering functions for SIDIS are not yet available (a few first steps toward their calculation have been taken in [10][11][12]), so that computations of the SIDIS cross section are currently restricted to next-to-leading order (NLO).
The Electron Ion Collider (EIC) is now firmly on its path toward construction [13]. The past few years have seen tremendous progress on the development of the theoretical framework for describing reactions relevant at the EIC. Further improvements will likely occur in the near term. Ideally, by the time the EIC will turn on, it would be hoped that the precision of theoretical calculations should be on par with what has by now been achieved for the LHC, with NNLO corrections available for many observables, and extractions of parton distributions and fragmentation functions routinely at NNLO using numerically efficient tools. As part of this, it is expected that also full calculations of the NNLO corrections to SIDIS will become available at some point. Until this is the case, it is useful to provide accurate approximations of the NNLO corrections for SIDIS. This is the main goal of this paper. The results we obtain may be used to carry out analyses of parton distributions and/or fragmentation functions using SIDIS data at (approximate) NNLO already now.
The strategy we will follow to derive approximate NNLO corrections to SIDIS is to use QCD threshold resummation. The partonic SIDIS process is characterized by two "scaling" variables, x = −q 2 /2p · q ≡ Q 2 /2p · q andẑ = p · p c /p · q, with q, p, p c the momenta of the virtual photon, the incoming parton, and the fragmenting parton, respectively. Whenx andẑ get close to 1, the partonic hard-scattering functions develop large double-logarithmic terms. These logarithms arise since largex,ẑ corresponds to scattering near a phase space boundary, where real-gluon emission is suppressed. At the kth order of perturbation theory, the SIDIS quark hard-scattering function contains terms of the form with m + n ≤ 2k − 2. Here the subscript "+" indicates the usual distribution. Threshold resummation addresses these large logarithmic terms to all orders in the strong coupling. The resummation for the case of SIDIS was discussed in Refs. [14][15][16][17] to next-to-leading logarithm (NLL), which amounts to the cases m = 2k − 1, 2k − 2, 2k − 3 and m + n = 2k − 2, 2k − 3 above, respectively. To NLO, this reproduces all double distributions, but only the three leading towers of logarithms at NNLO and beyond.
In the present paper we take a significant step further and extend the work of [15,16] to nextto-next-to-leading logarithm (NNLL). The close correspondence of SIDIS with the Drell-Yan cross section is particularly useful in this context [17], and so is the close correspondence between the totally inclusive Drell-Yan cross section and the cross section differential in rapidity [15,18,19]. We use results available in the literature [20][21][22][23] to determine the two-loop hard virtual contribution to the resummed expression for SIDIS. The NNLL results may then be expanded to fixed order, NNLO. The main important new result of our paper is that we derive all double distributions inx andẑ in the NNLO SIDIS quark coefficient function. We further improve our results by deriving the dominant part of NNLO contributions that are suppressed near threshold. These terms are of the form ln m (1 −x) ln n (1 −ẑ), with m + n = 3. We also show that the NNLO contributions near threshold are the same for the spin-averaged and spin-dependent cases. This is indeed expected for the terms with + -distributions, since these terms are associated with emission of soft gluons which does not care about spin, but it extends even to the threshold-suppressed contributions that we derive.
Our results are readily suited for phenomenology for the SIDIS cross section and spin asymmetry at (nearly) NNLO. At the very least, they provide important benchmarks for future full calculations of SIDIS at NNLO.
Our paper is organized as follows: Section 2 sets the stage by addressing the perturbative SIDIS cross section. In Sec. 3 we determine the threshold resummation for SIDIS to NNLL. Special emphasis is put on the derivation of the hard factor at two loops. Section 4 addresses the dominant threshold-suppressed contributions. Having determined all ingredients, we finally present the NNLO expansions in Sec. 5. Section 6 rounds off the paper by presenting some basic phenomenological results at approximate NNLO.

Perturbative SIDIS cross section
We consider the semi-inclusive deep-inelastic scattering (SIDIS) process (k) p(P ) → (k ) h(P h ) X with the momentum transfer q = k − k . It is described by the variables We have Q 2 = xys, with √ s the center-of-mass (c.m.) energy for the incoming electron and proton. We may write the spin-averaged SIDIS cross section as (see, for example [15]) where α is the fine structure constant and F h T ≡ 2F h 1 and F h L ≡ F h L /x are the transverse and longitudinal structure functions. For collisions of longitudinally polarized leptons and protons we obtain the helicity-dependent cross section as 1 2 with G h 1 = 2g h 1 in the more conventional notation of Ref. [24]. The subscripts in the first expression denote the helicities of the incoming lepton and proton.
Using factorization, the unpolarized structure functions may be written as In the following, it is convenient to take Mellin moments of the SIDIS cross section, for which the convolutions in Eqs. (4),(5) turn into ordinary products. We definẽ and in the same way for G h 1 . One readily finds from (4) wheref , We observe that the Mellin moments of the structure functions are obtained from the moments of the parton distribution functions and fragmentation functions, and the double-Mellin moments of the partonic hard-scattering functions. For the spin-dependent case we have in the same waỹ with the corresponding moments ∆f (N, µ F ) and ∆ω f f (N, M, . . .) of the polarized parton distributions and hard-scattering functions, respectively.
For the perturbative expansions given in Eqs. (6), (7), we have at lowest order according to (8) The corresponding moments of the next-to-leading order (NLO) terms ω (1) f f may be found in Refs. [15,30]. In the following, we address higher-order corrections to the hard-scattering functions that arise at large values ofx andẑ or, equivalently, at large N and M . As has been discussed in [15,16] (and as is familiar from numerous other situations in perturbative calculations of cross sections), in the "threshold limit" of large N and M the perturbative QCD corrections forω T f f and ∆ω f f develop large double-logarithmic corrections in ln(N ) and ln(M ). These corrections exponentiate and may thus be controlled to all orders in the strong coupling, amounting to a resummation of the logarithmic corrections. The exponentiated result may be used to obtain approximate fixed-order corrections to the SIDIS cross sections.
To achieve the resummation of the threshold logarithms for SIDIS, we will use the methods developed in Refs. [15,17,31,32]. Technically, the resummation for SIDIS with its two Mellin variables N and M bears much resemblance with that for the Drell-Yan or Higgs cross sections at measured rapidity, which are also described by two separate moments [14,18,19,33]. This is in contrast to observables characterized by a single moment variable N , such as the totally inclusive Drell-Yan cross section. However, as was shown in Ref. [15] to NLL, there is a simple correspondence between the threshold-resummed expressions for the case with two Mellin moments, and those with only a single moment. To state this correspondence, let us consider the resummed qq hard-scattering function for the Drell-Yan process as an example. For the totally inclusive cross section, we denote the function byω DY,incl √ z e ±y , respectively, with y the lepton pair's rapidity. Near threshold one then has This correspondence is a consequence of kinematics in the exponentiation of eikonal diagrams as discussed in Refs. [15,17,34]. It applies to all color-singlet processes and may therefore also be exploited for the SIDIS process § . As a result, we may obtain resummed expressions for SIDIS by considering those for the inclusive Drell-Yan process, "rescaling" N to √ N M appropriately, and "crossing" from timelike (Drell-Yan) kinematics to spacelike (SIDIS) kinematics. This is the strategy we will pursue in this paper.
There are various (of course, equivalent) ways of writing the all-order expression for the resummed inclusive Drell-Yan hard-scattering function near threshold. Here we will follow the approaches developed in Refs. [31,32]. We have, in the MS scheme, § As discussed in Ref. [17], one may actually define a simplified variant of SIDIS that is characterized by only a single Mellin variable, conjugate to τ SIDIS = xz, with x, z defined in Eq. (1).
with the Euler constant γ E . In Eq. (15) each of the functions H DY qq , C qq , A q , D q is a perturbative series in the strong coupling with expansion coefficients that are collected in Appendix A to the order required for resummation at next-to-next-to-leading logarithmic (NNLL) accuracy. The factor ∆ q in the first line contains all soft-gluon radiation near threshold (both collinear and wideangle), while the coefficient H DY qq collects hard virtual corrections to the underlying lowest-order (LO) process (here, qq → γ * ), which are independent of the moment variable. In the second line we have followed Refs. [31,32] to split up the soft-gluon factor ∆ q into the term C qq that is again independent of N , and an exponential that contains all N -dependence. The latter is in fact entirely a function of ln(N ) and contains no further N -independent terms.
The NNLL resummation formula for the SIDIS transverse structure function may now be written as follows: (17) As anticipated, we have "rescaled" N to √ N M in the moment-dependent part of the expression. The function ∆ q is otherwise identical to that for the Drell-Yan case in Eq. (15), including the function C qq . The hard coefficient H SIDIS qq is, however, different from H DY qq , owing to the different kinematics of the two processes. It will be derived in the next subsection. Inserting ∆ q from (15) into Eq. (17) we obtaiñ where (see (16))M = M e γ E . We note that the same resummation formula applies to the spindependent case: ∆ω res qq =ω T,res qq .

The hard factor H SIDIS qq
As already mentioned, the factor H SIDIS qq is derived from the finite part of the virtual corrections to the LO process, which for SIDIS is qγ * → q. Since we want to derive the resummed formula to NNLL (and ultimately the near-threshold NNLO corrections to SIDIS), we need H SIDIS qq to two loops. The relevant two-loop virtual corrections are known in terms of the "quark form factor" computed to two and even three loops in Refs. [22,23,35]. In case of the space-like kinematics (q 2 < 0) relevant for SIDIS the renormalized spacelike quark form factor is given to two loops in dimensional regularization with d = 4 − 2 space-time dimensions as [22,23] where with N f the number of flavors and C F = 4/3, C A = 3. In these expressions we have kept terms of order and 2 in the one-loop result since these turn out to make finite contributions in the end.
As shown in Refs. [20,21], the hard coefficient may be extracted from the form factor in the following way. Applied to the case of SIDIS we have from [21] whereĨ q is an operator that removes the poles of the form factor and makes the necessary soft and collinear adjustments needed to extract the hard coefficient. It is given in [21] in terms of a convenient all-order form: with functions R q and Φ q that each are perturbative series. The phase Φ q does not contribute in our case since we take the absolute square in Eq. (22). The function R q effects the cancelation of infrared divergences from the quark form factor. It can be expressed in terms of a soft and a collinear part: where for NNLL accuracy with The coefficient b 0 can be found in Appendix A. Inserting all terms into Eq. (22) and expanding in α s , all poles in powers of 1/ cancel, and we find for an arbitrary renormalization scale µ R , but for µ F = Q: with The factorization scale dependence of H SIDIS qq is trivially determined by the DGLAP evolution kernels of the parton distributions and fragmentation functions and will be addressed later.
With all ingredients to NNLL resummation at hand we are now also in the position to expand the hard-scattering function in (18) to NNLO (that is, O(α 2 s )) accuracy. This expansion will be carried out in Sec. 5. Before turning to it, we will discuss another class of corrections near threshold that are suppressed with respect to the terms addressed by resummation, but that can be significant as well in phenomenological studies.

Subleading contributions near threshold
All contributions contained in Eq. (18) are leading near threshold in the sense that they carry powers of ln(N ) or ln(M ), never accompanied by any suppression by 1/N or 1/M . Such terms are therefore often referred to as leading-power (LP) contributions. For the NNLL resummed cross section the LP terms contain the five "towers" α n s L m , with m ∈ {2n, . . . , 2n − 4}, where L m can be any product of (in total) m logarithms in N or M . The LP terms correspond to distributions ("+"-distributions and δ-functions) inx,ẑ space. In the full cross section there are, of course, also terms that are suppressed near threshold. The most important among these are terms still containing logarithms, but suppressed by a single power in 1/N or 1/M . Such terms are known as next-to-leading power (NLP) corrections. Their structure is α n s L m /N or α n s L m /M , with m ∈ {2n − 1, . . . , 2n − 3}, corresponding to terms of the form α n s m inx,ẑ space, where m is a product of ln(1 −x) and ln(1 −ẑ) with total power m.
The role of NLP terms in color-singlet hard-scattering cross sections has been addressed early on in Refs. [31,[36][37][38][39]. In recent years, the understanding of such corrections has further advanced, and numerous studies have been carried out  that address the NLP contributions from various angles, such as corrections to the eikonal approximation, resummations of NLP terms to leading logarithm and beyond, and generalized factorization theorems at NLP. For especially simple processes such as the fully inclusive Drell-Yan process, the results of these studies are quite mature. For processes described by two scaling variables (or, two Mellin moments), as relevant for SIDIS, comparably fewer studies are available [58,63]. In the present study we will derive the dominant NLP contributions at NNLL which, as described above, are of the form α n s L 2n−1 /N or α n s L 2n−1 /M . In terms of the NNLO expansion, these are the terms As discussed in [31,36,38], these dominant NLP terms may be incorporated to all orders via a particular treatment of the evolution of the parton distributions and fragmentation functions ¶ . To this end, we consider a specific SIDIS quark channel in the spin-averaged case and include the parton distribution and fragmentation function. From Eqs. (10),(18) the corresponding resummed contribution to the transverse SIDIS structure function in moment space may be written as where the function P q,δ corresponds to the coefficient of δ(1 − x) in the quark DGLAP splitting function and is also given in Appendix A.
We make the following observations concerning Eq. (29). We obviously have simply added and subtracted the terms involving P q,δ in the exponent, so that they cancel. However, each of the individual terms serves a separate purpose. The P q,δ term in the second line, when combined with H SIDIS qq α s (µ R ), µ R /Q, µ F /Q , removes the factorization scale dependence of the SIDIS hard function, so that we end up with H SIDIS qq α s (µ R ), µ R /Q, 1 , precisely as given in Eq. (28). Thanks to factorization, this must hold true to all orders of perturbation theory. The other two P q,δ terms in Eq. (29) combine with the terms A q lnN or A q lnM to reproduce the quark-to-quark splitting function in the large-N or large-M limit, at leading power. As a result, the last two exponential factors simply represent the DGLAP evolutions of the quark parton distribution function and the fragmentation function, respectively, from scale µ F to scale Q/ √NM . At leading power, this evolution is entirely diagonal, and evolution of parton distributions (spacelike) and of fragmentation functions (timelike) is identical. We can therefore carry out this evolution and write Eq. (29) as Again, this is correct to all orders. The trick now to obtain the dominant NLP corrections is to evolve the parton distributions and fragmentation functions from scale Q/ √NM back to scale µ F , but now using the DGLAP evolution including NLP corrections [31,36,38]. The latter are readily obtained from the 1/N or 1/M terms in the spacelike or timelike splitting functions, respectively. As it turns out, for the dominant NLP terms, only the 1/N (or 1/M ) terms in the leading-order splitting kernels need to be taken into account. The related terms in the higherorder splitting functions lead to contributions that have fewer logarithms. Let us for the moment continue to consider only diagonal evolution, corresponding to the SIDIS quark channel. We write the standard LO quark-to-quark splitting function at large values of the moment variable as The term proportional to Q (1) q is the NLP correction. At this order, the spacelike and timelike quark-to-quark splitting functions are identical so that also their NLP corrections are the same. The relation betweenq N, Q/ √NM andq(N, µ F ) including the dominant NLP correction is now given bỹ and in the same way for the quark fragmentation functions. As we discussed, only the LO term with Q (1) q is relevant for the dominant NLP corrections. The terms with A q and P q,δ remain, of course, the full all-order functions, needed to second order (NLO) for our purpose of obtaining NNLL/NNLO accuracy. We see from Eq. (32) that the dominant NLP corrections in the quark channel are obtained by multiplying the full resummed expression in Eq. (29) by the two factors exp − corresponding to the NLP terms related to diagonal evolution of the parton distribution and the fragmentation function.
As is well known, once the NLP terms are included, the evolution of parton distributions and fragmentation functions also involves quark-gluon mixing and hence is no longer diagonal, taking instead a matrix form. Transitions among quarks of different flavor turn out to be suppressed as 1/N 2 or higher, at least through NLO in the evolution kernels which is all we need here. Including the dominant NLP corrections, the full evolution equations for the parton distributions may be cast into the form to all orders, where P N s (α s ) denotes the NLO matrix of spacelike splitting functions in moment space, which may be found in [64]. A corresponding equation holds for the fragmentation functions, with however the timelike splitting functions P M t [64].
It is interesting to explore the implications of the singlet mixing and to see what NLP effects it generates beyond the quark-to-quark channel. We will do this as part of the NNLO expansion to be discussed in the next section. For this expansion we do not need to fully solve the evolution equation (although this could be done using the techniques of Ref. [27]). Instead, it suffices to just solve the equation to second order in the strong coupling, which may be achieved by iterating the kernel: and similarly for the fragmentation functions. This expression may then straightforwardly be expanded further in α s (µ R ). If we keep just the diagonal (quark-to-quark) contributions and their LP and lowest-order NLP parts, we recover the NLO and NNLO terms already contained in Eq. (32).
In the spin-dependent case the spacelike matrix in Eq. (34) is to be replaced by the polarized one, ∆P N s (α s ), given to NLO in [65][66][67]. The helicity evolution kernels ∆P N s (α s ) are identical to the unpolarized ones in the large-N limit at LP. This equality extends even to the first NLP (1/N ) corrections, except for a difference ∝ ln(N )/N in the NLO gq splitting function [68]. This difference, however, does not affect the dominant NLP corrections for SIDIS at NNLO. We thus conclude that the approximate NNLO corrections to be presented next apply to both the spinaveraged and the spin-dependent hard-scattering functions.

Expansion to NNLO
We are now ready to present the NNLO (O(α 2 s )) expansion for the SIDIS quark hard-scattering function near threshold, which is the main result of this paper. We insert the NLP evolved parton distributions and fragmentation functions of Eq. (32) into Eq. (30) and expand. To write our formulas compactly, we introduce We then find for the transverse hard-scattering function in the quark channel: T,(2) qq wherẽ The last term is the NLP contribution. We have kept "mixed" NLP corrections of the form ln(N )/M and ln(M )/N . Equation (38) reproduces the dominant part of the full NLO results given in [15,30], including the NLP terms. Its LP part is consistent with the results based on NLL threshold resummation presented in [15].
For the approximate NNLO terms we find Again, the last term is the dominant NLP correction. Here, two of the three powers of L arise from the LP part in the first line of Eq. (38), which then multiplies the NLO expansion of the NLP factor given in Eq. (33).
The results for the spin-dependent quark hard-scattering function near threshold are identical: for k = 0, 1, 2 and including NLP corrections. In fact, this will arguably hold to all orders of perturbation theory.
So far we have only addressed the q → q channel. As we discussed in the previous section, the off-diagonal evolution of parton distributions and fragmentation functions to NLP induces quark-gluon mixing. As a result, once we insert the NLP singlet evolution in (35) into the cross section (30), we also obtain terms withq(N, µ F )D g (M, µ F ) org(N, µ F )D q (M, µ F ). Evidently, these approximate the quark-to-gluon and gluon-to-quark channel contributions to SIDIS. The terms are of course suppressed by 1/N or 1/M , but they also carry logarithmic enhancement. We find, at NLO:ω with T R = 1/2. These expressions reproduce the corresponding full NLO transverse hardscattering functions of Refs. [15,30] at large moment variable. Again the contributions to the respective spin-dependent hard-scattering functions ∆ω (1) gq , ∆ω (1) qg are identical to the ones given in (41).
Unfortunately, the evolution method that we have used here to obtain the NLP corrections fails for the q → g and g → q channels beyond NLO. We have found this by inspecting related results for the Drell-Yan process at measured rapidity. Here, evolution gives the approximate resultω whereas the correct result is known to be [69][70][71][72][73][74][75] ω DY, (2) The difference of the two results is Interestingly, it depends only on one of the two Mellin variables. In the inclusive case, where N = M , this difference may be understood from Ref. [48] where the all-order resummation of the leading large-N contributions to the quark-gluon contribution to inclusive Drell-Yan was derived.
In the light of this, it is clear that evolution cannot correctly produce the leading NNLO terms for the SIDIS q → g and g → q channels. When expanding our corresponding results, we obtaiñ We note that the 1/N and 1/M terms in the NLO splitting functions contribute here. As already stated, the results in Eq. (45) are not expected to be complete, although it appears likely that a term identical to the one given in Eq. (44) (or with M → N ) would need to be added. It would be highly desirable to extend the work of [48] to the Drell-Yan process at measured rapidity and to SIDIS. This is of course beyond the scope of the present work. For now we therefore refrain from encouraging use of Eq. (45) in any phenomenological analysis. Our NNLO approximations given in this paper therefore only apply to the quark channel.
Appendix B presents our NLO and NNLO near-threshold results as functions ofx andẑ. These are obtained by a straightforward inverse transform of the above Mellin-space results.

Phenomenological predictions
We now turn to a few illustrative phenomenological applications of our approximate NNLO results. Here we only consider the unpolarized transverse structure function. We reserve a more detailed numerical analysis to future work [76], in which we will also investigate the phenomenology of NNLL resummation.
We first need to go back from Mellin space to x, z space. This is achieved by an inverse double-Mellin transform. The structure function F h i (x, z, Q 2 ) can be recovered from its moments F h i (N, M, Q 2 ) given in Eq. (9) in the following way: where C N and C M denote integration contours in the complex plane, one for each Mellin inverse. To be consistent with the NNLO approximation of the hard-scattering functions that we make, we also need to use NNLO parton distribution functions and fragmentation functions. For the former, we choose the CT18 NNLO set of Ref. [77], from which we also adopt the NNLO strong coupling. NNLO analyses of fragmentation functions are still scarce [2,78], partly because only the process e + e − → h + X is available at NNLO. For the present study we use the set of Ref. [2]. In order to be able to examine the sizes of the various corrections to the cross section, we stick to the NNLO sets of parton distributions and fragmentation functions also when computing LO or NLO results. Unless stated otherwise, we choose the renormalization and factorization scales as µ R = µ F = Q. Technically, in order to obtain Mellin moments of the parton distributions and fragmentation functions as needed for Eq. (10) in (46), we perform fits of a functional form P (x) to them, so that the Mellin moments of P (x) can be taken analytically. We have checked that our fits are accurate to better than 1% over the kinematic domain we are interested in.
We present results appropriate for the COMPASS experiment at CERN with c.m. energy √ s = 17.3 GeV, and for the EIC with √ s = 100 GeV. For both, we consider the process p → π + X. We compute the contribution by the transverse structure function to the SIDIS cross section, using Eq. (2) and dropping the longitudinal part. We focus on the z-dependence of the cross section and integrate over y ∈ [0.1, 0.9] and x ∈ [0.1, 0.8]. Note that we choose both x and z to be rather large so that we are safely in the threshold regime. Because of the relation Q 2 = xys, our choice of kinematics implies Q 2 > 3 GeV 2 for COMPASS, and Q 2 > 100 GeV 2 for the EIC. We furthermore require W > 7 GeV, where We note in passing that SIDIS experiments typically quote hadron multiplicities, which are ratios of the SIDIS cross section over the fully inclusive DIS one, for given kinematics. For the present paper we are interested in the actual NNLO corrections to SIDIS, so we do not compute multiplicities here. It would be straightforward to do this by computing DIS to full NNLO. We start by examining NLO, where the exact answer is of course known. In the following we normalize all results by the LO cross section. The left part of Fig. 1 presents results for COMPASS kinematics. The black line shows the ratio of the full transverse NLO cross section for the q → q channel to the LO one. As one can see, the NLO corrections show the expected strong increase toward large values of z. The dashed blue line shows the LP approximation to the NLO cross section, based on Eq. (38) but without the NLP term in the second line. The result shows overall good agreement with full NLO, indicating the dominance of the threshold regime, but has a nearly constant difference to the exact result. The agreement with full NLO becomes even much better when the dominant NLP corrections in the second line of Eq. (38) are included, as shown by the solid blue line. Clearly the full NLO is excellently approximated by this near-threshold result over the whole range in z, and especially so toward large z.
It is interesting to compare the NLO approximations based on the Mellin-space calculation (as shown so far) and on Eq. (60) inx,ẑ space. The two approximations differ by terms that are even more suppressed than the NLP terms. Nevertheless, their numerical difference is quite large, with the Mellin result yielding a far better approximation to the exact NLO result than the approximatex,ẑ space result. We thus conclude that Mellin space appears better suited for obtaining accurate approximations to the full result. Similar conclusions were obtained for other processes, such as for Higgs boson production [31].
The right part of Fig. 1 shows corresponding results for the EIC. They have a very similar trend as our COMPASS results, with a slightly reduced size of the corrections near threshold. This is expected due to the larger Q 2 relevant at the EIC. Again, the NLO corrections are extremely well reproduced by the approximate ones generated by Mellin-space LP+NLP resummation. The findings in Fig. 1 provide confidence that our Mellin-space NNLO expansions based on resummation also provide an accurate approximation to the full NNLO corrections for the q → q channel. Figure 2 presents our NNLO results, again normalized to LO. Here we have included the exact NLO part of the cross section, so that the approximation only applies to the NNLO terms. The dashed line shows the result based on the LP terms at NNLO, while for the solid one we have included the dominant NLP terms as well. We also display again the curves for full NLO that were already shown in Fig. 1. One can see that the NNLO corrections become sizable as z → 1, where the threshold logarithms grow in size. As in the NLO case, there is a rather significant positive contribution to the cross section by the NLP terms, both for COMPASS and the EIC.
Inclusion of the dominant NNLO terms is expected to reduce the dependence of the cross section on the renormalization and factorization scales. Figure 3 shows the variation of the LO, full NLO and the (approximate) NNLO cross sections with scale. Here we vary independently µ F = Q/2, Q, 2Q and µ R = Q/2, Q, 2Q. Among the nine combinations this results in, we discard the two with very disparate values, that is, µ F = Q/2, µ R = 2Q and µ F = 2Q, µ R = Q/2. We then take the envelope of the remaining seven results. The figure shows the resulting bands. We present them in terms of the ratio dσ(µ F , µ R ) − dσ(µ F = µ R = Q) /dσ(µ F = µ R = Q), so that the cross section with µ F = µ R = Q always produces the zero line in the plot. The result for COMPASS ‖ (left figure) shows that around z = 0.1 the NNLO scale uncertainty is large, but does improve significantly toward higher z where it becomes better than the NLO one. It does, however, remain non-negligible even at large z. The main patterns are reproduced also for EIC kinematics (right part of the figure); however, here the scale uncertainty is overall very small at NNLO, showing a band that is much narrower than the NLO one at medium to large z. We note We have varied µ F and µ R as described in the text. We have used Q 2 > 5 GeV 2 here. Right: Same for EIC kinematics.
that we have included the NLO contributions by the qg subprocesses in the results shown in the figure, whose effects are however relatively small.
We finally note that we do not consider the SIDIS spin asymmetry here. Since the approximate NNLO corrections are identical for the spin-averaged and spin-dependent cross sections (even for the dominant NLP terms), the asymmetry is expected to be affected very little by the corrections. This was indeed already observed in the NLL study [16].

Conclusions and outlook
We have presented approximate next-to-next-to-leading order corrections to semi-inclusive DIS, p → hX. These corrections apply to the quark channel and are based on the threshold resummation formalism. We have first determined all ingredients for threshold resummation for SIDIS at next-to-next-to-leading logarithmic accuracy, extending previous work by one logarithmic order. As SIDIS is characterized by two "scaling" variables,x = Q 2 /2p · q andẑ = p · p c /P · q, the moment-space resummation is naturally formulated in terms of two Mellin moments N and M . Although these are separate variables, the SIDIS resummation formula may be obtained by that for the inclusive Drell-Yan cross section by a simple rescaling N → √ N M , up to differences associated with the fact that the virtual photon in SIDIS is spacelike. These differences are accounted for by the hard factor in the resummed SIDIS cross section, which is related to the spacelike quark form factor rather than the timelike one contributing to Drell-Yan. We have subsequently expanded the resummed expressions to O(α 2 s ), to obtain the NNLO corrections.
To further improve the accuracy of the near-threshold NNLL and NNLO approximation we have determined also the dominant subleading terms which are suppressed by 1/N or 1/M near threshold, but still enhanced by logarithms. These "next-to-leading power" terms may be obtained by use of the DGLAP evolution of the parton distribution functions between scales Q and Q/ √NM . We have found that the approximate NNLO corrections are identical for the spinaveraged (transverse) cross section and the longitudinally polarized one, even including the NLP corrections. This has important ramifications for phenomenology as it means that the SIDIS spin asymmetry will be largely unaffected even by NNLO corrections.
We have presented a few basic phenomenological results at approximate NNLO. These indicate a significant increase of the cross section at large z, as well as a still sizable contribution of the NLP corrections. Our results are readily suited for initial studies of SIDIS at NNLO in "global" fitting frameworks for fragmentation functions and/or parton distributions, especially polarized ones. Furthermore, the corrections we have derived will provide important benchmarks for future full NNLO calculations of SIDIS.
There are several avenues for future improvements on our work. Extension to approximate N 3 LO near threshold and to N 3 LL resummation would be quite straightforward [76]. As already mentioned earlier, it will also be important to address the quark-gluon channels to SIDIS and to determine their dominant NLP corrections, following the lines in Ref. [48]. In the same vein, the longitudinal SIDIS structure function should be addressed at higher orders. For inclusive DIS, F L receives corrections as large as α 2 s ln 2 (1 − x) at high x, which were derived and extended to all orders in Ref. [45]. Although these are again NLP corrections, it will be relevant to investigate the corresponding logarithmic structure of F L in SIDIS. Finally, we note that the corrections we have derived here are really valid when both x and z are large. The recent study [63] considers the Drell-Yan cross section at measured rapidity and derives a factorization theorem that is valid when only one of the two kinematic variables √ z e ±y is large, while the other can have an arbitrary value. Extension of such a theorem to the SIDIS case when only x or z is large would be quite valuable as it would extend the validity of the threshold approximation for SIDIS.

A Appendix: Coefficients for resummation to NNLL
We use the following expansion of the running strong coupling [31,79] where with N f the number of flavors and The various functions we use in the main text have the following perturbative expansions: A q (α s ) = α s π A (1) q + α s π where to NNLL we use the coefficients of A q from [79][80][81][82][83]: Furthermore [31,32,79,84], and [85] P (1) q,δ Finally, the coefficient C qq in Eq. (18) is expanded as with [31] C (1) B Appendix: Near-threshold results in x, z-space Performing a double-inverse Mellin transform of the results given in Eqs. (38) and (39) we obtain approximate results for the quark-to-quark hard-scattering function ω T qq (x,ẑ, α s (µ R ), µ R /Q, µ F /Q) in Eq. (4), valid near threshold. To obtain compact expressions, we introduce the following abbreviations: We write and ∆ (2),N f qq We stress that the results in Eqs. (62), (64), (65) collect all double distributions inx andẑ that arise at NNLO. We also note that the Mellin-space and the x, z-space expressions near threshold are not strictly identical, but differ by terms that are suppressed near threshold. These terms are generically of the form ln m (N )/N 2 or 1/N (without logarithms) in Mellin space (and likewise with N replaced by M ), and of the form (1 −x) ln m (1 −x) or constant (and also withx replaced byẑ) in x, z-space.