Soft corrections to inclusive deep-inelastic scattering at four loops and beyond

We study the threshold corrections for inclusive deep-inelastic scattering (DIS) and their all-order resummation. Using recent results for the QCD form factor, related anomalous dimensions and Mellin moments of DIS structure functions at four loops we derive the complete soft and collinear contributions to the DIS Wilson coefficients at four loops. For a general SU(nc) gauge group the results are exact in the large-nc approximation and for QCD with nc = 3 we present precise approximations. We extend the threshold resummation exponent GN in Mellin-N space to the fifth logarithmic (N4LL) order collecting the terms αs3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\alpha}_{\mathrm{s}}^3 $$\end{document}(αs ln N)n to all orders in the strong coupling constant αs. We study the numerical effect of the N4LL corrections using both the fully exponentiated form and the expansion of the coefficient function in towers of logarithms. As a byproduct, we derive a numerical result for the complete pole structure of the QCD form factor in the parameter of dimensional regularization ε at four loops.


Introduction
The cross section for inclusive deep-inelastic scattering (DIS) of charged leptons or neutrinos off a nucleon target can be expressed in terms of structure functions, which factorize into the product of operator matrix elements for the nucleon under consideration and the DIS Wilson coefficients. The latter parametrize the hard partonic scattering cross section driven by short-distance physics and are calculable in perturbative Quantum Chromodynamics (QCD) order by order in the strong coupling constant α s . However, only a few terms in this expansion can be calculated completely. Currently, the QCD corrections at next-to-next-to-next-to-leading order (N 3 LO) for the neutral-and charged-current DIS structure functions F 1 , F 2 and F 3 are known [1,2]. In the exceptional region of phase space, on the other hand, when the Bjorken scaling variable x is close to unity, x → 1, the DIS Wilson coefficients develop large Sudakov logarithms at all orders in the form α n s ln l (1 − x)/(1 − x) + , with l = 2n − 1, . . . , 0. These are a consequence of the constrained phase space available for the real emission of soft and collinear gluons.
The DIS Wilson coefficients are subject to an additional factorization for x → 1 (or N → ∞ in Mellin-N space) which allows for resummation of those Sudakov double logarithms [3][4][5][6][7]. Based on the fixed-order QCD results up to N 3 LO for the DIS structure JHEP03(2020)116 functions [1,2] this resummation has been performed to next-to-next-to-next-to-leading logarithmic (N 3 LL) accuracy in Mellin-N space [8], where the resummation takes the form of an exponentiation which organizes the respective logarithms α n s ln l N , with l = 2n, . . . , 1 (see also [9][10][11][12] for related work in the Soft-Collinear Effective Theory).
From a phenomenological point of view, resummation in DIS is essential to obtain reliable QCD predictions in the (small) region of phase space dominated by such large threshold logarithms. This is of particular relevance since DIS structure functions also receive power corrections beyond the QCD factorization at leading twist and such highertwist contributions are supposed to be significant at large x. Therefore, good control of the perturbative series at leading twist is important to distinguish those effects in the analyses of DIS data, see for instance [13].
From a field-theoretic point of view, on the other hand, the threshold limit in DIS is of interest, since the factorization of the DIS Wilson coefficients for x → 1 links a number of fundamental quantities in QCD with each other due to the universality of the soft and collinear limit. In particular, this comprises the QCD form factor which features an all-order exponentiation in dimensional regularization with d = 4 − 2ε dimensions [14][15][16][17][18][19] and the QCD splitting functions, currently fully known to next-to-next-to-leading order (NNLO) [20,21]. The necessary cancellation of soft and collinear singularities in ε in inclusive observables due to the Kinoshita-Lee-Nauenberg theorem [22,23] offers a constructive approach to the DIS Wilson coefficients near threshold [24], see also [25][26][27].
Recent progress in the computation of QCD corrections has been obtained in particular for the quark form factor [28] and the non-singlet quark splitting functions [29], which have been calculated at four-loop order in the limit of large-n c for a general SU(n c ) gauge group. In addition, a low number of fixed Mellin moments for the DIS structure functions are available at four loops [30], as well as information on the n f -dependence of the QCD form factor [30][31][32][33][34] and on new colour factors (quartic Casimirs) [33,[35][36][37][38][39]. Taken together [40], these results allow us to extend the threshold resummation exponent G N for DIS Wilson coefficients in Mellin-N space to the fifth logarithmic (N 4 LL) order. At this accuracy we collect and resum in G N the terms α 3 s (α s ln N ) n to all orders in α s . To that end, we extract the resummation exponent B DIS of the quark jet function collecting final-state collinear emissions at the fourth order in α s , and we address the term α 4 s /(1−x) + in the fixed-order expansion of the DIS Wilson coefficients at four loops. With the help of the available four-loop QCD results, we can provide exact expressions in the large-n c limit and adequate numerical approximations for full QCD with n c = 3. As a byproduct, we also derive a numerical result including full colour dependence for the single pole 1/ε of the dimensional regulated QCD form factor at four loops.
This article is organized as follows. In section 2 we recall the general structure of the threshold resummation for DIS and provide the required formulae to perform resummation to N 4 LL accuracy. In section 3 we describe the extraction of the resummation coefficients and the calculation of the DIS Wilson coefficient in the soft and collinear limit. We present a numerical study for the exponentiated resummed DIS Wilson coefficients and their tower expansion in section 4 and summarize in section 5. Exact expressions for the DIS Wilson coefficients are given in appendix A.

JHEP03(2020)116 2 Theoretical framework
The general structure of the resummed DIS Wilson coefficients takes the following form in the Mellin N -space [3,4], for the neutral-or charged-current DIS structure functions F a , where a = 1, 2, 3. In Mellin space they factorize in a sum over all partons i as F N a (Q 2 ) = C N a,i (Q 2 )f N i with f N i denoting the Mellin transformed parton distributions f i (x) for the parton i, and the Born contributions to the are normalized as C N a,q | LO = 1. Here the resummation exponent G N contains all terms of the form ln k N to all orders in α s , and the pre-factor g 0 encompasses the N -independent terms. G N and g 0 depend on the physical hard scale Q 2 ( = −q 2 in DIS, with q the four-momentum of the exchanged gauge boson), and on the renormalization and factorization scales µ and µ f , suppressed for brevity. For inclusive DIS G N reads in the notation of [8,41] G N = ln ∆ q + ln J q + ln ∆ int DIS , (2.2) where the radiation factors (∆ q , J q , ∆ int DIS ) are given by well-known integrals over functions of the running coupling. The collinear soft-gluon radiation off an initial-state quark is collected by with the light-like cusp anomalous dimension A q (α s ) addressed below. The collinear emissions from an 'unobserved' final-state quark are summarized in the jet function, including the additional function B DIS (α s ). All expansions in terms of the strong coupling α s are normalized as etc. Any process-dependent contributions from large-angle soft gluons emissions are contained in the function ∆ int DIS which, however, evaluates to unity, i.e., ∆ int DIS = 1, since the corresponding evolution kernels vanish at all order in α s for inclusive DIS [42,43]. Thus, the last term in eq. (2.2) is absent for inclusive DIS. The running of the strong coupling is governed by the QCD beta function with β 0 = (11/3)C A − (2/3)n f and so on.

JHEP03(2020)116
The evaluation of the integrals in eqs. (2.3) and (2.4) proceeds in a standard manner, see for instance [8,44], leading to the following expression of the resummation exponent G N , (λ) collects all large logarithms up to α 3 s (α s ln N ) n in G N . The complete results for g DIS i with i ≤ 5, including for consistency the well-known lower order results [3,4,8,44,45], read For brevity, we use here the short-hand notations λ 1 = 1 − λ, L λ 1 = ln(1 − λ), L qr = ln(Q 2 /µ 2 ) and L fr = ln(µ 2 f /µ 2 ) . Also we suppress factors of β 0 which can be restored easily by the substitutions In addition, g 3 needs to be multiplied by β 0 and g 4 by β 2 0 .

JHEP03(2020)116
The new function g 5 (computed with methods of [8] and to be multiplied by β 3 0 ) is given by Note that for resummation to N 4 LL accuracy the function g DIS 5 (λ) needs the five-loop coefficient β 4 of the QCD beta function in eq. (2.6) which is available due to [46][47][48][49], as well as the cusp anomalous dimension up to five loops, i.e., the coefficient A 5 which has recently been estimated [50] (see below). In addition, one also needs the evolution kernel B DIS (α s ) of the jet function in eq. (2.4) to four-loop order, i.e., the term B DIS 4 , which will be addressed below. We will collect and discuss all necessary resummation coefficients in the next section.

JHEP03(2020)116 3 Resummation coefficients
We present results for a general SU(n c ) gauge group with n c colours and n f massless fermions; the QCD expressions can always be recovered by setting n c = 3. In this way, the resummation coefficients are expressible in terms of the usual SU(n c ) quadratic Casimir factors C A = n c and C F = (n 2 c − 1)/(2n c ), i.e., C A = 3, C F = 4/3 in QCD. Starting from three-loop order in perturbation theory, higher group invariants enter, such as the square of the symmetric part of the trace of three SU(n c ) generators T a F in the fundamental representation, , where x, y labels the representations with generators T a r and In QCD these factors evaluate to d FA /n c = 5/2 and d (4) FF /n c = 5/36.

Splitting functions at large x
The coefficients of the cusp anomalous dimension A n appear in the large-x expansion of the diagonal parts of the splitting functions. For a quark field, we consider the large-x behaviour of the non-singlet splitting functions P (n−1) ns (x) in the MS scheme at n-loops, recall eq. (2.5) for the normalization of the expansion in powers of α s . Disregarding terms that vanish for x → 0, these can be written [52,53] where 1/(1 − x) + on the right hand side represents the usual plus-distribution; B q n is sometimes referred to as virtual anomalous dimension.
The cusp anomalous dimension of a quark field A q is known up to third order for quite a long time [20] and the perturbative expansion reads according to eq. (2.5),

JHEP03(2020)116
The four-loop contribution A 4 has been of subject to intensive recent studies. Combining all available results it reads for quark fields, This expression combines available exact results in the large-n c limit of QCD [28,29], as well as for the terms proportional to n f [32,33], to n 2 f [31,54] and to n 3 f [55,56], the latter have been known for a long time. The quartic colour factors have been completed in exact form recently [33,36,38], and found to be in agreement with numerical estimates [29,35].

JHEP03(2020)116
The δ(1 − x) parts of the n-loop non-singlet splitting function P (n−1) ns (x), i.e. the coefficients B q n in eq. (3.4) are also known exactly to three loops [20], At four loops the exact result in the large-n c limit of QCD with an overall factor of C F reads [29] B q 4 Lnc The full colour dependence can be parametrized with coefficients b q  [29], thanks to the computation of more Mellin moments of the corresponding anomalous dimension and the exact result for A 4 in eq. (3.6). All parts proportional to n 2 f and n 3 f are known in analytic form [54,55]. The result for B q 4 can then be written as 51.03056 2.261237 Table 1. Numerical values for the colour coefficients of the δ(1 − x) part B q 4 in eq. (3.4) at fourth order. All exact values for the coefficients have been rounded to seven digits and the errors are correlated due to the known exact results at n f 0 and n f 1 in the large-n c limit given in the third row. in eq. (2.9) relevant at N 4 LL accuracy for the jet function in eq. (2.4) we have to combine the results summarized above with those for the QCD form factor and the DIS cross sections in the soft and collinear limit at four-loop order; this will be done next.

Quark form factor
The quark form factor summarizes the QCD corrections to the vertex of a photon with virtuality Q 2 and a massless external quark / anti-quark, the relevant amplitude being

JHEP03(2020)116
where e q denotes the quark's electric charge and the scalar function F the quark form factor. As we are interested in neutral-and charged-current DIS, we note that analogous expressions hold for general vector or axial-vector currents with the appropriate replacements of the couplings e q → v q and e q → a q , i.e., for DIS including Z-boson or W ± -boson exchange. The latter case is distinguished by certain diagram classes involving higher group theory invariants which will be addressed below. F is gauge invariant, but divergent. It has been computed in dimensional regularization with d = 4 − 2ε at four-loop order in the large-n c limit of QCD [28]. In addition, all terms proportional to n 2 f [31] and the quartic colour factor d FF /n c at four loops are also known [36,37]. For lower order results to sufficient depth in ε, see, e.g., [19,[57][58][59]. The exponentiation of F is achieved by solving the well-known evolution equations [14][15][16][17][18] where µ again represents the renormalization scale. Here K is a counter term containing all the poles in ε whereas the function G is finite in the limit ε → 0. The functions G and K follow the renormalization group equations [14], with the standard cusp anomalous dimension A q discussed above. The solution of eqs. (3.14) and (3.15) order by order in perturbation theory is straightforward, see, e.g., [18,19,26,27], and provides a perturbative expansion of the bare (unrenormalized) quark form factor in terms of the bare strong coupling a b s (recall the normalization a s = α s /(4π) in eq. (2.5)), as In terms of the n-th order coefficients A n in eqs. (3.5) and (3.6) and the coefficients G n (ε) of the function G in eq. (3.15), which are still functions of the parameter ε, the expansion of eq. (3.16) up to four loops leads to [19,26],

19)
JHEP03(2020)116 As observed at lower fixed orders [60,61] and generalized in [62] the function G in eq. (3.15) generates the subleading poles in ε at each order. G is the sum of three terms: twice the coefficient B q in eq. (3.4) of the δ(1 − x) part in the relevant parton splitting function, the single-logarithmic anomalous dimension of the eikonal form factor, and a term associated with the QCD beta function. Thus, the perturbative coefficients G n (ε) satisfy the following relations where the functions f q 0n (ε) at n loops are polynomials in ε. For consistency, their lower order expressions are needed to sufficient depth in ε in eqs. (3.17)- (3.20).
The finite coefficients f q n in the eq. (3.21) are related to the anomalous dimension of the eikonal form factor [62], see also the recent work [63]. In analogy to the cusp anomalous dimensions A q , they exhibit a maximally non-Abelian colour structure. This implies that the n f -independent colour factor for quarks in f q n is proportional to C F C A (n−1) up to three loops, where the explicit results read [61] f q 1 = 0 ,

