A novel determination of non-perturbative contributions to Bjorken sum rule
aa r X i v : . [ h e p - ph ] F e b A novel determination of non-perturbative contributions to Bjorken sum rule
Qing Yu, ∗ Xing-Gang Wu, † Hua Zhou, ‡ and Xu-Dong Huang § Department of Physics, Chongqing University, Chongqing 401331, People’s Republic of China (Dated: February 26, 2021)Based on the operator product expansion, the perturbative and nonperturbative contributionsto the polarized Bjorken sum rule (BSR) can be separated conveniently, and the nonperturbativeone can be fitted via a proper comparison with the experimental data. In the present paper,we first give a detailed study on the pQCD corrections to the leading-twist part of BSR, andthen give a novel fit of the non-perturbative high-twist contributions by comparing with the JLabdata. Previous pQCD corrections to the leading-twist part derived under conventional scale-settingapproach up to O ( α s )-level still show strong renormalization scale dependence. The principle ofmaximum conformality (PMC) provides a systematic way to eliminate conventional renormalizationscale-setting ambiguity by determining the accurate α s -running behavior of the process with thehelp of renormalization group equation. Our calculation confirms the PMC prediction satisfies thestandard renormalization group invariance, e.g. its fixed-order prediction does scheme-and-scaleindependent. In low Q -region, the effective momentum of the process is small and to have areliable prediction, we adopt four low-energy α s models to do the analysis, i.e. the model basedon the analytic perturbative theory (APT), the Webber model (WEB), the massive pQCD model(MPT) and the model under continuum QCD theory (CON). Our predictions show that even thoughthe high-twist terms are generally power suppressed in high Q -region, they shall have sizablecontributions in low and intermediate Q domain. Based on the more accurate scheme-and-scaleindependent pQCD prediction, our newly fitted results for the high-twist corrections at Q = 1 GeV are, f p − n | APT = − . ± . f p − n | WEB = − . ± . f p − n | MPT = − . ± .
013 and f p − n | CON = − . ± . µ | APT = 0 . ± . µ | WEB = 0 . ± . µ | MPT = 0 . ± . µ | CON = 0 . ± . I. INTRODUCTION
The Bjorken Sum Rule (BSR) [1, 2], which describesthe polarized spin structure of nucleon, has been mea-sured via polarized deep inelastic scattering (DIS) byvarious experimental collaborations [3–31]. Using the op-erator product expansion (OPE), the BSR of the spinstructure function can be calculated by separating theperturbative contribution of the matrix elements of lo-cal product operators from its non-perturbative contri-butions [1, 2], e.g.Γ p − n ( Q ) = Z dx [ g p ( x, Q ) − g n ( x, Q )]= g A (cid:2) − E ns ( Q ) (cid:3) + ∞ X i =2 µ p − n i ( Q )( Q ) i − , (1)where g p,n ( x, Q ) is the spin-dependent proton or neu-tron structure function with Bjorken scaling variable x ,and g A is the nucleon axial charge. The BSR relatesthe difference of the proton and the neutron structurefunctions Γ p and Γ n , and only the flavor non-singletquark operators appear in perturbative part, resulting asthe perturbative non-singlet leading-twist contributions ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] E ns ( Q ). The non-perturbative contribution is gener-ally power suppressed in comparison to the leading-twistterms, which has been written as a power series over1 /Q . Contributions from the high-twist terms could besizable in low and intermediate Q regions, and then theBSR provides a good platform for testing the perturba-tive and non-perturbative QCD contributions.Analyses of E ns ( Q ) under the MS-scheme have beengiven in the literature, such as Refs.[32–35]. Additionaltreatment on extending the pQCD prediction to low Q -region has been done by using low-energy models for thestrong coupling constant ( α s ) such as the analytic per-turbation theory (APT), the “massive analytic pQCDtheory” (MPT), the 2 δ - or 3 δ -analytic QCD variants [36–41]. In all those treatments, there are large renormaliza-tion scale ( µ r ) dependence for the perturbative part dueto the using of “guessed” µ r ; that is, in those analyses,the central (“optimal”) value of E ns ( Q ) is usually de-rived by setting µ r = Q , and then by varying it withinan arbitrary range such as [ Q/ , Q ] to estimate its un-certainty. Such guessing choice breaks the renormaliza-tion group invariance [42, 43] and leads to conventionalrenormalization scale-and-scheme ambiguities due to themismatching of the perturbative coefficients and the α s at each order. In the literature, the principle of maxi-mum conformality (PMC) [44–47] has been suggested toeliminate such renormalization scale-and-scheme ambi-guities. It is well known that the α s -running behavior isgoverned by the renormalization group equation (RGE).The existence of the { β i } -terms emerged in the pertur-bative series is thus helpful for fixing exact α s -value ofthe pQCD approximant of a physical observable. Andinstead of choosing an optimal µ r , the PMC fixes thecorrect magnitude of α s by using RGE, whose argumentis called as the PMC scale, which is independent to anychoice of µ r . The PMC prediction is scale-and-schemeindependent, more detail and applications of the PMCcan be found in the reviews [48–50].To achieve a reliable prediction for the BSR high-twistcontributions, it is important to have an accurate pQCDprediction on E ns ( Q ). In the present paper, we shall firstadopt the PMC single-scale approach [51] to deal withthe perturbative part of the BSR, and then give a newdetermination of the non-perturbative high-twist contri-butions by comparing with the JLab data. The PMCsinglet-scale approach follows the same idea of the orig-inal multi-scale approach [44–47], which determines anoverall effective momentum flow of the process by usingthe RGE, whose magnitude corresponds to the weightedaverage of the multi-scales of the multi-scale approach ateach order. It has also been demonstrated that the pre-diction under the PMC singlet-scale approach is scheme-and-scale independent up to any fixed order [52]. Thoughdifferent from conventional scale ambiguity, there is resid-ual scale dependence for fixed-order prediction due to un-known perturbative terms [53]. Such residual scale de-pendence can be greatly suppressed due to both α s -powersuppression and exponential suppression. A detailed dis-cussion on the residual scale dependence can be found inthe recent review [50].The remaining parts of the paper are organized as fol-lows. In Sec.II, we present the calculation technology forthe polarized Bjorken sum rule Γ p − n . The PMC treat-ment of the pQCD contributions to the leading-twist partand the non-perturbative high-twist contributions shallbe given. In Sec.III, we give the numerical results anddiscussions. Sec.V is reserved for a summary. II. CALCULATION TECHNOLOGY
In large Q -region, contributions from the leading-twist terms are dominant and those of the non-perturbative high-twist terms are generally power sup-pressed. In low and intermediate Q -region, contribu-tions from the high-twist terms may have large contri-butions. In the following, we shall analyze the pQCDcontributions to the leading-twist terms by using thePMC single-scale approach, and then give an estima-tion of the contributions from the non-perturbative high-twist terms. In low Q -region, the low-energy α s modelsshould be used; and for clarity, we shall adopt four low-energy α s models to do our discussion. A. Perturbative series of the leading-twist terms
The perturbative expansion over α s for the hard partof the leading-twist terms E ns ( Q ) has been calculated up to next-to-next-to-next-to leading order (N LO). Werewrite E ns ( Q ) in the following form, E ns ( Q ) = r , a ( µ r ) + ( r , + β r , ) a ( µ r )+ ( r , + β r , + 2 β r , + β r , ) a ( µ r )+ ( r , + β r , + 2 β r , + 52 β β r , + 3 β r , + 3 β r , + β r , ) a ( µ r ) + · · · , (2)where a ( µ r ) = α s ( µ r ) /π . We have implicitly writtenthe perturbative coefficients as the { β i } -series by us-ing general QCD degeneracy relations among differentorders [54]. The coefficients r i ∈ [1 , ,j ∈ [0 , can be ob-tained by transforming the original n f -series given inRefs.[55, 56]. To apply the PMC, we have also dividedthe coefficients r i,j into the conformal ones r i, and non-conformal ones r i,j =0 by using the standard PMC pro-cedures given in Refs.[46, 47]. Generally, the coefficients r i,j =0 are functions of the logarithm ln( µ r /Q ). If set-ting µ r = Q , all those log-terms becomes zero, leadingto a renormalon-free more convergent pQCD series; thisexplains why people usually choose µ r = Q as the opti-mal scale for conventional scale-setting approach. Thosecoefficients can be reexpressed as r i,j = j X k =0 C kj ln k ( µ r /Q )ˆ r i − k,j − k , (3)where the combination coefficients C kj = j ! /k !( j − k )!,and the coefficients ˆ r i,j = r i,j | µ r = Q . For convenience, weput the reduced coefficients ˆ r i,j in the Appendix A.The RGE, or the β -function, is defined as β ( a ( µ r )) = − ∞ X i =0 β i a i +2 ( µ r ) . (4)Following the decoupling theorem [57], the first two { β i ≥ } -functions β and β are scheme-independent, andwe have β = (11 − n f ) and β = (102 − n f ) forthe SU C (3)-color group with n f being the active flavornumbers of the process. The scheme dependent { β i ≥ } -functions have been calculated up to five-loop level un-der the MS-scheme [58–66]. A collection of all the known { β i } -functions can be found in Ref.[48].In Eq.(2), the { β i } -terms at each perturbative ordergovern the correct α s -running behavior, which inverselycan be used to determine the effective magnitude of α s . Practically, by requiring all the RGE-involved non-conformal { β i } -terms to be zero, one can achieve an over-all effective α s and hence the PMC scale Q ⋆ , and thenthe resultant pQCD series becomes the following scheme-independent conformal series: E ns | PMC ( Q ) = X i ≥ ˆ r i, a i ( Q ⋆ ) , (5)where the PMC scale Q ⋆ can be fixed up to next-to-next-to-leading-log (NNLL) accuracy by using the N LOperturbative series, e.g.ln Q ⋆ Q = T + T α s ( Q ) π + T α s ( Q ) π , (6) where T = − ˆ r , ˆ r , , T = 2(ˆ r , ˆ r , − ˆ r , ˆ r , )ˆ r , + ˆ r , − ˆ r , ˆ r , ˆ r , β , (7) T = 4(ˆ r , ˆ r , ˆ r , − ˆ r , ˆ r , ) + 3(ˆ r , ˆ r , ˆ r , − ˆ r , ˆ r , )ˆ r , + 3(ˆ r , − ˆ r , ˆ r , )2ˆ r , β − r , ˆ r , − r , ˆ r , ˆ r , − r , ˆ r , ˆ r , + 3ˆ r , ˆ r , ˆ r , β + 2ˆ r , ˆ r , ˆ r , − ˆ r , ˆ r , − ˆ r , ˆ r , β . (8)Those equations show the PMC scale Q ⋆ is exactly freeof µ r , together with the µ r -independent conformal coef-ficients ˆ r i, , the PMC prediction is exactly independentto any choice of µ r . Thus the conventional scale-settingambiguity can be eliminated at any fixed-order by ap-plying the PMC [52]. As a byproduct, due to the elim-ination of divergent renormalon terms in the resultantPMC perturbative series (5), the pQCD convergence canbe naturally improved. Those properties greatly improvethe precision of the pQCD theory.For a perturbative theory, it is important to have a re-liable way to estimate the magnitude of the uncalculatedhigher-order terms. The scale-invariant and scheme-invariant PMC conformal series, which is also more con-vergent than the conventional series, is quite suitable forsuch purpose. A way of using the PMC series togetherwith the Pad´ e approximation approach (PAA) [67–69]has been suggested in Ref.[70]. Some successful appli-cations of this method can be found in Refs.[71–74]. Weshall adopt this method to estimate the magnitude of theunknown O ( α s )-terms of E ns ( Q ), and in the following,we give a brief introduction of PAA.The PAA offers a feasible conjecture that yields the( n + 1) th -order coefficient by using a given n th -order per-turbative series. For the purpose, people usually adoptsa fractional function as the generating function. More ex-plicitly, the [ N/M ]-type generating function of a pQCDapproximant ρ n ( Q ) = n P i =1 ˆ r i, a i is defined as ρ [ N/M ] n ( Q ) = a × b + b a + · · · + b N a N c a + · · · + c M a M (9)= n X i =1 C i a i + C n +1 a n +1 + · · · , (10)where M ≥ N + M + 1 = n . The perturbativecoefficients C i in Eq.(10) can be expressed by the knowncoefficients b i ∈ [0 ,N ] and c j ∈ [1 ,M ] . Inversely, if we haveknown the coefficients C i ’s up to n th -order level, one candetermine the coefficients b i ∈ [0 ,N ] and c j ∈ [1 ,M ] , and then achieve a prediction for the uncalculated ( n + 1) th -ordercoefficient C n +1 .At the present, the leading-twist term E ns | PMC hasbeen known up to N LO-level, and the four coefficientsare known, C i = ˆ r i, for i ∈ [1 , LO-coefficient becomesˆ r , = ˆ r , − r , ˆ r , ˆ r , + ˆ r , ˆ r , + 2ˆ r , ˆ r , ˆ r , ˆ r , , (11)where the [0 /n − B. Contributions from the non-perturbativehigh-twist terms
The non-perturbative contributions to the BSR can beexpanded in 1 /Q -power series as Eq.(1). The O ( Q − )-term µ p − n can be written as [75–77] µ p − n = M a p − n + 4 d p − n + 4 f p − n ) , (12)where M ≈ .
94 GeV is the nucleon mass. The leading-twist target mass correction a p − n can be calculated byusing the leading-twist part of g p − n , which is kinemati-cally of high-twist [78] and its magnitude at Q = 1 GeV is 0 . ± .
010 [35]. The twist-3 matrix element d p − n isgiven by d p − n = Z dxx (2 g p − n + 3 g p − n ) , (13)whose magnitude at Q = 1 GeV is 0 . ± . x -weighted moment of the structure function in ordersof M /Q , thus a p − n and d p − n change logarithmically,and we shall fix their values to be the above ones at Q = 1 GeV . Then the remaining undetermined termin µ p − n is f p − n . The twist-4 term f p − n , which is relatedto the color electric and magnetic polarizabilities of nu-cleon, plays a pivotal role in phenomenological studies ofthe high-twist contributions. f p − n is sensitive to Q andits Q -evolution satisfies [76, 77] f p − n ( Q ) = f p − n (1) (cid:18) a ( Q ) a (1) (cid:19) γ / β , (14)where γ / β = 32 /
81 with n f = 3. The magnitude of f p − n (1) shall be fit by comparing with the data. More-over, it has been argued that the O ( Q − )-term µ p − n mayalso have sizable contribution, so we take µ /Q -terminto consideration to have a better fit of the data. C. The strong coupling constant α s The α s -running behavior in perturbative region is gov-erned by the RGE (4). Its solution can be written asan expansion over the inverse powers of the logarithm L = ln µ r / Λ ; and up to four-loop level, we have [79] α s ( µ r ) = πβ L (cid:26) − β β ln LL + 1 β L (cid:20) β β (ln L − ln L −
1) + β β (cid:21) + 1 β L (cid:20) β β ( − ln L + 52 ln L +2 ln L −
12 ) − β β β ln L + β β (cid:21)(cid:27) , (15)where Λ is the scheme-dependent asymptotic scale, whichcould be fixed by matching the measured value of α s ata reference scale such as M Z or m τ to its predicted valueunder a specific scheme.In infrared region, when the scale is close to Λ or evensmaller, α s becomes large whose magnitude cannot bewell described by the RGE. To make the QCD predictionmore reliable, we shall adopt four low-energy models forthe α s to do our calculation.The first low-energy model is based on the analyticalperturbation theory (APT) [80, 81], and we call it as theAPT model. In APT model, its strong coupling constant α APT s is described by applying the perturbation theorydirectly to the spectral function, which takes the follow-ing form, α APT s ( µ ) = πβ (cid:18) − y (cid:19) , (16)where µ is the energy scale, y = µ / Λ withΛ = µ exp[ − φ ( β α s ( µ ) /π )] , (17)where φ ( z ) satisfies 1 /φ ( z ) + 1 / (1 − exp[ φ ( z )]) = z . Itsfreezing value is close to α APT s (10 − ) /π ≈ .
43. The second low-energy model is an alteration ofEq.(16), we call it as the WEB model [82], which is sug-gested to suppress the nonperturbative power correctionsof the APT model, and it takes the following form α WEB s ( µ ) = πβ (cid:20) b (1 − y)(1 + b ) ( 1 + c y + c ) p (cid:21) , (18)where these phenomenological parameters b = 1 / p = c = 4. The obtained corresponding approximatefreezing value ∼ α WEB (10 − ) /π ≈ . m gl = √ ξ Λ as the in-frared regulator, and we call it as the MPT model. Ittakes the following form α MPT s ( µ ) = a cr (cid:26) a cr β π ln µ m gl ! + a cr β πβ × ln " a cr β π ln µ m gl ! + ... (cid:27) − , (19)whose freezing value at the origin satisfies a cr = π/ ( β ln ξ ). Under the Landau gauge, we have a cr | ξ =10 ± = 0 . ∓ .
05, which leads to the freezing point α MPT s (0) /π = 0 . − . . .The fourth low-energy model is based on the contin-uum theory [85] and we call it as the CON model, wherethe exchanging gluons with effective dynamical mass m g is adopted and the non-perturbative dynamics of gluonsis governed by the corresponding Schwinger-Dyson equa-tion. It takes the following from α CON s ( µ ) = πβ ln (cid:16) M g + µ Λ (cid:17) , (20)whose M g = m g [ln(y + 4 m g / Λ ) / ln(4 m g / Λ )] − / and m g = 500 ±
200 MeV [85, 86], which leads to the freezingpoint α CON s (0) /π = 0 . − . . . III. NUMERICAL RESULTS
To do the numerical analysis, we take the nucleon axialcharge ratio g A = 1 . ± . α s -value at the ref-erence point such as α MS s ( m τ ) = 0 . ± .
016 [87], whichgives Λ MS | n f =3 = 0 . +0 . − . GeV by using the four-loopRGE. Using the relation (17), we obtain Λ
APT | n f =3 =0 . +0 . − . GeV. In Fig. 1, we present the typical run-ning behaviors of α s /π under four low-energy models,where the parameters are set to be ξ = 10 for MPT and m g = 700 MeV for CON, respectively. The α s -runningbehavior derived from the RGE under MS-scheme isgiven as a comparison. Fig. 1 shows the importance of theusing of low-energy models in the region of small energy-scale. Using the criteria suggested in Ref.[88] for the ana-lytic matching of α s in perturbative and nonperturbative FIG. 1: Typical α s -running behavior in low-energy scales forfour typical low-energy models, APT, WEB, MPT, and CON,respectively. The α s -running behavior derived from RGE un-der MS-scheme is given as a comparison. regimes, we obtain the transition scales ( Q ) for variouslow-energy models, which are ∼ .
77 GeV, ∼ .
78 GeV, ∼ .
78 GeV, ∼ .
19 GeV for APT, WEB, MPT andCON models, respectively. As a subtle point, becausethe transition scales Q for the cases of WEB and MPTare slightly bigger than m τ , and for self-consistency, weuse the low-energy α WEB / MPT s ( m τ ) = 0 . ± .
016 tofix Λ, which is 0 . ± .
022 GeV or 0 . +0 . − . GeV,respecitvely.For later convenience, in the following discussions, wesimply use α MS s to stand for the case of using MS-scheme α s in all Q -region, α APT s to stand for the case of usingAPT model in low-energy region ( Q < Q , as mentionedabove, Q is different for different low-energy model) andMS-scheme α s in large Q -region, α WEB s to stand for thecase of using WEB model in low-energy model and MS-scheme α s in large Q -region, α MPT s to stand for the caseof using MPT model in low-energy model and MS-scheme α s in large Q -region, and α CON s to stand for the case ofusing CON model in low-energy model and MS-scheme α s in large Q -region. A. Perturbative contributions to the leading-twistpart of BSR up to N LO level The perturbative contributions to the leading-twistpart E ns ( Q ) has been known up to N LO. Under con-ventional scale-setting approach, the pQCD series is scaledependent, and by setting µ r = Q , we obtain E ns ( Q ) | Conv . = a ( Q ) + 3 . a ( Q ) + 20 . a ( Q )+175 . a ( Q ) . (21)On the other hand, the pQCD series becomes scale in-variant by applying the PMC, and we obtain E ns ( Q ) | PMC = a ( Q ⋆ ) + 1 . a ( Q ⋆ ) + 0 . a ( Q ⋆ )+0 . a ( Q ⋆ ) (22) for any choice of renormalization scale, where Q ⋆ is ofperturbative nature, which can be determined up toNNLL accuracyln Q ∗ Q = − . − . a ( Q ) − . a ( Q ) . (23)One may observe that the perturbative coefficients inPMC series (22) are much smaller than those of con-ventional series (21), especially for those of high-orders,which are due to the elimination of divergent renormalonterms as n ! β n a ns . This indicates that a much more con-vergent perturbative series can be achieved by applyingthe PMC. At the same time, the PMC scale Q ⋆ also showsa fast convergent at high Q -range, e.g. the relative abso-lute values of the LL, the NLL and the NNLL terms are1: 0.064 : 0.030 for Q = 100 GeV. Thus the residual scaledependence due to unknown even higher-order terms canbe greatly suppressed.Using the convergent PMC perturbative series, one canobtain a reliable prediction of unknown O ( a )-term byusing the PAA, e.g. by using Eq.(11), we obtain E ns ( Q ) | N LOPAA = 2 . a ( Q ⋆ ) . (24)We present the predicted leading-twist part of the spinstructure function Γ p − n ( Q ) under four low-energy mod-els in Fig. 2, where the results under conventional andPMC scale-setting approaches are presented. The ex-perimental data are from SLAC [4–7, 10], DESY [20–24], CREN [25–27] and JLab [32, 34, 35]. The PMCpredictions are independent to any choice of µ r , andthe shaded band shows the conventional renormalizationscale uncertainty by varying µ r ∈ [ Q/ , Q ]. Under con-ventional scale-setting approach, the spin structure func-tion Γ p − n ( Q ) shows large scale dependence, especiallyin low-energy region. In low-energy region, the resultsby using the IR-fixed couplings are much more reliable.And since couplings behaves differently in low-energy re-gion, the spin structure function Γ p − n ( Q ) behaves quitedifferently for Q →
0. When the energy scale is largeenough, such as
Q > . − . B. Analysis of high-twist contributions undervarious low-energy models
Following the discussions of Sec.II.B, we need to fittwo parameters, f p − n (1 GeV ) and µ , so as to deter-mine the high-twist contributions. We adopt the mostrecent data listed in Refs.[34, 35] to do the fitting, whosemomentum transfer lies in the range of 0 .
054 GeV ≤ (a)Leading-twist contributions using α APT s . (b)Leading-twist contributions using α WEB s . (c)Leading-twist contributions using α MPT s . (d)Leading-twist contributions using α CON s . FIG. 2: Perturbative leading-twist contributions to the spin structure function Γ p − n ( Q ) up to N LO versus momentum Q,under four α s models: (a) the APT model; (b) The WEB model; (c) the MPT model; and (d) the CON model. The solid line isfor conventional scale setting approach with µ r = Q and the shaded band shows its scale uncertainty by varying µ r ∈ [ Q/ , Q ].The dot-dashed line is the prediction Γ p − n ( Q ) up to N LO for PMC scale-setting approach, which is free of renormalizationscale dependence. α s models f p − n (1) µ χ /d.o.fµ r = Q/ − . ± . ± .
013 0 . ± . ± .
000 149APT | Conv µ r = Q − . ± . ± .
013 0 . ± . ± .
000 62 µ r = 2 Q − . ± . ± .
013 0 . ± . ± .
000 117APT | PMC µ r ∈ [ Q/ , Q ] − . ± . ± .
013 0 . ± . ± .
000 62 µ r = Q/ − . ± . ± .
013 0 . ± . ± .
000 193WEB | Conv µ r = Q − . ± . ± .
013 0 . ± . ± .
000 168 µ r = 2 Q − . ± . ± .
013 0 . ± . ± .
000 160WEB | PMC µ r ∈ [ Q/ , Q ] − . ± . ± .
013 0 . ± . ± .
000 45 µ r = Q/ − . ± . ± .
013 0 . ± . ± .
000 151MPT | Conv µ r = Q − . ± . ± .
013 0 . ± . ± .
000 56 µ r = 2 Q − . ± . ± .
013 0 . ± . ± .
000 126MPT | PMC µ r ∈ [ Q/ , Q ] − . ± . ± .
013 0 . ± . ± .
000 50 µ r = Q/ − . ± . ± .
013 0 . ± . ± .
000 125CON | Conv µ r = Q − . ± . ± .
013 0 . ± . ± .
000 60 µ r = 2 Q − . ± . ± .
013 0 . ± . ± .
000 138CON | PMC µ r ∈ [ Q/ , Q ] − . ± . ± .
013 0 . ± . ± .
000 49TABLE I: The fitted parameters f p − n ( Q = 1 GeV ) and µ and their corresponding quality of fit χ /d.o.f under four α s models before and after applying the PMC, where the first and the second errors are caused by the statistical and systematicerrors of the data [34, 35]. The twist-6 coefficient µ is almost independent to the choices of statistical and systematic errors. Q ≤ .
739 GeV . We adopt the APT, WEB, MPT, andthe CON couplings in doing the fitting. The quality of fit is measured by the parameter of χ /d.o.f , e.g. χ /d.o.f = 1 N − d N X j =1 (Γ p − n , the . ( Q j ) − Γ p − n , exp . ( Q j )) σ j, stat . , (25)where the symbol “ d.o.f ” (short notation of the degree offreedom) is equal to N − d with N = 31 being the numberof data points and d = 2 being the number of fitted pa-rameters, “the.” stands for theoretical prediction, “exp.”stands for measured value, and “ σ j, stat . ” is the statisticalerror at each point Q j . Comparing theoretical predictionΓ p − n , the . ( Q j ) with the measured value Γ p − n , exp . ( Q j ) at all thedata points Q j ∈ [1 ,N ] , we can derive the preferable f p − n and µ by requiring them to achieve the minimum valueof χ /d.o.f . To do the fitting, we also take into accountthe systematic error σ j, sys . at each point Q j , which hassizable contributions to the fitted values of f p − n and µ .For convenience, we put the detailed calculation technol-ogy in Appendix B.Our results for the two parameters f p − n (1 GeV ) and µ are presented in Table. I. The right-most columnshows the smallest χ /d.o.f for the predictions beforeand after applying the PMC under four α s models. Themagnitudes of those two parameters are small, whichagree with the usual consideration that at large Q -region, the high-twist terms are power suppressed and arenegligible. However in low Q -region, they will have siz-able contributions; especially f p − n (1 GeV ) is importantfor a reliable theoretical prediction on Γ p − n , the . ( Q ) in low Q -region. Table. I shows that the fitted parameters un-der conventional scale-setting approach have strong scaledependence, whose quality of fit χ /d.o.f varies from tensto hundreds, and the optimal fit are achieved for the caseof µ r ∼ Q . This, together with a better pQCD con-vergence due to the elimination of divergent log-termsln µ r /Q , in some sense explain why µ r = Q is usuallytaken as the preferable renormalization scale for conven-tional scale-setting approach. On the other hand, thefitted parameters for the PMC scale-setting approach isindependent for any choice of renormalization scale, thusa more reliable and accurate prediction is achieved. -0.3-0.2-0.100.10.2 FIG. 3: The twist-4 coefficient f p − n (1 GeV ) obtained fromthe PMC predictions under four α s low-energy models, inwhich the predictions using JLab data [32, 34, 35], the QCDsum rule predictions [90, 91], and the predictions using themodel of the instanton-based QCD vacuum [92, 93] and theBag model prediction [75] are also presented. We present a comparison of the twist-4 coefficient f p − n (1 GeV ) under various approaches [32, 34, 35, 75,90–93] in Fig. 3. The results of Refs.[32, 34, 35] are fittedby using conventional pQCD series for the leading-twistpart with fixing µ r = Q and the JLab data within differ-ent ranges, 0 . < Q j <
10 GeV [32], 0 .
66 GeV 10 GeV [34] and 0 . 84 GeV < Q j < 10 GeV [35].By using f p − n , we can evaluate the color polarizability, χ p − nE = (2 d p − n + f p − n ) and χ p − nB = (4 d p − n − f p − n ),which describes the response of the color magnetic andelectric fields to the spin of the nucleon [89, 90]. Usingthe PMC predictions for the hard-part of the leading-twist contributions, we obtain χ p − nB | APT = 0 . ± . , (26) χ p − nB | WEB = 0 . ± . , (27) χ p − nB | MPT = 0 . ± . , (28) χ p − nB | CON = 0 . ± . , (29) χ p − nE | APT = − . ± . , (30) χ p − nE | WEB = − . ± . , (31) χ p − nE | MPT = − . ± . , (32) χ p − nE | CON = − . ± . , (33)where the errors are squared average of those from∆ d p − n = ± . f p − n for the four low-energy α s models (e.g. Table. I).We present the prediction of Γ p − n ( Q ) with bothleading-twist and high-twist contributions in Fig. 4.Comparing with Fig. 2, Fig. 4 shows that a more rea-sonable prediction can be achieved by including high-twist contributions. Under conventional scale-setting ap-proach, the large scale dependence for the leading-twistprediction of Γ p − n , Conv . ( Q ) can be greatly suppressed byincluding high-twist terms due to the cancellation of scaledependence among different twist-terms. Under PMCscale-setting approach, the scale-invariant Γ p − n , PMC ( Q )under APT, MPT and CON α s models are close in shape,which as shown by Table. I also have close quality of fit χ /d.o.f ; while the PMC prediction under WEB modelis slightly different from those of other α s models.As a final remark, to improve the quality of fit, assuggested by Ref.[36], we use the JLab data points with Q > . 268 GeV to do fit . By using the scale-invariantPMC pQCD series, the quality of fit χ /d.o.f improvesto be ∼ 34 for APT model, ∼ 52 for WEB model, ∼ ∼ 38 for CON model, respectively,which correspond to the p -value around 95% − 99% [87]. As shown by Fig. 4, the predictions drops down quickly in verysmall Q -region, and the quality of fit is greatly affected by thedata within this Q -region, indicating the twist-expansion couldbe failed in very small Q -region. (a)APT model (b)WEB model (c)MPT model (d)CON model FIG. 4: The spin structure function Γ p − n ( Q ) with both leading-twist and high-twist contributions under four α s models: (a)the APT model; (b) the WEB model; (c) the MPT model; and (d) the CON model. The leading-twist perturbative contributionshave been calculated up to N LO level and N LO level before and after applying the PMC scale-setting approach, respectively.The shaded band shows the prediction under conventional scale-setting approach by varying µ r ∈ [ Q/ , Q ]. The solid line isthe scale-invariant PMC prediction. IV. SUMMARY In the paper, we have applied the PMC single-scale ap-proach to deal with the perturbative series of the leading-twist part of Γ p − n ( Q ) up to N LO level. The pQCDseries for both Γ p − n ( Q ) and the PMC scale Q ∗ are con-vergent in large Q -region. We have also provided aprediction on the uncalculated N LO by using the moreconvergent and scheme-and-scale invariant PMC confor-mal series. Thus a more accurate pQCD prediction onΓ p − n ( Q ) can be achieved by applying the PMC.Basing on the PMC predictions on the perturbativepart, we then provide a novel determination of the high-twist contributions by using the JLab data, whose mo-mentum transfer lies in the range of 0 . 054 GeV ≤ Q ≤ . 739 GeV . In large Q -region, the high-twist contribu-tions to Γ p − n ( Q ) are power suppressed and negligible,which are however sizable in low and intermediate Q -region; Fig. 2 shows that in low Q -region, the leading-twist terms alone cannot explain the JLab data. Takingthe high-twist contributions up to twist-6 accuracy, wehave fixed the twist-4 coefficient f p − n and the twist-6coefficient µ by using four typical α s -models, which give f p − n | APT = − . ± . 013 and µ | APT = 0 . ± . f p − n | WEB = − . ± . 013 and µ | WEB = 0 . ± . f p − n | MPT = − . ± . 013 and µ | MPT = 0 . ± . f p − n | CON = − . ± . 013 and µ | CON = 0 . ± . Acknowledgments: This work was supported in partby the Natural Science Foundation of China under GrantNo.11625520 and No.12047564, by the Fundamental Re-search Funds for the Central Universities under GrantNo.2020CQJQY-Z003, and by the Chongqing Gradu-ate Research and Innovation Foundation under GrantNo.ydstd1912. Appendix A: the reduced perturbative coefficients ˆ r i,j In this appendix, we give the required reduced coeffi-cients ˆ r i,j for the perturbative series of the leading-twistpart of Γ p − n ( Q , µ r ) up to four-loop level, i.e.,ˆ r , = 34 γ ns1 , ˆ r , = 34 γ ns2 − (cid:0) γ ns1 (cid:1) , ˆ r , = 34 Π ns1 + K ns1 , ˆ r , = 34 γ ns3 − γ ns2 γ ns1 + 2764 (cid:0) γ ns1 (cid:1) , ˆ r , = 34 Π ns2 + 12 K ns2 − γ ns1 (cid:18) K ns1 + 94 Π ns1 (cid:19) , ˆ r , = 0 , ˆ r , = 34 γ ns4 − γ ns3 γ ns1 − (cid:0) γ ns2 (cid:1) + 8164 γ ns2 (cid:0) γ ns1 (cid:1) − (cid:0) γ ns1 (cid:1) , ˆ r , = 34 Π ns3 + 13 K ns3 − γ ns1 ( K ns2 + 3Π ns2 ) − γ ns2 (cid:18) K ns1 + 32 Π ns1 (cid:19) + (cid:0) γ ns1 (cid:1) (cid:18) K ns1 + 274 Π ns1 (cid:19) , ˆ r , = − (cid:0) Π ns1 (cid:1) − K ns1 Π ns1 , ˆ r , = 0 , where γ ns i , Π ns i and K ns i can be found in Refs.[55, 56]. Appendix B: Derivation of the parameters f p − n (1 GeV ) and µ According to Ref.[36], it is straightforward to deducethe squares of the standard deviation of f p − n (1 GeV )and µ , based on the minimization of χ /d.o.f . Com-paring the experimental data Γ p − n , exp ( Q ) and theoreticalprediction Γ p − n , the ( Q ), we describe this difference at a spe-cific point Q j by using the following symbol: y j ≡ Γ p − n , exp ( Q j ) − Γ p − n , the ( Q j ) . (B1)The quality parameter χ is rewritten as χ ( µ , µ ) = X j w j ( y j − µ z j − µ z j ) . (B2)where z j ≡ /Q j and w j ≡ /σ j, stat from the squaredstatistical uncertainties of experimental values. The val-ues ˆ µ and ˆ µ can be obtained by the condition: thesimultaneous minimization of χ ( µ , µ ). Here the re-duced values of ˆ µ and ˆ µ are defined asˆ µ = − yz z + yz z z z − z z , ˆ µ = yz z − yz z z z − z z (B3)with the unnormalized “average” A ≡ X j w j A ( z j ) . (B4)Simplifying f p − n and µ as constants, there are thefollowing approximations: the statistical uncertainties atdifferent points are considered as uncorrelated; the sys-tematical uncertainties at different point are consideredas uncorrelated between different experiments but corre-lated under the same experiment. After using Taylor expansion of χ ( µ , µ ) around thepoint (ˆ µ , ˆ µ ) up to the terms quadratic in the deviations,the approximate relations are [36] χ (ˆ µ + σ (ˆ µ stat4 ) , ˆ µ ) = χ + z z z z − z z , (B5) χ (ˆ µ , ˆ µ + σ (ˆ µ stat6 )) = χ + z z z z − z z . (B6)Thus, the statistical uncertainties of fits parameters µ and µ can be obtained from Eqs.(B5,B6), and σ ( ˆ f stat ) = M σ (ˆ µ stat ) from Eq.(12).Moreover, the calculation of systematical uncertaintiesof the parameters µ and µ at different point should con-sider the weighted mean values of different experiments.From this, we firstly express the quantities yz and yz inform of ˆ µ and ˆ µ yz = ˆ µ z + ˆ µ z , yz = ˆ µ z + ˆ µ z , (B7)Considering the two different experimental group fromRefs.[34, 35], the ˆ µ and ˆ µ redefined byˆ µ = ˜ µ (1)4 + ˜ µ (2)4 , (B8)ˆ µ = ˜ µ (1)6 + ˜ µ (2)6 , (B9)and ˜ µ ( i )4 = ˜ α i ˆ µ ( i )4 − ˜ k i ˆ µ ( i )6 , (B10)˜ µ ( i )6 = ˜ β i ˆ µ ( i )6 + ˜ h i ˆ µ ( i )4 . (B11)where i = 1 , α i , ˜ β i , ˜ k i and ˜ h i satisfying that˜ α i = 1 D all ( X j =1 D ( ij ) ) , (B12)˜ β j = 1 D all ( X i =1 D ( ij ) ) , (B13)˜ k i = 1 D all 2 X j =1: j = i ( − z i ) z j ) + z j ) z i ) ) , (B14)˜ h i = 1 D all 2 X j =1: j = i ( − z i ) z j ) + z j ) z i ) ) , (B15) D ( ij ) = z i ) z j ) − z i ) z j ) ( i, j = 1 , , (B16) D all = X j =1 2 X i =1 D ( ij ) = z z − z z . (B17)The D all is the unnormalized averages over two experi-ments. Then, we estimate the systematical uncertaintyby averaging the deviations∆˜ µ ( i ) , sys N ≡ σ (˜ µ ( i ) ,sysN ) ≈ 12 ( | ˜ µ ( i ) N (UP) − ˜ µ ( i ) N | + | ˜ µ ( i ) N (DO) − ˜ µ ( i ) N | ) , (B18)0where N = 4 , 6, symbols “UP” and “DO” refer to thevalues ˜ µ ( i ) N extracted from the experimental data plus orminus the uncertainty σ j, sys at momentum Q j , respec-tively. Finally, the systematical uncertainty ∆ˆ µ sys N ex-tracted from independent experiments are∆ˆ µ sys4 ≡ σ (˜ µ sys4 ) = (cid:20) X i =1 σ (˜ µ ( i ) ,sys ) (cid:21) / , (B19)∆ˆ µ sys6 ≡ σ (˜ µ sys6 ) = (cid:20) X i =1 σ (˜ µ ( i ) ,sys ) (cid:21) / , (B20) and the corresponding uncertainties for the parameter f p − n (1 GeV ) is∆ ˆ f sys2 = σ ( ˆ f sys2 ) = 94 M σ (ˆ µ sys4 ) . (B21) [1] J. D. Bjorken, “Applications of the chiral U (6) × (6) al-gebra of current densities,” Phys. Rev. , 1467 (1966).[2] J. D. Bjorken, “Inelastic Scattering of Polarized Leptonsfrom Polarized Nucleons,” Phys. Rev. D , 1376 (1970).[3] P. L. Anthony et al. [E142 Collaboration], “Deep inelasticscattering of polarized electrons by polarized He and thestudy of the neutron spin structure,” Phys. Rev. D ,6620 (1996).[4] K. Abe et al. [E143 Collaboration], “Precision measure-ment of the proton spin structure function g p ,” Phys.Rev. Lett. , 346 (1995).[5] K. Abe et al. [E143 Collaboration], “Precision measure-ment of the deuteron spin structure function g d ,” Phys.Rev. Lett. , 25 (1995).[6] K. Abe et al. [E143 Collaboration], “Measurements ofthe proton and deuteron spin structure function g andasymmetry A ,” Phys. Rev. Lett. , 587 (1996).[7] K. Abe et al. [E143 Collaboration], “Measurements of the Q dependence of the proton and deuteron spin structurefunctions g p and g d ,” Phys. Lett. B , 61 (1995)[8] K. Abe et al. [E154 Collaboration], “Measurement of theneutron spin structure function g n and asymmetry A n ,”Phys. Lett. B , 377 (1997).[9] K. Abe et al. [E154 Collaboration], “Next-to-leading or-der QCD analysis of polarized deep inelastic scatteringdata,” Phys. Lett. B , 180 (1997).[10] K. Abe et al. [E143 Collaboration], “Measurements of theproton and deuteron spin structure functions g and g ,”Phys. Rev. D , 112003 (1998)[11] P. L. Anthony et al. [E155 Collaboration], “Measurementof the proton and deuteron spin structure functions g and asymmetry A ,” Phys. Lett. B , 529 (1999).[12] P. L. Anthony et al. [E155 Collaboration], “Measure-ment of the deuteron spin structure function g d ( x ) for1(GeV / c) < Q < / c) ,” Phys. Lett. B ,339 (1999).[13] P. L. Anthony et al. [E155 Collaboration], “Measure-ments of the Q dependence of the proton and neutronspin structure functions g p and g n ,” Phys. Lett. B ,19 (2000).[14] P. L. Anthony et al. [E155 Collaboration], “Precisionmeasurement of the proton and deuteron spin structurefunctions g and asymmetries A ,” Phys. Lett. B ,18 (2003).[15] D. Adams et al. [SMC Collaboration], “Measurement ofthe spin dependent structure function g ( x ) of the pro-ton,” Phys. Lett. B , 399 (1994).[16] D. Adams et al. [SMC Collaboration], “Spin asymmetry in muon - proton deep inelastic scattering on a trans-versely polarized target,” Phys. Lett. B , 125 (1994).[17] D. Adams et al. [SMC Collaboration], “A New measure-ment of the spin dependent structure function g ( x )1 of thedeuteron,” Phys. Lett. B , 248 (1995).[18] D. Adams et al. [SMC Collaboration], “The Spin depen-dent structure function g ( x )1 of the deuteron from polar-ized deep inelastic muon scattering,” Phys. Lett. B ,338 (1997).[19] D. Adams et al. [SMC Collaboration], “Spin structure ofthe proton from polarized inclusive deep inelastic muon- proton scattering,” Phys. Rev. D , 5330 (1997).[20] K. Ackerstaff et al. [HERMES Collaboration], “Measure-ment of the neutron spin structure function g n with apolarized He-3 internal target,” Phys. Lett. B , 383(1997).[21] K. Ackerstaff et al. [HERMES Collaboration], “Determi-nation of the deep inelastic contribution to the general-ized Gerasimov-Drell-Hearn integral for the proton andneutron,” Phys. Lett. B , 531 (1998).[22] A. Airapetian et al. [HERMES Collaboration], “Measure-ment of the proton spin structure function g p with a purehydrogen target,” Phys. Lett. B , 484 (1998).[23] A. Airapetian et al. [HERMES Collaboration], “Evidencefor quark hadron duality in the proton spin asymmetry A ,” Phys. Rev. Lett. , 092002 (2003).[24] A. Airapetian et al. [HERMES Collaboration], “Precisedetermination of the spin structure function g of theproton, deuteron and neutron,” Phys. Rev. D , 012007(2007).[25] V. Y. Alexakhin et al. [COMPASS Collaboration], “TheDeuteron Spin-dependent Structure Function g d and itsFirst Moment,” Phys. Lett. B , 8 (2007).[26] M. G. Alekseev et al. [COMPASS Collaboration], “TheSpin-dependent Structure Function of the Proton g p anda Test of the Bjorken Sum Rule,” Phys. Lett. B , 466(2010).[27] C. Adolph et al. [COMPASS Collaboration], “The spinstructure function g p1 of the proton and a test of theBjorken sum rule,” Phys. Lett. B , 18 (2016).[28] C. Adolph et al. [COMPASS Collaboration], “FinalCOMPASS results on the deuteron spin-dependent struc-ture function g d1 and the Bjorken sum rule,” Phys. Lett.B , 34 (2017).[29] F. R. Wesselmann et al. [RSS Collaboration], “Protonspin structure in the resonance region,” Phys. Rev. Lett. , 132003 (2007).[30] K. Slifer et al. [RSS Collaboration], “Probing Quark- Gluon Interactions with Transverse Polarized Scatter-ing,” Phys. Rev. Lett. , 101601 (2010).[31] Y. Prok et al. [CLAS Collaboration], “Precision measure-ments of g of the proton and the deuteron with 6 GeVelectrons,” Phys. Rev. C , 025212 (2014).[32] A. Deur et al. , “Experimental determination of the evolu-tion of the Bjorken integral at low Q ,” Phys. Rev. Lett. , 212001 (2004).[33] J. P. Chen, A. Deur and Z. E. Meziani, “Sum rules andmoments of the nucleon spin structure functions,” Mod.Phys. Lett. A , 2745 (2005).[34] A. Deur et al. , “Experimental study of isovector spin sumrules,” Phys. Rev. D , 032001 (2008).[35] A. Deur et al. , “High precision determination of the Q evolution of the Bjorken Sum,” Phys. Rev. D , 012009(2014).[36] C. Ayala, G. Cveti, A. V. Kotikov and B. G. Shaikhatde-nov, “Bjorken polarized sum rule and infrared-safe QCDcouplings,” Eur. Phys. J. C , 1002 (2018).[37] C. Ayala, G. Cveti, A. V. Kotikov and B. G. Shaikhat-denov, “Bjorken sum rule with analytic QCD coupling,”J. Phys. Conf. Ser. , 012016 (2020).[38] R. S. Pasechnik, D. V. Shirkov and O. V. Teryaev,“Bjorken Sum Rule and pQCD frontier on the move,”Phys. Rev. D , 071902 (2008).[39] R. S. Pasechnik, D. V. Shirkov, O. V. Teryaev,O. P. Solovtsova and V. L. Khandramai, “Nucleon spinstructure and pQCD frontier on the move,” Phys. Rev.D , 016010 (2010).[40] V. L. Khandramai, R. S. Pasechnik, D. V. Shirkov,O. P. Solovtsova and O. V. Teryaev, “Four-loop QCDanalysis of the Bjorken sum rule vs data,” Phys. Lett. B , 340 (2012).[41] V. L. Khandramai, O. P. Solovtsova and O. V. Teryaev,“Polarized Bjorken Sum Rule Analysis: Revised,” Non-lin. Phenom. Complex Syst. , 93 (2013).[42] S. J. Brodsky and X. G. Wu, “Self-Consistency Re-quirements of the Renormalization Group for Setting theRenormalization Scale,” Phys. Rev. D , 054018 (2012).[43] X. G. Wu, Y. Ma, S. Q. Wang, H. B. Fu, H. H. Ma,S. J. Brodsky and M. Mojaza, “Renormalization GroupInvariance and Optimal QCD Renormalization Scale-Setting,” Rep. Prog. Phys. , 126201 (2015).[44] S. J. Brodsky and X. G. Wu, “Scale Setting Using theExtended Renormalization Group and the Principle ofMaximum Conformality: the QCD Coupling Constantat Four Loops,” Phys. Rev. D , 034038 (2012).[45] S. J. Brodsky and X. G. Wu, “Eliminating the Renor-malization Scale Ambiguity for Top-Pair Production Us-ing the Principle of Maximum Conformality,” Phys. Rev.Lett. , 042002 (2012).[46] M. Mojaza, S. J. Brodsky and X. G. Wu, “SystematicAll-Orders Method to Eliminate Renormalization-Scaleand Scheme Ambiguities in Perturbative QCD,” Phys.Rev. Lett. , 192001 (2013).[47] S. J. Brodsky, M. Mojaza and X. G. Wu, “SystematicScale-Setting to All Orders: The Principle of MaximumConformality and Commensurate Scale Relations,” Phys.Rev. D , 014027 (2014).[48] X. G. Wu, S. J. Brodsky and M. Mojaza, “The Renor-malization Scale-Setting Problem in QCD,” Prog. Part.Nucl. Phys. , 44 (2013).[49] X. G. Wu, S. Q. Wang and S. J. Brodsky, “Importanceof proper renormalization scale-setting for QCD testing at colliders,” Front. Phys. , 111201 (2016).[50] X. G. Wu, J. M. Shen, B. L. Du, X. D. Huang, S. Q. Wangand S. J. Brodsky, “The QCD Renormalization GroupEquation and the Elimination of Fixed-Order Scheme-and-Scale Ambiguities Using the Principle of MaximumConformality,” Prog. Part. Nucl. Phys. , 103706(2019).[51] J. M. Shen, X. G. Wu, B. L. Du and S. J. Brodsky, “NovelAll-Orders Single-Scale Approach to QCD Renormaliza-tion Scale-Setting,” Phys. Rev. D , 094006 (2017).[52] X. G. Wu, J. M. Shen, B. L. Du and S. J. Brodsky, “Noveldemonstration of the renormalization group invariance ofthe fixed-order predictions using the principle of max-imum conformality and the C -scheme coupling,” Phys.Rev. D , 094030 (2018).[53] X. C. Zheng, X. G. Wu, S. Q. Wang, J. M. Shen andQ. L. Zhang, “Reanalysis of the BFKL Pomeron at thenext-to-leading logarithmic accuracy,” JHEP , 117(2013).[54] H. Y. Bi, X. G. Wu, Y. Ma, H. H. Ma, S. J. Brodskyand M. Mojaza, “Degeneracy Relations in QCD and theEquivalence of Two Systematic All-Orders Methods forSetting the Renormalization Scale,” Phys. Lett. B ,13 (2015).[55] P. A. Baikov, K. G. Chetyrkin and J. H. Kuhn, “AdlerFunction, Bjorken Sum Rule, and the Crewther Relationto Order α s in a General Gauge Theory,” Phys. Rev.Lett. , 132004 (2010)[56] P. A. Baikov, K. G. Chetyrkin, J. H. Kuhn and J. Rit-tinger, “Vector Correlator in Massless QCD at Order O ( α s ) and the QED beta-function at Five Loop,” JHEP , 017 (2012).[57] T. Appelquist and J. Carazzone, “Infrared Singularitiesand Massive Fields,” Phys. Rev. D , 2856 (1975).[58] D. J. Gross and F. Wilczek, “Ultraviolet Behavior ofNonabelian Gauge Theories,” Phys. Rev. Lett. , 1343(1973).[59] H. D. Politzer, “Reliable Perturbative Results for StrongInteractions?,” Phys. Rev. Lett. , 1346 (1973).[60] W. E. Caswell, “Asymptotic Behavior of NonabelianGauge Theories to Two Loop Order,” Phys. Rev. Lett. , 244 (1974).[61] O. V. Tarasov, A. A. Vladimirov and A. Y. Zharkov,“The Gell-Mann-Low Function of QCD in the ThreeLoop Approximation,” Phys. Lett. B , 429 (1980).[62] S. A. Larin and J. A. M. Vermaseren, “The Three loopQCD Beta function and anomalous dimensions,” Phys.Lett. B , 334 (1993).[63] T. van Ritbergen, J. A. M. Vermaseren and S. A. Larin,“The Four loop beta function in quantum chromodynam-ics,” Phys. Lett. B , 379 (1997).[64] K. G. Chetyrkin, “Four-loop renormalization of QCD:Full set of renormalization constants and anomalous di-mensions,” Nucl. Phys. B , 499 (2005).[65] M. Czakon, “The Four-loop QCD beta-function andanomalous dimensions,” Nucl. Phys. B , 485 (2005).[66] P. A. Baikov, K. G. Chetyrkin and J. H. Kuhn, “Five-Loop Running of the QCD coupling constant,” Phys.Rev. Lett. , 082002 (2017).[67] J. L. Basdevant, “The Pad´e approximation and its phys-ical applications,” Fortsch. Phys. , 283 (1972).[68] M. A. Samuel, G. Li and E. Steinfelds, “Estimating per-turbative coefficients in quantum field theory using Pad´eapproximants,” Phys. Lett. B , 188 (1994). [69] M. A. Samuel, J. R. Ellis and M. Karliner, “Comparisonof the Pad´e approximation method to perturbative QCDcalculations,” Phys. Rev. Lett. , 4380 (1995).[70] B. L. Du, X. G. Wu, J. M. Shen and S. J. Brodsky, “Ex-tending the Predictive Power of Perturbative QCD,” Eur.Phys. J. C , 182 (2019).[71] H. M. Yu, W. L. Sang, X. D. Huang, J. Zeng, X. G. Wuand S. J. Brodsky, “Scale-Fixed Predictions for γ + η c pro-duction in electron-positron collisions at NNLO in per-turbative QCD,” arXiv:2007.14553 [hep-ph].[72] X. D. Huang, X. G. Wu, J. Zeng, Q. Yu, X. C. Zhengand S. Xu, “Determination of the top-quark MS runningmass via its perturbative relation to the on-shell masswith the help of the principle of maximum conformality,”Phys. Rev. D , 114024 (2020).[73] Q. Yu, X. G. Wu, J. Zeng, X. D. Huang and H. M. Yu,“The heavy quarkonium inclusive decays using the prin-ciple of maximum conformality,” Eur. Phys. J. C , 362(2020).[74] Q. Yu, X. G. Wu, S. Q. Wang, X. D. Huang, J. M. Shenand J. Zeng, “Properties of the decay H → γγ using theapproximate α s corrections and the principle of maxi-mum conformality,” Chin. Phys. C , 093102 (2019).[75] X. D. Ji and P. Unrau, “ Q dependence of the proton’s G structure function sum rule,” Phys. Lett. B , 228(1994).[76] E. V. Shuryak and A. I. Vainshtein, “Theory of PowerCorrections to Deep Inelastic Scattering in QuantumChromodynamics. II. Q − Effects: Polarized Target,”Nucl. Phys. B , 141 (1982).[77] H. Kawamura, T. Uematsu, J. Kodaira and Y. Yasui,“Renormalization of twist four operators in QCD Bjorkenand Ellis-Jaffe sum rules,” Mod. Phys. Lett. A , 135(1997).[78] J. Blumlein and A. Tkabladze, “Target mass correctionsfor polarized structure functions and new sum rules,”Nucl. Phys. B , 427 (1999).[79] K. G. Chetyrkin, B. A. Kniehl and M. Steinhauser,“Strong coupling constant with flavor thresholds at fourloops in the Modified Minimal-Subtraction scheme,”Phys. Rev. Lett. , 2184 (1997).[80] D. V. Shirkov and I. L. Solovtsov, “Analytic model for the QCD running coupling with universal α s (0) value,”Phys. Rev. Lett. , 1209 (1997).[81] D. V. Shirkov, “On the analytic ’causal’ model for theQCD running coupling,” Nucl. Phys. Proc. Suppl. ,106 (1998).[82] B. R. Webber, “QCD power corrections from a sim-ple model for the running coupling,” JHEP , 012(1998).[83] D. V. Shirkov, “The Unitary mechanism of infrared freez-ing in QCD with massive gluons,” Phys. Atom. Nucl. ,1928 (1999).[84] D. V. Shirkov, “’Massive’ Perturbative QCD, regular inthe IR limit,” Phys. Part. Nucl. Lett. , 186 (2013).[85] F. Halzen, G. I. Krein and A. A. Natale, “Relating theQCD pomeron to an effective gluon mass,” Phys. Rev. D , 295 (1993).[86] Cornwall, M. John, “Dynamical mass generation in con-tinuum quantum chromodynamics,” Phys. Rev. D ,1453 (1982).[87] P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp.Phys. 2020, 083C01 (2020).[88] A. Deur, S. J. Brodsky and G. F. de Teramond, “Con-necting the Hadron Mass Scale to the Fundamental MassScale of Quantum Chromodynamics,” Phys. Lett. B ,528 (2015).[89] X. D. Ji, “Spin structure functions of the nucleon,” [hep-ph/9510362].[90] E. Stein, P. Gornicki, L. Mankiewicz and A. Schafer,“QCD sum rule calculation of twist four corrections toBjorken and Ellis-Jaffe sum rules,” Phys. Lett. B ,107 (1995).[91] I. I. Balitsky, V. M. Braun and A. V. Kolesnichenko,“Power corrections 1 /Q to parton sum rules for deepinelastic scattering from polarized targets,” Phys. Lett.B , 245 (1990).[92] N. Y. Lee, K. Goeke and C. Weiss, “Spin dependent twistfour matrix elements from the instanton vacuum: Flavorsinglet and nonsinglet,” Phys. Rev. D , 054008 (2002).[93] A. V. Sidorov and C. Weiss, “Higher twists in polarizedDIS and the size of the constituent quark,” Phys. Rev.D73