Nonsinglet polarized nucleon structure function in infrared-safe QCD

The polarized nucleon structure function in the nonsinglet case is investigated here by a new insight rather than conventional perturbative QCD (pQCD). For this purpose we note that the solution of the evolution equations in moment space involves noninteger powers of the coupling constant. Therefore it is possible to employ a new approach which is called fractional analytical perturbation theory Consequently, it is possible to remove the Landau singularities of the renormalized coupling, i.e., at the scales $Q \sim \Lambda$, using this approach. This provides an opportunity to continue the desired calculations toward small values of energy scales even less than the $\Lambda$ scale. To modify the analytical perturbation theory, a newer approach is introduced, called 2$\delta$anQCD, in which the spectral function of the holomorphic coupling is parameterized in the low-energy region by two delta functions. This model gives us more reliable results for the considered QCD observables, even in the deep infrared region. We calculate the nonsinglet part of the polarized nucleon structure function, using the 2$\delta$anQCD model, and compare it with the result from the underlying pQCD where both are in a new defined scheme, called the Lambert scheme. For this purpose we employ the anQCD package in the \textit{Mathematica} environment to establish the analytic (holomorphic) coupling constant. The results at various energy scales are also compared with the available experimental data, and it turns out that there is a good consistency between them. The results show that the obtained nucleon structure function at small energy scales has smoother behavior when using the 2$\delta$anQCD model than the underlying pQCD. In fact the coupling constant in analytic QCD behaves moderately and it makes the result approach the available data in a better way.


I. INTRODUCTION
QCD analysis of deep-inelastic scattering (DIS) data provides us new insight into the hadron physics. The reliability of theoretical results can be increased by testing the analysis, for instance, of hard lepton scattering off the hadrons. A reliable description of DIS at large momentum transfer Q 2 ≫ 1 GeV 2 can be obtained, considering the twist expansion and factorization theorem. Following that one can use the DIS processes to better understand the nucleon structure function. In this regard the polarized and unpolarized structure functions can be first modeled, using different parametrization methods. Then they can be evolved to high-energy scales, employing the QCD evolution equations. During the evolution process the numerical values of structure functions are affected by the energy scale, at which the renormalized coupling constant as the perturbation expansion parameter is evaluated. On the other hand the analysis of the structure functions should lead to consistent results with the related experimental data. To achieve this aim it is needed to model the parton densities, as ingredients of nucleon structure, at an initial energy scale. A further step is the evolution of these densities to high-energy scales. Therefore, the role of the renormalized coupling constant α s (Q 2 ) is important and should be considered in detail.
At low-Q 2 energy scales, below 1 or 2 GeV 2 , the qualification of DIS as a hard process is not appropriate. Its analysis is faced with two essential difficulties in the low-Q 2 region. One of them is related to the higher-twist corrections which are dominantly affecting the contribution of the leading twist at low values of the energy scale. Therefore, the twist corrections at higher orders in α s but at a low-energy scale causes the first difficulty. In addition, since α s (Q 2 ) as the QCD running coupling grows rapidly at the energy scales near the Landau singularities, the result of pQCD would not be reliable and this represents the second difficulty.
A new approach that is called the analytic perturbation theory (APT) can be applied to resolve these problems. This approach has been developed by Shirkov, Solovtsov et al. [1][2][3][4][5][6]. There, the running QCD coupling of conventional pQCD is transformed into an analytic (holomorphic) function of Q 2 which is accordingly called the APT coupling constant. In this regard, one needs a dispersion relation, involving a spectral or discontinuity function which is the imaginary part of the coupling along the timelike semiaxis (Q 2 < 0) in the complex Q 2 plane. In the framework of APT, the images (analogs) which correspond to integer powers of the original perturbative QCD (pQCD) coupling constant can also be constructed. Later on, an extension to the analogs of noninteger powers was performed, which is called fractional APT (FAPT). Some newer constructions of the coupling constant in analytic QCD (anQCD) are based on the modifications of the spectral function at low-energy scales. In the case of this modification, the spectral function is constructed in a specific manner, by parameterizing its behavior in the low-energy region by two positive delta functions [7,8]. The presented results of this paper are based on employing this model which is called 2δanQCD.
We remark that anQCD is not the only approach to obtain a proper behavior of the coupling constant at low energy, i.e., at infrared (IR) scales. To achieve an IRsafe QCD description one can also refer specifically to the AdS/CFT coupling of Brodsky, de Teramond, and Deur [9] or the dispersive approach of Dokshitzer et al. [10,11]. A new review of the QCD running coupling in various approaches can be found in Ref. [12].
The remainder of this paper consists of the following sections. In Sec. II we provide a brief discussion of the theoretical formalism for different aspects of anQCD, including the 2δanQCD model. In Sec. III we deal with the theoretical formalism of pQCD in DIS processes where the Jacobi polynomials and the rational expansion for the moment of polarized nucleon structure function are also presented. In Sec. IV we employ the 2δanQCD model as the approach to evaluate the perturbation expansions and to extract the moments of the polarized structure functions at the next-to-leading-orde (NLO) approximation. In Sec. V we also utilize this model to evaluate the Bjorken sum rule for polarized nucleon structure functions. Finally we give our summary and conclusions in Sec. VI.

