The N$^3$LO Scheme-invariant QCD Evolution of the Non-singlet Structure Functions \boldmath $F^{\rm NS}_2(x,Q^2)$ and $g_1^{\rm NS}(x,Q^2)$

We present the scheme-invariant unpolarized and polarized flavor non-singlet evolution equation to N$^3$LO for the structure functions $F_2(x,Q^2)$ and $g_1(x,Q^2)$ including the charm- and bottom quark effects in the asymptotic representation. The corresponding evolution is based on the experimental measurement of the non-singlet structure functions at a starting scale $Q_0^2$. In this way the evolution does only depend on the strong coupling constant $\alpha_s(M_Z)$ or the QCD scale $\Lambda_{\rm QCD}$ and the charm and bottom quark masses $m_c$ and $m_b$ and provides one of the cleanest ways to measure the strong coupling constant in future high luminosity deep-inelastic scattering experiments. The yet unknown parts of the 4-loop anomalous dimensions introduce only a marginal error in this analysis.


Introduction
The measurement of the strong coupling constant α s (M Z ) from precision data on deep-inelastic scattering is one of the cleanest ways to obtain this fundamental parameter of the Standard Model. The present world data allow measurements of δα s (M Z )/α s (M Z ) at the level of ∼ O(1%), [1][2][3]. In the standard analyses one is fitting the scaling violations of deep-inelastic structure functions and obtains α s (M Z ) together with the parameters of the parton distribution functions, observing the correlation of all parameters, cf. [4]. Compared to this, the measurement of α s (M Z ) from R e + e − , cf. [5], does not need to fit a larger amount of further parameters. The different precision determinations of α s (M Z ) using different processes do not yet agree and further precision measurements are therefore needed.
For deep-inelastic scattering a direct determination of α s (M Z ) is also possible in the case the starting distributions of the evolution are measured experimentally. In the flavor non-singlet case the corresponding non-singlet structure functions at a scale Q 2 0 provide the input and the data at Q 2 > Q 2 0 are fitted to a physical evolution function, which only depends on α s (M Z ), provided that the input for the values of the heavy quark masses m c and m b is known at high precision, as also required in the measurement of α s (M Z ) from R e + e − . Here Q 2 denotes the virtuality of the exchanged gauge boson. Therefore scheme-invariant evolution equations allow for a direct determination of the strong coupling constant using measured one dimensional input distributions which are here functions of the Bjorken variable x. An analysis of this kind has been performed in Ref. [6] in the unpolarized case, effectively up to N 3 LO, however, only considering massless quarks with some further assumptions, see also Ref. [7]. The sensitivity of the structure function on α s (M Z ) or Λ QCD is traditionally illustrated by the slope ∂ ln(F i (x, Q 2 ))/∂ ln(Q 2 ) determined experimentally [8][9][10].
In the present paper we complete the formalism by the single-and double-mass heavy flavor corrections for the Wilson coefficients up to three-loop order and provide numerical results for these effects, using the results given in Refs. [11][12][13][14][15][16][17] 1 , as well as the massless contributions. Scheme-invariant evolution is usually performed in Mellin N space for technical reasons, because it is simpler in analytic form compared to the corresponding z-space representation. In the polarized case the reference to factorization scheme-invariant quantities will also solve the associated γ 5 problem, which usually arises performing the loop calculation in D = 4 + ε spacetime dimensions. The strong coupling constant is used in the MS scheme. A main goal of the present study is to outline the different evolution effects for high precision measurements.
The measurement of scaling violations concerns that of massless parton evolution and is there implied by the ultraviolet singularity of the respective twist-2 local operators. Referring to the scale evolution using observables necessarily requests to map to the massless partons, which implies in N space the algebraic use of the Wilson coefficients, already at the hard scale Q 2 , as an obstacle, which are different for the various hard processes. For the latter ones it has to be guaranteed that only twist-2 terms contribute, cf. [7,19,20], by putting the cuts W 2 > 15 GeV 2 , Q 2 > 10 GeV 2 , where W is the hadronic mass. The natural scale of evolution in deep-inelastic scattering is Q 2 .
The evolution equation is solved analytically and one final numerical contour integral around the singularities of the corresponding problem delivers the evolved flavor non-singlet structure function in x space. We first consider the flavor decomposition to provide realistic input distributions and derive the corresponding evolution operators for the structure function. Here we also discuss the remainder systematics of the yet not completely known four-loop non-singlet splitting functions P NS,(3),± . Finally we are providing numerical results in the unpolarized and polarized case. Precision analyses of this kind can be performed in the physics programme at the future collider experiments at the EIC [21] and LHeC [22].

