Threshold resummation at N$^{3}$LL accuracy and approximate N$^{3}$LO corrections to semi-inclusive DIS

We advance the threshold resummation formalism for semi-inclusive deep-inelastic scattering (SIDIS) to next-to-next-to-next-to-leading logarithmic (N$^{3}$LL) order, including the three-loop hard factor. We expand the results in the strong coupling to obtain approximate next-to-next-to-next-to-leading order (N$^{3}$LO) corrections for the SIDIS cross section. In Mellin moment space, these corrections include all terms that are logarithmically enhanced at threshold, or that are constant. We also consider a set of corrections that are suppressed near threshold. Our numerical estimates show modest changes of the cross section by the approximate N$^{3}$LO terms, suggesting a very good perturbative stability of the SIDIS process.


Introduction
The semi-inclusive deep-inelastic scattering (SIDIS) process p → hX has become a widely used probe of hadronic structure and hadronization phenomena. Its main uses are extractions of (polarized) parton distribution and fragmentation functions or combinations thereof [1][2][3][4][5][6][7][8][9][10]. In global analyses of these quantities SIDIS data can add useful information on, for example, the flavor structure of the sea quarks. The future Electron Ion Collider (EIC) will allow precise measurements of SIDIS observables over wide kinematic regimes [11].
In a recent paper [12], we have studied higher-order QCD corrections to the SIDIS cross section. Our approach was to use the threshold resummation formalism for SIDIS and carry out fixed-order expansions of the resummed cross section. Threshold resummation for SIDIS was originally discussed in Ref. [13] and then further developed in more general terms in [14] and [15]. These papers formulated the resummation at next-to-leading logarithmic (NLL) accuracy. In [12] we extended the resummation to next-to-next-to-leading logarithm (NNLL), which also allowed us to obtain approximate fixed-order corrections to the hard scattering cross section for SIDIS at next-to-next-to-leading order (NNLO) level. These results were used recently to obtain the first NNLO set of fragmentation functions fit "globally" to SIDIS and electron-positron annihilation data [10].
The purpose of the present paper is to advance our previous study to N 3 LL and to again use the resummed cross section to derive approximate fixed-order corrections to the SIDIS cross section, in this case at N 3 LO. Our main motivation for this analysis is to further improve the perturbative framework for SIDIS and to set the stage for precision analyses of SIDIS data from the future EIC in terms of parton distributions or fragmentation functions at high perturbative order. While such analyses at N 3 LO may presently still seem far off, the study of the perturbative stability of the SIDIS cross section and its associated threshold resummation is in any case valuable. This becomes indeed possible by going to N 3 LL and N 3 LO and carrying out comparisons with lower orders. We also note that in our previous paper [12] we presented phenomenological results only for the fixed-order (NNLO) corrections. Here we wish to carry out numerical studies also for the resummed case, which provides another motivation for this study.
In Sec. 2 we give an overview of the kinematics of the process, introducing Mellin moments. Section 3 describes the threshold resummation framework. Section 4 is dedicated to the derivation of the three-loop hard factor to be used for obtaining N 3 LL or N 3 LO results. In Sec. 5 we carry out the expansion of the resummed results to N 3 LO. Finally, Section 6 presents some numerical studies in the EIC kinematical regime.

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 x = Q 2 2P · q , y = P · q P · k , z = P · P h P · q . (1) We have Q 2 = xys, with √ s the center-of-mass energy of the incoming electron and proton. We follow Ref. [14] to write the spin-averaged SIDIS cross section as 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. In what follows we will only treat the transverse structure function in the q → q orq →q channels, which is the only channel that appears already at the lowest order (LO) of perturbation theory. We write all equations for the spin-averaged case, although they will equally apply to the helicity-dependent one [12,15].
Using factorization, the unpolarized structure functions may be written as double convolutions. For example, for the transverse one we have wheref , As a result the structure functions can be obtained from the moments of the parton distribution functions and fragmentation functions, and the double-Mellin moments of the partonic hardscattering functions.
For the perturbative expansion given in Eq. (4) we have in moment space at lowest order according to Eq.
The corresponding moments of the next-to-leading order (NLO) terms ω T,(1) f f may be found in Refs. [14,16]. In the following, we consider logarithmic higher-order corrections to the hardscattering functions that arise at large values ofx andẑ or, equivalently, at large N and M .