II. AN OVERVIEW OF ANALYTIC QCD APPROACHES
The anQCD approaches are intended to address the problems which were mentioned in the introduction, especially the second problem, i.e., the presence of the Landau singularities in pQCD couplings α s (Q 2 ). FAPT [1][2][3][4][5][6][13][14][15][16][17] is one such approach. Other approaches [9][10][11] could be applied, but we do not follow this line here. The calculations in the present work will be performed using the 2δanQCD approach [7,8], which has significant modifications with respect to FAPT.
For this purpose, we first give a brief description of (F)APT, a description essentially using Ref. [18]. After that we give a summarized description of the 2δanQCD model.
In the FAPT approach the running pQCD coupling a s (Q 2 ) ≡ α s (Q 2 )/4π of pQCD is transformed into a holomorphic function of Q 2 in which first the conversion a(Q 2 ) → A 1 (Q 2 ) is performed for the original pQCD coupling. The new coupling A 1 (Q 2 ) is an analytic (holomorphic) function of Q 2 in the complex Q 2 plane with the only exception of the timelike semiaxis (Q 2 ≤ 0), and it is called the APT coupling constant.
As it will be seen later on, the APT coupling is given by a dispersion relation. In this relation there is the pQCD spectral density, denoted by ρ (pt) (σ) ≡ Im α s (Q 2 = −σ− iǫ)/π which in APT is not changed in the complex Q 2 plane for the whole negative axis, given by the condition σ ≥ 0. For the unphysical cut, i.e. 0 < Q 2 < Λ 2 , this spectral function is set equal to zero [18]. Following the same dispersion relation, based on the APT framework, the images (analogs) A n (Q 2 ) of integer powers a n (Q 2 ) can also be constructed.
In contrast with the original a n (Q 2 ) the couplings A n (Q 2 ) are changing slowly at low Q 2 values. On the other hand, at large Q 2 values these two couplings are close to each other . The FAPT approach, which was summarized here, is an extension to noninteger powers ν in which a ν maps to A ν [13][14][15][16][17]. This case has many applications to DIS processes. One can refer to Ref. [19] which contains reasonable results in analyzing with the FAPT approach the DIS data for the hadron characteristics. Alongside the FAPT framework, there are various analytic QCD models which can be found in Refs. [18,[20][21][22][23][24][25][26][27][28]. One of them is called the 2δanQCD model [7,8], which will be discussed later on in this section.
(2) and then insert in Eq. (1) the result for ρ 1 . The situation becomes more complicated when the two-loop approximation is considered. In this case, the solution of the QCD-β function will yield the following result for a s (2) in terms of the Lambert-W function [18]: .
In this equation z W (L) is given by where in which N f is denoting the number of active quark flavors, and Considering Eq.(3), the required expressions for R (2) and ϕ (2) would be, respectively, Here W −1 (z) represents the appropriate branch of the Lambert-W function. Extension of the calculations up to the four-loop level can be found in Ref. [29]. Because of the linearity of the transforms in Eq. (1) there is a one-to-one correspondence between the pQCD and FAPT expansion [5]. This linearity can be exposed, considering single-scale quantity D(Q 2 , µ 2 R ) at the renormalization scale µ 2 R = Q 2 within the minimal subtraction scheme. In this regard the expansions for D and its corresponding image, i.e., D, have the following expansions [18]: We note that d i coefficients are identical in both expansions, and the renormalization scale µ 2 R = Q 2 is used.
After the summarized description of the FAPT approach, we now present briefly the other anQCD model, which at low-energy scales (where Landau singularities appear in a) has a modified spectral function. It is therefore expected that, by applying this model, more reliable results are obtained, arising out of the anQCD approach. This model, as we mentioned earlier, is called 2δanQCD. Our calculations in this paper are entirely based on utilizing this model.
In 2δanQCD , in the low-σ regime, i.e. 0 < σ < M 2 0 , the behavior of the spectral function ρ(σ) ≡ ImA 1 (Q 2 = −σ − iǫ) is controlled by two positive delta functions. The coupling A 1 (Q 2 ) is presented by a dispersion integral which, in addition to the two delta functions, contains as well the ρ (pt) (σ) function. The unknown parameters of the delta functions, and the M 0 scale, are determined by two requirements [7,8,30]. The first one is to effectively match the model with the pQCD at large |Q 2 | > Λ 2 . Second, the model should reproduce the reported experimental rate for the τ lepton semihadronic nonstrange V + A decay ratio, i.e. r (D=0) τ ≈ 0.203 . The physical spacelike QCD quantities D(Q 2 ), to be evaluated in anQCD with A 1 , are in practice, in general, represented in pQCD as a truncated perturbations series of Eq. (9), and using the renormalization scale µ 2 = κQ 2 (0 < κ 1) We point out that in anQCD the simple replacement a(κQ 2 ) ν0+m → A 1 (κQ 2 ) ν0+m in Eq. (11) is not applicable, since it makes the perturbation series diverge rapidly by increasing the power index N (see [31]), and the nonperturbative contributions generated in this way become erratic. Therefore, the formalism introduced in Refs. [25,26] for the case of integer ν 0 and the formalism in Ref. [30] for the case of general real ν 0 , are needed to overcome the difficulty, in any anQCD. Consequently, the following replacement takes place in Eq. (11) in anQCD: The holomorphic power analog A ν0+m (Q 2 ) can be obtained from A 1 (Q 2 ) by the following construction. First one needs the logarithmic derivatives of A 1 (Q 2 ) which are given by [30] An+1 where β 0 = (11 − 2N f /3)/4, and it is obvious that Considering anQCD coupling A 1 , Eq. (13) can be written in terms of the discontinuity function ρ(σ) ≡ ImA 1 (−σ − iǫ) where the Cauchy theorem is used along the cut line. Therefore, Eq. (13) can be written as follows [25,26,30]: Here Γ(n) is denoting the usual Gamma function and Li −n is defined below in the general case by Eq. (17). The generalization to the noninteger n → ν is given by Instead of ρ the above quantity can be expressed in terms Here in which m = 0, 1, 2, . . . and 0 ≤ δ < 1. Equation (16) originates from Eq. (15), using the following expression for the Li −ν (z) function [32]: Combining various generalized logarithmic derivatives, the analytic analogs A ν of powers a ν are obtained [30] A where the coefficients k m (ν) were obtained in Ref. [30] for general ν. It can be shown that this approach, in the specific case of APT (ρ(σ) = ρ (pt) (σ); σ ≥ 0), gives the same result as the expression A (1). However, the approach Eqs. (15)-(18) is applicable in any anQCD. We recall that in the 2δanQCD model at high momenta, σ ≥ M 2 0 , the discontinuity function ρ(σ) ≡ In the low-momentum region it is parametrized by two delta functions, such that [7,8] where Θ is the step (Heaviside) function. Therefore the 2δanQCD coupling is represented by In Eq. (20) the parameters F 2 j and M j (j = 1, 2) are determined by the requirement that the deviation from the conventional pQCD result at high . These are, in fact, four requirements. The pQCD-onset scale M 0 is fixed by the condition that the model should reproduce the measured semihadronic decay ratio r τ ≈ 0.203 for the strangeless and massless V + A tau lepton [7,8,33].
In the 2δanQCD model the coupling a of the underlying pQCD is chosen for convenience in terms of the Lambert-W function (which Mathematica can evaluate without problems) such that As defined before, in this equation where φ is the required argument in which one can follow the calculations in the complex Q 2 plane. The two universal beta coefficients The upper and lower sign refer to φ ≥ 0 and φ < 0, respectively, and the z variable is The coupling (21) is in a renormalization scheme called the Lambert scheme. In these schemes, (n=3,4,· · · ), i.e., all the higher β n scheme parameters (n=3,4,· · · ) are fixed by c 2 [8]. The acceptable values of the scheme parameter c 2 can be in the interval −5.6 < c 2 < −2, with c 2 = −4.9 as the preferred (and chosen) value [8]. The Lambert scale Λ L = 0.255 GeV is determined in such a way that it corresponds in the MS scheme at µ 2 = k, Q 2 = M 2 Z (with N f = 5 and k = 1, see the text below Eq.(51)) to the value α s (M 2 Z ; MS) = 0.1184. The scaling in the pQCD coupling (21) corresponds to the MS scaling, i.e., In Fig. 1 the 2δanQCD running strong coupling is plotted as a function of Q 2 , and compared with the coupling of the underlying pQCD. As can be seen, the 2δanQCD coupling A 1 has a finite value when Q 2 → 0, while the coupling a from the underlying pQCD increases rapidly and goes to infinity near the Landau branching point Keeping in mind Eq. (18), any QCD observable can be calculated within 2δanQCD (or any anQCD). In the following sections we employ this model to analyze the polarized nucleon structure function and the related Bjorken sum rule.

III. NLO EVOLUTION OF POLARIZED PARTON DENSITIES, RATIONAL APPROXIMATION AND JACOBI POLYNOMIALS
We evaluate here the polarized nucleon structure function. The precision depends on the method to calculate the structure function. It is therefore necessary to characterize the employed method. One of the reliable methods is the use of the Jacobi polynomials which leads to sufficient precision in the evolution of the function with Q 2 . Alternatively, using the inverse Mellin transformation, the result in Bjorken x space can be directly obtained from the evolved moments of the structure functions in n space. Before employing the Jacobi method, we present some technical details about the evolution of the polarized parton densities in the moment n space. In this context, an approximate method (rational approximation) is introduced which yields required expansions for the considered moments. For this purpose, we need ∆P (0) N S (n) which is the nth moment of the nonsinglet splitting function at the leading order [34]: In this equation As we specify later on, our calculations throughout the whole manuscript are based on using the 2δanQCD model alongside the underlying pQCD as represented by the (pQCD) coupling Eq. (21). They are both in the specific Lambert (c 2 = −4.9) scheme, but in the rest of this section we use terms and phrases which are valid in general schemes. In the next section, which contains the main results for the nucleon structure function g 1 , the required expressions are used in the Lambert scheme, but here a general discussion is followed. In this regard we resort to a solution of the renormalization group equation for the pQCD coupling in NLO approximation that is given in Ref. [35]. Then the solution of the NLO evolution equation for the moments of the structure function leads to [36]: In the above equation the Wilson coefficient ∆C In Eq. (24), ∆M NS (n, Q 2 0 ) can be obtained from the moment of the polarized parton densities, ∆f NS (n, Q 2 0 ), as follows: ∆f NS (n, Q 2 0 ), the moment of the parton densities at the initial scale Q 2 0 , will be discussed later on in this section. The functions ∆d NS (n) and ∆p(n) in Eq. (24) are the moments of the splitting functions and are given by, respectively, by [36] ∆d NS (n) = −2∆P Here, the explicit expression for ∆P N S± is required which is presented in the Appendix.
One of the methods to obtain a solution for the moments of the structure function in Eq. (24) is to use the rational approximation, which is indispensable when we employ the analytic perturbation theory. On this basis, the following factor in Eq. (24): which is denoted by ∆m(n, Q 2 ) has the following expansions up to O(a) and O(a 2 ) accuracy, respectively [36]: The polarized splitting function ∆p(n) in Eq. (30) is obtained by Eq. (28).
To construct the moments of structure functions, one needs the polarized parton densities at initial Q 2 0 as the input densities. We take them from the data-based KATAO PDFs at Q 2 0 = 4 GeV 2 such that [37] One of the advantages of the value Q 2 0 = 4 GeV 2 here is that the 2δanQCD coupling at this scale practically coincides with the underlying pQCD coupling (21). The numerical values for η qv and α qv , β qv and γ qv parameters are listed in Table I. The normalization constants N qv , are chosen such that η qv is considered as the first mo- (32) is the Euler beta function. Details of computations to obtain numerical values for the first moments of ∆u v and ∆d v are presented in the following subsection. At the NLO approximation, the results of calculations are scheme independent (c 2 independent), and we can use the Lambert scheme.

A. First moments of ∆uv and ∆dv
The parameters η uv and η dv are the first moments of the polarized valence quark densities, ∆u v and ∆d v , respectively. These moments can be related to F and D quantities [38] as measured in neutron and hyperon βdecays In these equations a 3 and a 8 are the nonsinglet combinations of the first moments which are constructed from the polarized parton densities such that ∆q 8 = (∆u + ∆u) + (∆d + ∆d) − 2(∆s + ∆s) . (36) Doing a reanalysis for F and D with updated β-decay constants, one obtains F = 0.464 ± 0.008 and D = 0.806 ± 0.008 [38].
Based on Eq.(33) and considering the experimental values for F and D, the following numerical values are obtained for the first moments of the polarized valence densities: Full results for the moments of the nucleon polarized structure functions at the NLO approximation are now accessible. It is therefore necessary to obtain the structure functions in Bjorken-x space. This can be done using the Jacobi polynomial method, which we summarize in the next subsection.

B. The Jacobi polynomial method
Based on the Jacobi polynomial method,the polarized structure function, xg , can be expanded as follows [37,[39][40][41][42]: In the above equation, Θ α,β n (x) are Jacobi polynomials of order n, and N max is the maximum order of expansion. Jacobi polynomials provide a method to separate the main part of the x dependence of the structure function into the weight function x β (1 − x) α , while the Q 2 dependence is contained in the Jacobi moments a n (Q 2 ) [43].
Jacobi polynomials fulfill an orthogonality relation Using this, Eq. (39) can be inverted to yield the Jacobi moments a n (Q 2 ): To derive the last line in the above equation, it is needed to substitute Eq. (39) for xg N S 1 (x, Q 2 ) into the first line of Eq. (41) and to use the Mellin transform The polarized structure function xg N S 1 (x, Q 2 ) can now be related to Mellin moments as follows [37]: It is required to choose the set {N max , α, β} such that an optimal convergence of the series is achieved. This convergence should contain the whole kinematic region and cover the related experimental data. An improvement is achieved for α = 3.0, and β = 0.5 while N max varies between 7 and 9 [42].
As an alternative to the above method, the structure function xg N S 1 (x, Q 2 ) may be obtained by using the inverse Mellin transform for the moments of the structure function. On this basis, one obtains for xg N S 1 (x, Q 2 ) the following expression in which a convenient path of integration is chosen: In this integration, it is assumed c = 1.9 and φ = 3π 4 [44]. Now we are able to construct the polarized structure function in the 2δanQCD model, based on the Jacobi polynomials which is considered in the following section.

IV. EXTRACTING THE POLARIZED STRUCTURE FUNCTION, USING THE 2δANQCD MODEL
Here we are going to employ the 2δanQCD model to construct the polarized nucleon structure function. This model, like the underlying pQCD based on Eq.(21) as the coupling constant of series expansions, is considered in the Lambert scheme.
In the anQCD approach, we replace the power a ν s ≡ a ν /4 ν of the conventional pQCD by the analytic coupling A ν /4 ν . Note that A ν (Q 2 ) is the analog (image) of a(Q 2 ) ν ≡ (α s (Q 2 )/π) ν , and a s (Q 2 ) = a(Q 2 )/4. In this way, the moments ∆M NS as the analytic images of ∆M NS are obtained. We recall that ν in the analytic coupling A ν is an expansion order index rather than a power index. For simplicity, we will denote in this section i.e., A s,ν is the anQCD image (analog) of a ν s , while A ν is (always) the anQCD image (analog) of a ν . When the first rational approximation in Eq. ( 30) is employed, the mentioned replacement in Eq. (24) leads to the following result [18]: We then use this formalism, and evaluate A s,ν with the help of the Mathematica package that is called anQCD.m as introduced in Ref. [8].
Using the corresponding command for the analytic coupling A ν in the 2δanQCD model in the two-loop approximation, we write [8] Substituting Eq. (47) in Eq. (46), with Q2 = Q 2 > 0, leads to the numerical result for the moments of the structure function in the anQCD approach. We note that there are some other commands in Ref. [8] for the anQCD coupling in 2δanQCD model but the one in Eq. (47) is more appropriate for our numerical purpose with sufficient precision. Namely, A ν (Q 2 ) = A2d2l[3, 0, ν, Q2, 0], for positive Q 2 = Q2 > 0 and with N f = 3, is constructed as the truncated sum of two terms in Eq. (18), i.e., A ν = A ν + k 1 (ν) A ν+1 , which is sufficient for our NLO analysis. (x, Q 2 ) as a function of Q 2 at NLO. The dashed lines represent the 2δanQCD results, using ∆m (1) anQCD and ∆m (2) anQCD in Eqs. (48 and 49). The solid line is the underlying pQCD one.A comparison with the available experimental data [45][46][47]49] has also been done.
As more explanation, it should be said that A2dN l[N f , n, ν, |Q 2 |, φ] is representing the N -loop analytic 2δanQCD coupling A (2δ) n+ν (Q 2 , N f ) of fractional power n + ν where ν > −1 and index n is such that n = 0, 1, . . . , N − 1. Here the active quark flavor N f is fixed. In the Euclidean domain the energy scale Q 2 is given by Q 2 = |Q 2 | exp(iφ) ∈ C\[−M 2 thr. , −∞) in which M 2 thr. = M 2 2 that is applicable for the N n−1 LO truncation approach [8].
It is evident that the anQCD analog of Eq. (30) is ∆m (2) anQCD (n, Q 2 ) ≃ One can construct from Eqs. (30) a combined quantity such that ∆m (1,2) = |∆m (1) − ∆m (2) |/∆m (1) . As a result we obtain an accuracy better than 1% for any n ≤ 11 while the expansion of Jacobi polynomials contains nine terms to yield us a good approximation. This accuracy is obtained for both underlying pQCD and analytic perturbation theory based on the 2δanQCD model. The numerical results for this accuracy at the energy scale Q 2 ≈ 0.17 GeV 2 are collected in Table II.
∆m (1) in the underlying pQCD: ∆m (1,2) pQCD ; and for anQCD in 2δanQCD, ∆m (1,2) anQCD . The results are presented for Q 2 ≈ 0.17 GeV 2 . This confirms that the first rational approximation, ∆m (1) anQCD (n, Q 2 ) in Eq. (48), has sufficient precision, i.e., the corresponding Eq. (46) gives us the desired results for the moments of the polarized nucleon structure function in the 2δanQCD approach. The results, using these two approximations, are depicted in Fig. 2, i.e., the polarized structure function g (N S) 1 as a function of Q 2 at x = 0.03. The plots in this figure denoted by 2δanQCD (1) and 2δanQCD (2) , which are related to the two approximate expansions in Eqs. (48) and (49), are in good agreement with each other and also with the available experimental data. In this figure, some of data are from [45][46][47] which have been collected in Table 1 of Ref. [48]. The data at very low-energy scales, i.e. Q 2 < 0.5 GeV 2 , is from Ref. [49]. The plot for g (N S) 1 at x = 0.03 provides an opportunity to compare the available experimental data at low-energy scales with the theoretical prediction of 2δanQCD model, where the underlying pQCD does not give an acceptable behavior. This confirms the advantage and applicability of the anQCD approach in comparison to pQCD, especially at the low-energy scales. In Fig. 3 we present the results for g x, at Q 2 = 10, 3, 1.5 and 0.5 GeV 2 in 2δanQCD and in the underlying pQCD, and compare them with the available experimental data from SLAC [50], HERMES [51], SMC [46] and COMPASS [47] experimental groups. All these data have been collected in Table 1 of Ref. [48].
As expected, at high-and medium-energy scales, there is good agrement between the theoretical predictions and the experimental data. Figure 4 is like Fig. 3 but at low-energy scales Q 2 = 0.2, 0.12 and 0.07 GeV 2 . By (x, Q 2 ) as a function of x in the NLO approximation at low-energy scales. The solid line with smooth behavior represents the 2δanQCD result, and the dashed line the underlying pQCD one. The available experimental data [49] are included.
decreasing the energy scales, the difference between the 2δanQCD and the underlying pQCD results increases. At Q 2 = 0.07 GeV 2 , the underlying pQCD result for x ≥ 0.1 as a function of increasing x grows rapidly and also oscillates. As a consequence, we multiplied there the 2δanQCD results and the experimental results by a factor of 40 in order to facilitate the comparison with the pQCD results. Figure 4 indicates that the g (N S) 1 structure func-tion involves smooth behavior at very low-energy scales, using the 2δanQCD model. This can be considered as an advantage of this model in comparison with underlying pQCD.
There are not enough experimental data for g (N S) 1 at the very low-energy scales Q 2 . But there are individual experimental data for the polarized proton and neutron structure functions, g p 1 and g n 1 , respectively. Unfortunately, most of them are binned very differently with respect to the x-Bjorken variable and Q 2 scales and cannot be used to construct the data points even utilizing the linear interpolation. Nonetheless, there are limited data with the identical bins of x and Q 2 and few others which are binned not so much differently with respect to the x and Q 2 such that the linear interpolation method is applicable for them. These data are selected such as to be placed in the inelastic region. In fact, these data are selected such that they exclude the elastic and ∆ resonance particles. They are depicted in three different panels in Fig. 4. 1

V. Q 2 DEPENDENCE OF BJORKEN SUM RULE IN ANQCD APPROACH
Here we investigate the Bjorken sum rule (BSR) [52] in the 2δanQCD model. This sum rule is relating the spin dependence of quark densities to the axial charge.The BSR is important to understand the nucleon spin structure, that QCD can describe well the strong force in the polarized case; i.e. its Q 2 dependence reflects the strong force in the polarized case. BSR has been measured at SLAC, DESY and CERN [45,46,50,51,[53][54][55][56][57][58][59][60][61][62][63] via polarized deep inelastic lepton scattering process. It has also been measured at moderate values of Q 2 by Jefferson Lab (JLab) [64][65][66][67]. The Q 2 momentum that is probing the nucleon is related to the inverse of the space-time scale. BSR at high Q 2 was evaluated in conventional pQCD by Refs. [68][69][70][71] and has the form In this equation the spin-dependent proton and neutron structure functions are g p 1 and g n 1 , respectively. The strength of the neutron β decay is controlling by g A which is the nucleon axial charge. Soft gluon radiation in conventional pQCD causes the leading-twist term (known 1 We thank S. E Kühn for useful comments on this point. as µ 2 ), given by the first term in the right-hand side of Eq. (50), to depend smoothly on Q 2 . The d i coefficients up to α 4 s order can be found in Ref. [68]. The terms with µ 4 , µ 6 , etc., are the nonperturbative power corrections, usually known as higher twists (HT). The higher twists, in fact, reflect a correlation between quarks and gluons. For a good understanding of the nucleon structure at the low-energy scales, analyzing the HT is essential. This is why it is interesting to follow as well this subject in the anQCD approaches, as we mentioned in the introduction. We recall that in the introductory part, a brief review has been done independent of the specific scheme. Therefore, the references cited here are not also in the Lambert scheme (usually they are in the MS scheme), but the evaluations here are in the (c 2 = −4.9) Lambert scheme.
The anQCD modification of the BSR , using the 2δanQCD model, has the form where where d i (k) coefficients are in the Lambert scheme with c 2 = −4.9 (cf. Appendix A in Ref. [72]). The k parameter is the renormalization scale parameter k = µ 2 /Q 2 , the couplings are A n = A (2δ) n (µ 2 ), and we fixed the value of k to k = 1.
The higher-twist effects are included in Eq. (51). The term with dimension D = 2, i.e., µ p−n 4 /Q 2 , has the following coefficient [73]: In this equation the nucleon mass is M N ≈ 0.94 GeV.
The coefficient a p−n 2 represents the twist-2 target mass correction; d p−n 2 is related to the twist-3 matrix element [73]. These coefficients can be computed, using the following relations [73]: To calculate d p−n 2 we need the g 2 structure function. However, there is a relation which gives us g 2 in terms of the g 1 structure function in Ref. [74] but this relation includes just the twist-2 part of g 2 and it cannot be substituted in Eq.(53) since it is related to the twist-3 matrix element. Instead, one can resort to the quoted value for d p−n 2 in Refs. [48,67], a result from recent experimental analysis. The numerical value which is reported there, is d p−n 2 = 0.008 ± 0.0036.
Since one can compute the polarized structure functions g p 1 and g n 1 [see Eq. (43) it is required at first to achieve the polarized valence densities at initial scale Q 0 , given by Eq. (31). The unknown parameters of these densities can be obtained using the traditional global fit which has been done in Refs. [37,40]. It should be recalled that the polarized structure function g N S 1 can be calculated, based on the formalism which we applied in this article. But the numerical values of the unknown parameters for the required parton densities have been quoted from Refs. [37,40], listed in Table I. The f p−n 2 function in Eq.(52) is the expectation value of a well-defined operator with a specific physical meaning [75,76]. Here it is regarded as a parameter which parameterizes a power correction to the anQCD analysis. From this point of view, we take the anQCD evolution form (based on its pQCD form) such that as has been quoted in Ref. [72].
While f p−n 2 is evaluated in the 2δanQCD approach, the coefficients a 2 and d 2 , as twist-2 and twist-3 quantities, are obtained from the underlying pQCD and the analysis of related experimental data, respectively. In fact the first integral in Eq. (53) can be calculated from the computed polarized structure function g N S 1 . As we explained in the two previous paragraph, to compute it one needs to use the traditional global fit as in Refs. [37,40] to obtain the unknown parameters of the polarized valence densities in Eq.(31).
We extend our computations up to the fourth higher twist to include as well the µ 6 term. In our calculations we have in total two free parameters µ 6 and f p−n 2 (Q 0 ) which can be obtained by fitting to the available experimental data for the BSR. The fitted numerical values are µ 6 = 0.0007 ± 0.0001 and f p−n 2 (Q 0 ) = −0.020 ± 0.001 ,respectively. Based on Eq.(52) the numerical value for µ 4 would be µ 4 = −0.0031 ± 0.001.
Equation (51) is the analytical result for Γ p−n 1 in the 2δanQCD model. In Fig. 5, we present Γ p−n 1 in the 2δanQCD approach and compare with the result from the underlying pQCD where both are considered in the Lambert scheme, as well as with the E143 , E154 and E155 [50,53,63] and JLab [64][65][66][67] experimental data. The good agreement between the available experimental data and the analytical result from the 2δanQCD model, in contrast to the underlying pQCD result especially at low energies, confirms that the anQCD approach, based on the employed model, is working well. (Q 2 ), resulted from two models.The red line gives the result from 2δanQCD model and black line from the underlying pQCD. Other symbols show the data from E143 , E154 and E155 [50,53,63] and JLab [64][65][66][67] experiments.

VI. SUMMARY AND CONCLUSION
As a new theoretical approach to analyze the nonsinglet polarized structure function at the low Q 2 continuum region, the analytic QCD (anQCD), specifically the 2δanQCD model, is employed to evaluate this function and also the Q 2 dependence of Bjorken sum rule.
Using this approach, it is possible to analyze the polarized structure function in the whole Q 2 range at the leading-twist order. On this basis, one may conclude some outstanding characteristics arising from the anQCD approach. At first, it can be shown that the g (N S) 1 (x, Q 2 ) structure function at a fixed x value and in the whole range of Q 2 is slowly changing. The plot in Fig. 2 reveals this feature. Second, the results of the underlying pQCD and anQCD approaches at moderate and high Q 2 scales are consistent with each other as can be seen from the depicted plots in Fig. 3. Third, at low values of Q 2 , the evolution of g (N S) 1 (x, Q 2 ) as a function of x in the anQCD approach is slower and smoother in comparison with the underlying pQCD which are both in the Lambert scheme. This fact can be seen from the three plots in Fig. 4 at the low-energy scales Q 2 =0.2, 0.17 and 0.07 GeV 2 , respectively. Consequently, we conclude that the result from analytic series in the anQCD approach, using the 2δanQCD model, better reproduces the data and behaves more smoothly than in underlying pQCD, especially at low Q 2 values.
We also investigated the Q 2 dependence of the Bjorken sum rule in the anQCD approach, using the 2δanQCD model. The result is shown in Fig. 5 which indicates good agreement with the available experimental data in the whole interval of energy scales. The agreement between them especially at the low-energy scales confirms the validity of the employed anQCD approach. A comparison has also been done with the result from the underlying pQCD approach, suggesting an advantage of the anQCD over the pQCD approach in the evaluation of spacelike QCD observables at low energies .
All these features have their origin in the fact that the behavior of analytic (holomorphic) coupling is under control in the entire range of spacelike Q 2 values, including at Q 2 → 0. This outstanding feature has been presented in Fig. 1.
As a final point we mention that the plots in Fig. 2, toward small values of Q 2 , resulted from Eq. (46) in the 2δanQCD model. In this case the Jacobi polynomial expansion is used [see Eq. (43)] to convert the result from the moment space to x-Bjorken space. In Eq. (43) ,the notation ∆M [xg N S 1 , j + 2] for the polarized moment of structure function is equivalent to the ∆M NS (n, Q 2 ) in Eq. (46), with n = j + 2. It should be noted that in converting the structure function from the moment space to Bjorken-x space in underlying pQCD, we used the inverse Mellin technique based on Eq. (44), while for the anQCD approach, due to computational difficulties which arise in employing inverse Mellin technique, we resorted to employing the technique of Jacobi polynomial expansion.
A further research task and an interesting subject would be to extend the anQCD analysis to the nonsinglet case for the unpolarized structure function which contains heavy quark flavors.
In this case the effect of quark mass cannot be ignored, and it should be considered in the calculations. In this regard the related commands in the 2δanQCD model would be different with respect to what we used here for light quarks. Also of interest would be to consider the singlet case of structure functions in the anQCD approach, a subject on which we hope to report in the future.

ACKNOWLEDGMENTS
The authors are indebted to G. Cvetič and C. Ayala for reading the manuscript and providing crucial comments which helped us improve the manuscript. We are grateful to S. E. Kühn for the very essential discussions about the CLAS data. We also acknowledge G. Mallot for productive comments and providing us with the required COMPASS data. We are finally thankful to Charles C. Young for his useful comments in connection to the related CERN data. S.A.T is grateful to the School of Particles and Accelerators, Institute for Research in Fundamental Sciences (IPM) to make the required facilities to do this project.

APPENDIX A: THE SPLITTING FUNCTIONS IN THE NLO APPROXIMATION
The analytical expression for the splitting functions in the moment space for nonsinglet sectors and at NLO approximation are taken from [34] and presented below. The usual quadratic Casimir operators are fixed to their exact values, using C A = 3, T F = 1 2 and C F = 4 3 .