JHEP03(2020)116
With the currently available results for B q 4 in eq. (3.11) and G 4 in eq. (3.21) in the planar limit we can determine f q 4 for quarks completely in the large-n c limit as With these results the perturbative expansion of f q through four loops reads where we have used the full QCD coefficients f q n in eq. (3.22) and, as indicated, the largen c result of eq. (3.23). Based on experience with other anomalous dimensions and lower orders, the latter is expected to provide a very good approximation of full QCD for the coefficients of individual powers of n f , as well as the complete result for the four loop term, putting n f = 3, 4 or 5.
With the results of section 3.1 and using all available results for the quark form factor we can, in addition, also determine the complete colour decomposition of the four-loop coefficient f q 4 . Given its relation to the anomalous dimension of the eikonal form factor, we assume here that f q 4 exhibits the same generalized maximal non-Abelian property in the presence of quartic colour factors as the cusp anomalous dimension A 4 , see [35]. This assumption, supported also by recent studies of amplitude factorization in [63], implies the absence of the colour coefficients C 4 with three still unknown coefficients f q We do have, however, a low number of fixed Mellin moments for the DIS structure functions at four loops at our disposal [30]. As will be explained in the section 3.3 below, this information can be used to constrain the DIS Wilson coefficients at large x, specifically the term proportional to the plus-distribution [1/(1−x)] + , so that we can extract numerical values for the unknown colour coefficients of f q 4 in eq. (3.25). We obtain where the errors are correlated due to the known exact results in the large-n c limit. These are the best estimations for these coefficients given our current knowledge. The four-loop contributions in eq. (3.26) display still a significant uncertainty, but in the bare form factor F 4 at four loops in eq. (3.20) they are numerically not dominant. This leads to the expression for the full colour dependence of the single pole in ε at four loops where all exact values have been rounded to seven digits and the errors are inherited from the numerical results in table 1 and eq. (3.26) and therefore are correlated. Note that the form factor F in eq. (3.13) receives additional corrections due to a new flavour structure starting from three loops, where the photon couples to a closed quark loop [1]. This requires the summation over the charges of all quark flavours in the loop, leading to additional terms with a relative factor ( q e q )/e q = n f e /e q and proportional to the group invariant (d abc d abc )/n c in eq. (3.1). For the analogous expression of a chargedcurrent, i.e., the form factor with the coupling of a W ± -boson, such contributions are absent. For the photon form factor the respective terms proportional to (d abc d abc )/n c are JHEP03(2020)116 finite at three loops, thus appearing in f q 03 in G 3 in eq. (3.21). At four loops, those terms enter in G 4 through β 0 f q 03 in eq. (3.21) and generate a single pole in ε in the bare form factor, as confirmed by the exact result for the coefficient of n f (d abc d abc )/n c in [37]. This colour factor has been omitted in eq. (3.27).
Eqs. (3.23), (3.25) and (3.27) are new results of the present paper, with the large-n c results in eq. (3.23) being exact. Eq. (3.25) for f q 4 and, as a consequence, eq. (3.27) for the single pole in ε of F 4 are based on the (by now well-supported) conjecture that f q 4 has the same generalized maximal non-Abelian colour structure as the cusp anomalous dimension.