Flavor Decomposition
Since the evolution for the three different flavor non-singlet distributions and the singlet evolution are different, the experimental input at the starting scale Q 2 0 has to be projected correspondingly by combining different deep-inelastic structure functions. One may refer to the following wellknown decomposition, see also [23]. Let with q i the quark distributions and one has The nucleon structure functions for pure photon exchange in the case of three light flavors (u, d, s) are then given at leading order (LO) by A synonymous decomposition applies to the structure functions xg p 1 and xg d 1 . To project onto the singlet distribution Σ directly one usually needs charged current structure functions, as measured in neutrino experiments [24] or at facilities like the planned LHeC project [22] or the EIC [21] in the unpolarized and polarized case, at high luminosity, with [25] above the charm threshold [26]. Here the index ± denotes the exchange of a W + or a W − boson, respectively. Otherwise additional information on sea-quarks is necessary. The flavor non-singlet combinations we are considering are given by where the Wilson coefficients are C q and ∆C q contain light and heavy flavor contributions, and ⊗ denotes the Mellin convolution The Mellin transform Before forming the structure function difference in (10,11) one has to unfold the nuclear corrections of the deuteron structure functions. The lowest Mellin moment of (11) is given by the polarized Bjorken sum rule [27].

The Non-singlet Evolution
The evolution operator of scheme-invariant flavor non-singlet evolution, E NS , is obtained as follows. 2 The evolution equation for the non-singlet structure functions can be written as Its solution is given by The Wilson coefficient is given by Here c k denote the expansion coefficients of the massless Wilson coefficients and h k of the massive Wilson coefficient, with and m c,b are the on-shell charm and bottom quark masses.
In the non-singlet case the heavy flavor corrections contribute form O(a 2 s ) onward. One has 2 Here and below we will work in Mellin N space.
whereĥ i denote the single mass andĥ 3 the double mass contributions. One may rewrite the differential operator The evolution equation for the non-singlet quark density is given by where β k are expansion coefficients of the QCD-β function and P k are the splitting functions.
The anomalous dimensions are related to the splitting functions 3 by The solution of Eq. (15) to N 3 LO reads and a = a s (Q 2 ), a 0 = a s (Q 2 0 ) and P i (N ) ≡ P i denotes the Mellin transform of P i (z). The heavy quark contributions to the Wilson coefficients are given by [15,16,30,31] 3 Our normalizations are such that a factor of two has to be applied to those given in [28,29].
1,Q + The two-mass three-loop contributions [32] read The two-mass term is the same in the unpolarized and polarized case. In the r.h.s. of Eqs. (27)(28)(29) we definef In the MS scheme the iterative solution for a s (Q 2 ) is [33] with L = ln(Q 2 /Λ 2 QCD ). Here the integration constant for solving (23) is chosen by (β 1 /β 2 0 ) ln(β 0 ) [34]. The expansion coefficients of the β-function to N 3 LO were calculated in [35]. The flavor matching conditions were given in [33]. The expansion coefficients of the renormalized mass were given in [36,37]. The constant and O(ε) parts of the massive unrenormalized operator matrix elements at O(a k s ) are denoted by a ij , respectively, cf. [11,12,38,39]. The three-loop massless Wilson coefficients are expressed by effective representations. Otherwise we use the analytic Mellin space representations. After algebraic reduction [40], they depend on 32 harmonic sums [41,42] up to weight w=6 only, which is particular to the flavor non-singlet case to three-loop orders. Other massive Wilson coefficients have a more involved structure [43][44][45][46]. It is useful to represent harmonic sums with an alternating index by their Mellin transform, to eliminate spurious (−1) N terms [42] before performing the analytic continuation from even or odd integers to the complex plane. The singularities of the problem are located at the non-positive integers, around which the contour integral is performed, cf. [47][48][49][50]. For the analytic continuation we follow Refs. [51,52].
In the case of the structure function F 2 the non-singlet anomalous dimensions are P NS,+ and the expansion coefficients of the Wilson coefficient up to c 3 were calculated in [53]. In the case of the structure function g 1 the anomalous dimensions are P NS,− , cf. [28], and the Wilson coefficient has been given in [54]. 4 Below we will also use the combination The four-loop non-singlet anomalous dimensions are not yet completely known as a function of N . However, a series of moments has been calculated in [79][80][81][82][83][84]. We follow the earlier investigation in Ref. [6] and compare these moments with the Padé-approximation with the exact moments in [79][80][81][82][83][84] in Table 1. Furthermore, the leading N F terms for the even moments have been predicted in [85]. From the 2nd moment, which agrees within 21%, cf. also [6], the accuracy improves to 2.2 % for the known even moments at N = 16 and to 2.6% for the odd moments at N = 15. For the first moment for γ −,NS the Padé-approximation delivers even the correct result, after using the l'Hospital rule.  have been given in [86] after correcting some misprints there in Ref. [87], see also [88]. However, yet unknown sub-leading terms, having been studied at the known lower orders, are numerically dominant over the leading small x corrections and several of these sub-leading terms are needed, cf. [49,87], to be calculated to obtain reliable quantitative prediction in this region of x. Table 1 shows that (33) provides an excellent model. The earlier assumption of an error of 100% in [6] on this relation has been very conservative and still led to an error in Λ QCD of 2 MeV only, which could now be improved in principle. Yet the experimental accuracy at this level cannot be reached at present, since the current experimental error amounts to δΛ QCD = 26 MeV [6].