Threshold resummation at N LL accuracy
The resummation of threshold logarithms for SIDIS was extensively studied in Refs. [12][13][14]17]. The NNLL resummation formula for the unpolarized SIDIS transverse structure function was discussed in Ref. [12]. The resummed partonic transverse structure function takes the form which actually holds to any logarithmic order. As stated earlier, our goal is to set up the formalism for resummation to N 3 LL. In Eq. (9) we havē with the Euler constant γ E . Each of the functions A q , D q , H SIDIS qq , C qq is a perturbative series in the strong coupling. We write the corresponding expansions generically as where q = 0. To achieve N 3 LL accuracy, we need A q to order α 4 s and all other functions to order α 3 s . The corresponding coefficients are collected in Appendix A. The main new ingredient not directly known from the literature is thē N ,M -independent coefficient H SIDIS,(3) qq whose derivation will be presented below in Sec. 4. The other prefactor C qq in Eq. (9) collects all moment-independent terms of the resummed exponent; see [12,18,19]. The formulas needed for its derivation to order α 3 s may be found in Ref. [18].
In order to explicitly obtain the structure function resummed to N 3 LL we now expand the exponents in Eq. (9) appropriately. The operations are quite standard. We obtaiñ where The functions h (k) q impart resummation to N k−1 LL accuracy. The first three are well known in the literature: The function h (4) q , needed for N 3 LL resummation, is found to be This result is in agreement with that given in Ref. [20] for the Drell-Yan process. For simplicity, we have set the renormalization and factorization scales to Q. The results presented in Eq. (12) may be used to obtain N 3 LO (that is, O(α 3 s )) expansions of the hard-scattering functionω T qq . This expansion will be carried out in Sec. 5.
We stress that all terms generated by Eq. (12) are either logarithmic or constant near threshold. The full hard-scattering function in Mellin space will, at any order in perturbation theory, also contain terms that are suppressed by powers of 1/N and/or 1/M . Such terms are often referred to as "next-to-leading power (NLP)" corrections. As discussed in our previous paper [12] (see also references therein), one can straightforwardly account for the dominant NLP terms by multiplying the resummed cross section in Eq. (9) by the two factors exp − where the coefficients −C F /(2N ) and −C F /(2M ) in the exponents correspond to the NLP terms in the LO diagonal evolution kernels for the quark parton distributions and the quark fragmentation functions, respectively. At N 3 LO the two exponential factors, when combined with the resummed exponents in Eq. (9), will generate all terms of the form α 3 s ln n (N ) ln m (M )(1/N + 1/M ), with n + m = 5.

The hard factor at three loops
The factor H SIDIS qq is derived from the finite part of the virtual corrections to the process γ * q → q. The basic ingredient is the renormalized spacelike form quark factor, from which one needs to subtract the infrared divergencies via a suitable method developed in Refs. [21,22]. For our present purposes, we will need the renormalized three-loop form factor, which was derived in [23][24][25] ‡ and reads in dimensional regularization with d = 4 − 2 space-time dimensions: where 320 ‡ We note that recently even the four-loop results were published [26].
Here we have kept terms of higher order in in the one-loop and two-loop results since these turn out to make finite contributions in the end. In the above expressions, ζ(j) is the Riemann zeta function, N f is the number of flavors, and C F = 4/3, C A = 3. For purely electromagnetic interactions the factor N f,V =γ becomes [24] As shown in Refs. [21,22], the hard coefficient may be extracted from the form factor in the following way. Adapted to the case of SIDIS we have from [22] H SIDIS 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 [22] 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. (20). 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 N 3 LL accuracy with and The coefficients b 0 and b 1 can be found in Appendix A. Inserting all terms into Eq. (20) and expanding in α s , all poles in powers of 1/ cancel. The final expression for H SIDIS qq up to three loops can be found in Appendix A.