DIS Wilson coefficients at large x
The direct link between the QCD form factor and the DIS Wilson coefficients near threshold for x → 1 through factorization in the soft and collinear limit allows to relate the knowledge on either quantity at a given order in perturbation theory. To that end, we write the partonic coefficient function W b as a series in bare coupling a b s , see [19,24]. With the convention in eq. (3.16) and choosing the scale µ = Q, the expansion coefficients W b n have the following structure up to four loops, Here F n are the space-like form factors discussed in the previous section 3.2, whereas S n denote the real emission contributions. In the limit Bjorken x → 1, the singular terms are proportional to δ(1 − x) and to the usual plus-distributions, The dependence of the pure real-emission contributions S n on the scaling variable x is given by the d-dimensional plus-distributions up to f n,ε defined by The bare strong coupling in the coefficients W b n has to be renormalized (in the MS scheme) according to With the ingredients of sections 3.1 and 3.2, using eq. (3.28), the renormalized coupling in eq. (3.31), and the known results for the DIS Wilson coefficients up to third JHEP03(2020)116 order [1] we can then derive the coefficient functions C a,q for the DIS structure functions F a for a = 1, 2, 3 in the soft and collinear limit up to four-loop order. In xspace, with distributions f i for the parton i, the DIS structure functions F a factorize as F a (x, Q 2 ) = C a,i (α s , µ 2 /Q 2 ) ⊗ f i (µ 2 ) (x), cf. eq. (2.1) for the corresponding Mellin space expressions. With the standard normalization for a = 1, 2, 3, we determine the coefficient proportional to D 0 in c a,q . For the structure function F 2 it takes the following form, 2,q is given in eq. (A.4) in appendix A, where also the lower order soft-virtual expressions for C 2,q from [8] have been collected for convenience. The expressions for the DIS Wilson coefficients in eq. (3.32) for neutralor charged-current DIS are identical up to two loops. They start to differ from threeloop order onwards, because in the case of neutral-current DIS additional corrections need to be considered, where the exchanged virtual photon couples to a closed quark loop [1]. These contributions gives rise to a new flavour structure, depending on the quark flavour composition of the nucleon target through their charges e q . In the normalization of [1] this leads to a relative factor f l 11 = 3 e = (3/n f ) q e q in the non-singlet and f l 11 = (1/n f ) e 2 / e 2 = (1/n f )( q e q ) 2 /( q e 2 q ) in the singlet case, respectively. These f l 11 terms are implicitly contained in g DIS 03 in eq.
2,q . Following a well-established procedure (see, e.g., [29]) the unknown coefficients of a given functional form for c (4) 2,q can be approximately determined using the available Mellin moments for N ≤ 9 of the DIS Wilson coefficients at four loops [30]. In this way, fixing the individual colour coefficients in terms proportional to D 0 , the numerical constraints for the remaining unknown coefficients, i.e., f q have been derived. The terms presented in eq. (3.26) display still a significant uncertainty, whereas the results in large-n c limit are exact. Therefore, we determine the best estimate for the term proportional to D 0 in c 2,q in eq. (3.33) by using the large-n c limit of eq.
where references for all quantities have been already given below eq. (3.33). This leads to gives the following expression which can be evaluated numerically using the results in table 1 and eq. (3.26).
As discussed above, the four-loop terms in eq. (3.26) are not very well constrained and therefore we follow the same approach as in the derivation of eq. (3.34) to determine the present best estimate for B DIS 4 in eq. (3.38). We use the large-n c limit of eq. where the numerical uncertainties at four loops follow from varying the large-n c limit for the terms proportional to n 1 f in f q 4 of eq. (3.23) by 10%, as done also in eq. (3.34). With the knowledge of all resummation coefficients we are now ready to perform the numerical analysis up to N 4 LL accuracy.

