Bjorken polarized sum rule and infrared-safe QCD couplings

Experimental data obtained for the polarized Bjorken sum rule (BSR) $\Gamma_1^{p-n}(Q^2)$ are fitted by using predictions derived within a truncated operator product expansion (OPE) approach to QCD. Four QCD versions are considered: perturbative QCD (pQCD) in the MSbar scheme, Analytic Perturbation Theory (APT), and 2$\delta$ and 3$\delta$ analytic QCD versions. In contrast to pQCD, these QCD variants do not have Landau singularities at low positive $Q^2$, which facilitates the fitting procedure significantly. The fitting procedure is applied first to the experimental data of the inelastic part of BSR, and the known elastic contributions are added after the fitting. In general, when 2$\delta$ and 3$\delta$ QCD coupling is used the fitted curves give the best results, within the $Q^2$-range of the fit as well as in extended $Q^2$-intervals. When the fitting procedure is applied to the total BSR, i.e., to the sum of the experimental data and the elastic contribution, the quality of the results deteriorates significantly.


I. INTRODUCTION
The polarized Bjorken sum rule (BSR) Γ p−n 1 (Q 2 ) [1] is an important spacelike QCD observable for various reasons. It is a difference of the first moment of the spin-dependent structure functions of proton and neutron, therefore its isovector nature makes it easier to describe it theoretically, in pQCD in terms of OPE, than the separate integrals of the two nucleons. Further, high quality experimental results for this quantity, obtained in polarized deep inelastic scattering (DIS), are now available in a large range of spacelike squared momenta −q 2 ≡ Q 2 : 0.054 GeV 2 ≤ Q 2 < 5 GeV 2 [2][3][4][5][6][7][8][9][10]. In particular, the newer experimental results [4] with highly reduced statistical uncertainties, extracted mainly from the Jefferson Lab CLAS EG1-DVCS experiment [11] taken on polarized targets made of protons and deuterons, render BSR an attractive quantity to test on it various extensions of pQCD to low Q 2 1 GeV 2 .
In this work we fit the theoretical OPE expressions to the experimental BSR results in pQCD, in (F)APT, and two additional extensions of QCD to low Q 2 , namely the 2δ [22,23] and 3δ [24,25] AQCD. The latter two extensions have the coupling A(Q 2 ) [the analog of the pQCD coupling a(Q 2 )] which is free of Landau singularities and physically motivated in the entire relevant regime of Q 2 in the complex plane, Q 2 ∈ C\(−∞, −M 2 thr ], where M 2 thr 1 GeV 2 is a positive threshold scale. The present work is an extension of our previous work on BSR [26]; in comparison with the latter work, we now vary the Q 2 -range of the fit, include the consideration of the elastic contribution, and estimate the uncertainties of the values of the extracted fit parameters due to systematic (in addition to statistical) uncertainties of the experimental BSR data.
The mentioned QCD variants (F)APT, 2δ and 3δ AQCD were constructed directly by imposing certain physically motivated conditions on the QCD coupling. In this context, we point out that there exist several other approaches to the construction of the QCD coupling, among them those using Dyson-Schwinger equations which involve various versions of the dynamically generated gluon mass (the mentioned 2δ AQCD gives a coupling with similar properties). For a recent review of various approaches, we refer to [27].
In Sec. II we present the theoretical leading-twist and higher-twist OPE contributions for the considered quantity BSR, as well as the parametrization of the elastic contribution to BSR. In Sec. III we present the results of various fits of theoretical expressions to the experimental results. Finally, in Sec. IV we summarize our conclusions.
Most of the formal aspects of the calculations are relegated to Appendices: in Appendix A we present the form of the leading-twist perturbation coefficients for a general renormalization scale and scheme; in Appendix B we present construction of A n , the analogs of pQCD powers a n , in extensions of pQCD without Landau singularities; in Appendix C we summarize such extensions [(F)APT, 2δ and 3δ]; in Appendix D we explain how the statistical and systematic uncertainties of the experimental data are reflected in the corresponding uncertainties of the parameters extracted in the fits; and in Appendix E we estimate the effects of the finiteness of the charm quark mass in our evaluations.

II. BJORKEN SUM RULE: THEORETICAL EXPRESSIONS
The polarized Bjorken sum rule (BSR), Γ p−n 1 , is defined as the difference between proton and neutron polarized structure functions g 1 integrated over the whole x-Bjorken interval Based on the various measurements of these and the related structure functions, the inelatic part of the above quantity, Γ p−n 1 (Q 2 ) inel. , has been extracted at various values of squared momenta Q 2 j (0.054 GeV 2 ≤ Q 2 j < 5 GeV 2 ): [2][3][4] (Jefferson Lab), [7,8] (SLAC), [5,6] (COMPASS at CERN). 1 Theoretically, this quantity can be written in the Operator Product Expansion (OPE) form [1] Γ p−n,OPE Here, |g A /g V | is the ratio of the nucleon axial charge, (1 − D BS ) is the perturbation expansion for the leading-twist contribution, and µ 2i /Q 2i−2 are the higher-twist contributions. The value obtained from neutron β decay is known to a high accuracy, |g A /g V | = 1.2723 ± 0.0023 [28], and we will use the central value. In the higher-twist terms, we will take only the terms ∼ 1/Q 2 and 1/(Q 2 ) 2 .

A. Perturbation expansion of the leading-twist
The leading-twist term has the canonical part D BS (Q 2 ) whose perturbation expansion in a ≡ α s /π is known up to N 3 LO (∼ a 4 ) where the bar indicates that the expansion is in the MS scheme, and the renormalization scale µ 2 is implicitly understood to be equal to the physical scale Q 2 . The NLO, N 2 LO and N 3 LO coefficientsd j (j = 1, 2, 3) were obtained in [29][30][31], respectively. In the considered range of momentum transfer 0 < Q 2 < 5 GeV 2 , we will assume that the effective number of active quark flavors is N f = 3, and therefore only the nonsinglet (NS) contributions appear. 2 When the renormalization scale is changed from µ 2 = Q 2 to a general value µ 2 = kQ 2 (0 < k ∼ 1), and when the renormalization scheme parameters are changed from the MS valuesc j ≡β j /β 0 to general scheme parameter values c j (j ≥ 2), the perturbation expansion changes accordingly D BS (Q 2 ) pt = a(kQ 2 ) + d 1 (k)a(kQ 2 ) 2 + d 2 (k; c 2 )a(kQ 2 ) 3 + d 3 (k; c 2 , c 3 )a(kQ 2 ) 4 + O(a 5 ).
The expressions for the new coefficients d 1 (k), d 2 (k; c 2 ) and d 3 (k; c 2 , c 3 ) are obtained on the basis of independence of the observable D BS (Q 2 ) pt from k and c j (j ≥ 2), and are given in Appendix A.
In those versions of AQCD where the coupling is a holomorphic function A(Q 2 ) [instead of the nonholomorphic a(Q 2 )] with nonperturbative contributions, the power expansion (4) becomes a nonpower expansion where a n get replaced by A n ( = A n ) The construction of the power analogs A n of a n were obtained in Ref. [34,35] for integer n and in Ref. [36] for general real n > −1. A brief description for the construction of A n is given in Appendix B. These expressions are based on the renormalization group equation (RGE), in close analogy with the RGE in the perturbation theory. The couplings A n (Q 2 ) can be obtained once the coupling A(Q 2 ) is known. The construction of A(Q 2 ) coupling is summarized in Appendix C for various variants of QCD with holomorphic coupling: (F)APT, 2δ and 3δ AQCD, and we refer to that Appendix for more details.