Expansion to N 3 LO
We are now ready to present the N 3 LO (O(α 3 s )) expansion for the SIDIS q → q hard-scattering function near threshold. To write our formulas compactly we introduce The coefficientsω T,(1) qq andω T,(2) qq in Eq. (4) were already given in our previous paper [12]; for completeness, we recall them in Appendix B. For the approximate N 3 LO terms we find: As before, we have set µ R = µ F = Q for simplicity. We stress that the corrections given by this expression include all terms that are logarithmically enhanced at threshold, or that are constant.
In physical space these are terms with double distributions (that is, "plus" distributions and δ-functions) inx andẑ.
The last term in Eq. (27) represents the dominant NLP contributions. Note that upon expansion beyond NLO the exponential factors in (16) will also generate terms with inverse powers 1/N 2 , 1/M 2 and higher, which we have discarded for consistency since they are far beyond the approximations we make. We will see later that these terms are numerically very small.

Phenomenological predictions
We will now present some phenomenological predictions for the transverse SIDIS cross section at NNLL and N 3 LL, as well as for the expansion to N 3 LO. We will also compare to our previous NNLO results of [12].
In order to obtain results for the transverse structure function F h T (x, z, Q 2 ) in physical x, z space we need to invert its Mellin momentsF h T (N, M, Q 2 ) in Eq. (6). This is achieved by the inverse double-Mellin transform where C N and C M denote integration contours in the complex plane, one for each Mellin inverse. We adopt the minimal prescription of Ref. [27] to treat the Landau pole present in the resummed exponents in Eqs. (12), (14) at λ N M = 1, or (see Eq. (13)), According to the minimal prescription, the two contours need to be chosen such that all singularities in the complex plane lie to their left, except for the Landau pole. We parameterize the two contours as with ζ, ξ ∈ [0, ∞] as contour parameters, where c N = 1.8 and c M = 3.3. We furthermore choose φ N = 3π/4; the two signs for the N -contour in (30) select the two branches in the complex plane. As N moves along its contour, the position of the Landau pole relevant for the M -integral will move as well, mapping out a trajectory in the plane. This implies that the angle φ M needs to be chosen as a function of N , so that during the M integration this trajectory is never crossed. A more detailed description of the inverse double-Mellin transform can be found in Ref. [14].
We note that we only consider the transverse structure function in (2) here. The longitudinal one is suppressed near threshold, even beyond the dominant NLP terms we have included for F h T . While it would be very interesting to also investigate higher-order corrections to F h L , this is beyond the scope of this work. In what follows, we also discard the contributions by the q → g and g → q channels to the structure function. These are fully known only to NLO. We could include the contributions at NLO level in our approximate NNLO, N 3 LO results to be presented below, but this would simply amount to a uniform shift of all results by a few per cent, which is not really relevant for our main goal of analyzing the structure of higher-order contributions in the q → q channel.
For the parton distribution functions and fragmentation functions we choose the NNLO sets of Ref. [28] and Ref. [29], respectively. Clearly, in order to present true N 3 LL or N 3 LO results, we would need PDFs and FFs evolved at those same orders, which are currently not yet available. We therefore keep the renormalization and factorization scales at µ R = µ F = Q and do not investigate the scale dependence of our results. In this sense we use the NNLO parton distributions and fragmentation functions as templates for the N 3 LO ones, which should be adequate for a first analysis of the beyond-NNLO effects we are interested in. We note that the scale dependence of the transverse SIDIS cross section was anyway found to be rather small already at NNLO in Ref. [12]. Note that we "match" all results to NLO, so that the NLO corrections for the q → q channel are always included exactly.
Our predictions will refer to the unpolarized p → π + X process appropriate for the future EIC with √ s = 100 GeV. We focus on the z-dependence of the cross section and integrate over y ∈ [0.1, 0.9] and x ∈ [0.1, 0.8]. We choose 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 > 100 GeV 2 for the EIC. We furthermore require W > 7 GeV, where W 2 = Q 2 (1 − x)/x + m 2 p , with m p the proton mass.
We begin by comparing fully resummed results obtained at various different levels of logarithmic accuracy. The upper left part of Fig. 1 shows the NLL, NNLL, and N 3 LL resummed cross sections as functions of z, normalized to the LO one. As one can see, the NNLL terms lead to a significant enhancement over NLL, while the additional terms arising at N 3 LL only lead to a modest further increase of the cross section. This result demonstrates that the resummed SIDIS cross section has excellent perturbative stability.
We can further investigate the improvements provided by going to NNLL and N 3 LL. To this end, we note that even at a given logarithmic order the resummation formula in Eq. (9) may actually be used in various ways that are all equivalent in terms of their perturbative content, but differ numerically. Let us refer to the corresponding choices as resummation schemes. We consider three such schemes: Scheme (a) Here we use Eq. (9) as written. That is, we keep the functions H SIDIS qq and C qq as separate factors, each its own perturbative series of the form (11). Also, we use the Mellin moments N and M precisely in the formN andM as defined in (10). This scheme has been used for the first plot in Fig. 1.