Resummation at N 4 LL
In the following we present numerical results for the resummed series up to N 4 LL accuracy. We use the approximate value for B DIS 4 as provided in eq. (3.44) and study the convergence of the resummed series and also the soft-virtual (SV) results at the fourth order in strong coupling and beyond. We set the renormalization and factorization scales to be equal to the momentum transfer squared (Q 2 = −q 2 ). Since we are interested in the region of (very) large x at moderate Q 2 , the charm and bottom quarks cannot be treated as effectively massless. Hence we will focus on the results for n f = 3 light quark flavours, which implies JHEP03(2020)116 f l 11 = 0 for the new flavour structure due to a vanishing combination of the quark charges. For n f = 3 this flavour structure contributes to neutral-current DIS. We will come back to the resulting numerical differences between neutral-and charged-current DIS at the end.
In the left panel of figure 1, we present the resummation exponent up to N 4 LL, normalized to the NLL results for a better visibility of the small higher-order effects. Recall that the resummation accuracy is defined by truncating the resummation exponent in eq. (2.7) to a particular order; i.e., the first term ln N g (1) defines the LL approximation, the first two terms ln N g (1) + g (2) define the NLL accuracy and so on. Using the order-independent value α s = 0.2 for the strong coupling, we observe a good perturbative convergence up to N 4 LL level. At the largest N -value shown, N = 40, the exponent increases from LL to NLL by around 66.7%, from NLL to NNLL by around 6.0%, from NNLL to N 3 LL by around 1.5% and from N 3 LL to N 4 LL by around 0.6%. The resummed series thus shows a good stabilization at this order. The uncertainty in B DIS 4 quoted in eq. (3.44) as a result of changing the respective large-n c terms by ±10% does not affect the cross section much. It amounts at most to 0.1%, which is small compared to the size of the N 4 LL correction.
In the right panel of figure 1 we plot exp(G N ) ⊗ f , the resummation exponential convoluted with a schematic but sufficiently typical shape for a quark distribution given by xf = x 0.5 (1 − x) 3 . For the Mellin inversion we have used the minimal prescription [6] with the standard contour choice as in QCD-Pegasus [64]. Also here we find a good perturbative convergence up to N 4 LL. At x = 0.9 the cross section changes from LL to NLL by 67.1%, from NLL to NNLL by 4.2%, from NNLL to N 3 LL by 0.09% and from N 3 LL to N 4 LL by −0.17%. At x = 0.95 the corresponding changes are much larger, yet the N 4 LL effect stays well within 0.5%. For n f = 4 light quark flavours on the other hand, we see an even faster perturbative convergence for the resummation exponent as well as for the convolution of exponential with the same input shape. For example, at N = 40 we observe a 61.4% increase from LL to NLL, another increase of 1.8% from NLL to NNLL, and changes of −1.1% from NNLL to N 3 LL and of −0.9% from N 3 LL to N 4 LL accuracy. For the case of the convolution of the exponential with the above input shape, the corresponding values read 61%, 0.3%, −1.7% and −0.8%, respectively, at x = 0.9.
The threshold limit in the Mellin-N space is based on the large-N limit, which we have discussed so far in terms of the variable N = N exp(γ E ), including the exponentiation of all ln N terms. One can, in principle, define the large Mellin variable as N only and collect only those terms in the resummed exponent which are enhanced as N → ∞. As far as the threshold limit is concerned, both, N → ∞ and N → ∞ are physically equivalent. However, numerically the two options differ since in the former case all γ E terms associated with N are also exponentiated to all orders. Due to this feature the results for both cases are different starting already at LL accuracy. The resummed exponent shows a faster convergence for the N exponentiation. However, note that now also the function g 0 in eq. (2.1) differs between the two cases. Whereas in the N exponentiation, the γ E terms are collected by definition in N as the large-N variable, for the N -exponentiation they are appear in the finite function g 0 and in the resummed exponent G. As mentioned already, we collect the explicit results for the perturbative expansion of g 0 in the appendix B. A useful comparison for these two approaches is therefore to check the convergence at the cross-section level itself. We provide this comparison in figure 2. For the Nexponentiation satisfactory convergence occurs in the threshold region only beyond NLL accuracy whereas the N -exponentiation shows a systematic behaviour for perturbative convergence for the successive orders of the resummation and good convergence in the threshold region is achieved already at NLL accuracy. At sufficiently high logarithmic order both approaches converge. However they start to differ away from the threshold. At x = 0.7, the LL result in N -exponentiation differs by as large as 15% compared to the N -exponentiation, whereas at the N 4 LL level this difference decreases to 0.1% showing that at higher logarithmic accuracy they tend to converge. We see similar observation in the higher x region as well. In the rest of the article, we apply the standard N exponentiation throughout.
It is instructive to study the quality of the specific large-n c approximation used to derive the approximate value for B DIS 4 in eq. (3.44) at each lower order. To that end, we have compared the large-n c behaviour of results at lower loops with the exact full-colour expressions available up to three loops. Note that in order to have a meaningful estimate of the resummed coefficient function, we take the large-n c approximation only for the process-dependent pieces appearing at that order, whereas all the lower order terms are kept exact. For example, at the third order, we take the large-n c approximation only for the combination −(f q 3 +B q 3 ) appearing in eq. (3.37) whereas all the other terms proportional to coefficients of the QCD beta function are kept exact. The LL result is always exact since JHEP03(2020)116 it only depends on the universal cusp anomalous dimension A 1 while the first processdependent coefficients enter as −(f q 1 + B q 1 ) at NLL order. Note that up to the third order the large-n c limit for the terms proportional to n 1 f and n 0 f in f q i and B q i coincides with the exact large-n c expressions of these coefficients when restoring the overall factor C F . However this is not true anymore at the N 4 LO level.
In figure 3, we compare the large-n c and the exact results for the resummation exponent (left) and the resummed series convoluted with the same quark distribution xf used above (right). In particular we have taken the ratio between the large-n c and exact results at lower orders. We have taken n f = 3 light quark flavours also for this study. Up to NLL the largen c result is exact. The difference in the cross section starts from NNLL and is about 0.4% in the large-N region (N = 40), whereas at N 3 LL it is about 0.2%; the large-n c approximation always overestimates the exact result in the region of interest. These differences are to be compared with the absolute size of the exact corrections. These are much larger with about 6.0% at N = 40, when increasing the logarithmic accuracy from NLL to NNLL, and 1.5% for NNLL to N 3 LL. Similar observations also hold at lower N values.
Using the approximate value for B DIS 4 in eq. (3.44) with the specific large-n c limit, the step from N 3 LL to N 4 LL accuracy amounts to a correction of about 0.6% at N = 40. Based on the lower order studies, we thus expect that the exact result at fourth order for the complete N 4 LL calculation should not differ by more that ±0.2% from our prediction for the resummed exponent at large-n c . In summary, the above procedure of taking the specific large-n c approximation provides robust estimates for the cross sections of interest.