B. Higher-twist
In the theoretical OPE expression (2), the term with the dimension D = 2 (i.e., ∝ 1/Q 2 ) has the coefficient where M N ≈ 0.94 GeV is the nucleon mass, a p−n 2 is the (twist-2) target mass correction, and d p−n 2 is a twist-3 matrix element At Q 2 = 1 GeV 2 , we have a p−n 2 = 0.031 ± 0.010 and d p−n 2 = 0.008 ± 0.0036 [4]. We will neglect Q 2 -dependence of these two quantities [as was done also in Ref. [4]], and will take the central values, i.e., a p−n 2 + 4d p−n 2 = 0.063. On the other hand, the coefficient f p−n 2 (Q 2 ) will be a parameter of the fit, and its Q 2 -dependence will not be neglected, its evolution is known [37,38] in pQCD where ν ≡ γ 0 /(8β 0 ) = 32/81 when N f = 3, and the reference scale will be taken Q 2 in = 1 GeV 2 . In QCD with holomorphic coupling A(Q 2 ), the power a γ0/8β0 gets replaced by A γ0/8β0 (which is in general not equal to the simple power A γ0/8β0 , as mentioned above, cf. also Appendix B, In addition, in some of the fits we will also include the D = 4 term in the theoretical OPE expression (2) µ 6 /(Q 2 ) 2 , where we will consider the coefficient µ 6 as Q 2 -independent. Thus, our theoretical expression for BSR will be the truncated (at D = 4, or D = 2) OPE expression where the leading-twist truncated expressions are given in Eqs. (4) and (5), and the Q 2 -dependent part of the D = 2 term (twist-4) in Eqs. (8) and (9), for the pQCD and AQCD version of QCD, respectively. We will regard the renormalization scale parameter k ≡ µ 2 /Q 2 and the higher-twist coefficients f p−n 2 (Q 2 in ) (with Q 2 in = 1 GeV 2 ) and µ 6 as the free parameters to be determined in the fitting procedure.
In our approach, we will include in the leading-twist part of the OPE all the known terms (i.e., up to order ∼ a 4 ∼ A 4 ). The order of the leading-twist terms in general affects the fitted higher-twist contributions (cf. [39][40][41] for the fit of truncated OPE to structure functions).

C. Elastic contribution
The OPE evaluation is, in principle, for inclusive observables; in the case of BSR, this means that the OPE fit should be applied to the experimental values of the sum of the inelastic and elastic contribution.
The elastic contribution to BSR can be expressed [42,43] in terms of the proton and neutron electromagnetic form factors F 1 and F 2 The most recent parametrization of these form factors was performed in [44] from light-front holographic QCD (LFH QCD). Namely, these form factors can be expressed in terms of the inverse power expressions which are products of τ − 1 poles along the vector meson Regge radial trajectory in terms of the ρ-vector mass M ρ = 0.755 GeV and its radial excitations. We have Here: χ p = µ p − 1 = 1.793 and χ n = µ n = −1.913 are the anomalous moments of p and n, respectively; γ n,p are the higher Fock probabilities for the L = 0 → L = 1 spin-flip electromagnetic form factors, and r is a phenomenological factor. The values of these three parameters are obtained by fitting to polarization data for the form factors [44]: γ p = 0.27; γ n = 0.38; r = 2.08. As a consequence, the elastic contribution (11) can be represented as a combination of inverse powers (1 + Q 2 /M 2 ρ /(2m − 1)), and for large Q 2 this can be expanded in inverse powers of Q 2 This means that at high Q 2 the elastic contribution starts with the dimension D = 8 term ∼ (1/Q 2 ) 4 . Theoretically, it is not included in the truncated OPE expression (10) in pQCD. 3,4 Further, in 3δ and 2δ AQCD, the higher dimensional terms up to (and including) the D = 8 term are not included in the leading-twist contribution (5) because in these approaches A(Q 2 ) − a(Q 2 ) ∼ (Λ 2 /Q 2 ) 5 as explained in Appendix C 2 Eq. (C7). Therefore, the truncated OPE expression (10) does not contain the dimension D = 8 term ∼ (1/Q 2 ) 4 also in 3δ and 2δ AQCD approaches. The only exception is the (F)APT where A(Q 2 ) − a(Q 2 ) ∼ (Λ 2 /Q 2 ) 1 , cf. Appendix C 1, and where the truncated OPE (10) could contain, in principle, all the elastic contributions, including the ∼ (1/Q 2 ) 4 term.
Therefore, it looks more natural to first fit the truncated OPE expression (10), with pQCD and (3δ and 2δ) AQCD, to the BSR experimental results [2-4, 6, 7] which are obtained for the inelastic contribution; and after such a fit, add the elastic contribution (11) [with the parametrization (12)-(13b)]. We note that the theoretical expression obtained for the total BSR in this way is the sum of the expressions (10) and (11), which is again an OPE expression as it should be for such an inclusive spacelike observable as the total BSR.
The other approach would be to fit the truncated OPE expression (10) with the experimental points for the total BSR, i.e., fit with the experimental points for inelastic BSR with the elastic contribution added to them. Such an approach seems less natural, at least in pQCD and 3δ and 2δ AQCD, because the truncated expression (10) in these cases in principle does not contain the leading elastic contribution ∼ 1/(Q 2 ) 4 , cf. Eq. (C7) in Appendix C 2.
Nonetheless, below we will apply, for completeness, both approaches in our numerical fitting procedures.

III. NUMERICAL FITS
In the numerical fits, we will consider the following parameters to be fitted in the truncated OPE expression (10): (i) the renormalization scale parameter k ≡ µ 2 /Q 2 of the (truncated) leading-twist contribution (4) [for pQCD] or (5) [for AQCD]; (ii) the parameter f p−n 2 (1 GeV 2 ) appearing in the D = 2 term of Eq. (10) [cf. Eqs. (8) and (9)]; (iii) and in some fits the parameter µ 6 in the D = 4 term of Eq. (10) will be included. We will take for the experimental data points for the inelastic contribution to BSR the data from [2-4, 6, 7]; 5 they are in the momentum interval 0.054 GeV 2 ≤ Q 2 j < 5 GeV 2 . Among these data points, we exclude from the fit the four points with Q 2 j ≥ 3 GeV 2 (three from [4], at Q 2 j = 3.223, 3.871 and 4.739 GeV 2 ; and one from [6] at Q 2 = 3 GeV 2 ) because they tend to decrease the quality parameter χ 2 /(d.o.f.) of the fit significantly. Nonetheless, we will show both quality parameters 3 However, the last coefficient in the truncated OPE, e.g. µ 6 , is sometimes in the literature regarded to include, in a certain effective way, the contributions from the higher dimensional terms D ≥ 6, cf. [4] for a discussion of these aspects. 4 The elastic contributions parametrized in a way different from that of Ref. [44] could in general contain terms of dimension D < 8.
Nonetheless, the fitted expressions of Ref. [43] give for the elastic contributions an expansion similar to Eq. (14), where the first nonzero term is a very small D = 6 contribution which is negligible in comparison with D = 8 term (the coefficient at D = 6 term is by about a factor of 10 −4 smaller than at D = 8 term). All these terms there are divided by a large weakly Q 2 -dependent factor [ln(27.81 + 1.72 Q 2 /M 2 ρ )] 4.296 . 5 We will not take into account the SLAC E155, SMC and HERMES points [8][9][10], at Q 2 = 5 GeV 2 , 10 GeV 2 , and 2.5 (and 5) GeV 2 , respectively, as they were obtained by NLO pQCD Q 2 -evolution from data distributed across a wide range of Q 2 . Among the points of SLAC E143 we exclude the point at Q 2 = 5 GeV 2 for the same reason.
obtained from the mentioned fits, i.e., the one with the four points excluded, and included. 6 The number of active flavors will be kept all the time at N f = 3. The fits will be performed by the least squares method, taking into account the statistical uncertainties σ j,stat in the points Q 2 j , we will consider them independent of each other. These uncertainties σ j,stat are in general significantly smaller than the systematic uncertainties σ j,sys . The latter are strongly correlated, and we will consider them as completely correlated. The uncertainties in the values of the extracted parameters k ≡ µ 2 /Q 2 , f p−n 2 (Q 2 in ) and µ 6 are then due to statistical (small) and systematic (larger) uncertainties of the data. We refer to Appendix D on how we obtained the uncertainties of the extracted values of the fit parameters.
Using specific values of the renormalization scale parameter k ≡ µ 2 /Q 2 may allow us to incorporate in our evaluations (5) at least a part of the contributions of higher orders of the series. It is expected that k ∼ 1, and usually it is taken in the literature in the range 1/2 < k < 2, sometimes 1/4 < k < 4. In the considered work, after replacing the powers of the pQCD coupling by their analytic counterparts, cf. Eqs. (4)-(5), a spacelike observable depends usually weakly on the contributions of higher orders. Still, in order not to miss the possibly relevant influence of higher orders, we decided to increase the range of possible k values to: 1/16 < k < 16. This choice still avoids very large or very small renormalization scales where the corresponding coefficients d n (k) at the powers a(kQ 2 ) n+1 , or at their analytic counterparts A n+1 (kQ 2 ), contain powers of the large terms ∼ ln k [cf. Eqs. (A3)] which may destroy the convergence of the series already at low n. We will see that the results of the fitting will rarely give the extreme values k = 16 or k = 1/16 ≈ 0.063, mostly in the cases of (MS) pQCD and (F)APT.
We present below the results of the fits of the inelastic contributions to the truncated OPE expression (10), first in Sec. III A with µ 6 = 0, and then in Sec. III B for µ 6 = 0 included as a fit parameter; the elastic contribution is added afterwards. In Sec. III C we present the corresponding fits for the case when the higher-twist term in the OPE has a mass parameter. In Sec. III D we include the fits for two different ansätze for BSR at very low Q 2 . Finally, in Sec. III E we present the fits of the truncated OPE expression (10), with µ 6 = 0 included as a fit parameter, applied to the total BSR, i.e., to the sum of the inelastic data points (at Q 2 j ) and the corresponding elastic contribution Γ p−n 1 (Q 2 j ) el. . In each case, the fits are performed for four variants of QCD: in pQCD (in MS), in (F)APT (in MS), and in 2δ [23] and 3δ AQCD [25]. Each of these fits is performed by excluding among the data points those with Q 2 < Q 2 min , where Q 2 min = 0.268 or 0.66 GeV 2 (and sometimes also: Q 2 min = 0.47 GeV 2 ). As mentioned earlier, the four points at high Q 2 ≥ 3 GeV 2 are excluded from the fit as well.
It is reasonable to assume that the number of active quark flavors in the considered interval 0.054 GeV 2 < Q 2 < 5 GeV 2 is N f = 3, cf. Appendix E. All the fits will be performed in MS pQCD approach, in the Fractional Analytic Perturbation Theory [(F)APT], and in 2δ and 3δ AQCD approach, cf. Appendix C where also the values of the parameters of these three QCD variants are presented. In MS pQCD and in 2δ and 3δ AQCD the (underlying) pQCD running coupling a(Q 2 ) is determined by the requirement πa(M 2 Z ; MS) = 0.1185 [28,45], and in (F)APT we use Λ 3 = 0.45 GeV (at the end of Sec. III B we also comment on the case Λ 3 = 0.40 GeV); we refer to Appendix C for more details.  Table I, for Q 2 min = 0.66, 0.47 and 0.268 GeV 2 . As mentioned earlier, the fit is performed for the data in the momentum interval Q 2 min ≤ Q 2 j ≤ 3 GeV 2 (with the COMPASS data point [6] excluded), and the corresponding quality parameter for that interval is denoted as χ 2 /d.o.f. In addition, χ 2 ext /d.o.f. is the quality parameter for the wider interval Q 2 min ≤ Q 2 j < 5 GeV 2 ; is the parameter when 0.268 GeV 2 ≤ Q 2 j ≤ 3 GeV 2 ; and χ 2 all /d.o.f. is the parameter when all the experimental points [2-4, 6, 7] are included, 0.054 GeV 2 ≤ Q 2 j < 5 GeV 2 . In these quantities χ 2 , the central values of k and f p−n 2 (Q 2 in ) (Q 2 in = 1 GeV 2 ) parameters obtained from the mentioned fit interval Q 2 min ≤ Q 2 j ≤ 3 GeV 2 are used; we take here "d.o.f" as (N − p) where N is the number of fitted data points entering the considered χ 2 and p   is the number of parameters of the fit (p = 2 here)

QCD variant
The values of these quantities are always large, in the best cases between 1 and 10. This is so because the statistical uncertainties of the newer JLAB data [4] are very small, σ j,stat 10 −3 , and simultaneously, our fit function (truncated OPE) is different from the ideal function which we do not know. Nonetheless, they decrease when analytic QCD variants are employed, especially 2δ and 3δ AQCD.
We wish to point out that the approach of (MS) QCD in the case of Q 2 min = 0.268 GeV 2 is, in principle, not applicable. This is so because the corresponding coupling a(Q 2 ) has a Landau branching point at Q 2 branch = 0.371 GeV 2 , which makes the D = 2 running coefficient f p−n 2 (Q 2 ), Eq. (8), undefined at Q 2 ≤ 0.371 GeV 2 . Nonetheless, in order to be able to present a curve, we applied in the fitting case Q 2 min = 0.268 GeV 2 in the MS pQCD approach a restriction on the leading-twist renormalization scale µ 2 = kQ 2 , namely k > 1.383; and this not just in the leadingtwist contribution, but we also made an ad hoc replacement in the D = 2 running coefficient f p−n 2 (Q 2 ), Eq. (8): . We applied this also in the case when µ 6 = 0 (Sec. III B, for MS pQCD approach with Q 2 min = 0.268 GeV 2 ). In other approaches (APT and AQCD's) this is not necessary, as there are no Landau singularities. Further, we note that in the last column in Table I we have χ 2 /d.o.f. = ∞ in the case of MS pQCD for Q 2 min = 0.66 and 0.47 GeV 2 . This is so because at the lowest available experimental point, Q 2 j=1 = 0.054 GeV 2 , the coupling is a(kQ 2 1 ; MS) = ∞ because kQ 2 1 < Q 2 branch = 0.371 GeV 2 and thus we hit Landau singularities there. We can deduce from Figs. 1 and Table I: (a) the best results in the considered approach (µ 6 = 0) are obtained in 3δ AQCD; (b) the quality of extrapolation of the obtained fitted curves from the fitting interval, Q 2 min ≤ Q 2 ≤ 3 GeV 2 , to the entire interval, 0.054 GeV 2 ≤ Q 2 < 5 GeV 2 , does not improve significantly when the fitted interval is extended (i.e., when Q 2 min is lowered), cf. also the last column in Table I. This indicates that the behavior of the curves in the extrapolated regions, 0.054 GeV 2 < Q 2 < 0.268 GeV 2 and 3 GeV 2 < Q 2 ≤ 5 GeV 2 , is of similar quality in the cases of different values of Q 2 min (fit), and gives the dominant part of χ 2 all /d.o.f.. The same can be observed for the values of f p−n 2 (1): they do not depend much on the value of Q 2 min (fit), but only on the QCD variant used in the fit. This appears to be related with the fact that only one parameter beyond the leading-twist contribution is used here (µ 4 ↔ f p−n 2 ), representing a truncated OPE ansatz with a significantly restricted freedom. , at the higher Q 2 and the lower Q 2 < 1 GeV 2 momenta, respectively. The corresponding results are given in Table  II, for Q 2 min = 0.66, 0.47 and 0.268 GeV 2 , with the same notations as in Table I. The various versions of     (1 GeV 2 ) and µ6 (in GeV 4 ), obtained by fitting the OPE expression (10), truncated at D = 4 (µ6 = 0). The notations are as in Table I are those as explained in the previous Sec. III A, except that now in the relation (15) the factor in front of the sum Comparing Table II and Figs. 2-3 with Table I Tables I and II). This means that the extrapolation down to the lowest experimental point Q 2 = Q 2 j=1 = 0.054 GeV 2 is better when µ 6 = 0 (with the exception of MS pQCD case where problems with Landau singularities appear). This occurs because the inclusion of the µ 6 /(Q 2 ) 2 term in the truncated OPE makes this expression less stable at very low Q 2 .
On the other hand, when Q 2 (fit) ≥ 0.268 GeV 2 , some of the fits (2δ and 3δ AQCD) give better extrapolation when µ 6 = 0. Table II and Figs. 3 also show that when µ 6 = 0 and Q 2 (fit) ≥ 0.268 GeV 2 , the best extrapolation to low Q 2 is obtained in 3δ AQCD, followed by 2δ AQCD. As Fig. 3(b) suggests, the fitted curve in pQCD MS extrapolated to low Q 2 appears to be almost as good; in this case, however, we should keep in mind that the renormalization scale is µ 2 = kQ 2 (k = 3.00) and that this scale was used also in f p−n 2 (Q 2 ), i.e., the ad hoc replacement f p−n 2 (Q 2 ) → f p−n 2 (kQ 2 ) was performed in order to avoid the Landau singularities in the D = 2 term at Q 2 ≥ 0.268 GeV 2 (cf. also the discussion about that point in Sec. III A). Despite this replacement, in χ 2 all /d.o.f. Landau singularities are hit, because the three lowest experimental points give kQ 2 j < Q 2 branch = 0.371 GeV 2 (j = 1, 2, 3) and are thus in the Landau singularity region. 7 In the case of (F)APT, in contrast to 2δ and 3δ AQCD, the coupling A(Q 2 ) differs from the underlying pQCD coupling a(Q 2 ) nonnegligibly at high |Q 2 | > 1 GeV 2 [the index in Eq. (B1) is N = 1 in (F)APT; N = 5 in 2δ and 3δ AQCD]. This implies that (F)APT has a certain ambiguity when "normalizing" the strength of A(Q 2 ). As mentioned, we fixed the strength of the coupling A(Q 2 ) in (F)APT to the value Λ N f =3 = 0.45 GeV, since in such a case (F)APT reproduces approximately the QCD phenomenology at high energies (cf. also [15]). A question appears whether the results of our fits in (F)APT depend significantly on this value. We repeated the analysis in (F)APT with the value Λ N f =3 = 0.40 GeV, and it turned out that the results of the fits did not change significantly.  Table II, the third line from below).
When we add the parametrized elastic contribution of BSR, Eqs. (11)- (13), to the experimental points and to the theoretical curves of Figs. 3, we obtain the results presented in Figs. 4(a),(b). In comparison with Figs. 3, the values of BSR are shifted to significantly higher values at low Q 2 . With this approach, the quality of fits (χ 2 /d.o.f) does not change, as we consider the elastic contribution as known (and parametrized) and added here simultaneously to the 7 We have Q 2 1 = 0.054 GeV 2 , Q 2 2 = 0.078 GeV 2 , and Q 2 3 = 0.101 GeV 2 . (inelastic) theoretical fitting curves and to the data points. This means that the results of Table II  (1) and µ 6 , are in the analytic variants of QCD smaller than in pQCD, this reduction being especially strong in the 3δ QCD variant. It has been noted in the literature that in pQCD OPE there is a duality between the order of truncation of the leadingtwist series and the higher-twist contribution [21,[39][40][41]: higher-twist contribution often significantly decreases with the inclusion of higher orders in the leading-twist part. This effect and ambiguity become stronger in the ranges where the perturbation theory becomes questionable (for example, at the large and low values of the Bjorken variable x, as it was shown in Refs. [46,47], respectively). It has been observed that the higher-twist contribution is smaller, but also more stable (under the inclusion of more terms in the leading-twist), in QCD variants with infrared modifications of the coupling (various modifications lead to quite similar results [48]). The latter probably incorporate a part of the higher-twist contributions (which are rather cumbersome [49]) into (formally) the leading-twist contribution for small x range at moderately small Q 2 values ( 1 GeV 2 ) (see Ref. [50] and more recent studies [51] of the precise combined H1 and ZEUS data [52] for the DIS structure function F 2 ).
Following the above observations in Tables II and I, we can conclude that the applications of the 2δ and especially 3δ AQCD are very appropriate frameworks for the BSR studies because they appear to resum effectively a large part of the perturbative contribution into the leading-twist part (5). In this context, we wish to recall that 3δ AQCD is significantly different from the other two AQCD variants [(F)APT and 2δ AQCD] in the infrared region, because its coupling is not just finite there but goes to zero, A(Q 2 ) ∼ Q 2 → 0, as motivated by large-volume lattice calculations, cf. [53][54][55][56] and Appendix C 2. Further, we wish to point out that the 2δ and 3δ AQCD couplings A (and thus A n and A n ) are at large Q 2 indistinguishable from their underlying pQCD couplings a ( a n , a n ), cf. Eq. (C7), in contrast to (F)APT which satisfies the relation Eq. (B1) with N = 1. Therefore, theoretically, neither the higher-twist contribution of order 1/(Q 2 ) N with N ≤ 4 (D ≤ 8), nor a part of it, is incorporated in the leading-twist contribution (5) in the 2δ and 3δ AQCD, in contrast to (F)APT. This indicates that the higher-twist terms extracted here with 2δ and 3δ AQCD (in truncated OPE) represent an effective form for the true higher-twist contribution with dimension D ≤ 8 [and a part of the other (D ≥ 10) presumably small contribution]. In pQCD this is definitely not so, because of the mentioned duality there between the order of truncation of the leading-twist series and the extracted higher-twist contribution. As a consequence, the extracted effective higher-twist contribution in pQCD represents a sum of the true higher-twist contribution and a significant part of the perturbative (leading-twist) contribution; this effective higher-twist contribution appears to be in general larger than the true higher-twist contribution. For comparison, we performed a similar fit, but now with a "massive" higher-twist term instead of the truncated OPE expression (10) where the squared mass M 2 in the denominator of the higher-twist part 8 is taken to be constant (not running), and is expected to be 0 < M 2 1 GeV 2 . Now, instead of f p−n

2
(1) and µ 6 , the fit parameters are f p−n

2
(1) and M 2 . The resulting curves, for Q 2 min = 0.268 GeV 2 , are given in Figs. 5(a),(b), at the higher Q 2 and the lower Q 2 < 1 GeV 2 momenta, respectively. The corresponding results are given in Table III. These curves are analogous to those in  Table III with the corresponding "nonmassive" case Figs. 3 and Table II, we see that the results and extrapolations in the case of 3δ AQCD are now comparably good in the "massive" and the µ 4 &µ 6 approaches. Stated differently, the corresponding χ 2 Tables III and II (with Q 2 fit = 0.268 GeV 2 ) are very similar. In the case of 2δ AQCD and (F)APT, the extrapolations are better in the "massive" than in the µ 4 &µ 6 approach, i.e., χ 2 all /d.o.f. is significantly reduced in the "massive" case. One reason for this lies perhaps in the fact that the massive higher-twist term is under control at very low Q 2 , unlike the separate µ 4 (Q 2 )/Q 2 and µ 6 /(Q 2 ) 2 terms. Further, the extracted values of f p−n

2
(1) are in general similar in the µ 6 = 0, µ 4 &µ 6 and the "massive" approaches, although the uncertainties of the extracted parameters are quite high in the "massive" approach. 9 The results of Table III show that the QCD variants with infrared-finite analytic coupling, and especially 3δ AQCD, give smaller values of higher-twist parameters f p−n (1) mean that at higher values of Q 2 the higher-twist contribution is significantly reduced; this can be seen also by expanding the massive higher-twist term of Eq. (16) in powers of M 2 /Q 2 . Such effect is in full agreement with one observed in Ref. [21]; the effect can be considered as a stabilization of the highertwist contribution. We also recall that 3δ AQCD is significantly different from (F)APT and 2δ AQCD in the infrared region, since its coupling goes to zero there, A(Q 2 ) ∼ Q 2 → 0.
It is possible to choose for the higher-twist term a massive form with a running mass, in the spirit of a dynamical effective gluon mass of the gluon propagator at low Q 2 [60]. Such masses appear in the literature often in definitions of QCD couplings at low Q 2 , and are responsible for the freezing (finiteness) of the coupling at Q 2 → 0. The couplings in (F)APT and 2δ AQCD have a freezing which could be described also via a running effective gluon mass. Our view is that the OPE higher-twist terms represent a new contribution not contained in the QCD coupling itself. The squared mass M 2 in such terms, Eq. (16), is considered constant, in the spirit of the approach of Ref. [57] (cf. also [58]), where the basic component (delta function) in the spectral function of the higher-twist contribution gives such a mass term. Nonetheless, we repeated the aforementioned analysis for the case of a running squared mass M 2 (Q 2 ) representative of a dynamical effective gluon mass, chosen with a simple parametrization of Ref.
where we chose for the parameters M and p values within the expected regions [61]. The adjustable squared mass scale was taken (instead of m 2 0 ) to be M 2 (1 GeV 2 ). The same analysis then gave the results presented in Table  IV. We can see that the results are qualitatively similar to those with the constant squared mass M 2 , Table III, except that the values of M 2 (1 GeV 2 ) are now lower. Nonetheless, a reasonable definition of the average squared mass value in the present analysis, for the considered fit interval 0.268 GeV 2 < Q 2 < 3 GeV 2 may be M 2 (Q 2 ) = (1/2) × (M 2 (0.268) + M 2 (3)) ≈ 1.52 × M 2 (1 GeV 2 ). Or, if we regard that the squared mass is most relevant only in the low-Q 2 part 0.268 GeV 2 < Q 2 < 1 GeV 2 of the fit interval, a reasonable definition of the average squared mass would be M 2 (Q 2 ) = (1/2) × (M 2 (0.268) + M 2 (1)) ≈ 1.85 × M 2 (1 GeV 2 ).

D. Testing low-Q 2 regime ansätze
At low Q 2 , the inelastic contribution to BSR behaves as ∼ Q 2 , according to Gerasimov-Drell-Hearn sum rule [62] as pointed out and used in [21,63,64]. Based on this, an expansion [3] motivated by chiral perturbation theory (χPT) can be constructed 9 The systematic uncertainties of the extracted parameters in Table III are large and should therefore be regarded as crude estimates only, cf. comments in Appendix D 3.
where, according to Gerasimov-Drell-Hearn sum rule [62], χ n and χ p are anomalous magnetic moments of nucleons [which appear also in the elastic BSR contributions, cf. Eqs. (13)]; the parameters A and B are determined in the fit. When we fit with this expression the inelastic BSR data [2][3][4] for Q 2 ≤ Q 2 max = 0.2 GeV 2 , we obtain A = 0.765 GeV −2 and B = 0.678 GeV −4 , with χ 2 (Q 2 ≤ Q 2 max )/d.o.f. = 0.720. On the other hand, if we take Q 2 max = 0.5 GeV 2 in the fitting, we obtain A = 0.744 GeV −2 and B = −1.033 GeV −4 , χ 2 (Q 2 ≤ Q 2 max )/d.o.f. = 1.313. Another possible ansatz for the inelastic contribution to BSR at low Q 2 is the form of the light-front holographic (LFH) effective charge A (LFH) in the BSR (g 1 ) scheme [65,66] When we fit with this expression the inelastic BSR data [2][3][4] for Q 2 < Q 2 max = 0.5 GeV 2 , we obtain κ = 0.479 GeV. 10 In Figs. 6(a),(b), we present these low-Q 2 expressions. The MS pQCD and 3δ AQCD curves (obtained from fit with Q 2 ≥ Q 2 min = 0.268 GeV 2 ) are also included for comparison. In Fig. 6(b) we can see that the theoretical fitted curves    (18) and (19) included. Fig. (b) is a zoomed-in version of Fig. (a), for Q 2 < 1 GeV 2 . The lower bound of the fitting interval, Q 2 min = 0.268 GeV 2 , is included as the thin dotted vertical line.
of the presented QCD variants connect smoothly with the nonperturbative low-Q 2 curves (χPT and LFH, both fitted to BSR up to Q 2 = 0.5 GeV 2 ) at the values of Q 2 around 0.3 GeV 2 . Nonetheless, we recall that the apparent success of the (MS) pQCD curve, down to about 0.3 GeV 2 , was achieved due to the ad hoc change of scale in f p−n 2 (Q 2 ) to kQ 2 (with k = 3), to avoid the Landau singularities in the D = 2 term at Q 2 ≥ 0.268 GeV 2 (cf. also the discussion about that point in Sec. III B). We further notice that the 3δ QCD curve agrees well with the mentioned χPT curve in a broader interval, 0.17 GeV 2 < Q 2 < 0.3 GeV 2 . Similar analyses with the goal of connecting the curves of pQCD (or of specific QCD variants) with nonperturbative curves at low Q 2 were performed in some of the references [21], in [67], and was discussed also in [68].
In Refs. [69,70], BSR at low Q 2 < 0.3 GeV 2 was calculated with baryon chiral perturbation theory (BχPT) up to NLO. The authors of [69,70] did not give the values of their parameters for BSR. However, careful visual comparison of their obtained Q 2 -dependence [cf. curves and bands in their Fig. 6(b) in Ref. [70]] with our χPT curve (18) with Q 2 fit ≤ 0.5 GeV 2 [the long-dashed curve in our Fig. 6(b)] shows very good agreement between them.

E. Fitting to total (inelastic + elastic) BSR data
We perform also the fitting of the truncated (at D = 4) OPE expression (10) to the data for the total BSR Γ p−n 1 (Q 2 ) inel.+el. . As argued at the end of Sec. II C [after Eq. (14)], such a fit has problematic aspects. These data are obtained by adding to the experimental data Γ p−n 1 (Q 2 j ) inel. of Refs. [2][3][4] the parametrized elastic contribution Eqs. (11)-(13) of Ref. [44]. The uncertainties of the parametrized elastic form factors are considered to be less than 10% [44]. Since the elastic contribution (11) is quadratic in the form factors (and is numerically dominated by the proton contribution), the relative uncertainties of the elastic contribution are less than 5%. Since we do not have more information about these uncertainties, we will neglect them in this analysis. This means that in the total BSR, we will consider that the statistical and the systematic uncertainties σ j,stat and σ j,sys are those of the inelastic contribution. Otherwise, the fitting is performed as in the previous Sec. III B.  The resulting curves are presented in Figs. 7(a),(b), for Q 2 ≥ Q 2 min with Q 2 min = 0.268 GeV 2 . The obtained results are given in Table V Table II (where the elastic part was not included in the fit procedure), we see that the inclusion of the elastic contribution in the fit procedure significantly deteriorates (increases) the values of the various fit quality parameters χ 2 , especially when the fit is performed in the larger interval (Q 2 min =)0.268 GeV 2 ≤ Q 2 < 3GeV 2 . In particular, the differences in quality are clearly visible when comparing Fig. 7(b) with Fig. 4(b) at low Q 2 . Only when Q 2 min is relatively high, Q 2 min = 0.66 GeV 2 , are some of the χ 2 parameters comparable in the two cases (but not the extrapolation quality parameters χ 2 0.268 and χ 2 all ). Further, the inclusion of the elastic contribution in the fit in general does not change significantly the extracted (mostly negative) values of f p−n 2 (1GeV 2 ), but increases significantly the (positive) values of µ 6 . For the fit analysis with the "massive" truncated OPE, Eq. (16), the results are presented in Table VI for Q 2 min = 0.268 GeV 2 . In general, the results are worse than in the corresponding "massive" OPE case fitted to the inelastic contribution (cf. Table III    part is not included in the fit); this is so because the elastic parts increase significantly BSR Γ p−n 1 (Q 2 ) at low Q 2 where they are represented mostly by the "massive" higher-twist term, and hence the relative uncertainties of BSR at low Q 2 become smaller.
In addition to the massive case with constant squared mass M 2 , we include also the case of Q 2 -dependent squared mass Eq. (17) in the analysis with the elastic contribution included in the fit. The results for this case are presented in Table VII. These results are qualitatively similar to those with constant squared mass, Table VI; and the comparison  of Table VII with its counterpart Table IV without the elastic contribution is similar to the above comparison when the squared mass is constant. We recall that an average mass in the Q 2 -independent case can be regarded to be

IV. SUMMARY
Experimental results for the polarized Bjorken sum rule (BSR) Γ p−n 1 (Q 2 ) were fitted, for various ranges of Q 2 , with OPE theoretical expressions using QCD couplings obtained in four different approaches: perturbative QCD (pQCD) in MS scheme; (Fractional) Analytic Perturbation Theory [(F)APT]; Two-delta AQCD (2δ); and Three-delta latticemotivated AQCD (3δ). The QCD running coupling A(Q 2 ) in the latter three QCD variants does not have Landau singularities, in contrast to the pQCD coupling a(Q 2 ) [≡ α s (Q 2 )/π].
In the fit of the inelastic experimental BSR results, up to two higher-twist terms [∼ 1/Q 2 , 1/(Q 2 ) 2 ] were added to the theoretical leading-twist contribution. The elastic contributions, which are ∼ 1/(Q 2 ) n with typically n ≥ 4 at Q 2 > 1 GeV 2 , were then added by using the parametrization obtained from the literature. The fits were performed for the ranges Q 2 min ≤ Q 2 ≤ 3 GeV 2 , where Q 2 min = 0.66, 0.47 and 0.268 GeV 2 . In general, the best curves were obtained when 2δ or 3δ-couplings were used. When only D = 2 (∼ 1/Q 2 ) higher-twist term was included in the fit, the quality of the fitted curves, in the range of the fit and in the extrapolated ranges of Q 2 , in general did not depend significantly on Q 2 min of the fit. On the other hand, when both D = 2 and D = 4 terms were included, the quality in the extrapolated ranges of Q 2 was in general better for the lowest Q 2 min value (0.268 GeV 2 ), i.e., when the Q 2 -range of the fit was the largest. Comparably good results were obtained when "massive" higher-twist term was used in the OPE and the QCD coupling was either from (F)APT or 2δ or 3δ AQCD.
When the range of fit had Q 2 min = 0.268 GeV 2 , the pQCD MS coupling approach worked and gave acceptable results only if the renormalization scale of the coupling was maintained everywhere at sufficiently high values, and the coefficient f p−n 2 (Q 2 ) [∼ a(Q 2 ) γ0/8β0 ] at the D = 2 term had an (ad hoc) increased scale Q 2 → kQ 2 , in order to avoid the problem of the Landau singularities.
When the fit procedure was performed by fitting the theoretical OPE, truncated at D = 4 (∼ 1/(Q 2 ) 2 ) terms, to the sum of (experimental) inelastic and (parametrized) elastic BSR, the quality of the results turned out to be significantly worse in all the cases of the theoretical curves, something expected by the arguments presented at the end of Sec. II C. Namely, the elastic contribution is dominated by terms which behave at high values of Q 2 as ∼ 1/(Q 2 ) (D/2) where usually D ≥ 8, and these terms are not contained in the theoretical expressions for BSR which are usually OPE series truncated at 1/(Q 2 ) 2 .
The results of this work can be interpreted as an additional indication of the following important property: the evaluation of the (truncated) leading-twist contribution of spacelike low-Q 2 QCD observables such as inelastic BSR, in QCD variants 2δ and in particular 3δ AQCD [both have infrared finite and holomorphic coupling A(Q 2 )], appear to resum effectively a large part of the perturbative contribution of the observables, and leads to reduced extracted values of the higher-twist terms (D = 2, 4) in the truncated OPE. This property was noted earlier, for different observables, in Refs. [24,25,71]. In this context, it appears to be important that in 2δ and 3δ AQCD the coupling practically merges with the underlying pQCD coupling a(Q 2 ) at higher values of Q 2 Λ 2 QCD . This property is not shared by the (F)APT holomorphic coupling where the leading-twist series contains parts of the higher-twist contribution of as low dimensionality as D = 2. The extracted parameters in the higher-twist contribution, including those in the "massive" OPE, are especially reduced in 3δ AQCD. This suggests the possibility that the true higher-twist contribution is small, including the (sum of) terms of high dimension; and that the (truncated) OPE with 3δ AQCD leading-twist gives, through fitting, an extracted value which is a good approximation to this true value of the higher-twist contribution. Numerically, the significantly reduced extracted value in (truncated) OPE with 3δ AQCD is probably partly related with the fact that 3δ AQCD differs from both 2δ and (F)APT AQCD variants in that its coupling goes to zero in the deep infrared regime, A (3δ) (Q 2 ) ∼ Q 2 → 0. The latter property, we recall, is suggested by the large-volume lattice calculations of the dressing functions of the Landau-gauge gluon and ghost propagators at low Q 2 values.
The dependence of the coupling a(µ 2 ) ≡ α s (µ 2 )/π (where µ 2 > 0 means the spacelike region) is governed by the perturbative renormalization group equation (RGE) where in the mass-independent schemes the coefficients β 0 = (11−2N f /3)/4 and β 1 = (102−38N f /3)/16 are universal (scheme-independent), while the coefficients β j (or equivalently, c j ≡ β j /β 0 ) for j ≥ 2 are the (arbitrary) parameters which characterize the renormalization scheme. 11 The running of the coupling a(µ 2 ; c 2 , c 3 , . . .) with these scheme parameters is governed by the following relations (cf. App. A of Ref. [72], and App. A of Ref. [73]): When we use the relations (A1)-(A2) in the perturbation expansion (4) and account for the fact that D BSR (Q 2 ) is a (spacelike) observable and thus independent of the scale µ 2 (i.e., independent of k ≡ µ 2 /Q 2 ) and of the scheme parameters c j (j ≥ 2), we obtain the following expressions for the perturbation coefficients d j in terms of the general renormalization scale and scheme parameters (k, c 2 , c 3 ): Here we used the bar symbol to denote the choice of the scheme MS with the renormalization scale µ 2 = Q 2 (k = 1). We recall that k ≡ µ 2 /Q 2 is the renormalization scale parameter (0 < k ∼ 1), and thatc j ≡β j /β 0 (j ≥ 2) are the MS scheme parameters.
More specifically, we have where, as mentioned in Appendix A, c j ≡ β j /β 0 . When we invert these relations, we obtain The linearity of analytization, Eq. (B2), then gives us the analogs A n of the powers a n We note that in general A n (Q 2 ) = A(Q 2 ) n . The described construction (for n = 1, 2, 3, . . .) was performed in [34,35]. The above approach was extended in Ref. [36] to the case of general real index n = ν where Li −ν+1 (z) is the polylogarithm function of order −ν+1, and ρ A (σ) = ImA(Q 2 = −σ−i ) is the cut discontinuity (spectral) function of A. The coupling A ν can also be presented in an alternative form applicable in an extended region −1 < ν (cf. [36] for details). The expression A ν , the analog of the power a ν , was then obtained in the form with the coefficients k m (ν) given in Appendix A of Ref. [36]. Eqs. (B7) are a special case of Eq. (B9). The perturbation expansion of the type Eq. (4) in Sec. II A, for any spacelike obervable D(Q 2 ) in pQCD, can be reexpressed in terms of derivatives a n of Eq. (B3a) where we denoted a ≡ a(kQ 2 ; c 2 , . . .), a n ≡ a n (kQ 2 ; c 2 , . . .), and the coefficients d n ≡ d n (k; c 2 , . . . , c n ) of this "modified" perturbation expansion (mpt) can be obtained by using the RGE-relations (B6) and the coefficients d j ≡ d j (k; c 2 , . . . , c j ) are those of Eqs. (A3). The expressions in AQCD, corresponding to the perturbation expansions (B10) and (4), are then The expansion (B12b) is written again, in a more detailed form, in Eq. (5) in Sec. II A. Both expressions (B12) are equivalent, but in practice it is more economical to do numerical evaluations using the expression (B12a).
We have A n (Q 2 ) = A(Q 2 ) n only when A(Q 2 ) is a perturbative coupling, i.e., when it has no NP terms (∼ 1/(Q 2 ) m ). It is important not to use the power expansion in A for the evaluation of spacelike observables. Namely, if we used power expansion in A, the truncated series for D(Q 2 ) would have increasingly large (out of control) NP contributions when more power terms were included, and renormalization scale invariance would be increasingly violated, as emphasized in [74]. It turns out that in practice the sequence A n (Q 2 ) (n = 0, 1, 2, . . .) is, in a general holomorphic AQCD, a sequence with decreasing absolute values, at any finite Q 2 : | A n (Q 2 )| > | A n+1 (Q 2 )| > . . .. In pQCD ( a n (Q 2 )) this is in general not valid at low values |Q 2 | 1 GeV 2 .
The analog A (APT) ν (Q 2 ) of the power a(Q 2 ) ν is then constructed in complete analogy, by replacing ρ (pt) [13,15] We use in this work this form of (F)APT couplings A ν for ν = 1, 2, . . ., in MS scheme. Specifically, we apply the underlying MS pQCD coupling a(Q 2 ; MS) with N f = 3 to evaluate ρ (pt) n and thus A (APT) n (Q 2 ). The authors of [16,17] obtained explicit form of A (APT) ν (Q 2 ) at the one-loop level of the underlying pQCD coupling, and extended it to higher loop level by a perturbative approach [18,19]. This theory has the name Fractional Analytic Perturbation Theory (FAPT).
Numerical programs were constructed to calculate A (APT) ν (Q 2 ) up to four-loop level of the underlying MS pQCD coupling a, in [75,76] in Maple and/or Fortran, and in [23,77] in Mathematica.
We note that the general approach presented in Appendix B for calculation of A ν gives in the APT case (where ρ A (σ) = ρ (pt) 1 (σ), for σ > 0) approximately the same numerical results as the approach of Eq. (C3). If in the approach described in Appendix B we take into account in Eq. (B9) [or (B7) when ν = n is integer] a large number of terms A ν+m , 12 it turns out that the obtained A ν (Q 2 ) numerically converges to that in Eq. (C3). We wish to stress that the approach presented in Appendix B works for any AQCD, i.e., QCD with any holomorphic A(Q 2 ), while the approach Eq. (C3) is applicable only in the (F)APT case, i.e., when ρ A (σ) = ρ ; the high energy QCD is approximately reproduced with N f = 3 (F)APT when Λ 3 = 0.45 GeV (cf. also Ref. [15]); we use this value in our analysis, but we also comment on the case Λ 3 = 0.40 GeV.

2δ and 3δ AQCD
This type of QCD variants [22,24,25] with coupling A(Q 2 ) holomorphic in Q 2 ∈ C\(−∞, −M 2 thr ] are constructed on the idea that: (a) the coupling A(Q 2 ) at high |Q 2 | > 1 GeV 2 should practically coincide with the underlying pQCD 13 coupling a(Q 2 ); (b) and at moderate and low |Q 2 | 1GeV 2 the coupling should reproduce the well measured semihadronic τ -lepton decay physics and possibly some other experimental indicators. We note that (F)APT does not fulfill these requirements.
The condition (a) then implies that the spectral function ρ A (σ) ≡ ImA(−σ − i ) is at large σ > 1 GeV 2 (approximately) equal to the spectral function of the underlying pQCD, ρ a (σ) ≡ ρ Implicitly, we expect M 2 1 < M 2 2 < · · · < M 2 n < M 2 0 , where M 2 1 = M 2 thr. is the mentioned threshold scale (expected to be ∼ m 2 π ∼ 10 −2 GeV 2 ), and M 2 0 (∼ 1 GeV 2 ) can be called the pQCD-onset scale. Application of the Cauchy theorem then gives for the coupling the expression The parametrization with n delta functions in Eq. (C4) means that the part ∆A of the QCD coupling originating from the unknown low-σ part of the spectrum is parametrized by a near-diagonal Padé [n − 1/n](Q 2 ) approximant where P n−1 and Q n are polynomials of degree n − 1 and n in Q 2 , respectively. If ∆A is a Stieltjes function 14 , a theorem in the Padé theory [78] (cf. also [79]) states that there exists a sequence of Padé approximants [n − 1/n](Q 2 ) which converges to ∆A(Q 2 ) for any Q 2 ∈ C\(−∞, −M 2 thr ] when n increases (if ∆A is not Stieltjes, it is not known whether such a convergence is guaranteed).
The underlying pQCD coupling a(Q 2 ) is determined, to a rather high degree of accuracy, by the world average value a(Q 2 = M 2 Z ; MS) = 0.1185/π [45]. We use this value, and we RGE-evolve a(MS) to lower values of |Q 2 |, by using the four-loop RGE (A1) and three-loop quark mass threshold conditions [80] in MS, into the region of interest where N f = 3; subsequently, we change the coupling a(Q 2 ; MS) to the considered renormalization scheme to obtain the underlying coupling a(Q 2 ); cf. [25] for more details. Thus we also obtain the perturbative spectral function ρ (pt) 1 (σ) = Ima(Q 2 = −σ − i ), in the considered scheme and for N f = 3. At that point, the considered AQCD, Eqs. (C4)-(C5) has altogether (2n + 1) parameters, namely F j and M 2 j (j = 1, . . . , n) and the pQCD-onset scale M 2 0 . These are to be fixed by the conditions (a) and (b) for A(Q 2 ) at high and low |Q 2 | mentioned at the beginning of this Section C 2. The condition (a) is implemented by requiring a large index value for the difference Eq. (B1); in our considered cases of 2δ and 3δ coupling we took N = 5 This represents altogether four conditions (for each increase of N by one, from N = 1, there is one condition). When we take n = 2, i.e., two delta functions in the spectral function (C4), we have five parameters to determine; therefore one additional condition is needed. This condition will be the reproduction of the measured values of the quantity r (D=0) τ , the semihadronic strangeless τ decay rate ratio (the leading-twist part, and with mass effects subtracted). 15 Its experimental value is approximately in the range 0.201 ± 0.002. Its theoretical expression can be represented as a weighted countour integral of the (massless) Adler function 16 We refer to [22,24,25] for details. In the n = 2 case (2δ AQCD), we still have the freedom of choosing the renormalization scheme. We took it as the Lambert scheme (c j = c j−1 2 /c j−2 1 , for j ≥ 3), with c 2 = −4.9 [23]. It is possible to vary the value of c 2 , but when it is different by several units from this value, either the pQCD-onset scale M 0 becomes appreciably higher than ≈ 1 GeV, or the value A(0) becomes larger than one, cf. Table 2 of [23].
When we take n = 3, i.e., three delta functions in the spectral function (C4), there are two additional parameters to be fixed. These two parameters are fixed by the condition A(Q 2 ) ∼ Q 2 when Q 2 → 0 and the local maximum of A(Q 2 ) achieved at Q 2 ≈ 0.135 GeV 2 . These conditions are motivated by the results of lattice calculations of the gluon and ghost dressing functions in the Landau gauge at low positive Q 2 [53] (cf. also [54]). 17 In our couplings we use throughout N f = 3, which makes them applicable in the regions |Q 2 | < 3 GeV 2 . The lattice calculations were performed in the MiniMOM renormalization scheme (MM) [83][84][85][86]. Our coupling A(Q 2 ) was constructed in the N f = 3 MM scheme, but rescaled from the MM scale convention (Λ MM ) to the usual scale convention (Λ MS ), the latter representing what we call Lambert MM (LMM) scheme. In [24] we used the three-loop LMM, and in [25] the four-loop LMM scheme. In the present work we use the latter (four-loop LMM) scheme, i.e., the coupling from [25].
In Table VIII we specify the parameters of the 2δ AQCD [22,23] and 3δ AQCD coupling [25], used in the present work. The values of the parameters of the 2δ coupling, with c 2 = −4.9, are slightly different from the corresponding values in [23] (Table 2 there Other QCD running couplings A(Q 2 ) without Landau singularities have been constructed in the literature, some of them having at Q 2 = 0 zero value A(0) = 0 [87][88][89], finite nonzero value (for reviews, cf. [27,90]), or infinite value [91]. The construction of such couplings is mostly based on the dispersive approaches similar to the ones described in this Appendix; however, such kind of dispersive approaches can also be applied to entire physical observables, cf. [14,89,[92][93][94][95][96][97][98]. Yet another approach is the Light Front Holographic (LFH) QCD [65,66], where the coupling A ((LFH )(Q 2 ) has the form exp(−Q 2 /2κ 2 ) at low positive Q 2 , cf. Eq. (19), and can be extended to higher positive Q 2 by matching of A(Q 2 ) and dA(Q 2 )/dQ 2 at a matching scale Q 0 ∼ 1 GeV 2 to pQCD [67]. The experimental data in the fit procedure are Γ p−n 1 (Q 2 j ; inel.) exp. . The leading-twist part of the (theoretical) OPE expression (10) has no fit parameters in our approach, apart from the renormalization scale parameter k ≡ µ 2 /Q 2 . The latter parameter will be kept fixed at its value obtained by the approach of minimization of χ 2 . We will use the following notation: In order to estimate the uncertainties of the extracted parameters f p−n

2
(1) and µ 6 , we will consider for simplicity that not just µ 6 but also f p−n 2 does not run (i.e., that µ 4 does not run) with Q 2 . Further, we will make the following approximations: the statistical uncertainties at different points will be considered as completely uncorrelated; the systematic uncertainties at different points will be considered as completely correlated within one experiment, and completely uncorrelated between two different experiments. The systematic uncertainties δΓ p−n 1 (Q 2 j ; inel.) sys are in general significantly larger than the statistical uncertainties δΓ p−n 1 (Q 2 j ; inel.) stat , cf. Fig. 8. 1. µ6 = 0 case First we will consider the case when µ 6 = 0 (and f p−n 2 is not running with Q 2 ). We will estimate the uncertainty of f p−n 2 due to the statistical uncertainties σ j,stat at Q 2 j 's and then due to the systematic uncertainties σ j,sys . The renormalization scheme parameter k will be kept fixed during these estimations.
In this approach, the considered χ 2 is where we will denote i.e., the squared statistical uncertainty at the measured point Q 2 j . The minimization of χ 2 then gives the extracted value of µ 4 where we will use throughout the notation for the unnormalized "average" Since the deviations of y j ≡ y(z j ) due to statistical uncertainties are considered independent at different z j , the rule of the sum of the squares of standard deviations is valid. It is then straightforward to deduce the square of the standard deviation forμ 4 .
Any of the two relations (D6)-(D7) can be used to evaluate the uncertainty (δμ 4 ) stat ≡ σ(μ 4 ) stat . We used the relation (D7), which gives us in practice somewhat higher values of the uncertainty. 18 Further, using the relation (6) at Q 2 = 1 GeV 2 , this then gives us the uncertainty (δf 2 ) stat The systematic uncertainties were estimated in the following way. For simplicity we consider for this only two experiments, namely EG1b [3] (experiment "1") and the newer JLAB EG1-DVCS results [4] (experiment"2"), in the interval 0 < Q 2 < 3 GeV 2 . These two experiments represent most of the available data points, and each of them covers most of the mentioned Q 2 -interval. In each experiment, the systematic uncertainties at different points will be considered as strongly correlated, in the sense that in each experiment we will estimate the systematic uncertainty of f 2 as the one obtained by averaging the deviations where N = 4 or N = 6. For explanations of the notation 'UP' and 'DO' we refer to the previous Sec.D 2. In our approximation, we consider that the systematic uncertainties in the experiments 1 and 2 are mutually independent; this, in conjunction with the relations (D21), then gives (δμ 4 ) sys ≡ σ(μ 4 ) sys = σ 2 ( µ The systematic uncertainty of the extracted value f p−n 2 (1GeV 2 ) is obtained then from Eq. (D26a) via the relation (6) (δf 2 ) sys ≡ σ(f 2 ) sys = 9 4M 2 N σ(μ 4 ) sys .
(D27) and m c ≈ 1.59 GeV is the pole mass. This function appears in the NLO coefficient (at a 2 ) when the perturbation expansion (4) is reexpressed in terms of the N f = 4 coupling a N f =4 In this limit (ξ c → ∞) this then gives in the expansion (E1) the massless expression for the NLO coefficient As noted, the convergence to the N f = 4 case (four massless quarks) is rather slow in BSR. Specifically, 3 + C mass.,(2) pBJ (Q 2 /m 2 c ) ≈ 3.13, 3.36, 3.73, 3.83 for Q 2 = 5, 10, 50, 100 GeV 2 , respectively. This indicates that at Q 2 = 5 GeV 2 (the highest considered experimental Q 2 in BSR) we are still rather far away from the N f = 4 case. Therefore, it appears reasonable to use N f = 3 (i.e., decoupled charm quark) in the polarized BSR for the interval 0 < Q 2 < 5 GeV 2 considered in the present work.