Scheme (b)
Here we expand the product H SIDIS qq × C qq in Eq. (9) strictly to the desired order.
That is, suppose we are at NLL where H SIDIS qq = 1 + αs π H SIDIS,(1) qq qq ) and drop terms of O(α 2 s ). We continue to use the variablesN andM .

Scheme (c)
Here we first use the expansion of H SIDIS qq × C qq as for scheme (b). In addition, we use (10), (13) to write The terms with the Euler constant lead to modifications of the functions h (k>1) q in Eq. (14) [18,20]. They evidently also generate non-logarithmic corrections in the resummed exponent.
These may be expanded out perturbatively, so that they migrate from the exponent to an N, M -independent prefactor. This prefactor is then expanded along with the factor H SIDIS qq × C qq into a single perturbative function that now multiplies the resummed exponent, the latter now being a function of N and M rather than ofN andM .
It is immediately clear that the three resummation schemes are indeed equivalent for a given logarithmic accuracy. The remaining three plots in Fig. 1 compare the three schemes at NLL (upper right), NNLL (lower left), and N 3 LL (lower right). It is striking to see how the difference among the three schemes is still rather large at NLL, then strongly decreases at NNLL, and finally becomes extremely small at N 3 LL. Of course, one does expect the details of how the expansions are performed to matter less and less with increasing logarithmic order. Nevertheless, the level at which the resummed predictions become independent of the resummation scheme at NNLL and especially at N 3 LL is truly remarkable.
Encouraged by these observations, we now turn to fixed-order expansions of our resummed results. Figure 2 (left) shows again the NNLL-resummed result for scheme (a), along with its expansion to NNLO as given by Eqs. (42) and (43) (black solid line) and already obtained in Ref. [12]. All results are again normalized to the LO cross section. We note that finite-order expansions are independent of the resummation scheme chosen. We observe that resummation within scheme (a) leads to a suppression of the cross section at lower z and to the expected enhancement at high z where the threshold logarithms become particularly important. In addition to these two results, we also expand the resummed cross section numerically to orders α 2 s and α 3 s . As expected, the result for the O(α 2 s ) expansion (dash-dotted line) is extremely close to the NNLO one. The only difference between these two results comes from the fact that the formal expansion of the NLP factors in (16) will produce also terms with higher inverse powers of N and M , as noted at the end of Sec. 5. These terms are not included in our explicit NNLO expansions, but do contribute to the numerical O(α 2 s ) expansion of the cross section. As one can see by comparing the two corresponding curves, they are of very small size. The dashed line in the left part of Fig. 2 shows the O(α 3 s ) expansion of the NNLL resummed cross section. We observe that this result is already very close to the full NNLL one, indicating that terms of order O(α 4 s ) or higher are small.
The right part of Fig. 2 presents the same analysis one order higher. We show the N 3 LLresummed cross section for scheme (a), and now the expansion to N 3 LO as given by Eqs. (42),(43) and (27). This time, we numerically expand the N 3 LL result to orders α 3 s and α 4 s . Again the numerical expansion to O(α 3 s ) essentially coincides with the approximate N 3 LO one, up to tiny corrections suppressed as 1/N 2 , 1/M 2 or higher. The result at O(α 3 s ) is almost indistinguishable from the full N 3 LL-resummed one, demonstrating again that corrections beyond third order are all but negligible. We note that the O(α 3 s ) expansion obtained from the N 3 LL-resummed result is more complete than the O(α 3 s ) expansion shown in the left part of Fig. 2: It contains all seven "towers" of threshold logarithms, that is, terms of the form α 3 s ln n (N ) ln m (M ) with 0 ≤ n+m ≤ 6, whereas NNLL resummation can only correctly reproduce the five towers with 2 ≤ n + m ≤ 6.