Soft-virtual cross section at the fourth order
The convergence of the threshold expansion in Bjorken-x or Mellin-N space can be illustrated by means of the successive addition of the plus-distributions D k , k = 2n − 1, . . . , 0 in the SV cross section in x-space compared to the successive addition of the logarithms ln i N , i = 2n, . . . , 1 in N -space. 1 The former procedure shows poor convergence, a fact that has been already observed in [6] for two (n = 2) and in [8] for three (n = 3) loops. In figure 4 we compare the convergence of the DIS Wilson coefficient at four loops c (4) 2,q constructed in this manner either in x-or in N -space and convoluted with the above input shape xf . As expected we see a very similar behaviour compared to the known lower orders. The x-space result converges poorly and only stabilizes after the addition of the six highest plus-distributions D k at this order. On the other hand in N -space, the same convergence is achieved much faster after adding only the four highest logarithms ln i N in the series. This confirms that N -space is better suited for studies of the threshold limit also at the fourth order.
The exact expression of D 0 at fourth order still contains the poorly constrained coefficients f q 4, d 2,q in eq. (3.33) convoluted with the input shape xf of figure 1 with the successive addition of the plus-distributions D k starting from the highest term. Right: the same with the successive addition of the N -space logarithms. rather small change considering the errors in these coefficients in eq. (3.26). On the other hand setting the f 4 coefficients in eq. (3.26) to zero, which is within the uncertainty range quoted, yields a value for the D 0 coefficient which differs only by 0.5% from its best prediction in eq. (3.34). These are small effects given that the best estimate for D 0 in eq. (3.34) changes the whole N 4 LO SV cross-section (up to the contribution proportional to δ(1 − x) in eq. (A.5)) by only 0.9% even at x = 0.95. Altogether, this demonstrates a small dependence of the total cross section on those, as of yet, still poorly known coefficients. At the same time, it corroborates the assumptions made in the derivation of the best estimate for D 0 in eq. (3.34).