Numerical Results
The measured input distributions contain correlated errors, which are parameterized by the respective fits. cf. [6,19,89]. Their evolution, including the corresponding correlation matrix, has to be performed to provide the error prediction of the structure function at every fixed value of Λ QCD , m c and m b respectively. Here one may use the corresponding formulae for the correlated error propagation given in Ref. [6] and extensions of them, which are straightforward. Within future data analyses one will probably import the values of the heavy quark masses from the world data analyses. For the charm quark mass it has already been shown that its value obtained from deep-inelastic scattering data fully agrees with other precision measurements [90]. In our illustrations we will use the values m c = 1.59 GeV [90] and m b = 4.78 GeV [91].
For the input distribution in the unpolarized case we refer to the one of Ref. [6] F NS 2 (x, at Q 2 0 = 4 GeV 2 . In the polarized case we use the fit to xg NS 1 (x, Q 2 0 ) of the structure function of Ref. [89] at Q 2 0 = 10 GeV 2 , In the numeric illustration our reference starting scale will be chosen to be Q 2 0 = 10 GeV 2 both in the unpolarized and polarized case. In Figures 1 we show     In Figures 4 we show the ratio of the results obtained at leading order (LO), next-to-leading order (NLO), and next-to-next-to leading order (NNLO) to the N 3 LO results at Q 2 = 100 GeV 2 . This illustrates the convergence of the perturbative series of the corrections and the necessity to include N 3 LO corrections at an accuracy of the data in the 1% region. and xg NS 1 . This rescaled correction is in the sub-percent range. Moreover, the impact on Λ QCD comes from the slope in Q 2 which is seen to be rather small.

Conclusions
We have calculated the evolution of flavor non-singlet structure functions in the unpolarized and polarized case, which become available in high luminosity experiments performed both on proton and deuteron data with comparable statistics. The evolution has been performed in Mellin N space using our code QCDEVO. Referring to the structure functions directly has the advantage that the input distributions are measured. The evolution does only depend on the QCD scale Λ QCD and the heavy quark masses m c and m b , the latter of which have been measured very precisely already [91]. Therefore, this method allows for a one parameter fit of Λ QCD only, and provides a method which is systematically very stable. The input, after corrections for the deuteron wave function, is given by the difference of two structure functions which both can be measured at very high statistics at future facilities [21,22]. The heavy flavor corrections amount to a contribution of ∼ 1% in the region x < 0.4. Let us finally mention that one may also consider the scheme-invariant singlet evolution [62,65,[92][93][94][95][96][97][98][99][100] in addition to the one in the flavor non-singlet case, provided the corresponding flavor decomposition can be carried out. This will usually require more than just the deepinelastic data, because of the sea-quark combinations. A possible way would consist in referring to the Drell-Yan cross section here. Already this part complicates the analysis. Furthermore, one needs to measure the slope ∂F i /∂ ln(Q 2 ) at the input scale Q 2 0 together with F i (Q 2 0 ) in an uncorrelated way. A second possibility would consist in using the pair F S 2 (Q 2 0 ), F S L (Q 2 0 ). However, the measurement of the structure function F S L is much more difficult than that of F S 2 and in the past not enough statistics has been collected in this case, cf. [101][102][103], see, however, [104]. To perform the necessary flavor decomposition, one needs to perform all these measurements both at proton and deuteron targets. By scheme-invariant evolution the fitting problem of input distributions does not exist and no further optimization by neural network techniques is necessary since the input distribution is fully determined by its measurement.