Conclusions and outlook
We have explored higher-order QCD corrections to the quark-to-quark hard-scattering cross section relevant for semi-inclusive DIS. We have developed the threshold resummation framework for SIDIS to N 3 LL accuracy, hereby extending previous work carried out at NNLL [12]. Among the main tasks to be completed for achieving N 3 LL resummation was the derivation of the three-loop hard factor from the spacelike form factor. We have used our N 3 LL results to derive approximate N 3 LO corrections for SIDIS. These corrections contain all seven "towers" of threshold logarithms that are present at this order. We have also included dominant subleading logarithmic terms that are suppressed near threshold.
We have presented phenomenological results for resummed and approximate fixed-order SIDIS cross sections for EIC kinematics. These show an excellent perturbative stability of the cross section in the sense that the N 3 LL cross section is only modestly enhanced over the NNLL one, and that generally corrections beyond O(α 3 s ) seem unimportant. A particularly striking result is that the actual treatment of resummation, in terms of how the relevant expansions are carried out in practice, matters less and less when the logarithmic accuracy of resummation increases, so that the N 3 LL result is essentially insensitive to the resummation scheme adopted. Clearly, our results show that the SIDIS cross section may serve as an excellent testbed for studies of higher orders in perturbation theory. We believe that our results are a valuable addition to the general "library" of QCD observables that are known to NNLO and beyond.
Future extensions of this work should also address non-perturbative power corrections to the SIDIS cross section, very little about which is currently known. It would be an interesting phenomenlogical study to confront experimental data with our perturbative results at various high orders ranging from NLO to N 3 LL, ascertaining how the size of phenomenlogically extracted power corrections depends on the order of perturbation theory that is employed.
We finally note that while we have focused our studies entirely on the spin-averaged SIDIS cross section, all our results equally apply to the helicity-dependent one. More precisely, the N 3 LL result and hence its approximate N 3 LO expansion are identical in the spin-averaged and the spindependent cases. This further corroborates the finding of Ref. [15] that the SIDIS spin asymmetry is insensitive to higher-order perturbative QCD corrections.

A Appendix: Coefficients for resummation to N 3 LL
We use the following expansion of the running strong coupling [18,30] where with N f the number of flavors and For the b 3 coefficient we have inserted the values of C F , C A ; the full expression can be found in [31]. The relevant expansion coefficients for A q in Eq. (11) read [30,[32][33][34][35][36][37][38]: with Furthermore for the expansion of the function D q we have [18,19,22,30,39], The expansion coefficients for C qq in Eq. (9) read [18] C (1) Finally, for H SIDIS qq we find for an arbitrary renormalization scale µ R , but for µ F = Q: B Appendix: NLO and NNLO hard-scattering coefficient functions The full NLO coefficient function can be found in [14,16]. Its near-threshold approximation was given in [12] and reads where L ≡ 1 2 ln(N ) + ln(M ) . The corresponding approximate NNLO result is given by