Tower expansion vs exponentiation
The resummed cross section in Mellin N -space can also be reorganized in a different manner [8,65] by re-expanding the exponential in eq. (2.1) as follows: where the coefficients c kl for DIS are given in table 2 for all known terms. By successive addition of terms in this expansion, one can also predict the resummed cross section up to a certain accuracy. The truncation of this series at any particular order will, of course, recover the fixed-order SV result in N -space. The coefficients c kl of the logarithms in the above JHEP03(2020)116 series vanish factorially at sufficiently higher orders in α s for a fixed logarithmic structure (i.e., fixed by l). The coefficients c kl are determined up to l = 1, 2, 3 by the complete one loop calculation with NLL resummation. The complete two-loop calculation along with the NNLL resummation fixes the next two towers (l = 4, 5), whereas the complete threeloop calculation along with the N 3 LL resummation determines again the next two towers, i.e., up to l = 7. The N 4 LL resummation derived here allows to completely fix the eighth tower, which corresponds to l = 8. Note that the quartic colour factors, i.e., d FA and d (4) FF in eq. (3.3), appear in this tower for the first time. To obtain all terms up to the ninth tower (l = 9), the complete four-loop calculation is needed, which is currently not available. However from the knowledge of higher moments for DIS, we can estimate approximate values for the unknown δ(1 − x) coefficient (see appendix A) at the four loops which in turn can be used to extract the coefficient g DIS 04 . For the standard N exponentiation, then one can have at the fourth order, are added in quadrature to get the final uncertainty in column c k9 . Although we note that this uncertainty is mostly dominated by the uncertainty in g DIS 04 estimate. In figure 5 on the left we illustrate the expansion in decreasing powers of ln N for the five-loop coefficient function c (5) 2,q in the same manner as for the four-loop coefficient function c (4) 2,q in the right part of figure 4 above. We confirm the general pattern, that the first l logarithms provide a good estimate up to the l-th order in α s , i.e., l = 5 for the case at hand. In the right panel of figure 5 we compare the resummed result (exponentiated) for the Wilson coefficient with the tower expansion given in table 2. Although they both converge in the region of x 0.5, they start to diverge at large x where exponentiation shows a better stability compared to the tower expansion. In fact the expansion with a fixed number of towers severely underestimates the effect of the coefficient functions of much higher orders, which become more important at the region of very large-x around x 0.9.
So far we have considered f l 11 = 0 for the flavour structure, which holds for n f = 3, due to the vanishing combination of the quark charges. Starting from three loops, the neutral-current DIS coefficient functions differ due to f l 11 = 0, which captures the effect of incoming and outgoing photons coupling to quark lines of different flavour, see e.g., [1].  and corresponding errors and the ninth column in addition requires the estimate for g DIS 04 in eq. (4.2).

JHEP03(2020)116
2,q by the large-N terms specified in table 2, illustrated by the convolution with the input shape xf of figure 1. Right: the corresponding results for the effect of higher terms beyond α 4 s as obtained from the tower expansion up to nine towers and from the exponentiation up to N 4 LL accuracy on the Wilson coefficient.

JHEP03(2020)116
The effect of this contribution is numerically insignificant as already pointed out in [8]. To estimate its numerical impact on the SV cross section as well as in the tower expansion, we use n f = 4 and keep the exact value f l 11 = 1/10. Note that a non-zero f l 11 will change the tower expansion starting from seventh tower, i.e., the coefficients c k7 , and from the lowest order for this towers (k = 3). In the seventh tower, the contribution proportional to f l 11 changes the finite piece of third order SV result, whereas at the fourth order in addition it also changes the eight tower. Keeping exact value for f l 11 we get a difference of 0.1% in the coefficient of the constant term at third order and a change of 0.09% in the last logarithmic term at fourth order.
The additional new term in the charged-current DIS coefficient function C 3,q proportional to the group invariant (d abc d abc )/n c in eq. (3.1) which corresponds to incoming and outgoing W ± -bosons coupling to a closed internal quark loop [2]. This term is power suppressed for x → 1, hence does not contribute to the threshold approximation considered here.

Summary
We have applied recent progress in the calculation of the cusp anomalous dimension, the quark (non-singlet) splitting functions and the quark form factor to obtain predictions for the DIS cross sections with charged-and neutral-current interactions at four loops in perturbation theory. We have performed the threshold resummation in the large-n c limit up to N 4 LL collecting the terms α 3 s (α s ln N ) k to all orders in α s in the threshold resummation exponent G N in Mellin-N space. To that end, the relevant quantity B DIS controlling the resummation of the quark jet function has been extracted at four-loop order.
We have also derived the full colour dependence of the DIS Wilson coefficients at four loops down to the term proportional to the plus-distribution [1/(1−x)] + using a conjecture about the colour structure of the anomalous dimension of the eikonal form factor, namely that it exhibits the same generalized maximal non-Abelian property regarding quadratic and quartic Casimir coefficients as the cusp anomalous dimension. As a by-product, also the single pole of QCD form factor in the parameter of dimensional regularization ε has been determined numerically to good accuracy.
Using the analytical expressions for the resummed DIS Wilson coefficient at the N 4 LL level, we have shown that the lower order corrections are sizable and that the threshold logarithms at N 4 LL resummed order amount to corrections well within 1% and improve the perturbative convergence. We have compared the effect of successive plus-distributions in x-space as well as the logarithms in N -space for the cross section up to the fourth order and we confirm known observations that the convergence of the logarithmic series in N -space is much smoother compared to the one in x-space. We have also reorganized the resummed series as a tower of logarithms up to the ninth tower. While the tower expansion is useful to study the effect of higher logarithms, it lacks contributions at very large x compared to the exponentiated result. The latter exhibits an even better stability even in the very large-x region. Near threshold the difference between DIS Wilson coefficients for charged-and JHEP03(2020)116 neutral-current exchange, which manifests itself in the flavour structure f l 11 at three loops and beyond, has only a minor numerical effect (for n f = 4) on the respective cross sections.
The current accuracy can be significantly improved with the knowledge of the complete QCD form factor at four loops or an exact computation of the corresponding eikonal anomalous dimension at this order. In particular analytic results for all n f dependent terms in those quantities would greatly reduce the present residual numerical uncertainties. At the accuracy reached with N 4 LL resummation also power corrections proportional to 1/N in Mellin-N space or logarithmic enhancement ln k (1 − x) in x-space become numerically relevant, see [66,67]. Finally, the effects from the massive quarks are important in the threshold region (see for example [68,69]) and a proper treatment of heavy quark mass effects becomes essential for phenomenological applications.
Note added in proof. During review of the present paper, the complete analytic fourloop anomalous dimensions in massless QCD from form factors have been obtained in [72]. This provides information on the missing four-loop contributions f q , which we can determine as 2) agrees very well with the best estimate in eq. (3.34). This also corroborates the arguments based on the large-n c limit.

JHEP03(2020)116
With the full colour dependence in eq. (5.2) at hand, the resummation coefficient B DIS

A DIS Wilson coefficients
For convenience of the reader we have collected all available information on the soft-virtual expressions for the DIS coefficient functions c (n) 2,q through four loops (n ≤ 4). The leading order is c (0) 2,q = δ(1−x) and the one-loop result c (1) 2,q in eq. (A.1) has been obtained in [73]. The two-and three-loop expressions c