A novel determination of non-perturbative contributions to Bjorken sum rule

In the present paper, we first give a detailed study on the pQCD corrections to the leading-twist part of BSR. Previous pQCD corrections to the leading-twist part derived under conventional scale-setting approach up to ${\cal O}(\alpha_s^4)$-level still show strong renormalization scale dependence. The principle of maximum conformality (PMC) provides a systematic way to eliminate conventional renormalization scale-setting ambiguity by determining the accurate $\alpha_s$-running behavior of the process with the help of renormalization group equation. Our calculation confirms the PMC prediction satisfies the standard renormalization group invariance, e.g. its fixed-order prediction does scheme-and-scale independent. In low $Q^2$-region, the effective momentum of the process is small and to have a reliable prediction, we adopt four low-energy $\alpha_s$ models to do the analysis. Our predictions show that even though the high-twist terms are generally power suppressed in high $Q^2$-region, they shall have sizable contributions in low and intermediate $Q^2$ domain. By using the more accurate scheme-and-scale independent pQCD prediction, we present a novel fit of the non-perturbative high-twist contributions by comparing with the JLab data.


I. INTRODUCTION
The Bjorken Sum Rule (BSR) [1,2], which describes the polarized spin structure of nucleon, has been measured via polarized deep inelastic scattering (DIS) by various experimental collaborations . Using the operator product expansion (OPE), the BSR of the spin structure function can be calculated by separating the perturbative contribution of the matrix elements of local product operators from its non-perturbative contributions [1,2], e.g.
where g p,n 1 (x, Q 2 ) is the spin-dependent proton or neutron structure function with Bjorken scaling variable x, and g A is the nucleon axial charge. The BSR relates the difference of the proton and the neutron structure functions Γ p 1 and Γ n 1 , and only the flavor non-singlet quark operators appear in perturbative part, resulting as * Electronic address: yuq@cqu.edu.cn † Electronic address: wuxg@cqu.edu.cn ‡ Electronic address: zhouhua@cqu.edu.cn § Electronic address: hxud@cqu.eud.cn the perturbative non-singlet leading-twist contributions E ns (Q 2 ). The non-perturbative contribution is generally power suppressed in comparison to the leading-twist terms, which has been written as a power series over 1/Q 2 . Contributions from the high-twist terms could be sizable in low and intermediate Q 2 regions, and then the BSR provides a good platform for testing the perturbative and non-perturbative QCD contributions.
Analyses of E ns (Q 2 ) under the MS-scheme have been given in the literature, such as Refs. [32][33][34][35][36]. Additional treatment on extending the pQCD prediction to low Q 2region has been done by using low-energy models for the strong coupling constant (α s ) such as the analytic perturbation theory (APT), the "massive analytic pQCD theory" (MPT), the 2δ-or 3δ-analytic QCD variants [37][38][39][40][41][42]. In all those treatments, there are large renormalization scale (µ r ) dependence for the perturbative part due to the using of "guessed" µ r ; that is, in those analyses, the central ("optimal") value of E ns (Q 2 ) is usually derived by setting µ r = Q, and then by varying it within an arbitrary range such as [Q/2, 2Q] to estimate its uncertainty. Such guessing choice breaks the renormalization group invariance [43,44] and leads to conventional renormalization scale-and-scheme ambiguities due to the mismatching of the perturbative coefficients and the α s at each order. In the literature, the principle of maximum conformality (PMC) [45][46][47][48] has been suggested to eliminate such renormalization scale-and-scheme ambiguities. It is well known that the α s -running behavior is governed by the renormalization group equation (RGE).
The existence of the {β i }-terms emerged in the perturbative series is thus helpful for fixing exact α s -value of the pQCD approximant of a physical observable. And instead of choosing an optimal µ r , the PMC fixes the correct magnitude of α s by using RGE, whose argument is called as the PMC scale, which is independent to any choice of µ r . The PMC prediction is scale-and-scheme independent, more detail and applications of the PMC can be found in the reviews [49][50][51].
To achieve a reliable prediction for the BSR high-twist contributions, it is important to have an accurate pQCD prediction on E ns (Q 2 ). In the present paper, we shall first adopt the PMC single-scale approach [52] to deal with the perturbative part of the BSR, and then give a new determination of the non-perturbative high-twist contributions by comparing with the JLab data. The PMC singlet-scale approach follows the same idea of the original multi-scale approach [45][46][47][48], which determines an overall effective momentum flow of the process by using the RGE, whose magnitude corresponds to the weighted average of the multi-scales of the multi-scale approach at each order. It has also been demonstrated that the prediction under the PMC singlet-scale approach is schemeand-scale independent up to any fixed order [53]. Though different from conventional scale ambiguity, there is residual scale dependence for fixed-order prediction due to unknown perturbative terms [54]. Such residual scale dependence can be greatly suppressed due to both α s -power suppression and exponential suppression. A detailed discussion on the residual scale dependence can be found in the recent review [51].
The remaining parts of the paper are organized as follows. In Sec.II, we present the calculation technology for the polarized Bjorken sum rule Γ p−n 1 . The PMC treatment of the pQCD contributions to the leading-twist part and the non-perturbative high-twist contributions shall be given. In Sec.III, we give the numerical results and discussions. Sec.V is reserved for a summary.

II. CALCULATION TECHNOLOGY
In large Q 2 -region, contributions from the leadingtwist terms are dominant and those of the nonperturbative high-twist terms are generally power suppressed. In low and intermediate Q 2 -region, contributions from the high-twist terms may have large contributions. In the following, we shall analyze the pQCD contributions to the leading-twist terms by using the PMC single-scale approach, and then give an estimation of the contributions from the non-perturbative hightwist terms. In low Q 2 -region, the low-energy α s models should be used; and for clarity, we shall adopt four lowenergy α s models to do our discussion.
A. Perturbative series of the leading-twist terms The perturbative expansion over α s for the hard part of the leading-twist terms E ns (Q 2 ) has been calculated up to next-to-next-to-next-to leading order (N 3 LO), which can be written as where a(µ r ) = α s (µ r )/π and the perturbative coefficients r i are power series of the active flavor numbers n f , The explicit expressions of the coefficients c i,j have been given in Refs. [55,56]. To apply the PMC, we need to use the general QCD degeneracy relations [57] among different orders to make the transformation of the n fseries to {β i }-series, i.e. we need to rewrite E ns (Q 2 ) in the following form, where the coefficients r i,j up to N 3 LO-order level are Generally, the coefficients r i,j =0 are functions of the logarithm ln(µ 2 r /Q 2 ). If setting µ r = Q, all those types of log-terms becomes zero, leading to a renormalon-free more convergent pQCD series; this explains why people usually choose µ r = Q as the optimal scale for conventional scale-setting approach. Those coefficients can be reexpressed as where the combination coefficients C k j = j!/k!(j − k)!, and the coefficientsr i,j = r i,j | µr =Q . For convenience, we put the reduced coefficientsr i,j in the Appendix A.
In Eq.(3), the {β i }-terms at each perturbative order govern the correct α s -running behavior, which inversely can be used to determine the effective magnitude of α s . Practically, by requiring all the RGE-involved nonconformal {β i }-terms to be zero, one can achieve an overall effective α s and hence the PMC scale Q ⋆ , and then the resultant pQCD series becomes the following schemeindependent conformal series: where the PMC scale Q ⋆ can be fixed up to next-tonext-to-leading-log (NNLL) accuracy by using the N 3 LO perturbative series, e.g. where T 2 = 4(r 1,0r2,0r3,1 −r 2 2,0r 2,1 ) + 3(r 1,0r2,1r3,0 −r 2 1,0r 4,1 ) r 3 Those equations show the PMC scale Q ⋆ is exactly free of µ r , together with the µ r -independent conformal coefficientsr i,0 , the PMC prediction is exactly independent to any choice of µ r . Thus the conventional scale-setting ambiguity can be eliminated at any fixed-order by applying the PMC [53]. As a byproduct, due to the elimination of divergent renormalon terms in the resultant PMC perturbative series (16), the pQCD convergence can be naturally improved. Those properties greatly improve the precision of the pQCD theory. For a perturbative theory, it is important to have a reliable way to estimate the magnitude of the uncalculated higher-order terms. The scale-invariant and schemeinvariant PMC conformal series, which is also more convergent than the conventional series, is quite suitable for such purpose. A way of using the PMC series together with the Padé approximation approach (PAA) [68][69][70] has been suggested in Ref. [71]. Some successful applications of this method can be found in Refs. [72][73][74][75]. We shall adopt this method to estimate the magnitude of the unknown O(α 5 s )-terms of E ns (Q 2 ), 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 perturbative series. For the purpose, people usually adopts a fractional function as the generating function. More explicitly, the [N/M ]-type generating function of a pQCD approximant ρ n (Q) = n i=1r i,0 a i is defined as where M ≥ 1 and N + M + 1 = n. The perturbative coefficients C i in Eq.(21) can be expressed by the known coefficients b i∈[0,N ] and c j∈ [1,M] . Inversely, if we have known the coefficients C i 's up to n th -order level, one can determine the coefficients b i∈[0,N ] and c j∈ [1,M] , and then achieve a prediction for the uncalculated (n + 1) th -order coefficient C n+1 . At the present, the leading-twist term E ns | PMC has been known up to N 3 LO-level, and the four coefficients are known, C i =r i,0 for i ∈ [1,4]. Then the predicted N 4 LO-coefficient becomeŝ where the [0/n − 1]-type PAA generating function has been implicitly adopted, which is the preferable type for the convergent PMC series [71].

B. Contributions from the non-perturbative high-twist terms
The non-perturbative contributions to the BSR can be expanded in 1/Q 2 -power series as Eq.
can be written as [76][77][78] where M ≈ 0.94 GeV is the nucleon mass. The leadingtwist target mass correction a p−n 2 can be calculated by using the leading-twist part of g p−n 1 , which is kinematically of high-twist [79] and its magnitude at Q 2 = 1 GeV 2 is 0.031 ± 0.010 [35]. The twist-3 matrix element d p−n 2 is given by whose magnitude at Q 2 = 1 GeV 2 is 0.008 ± 0.0036 [35].
The dynamical values of the twist-2 and twist-3 contributions can be measured by polarized lepton scattering off transversely and longitudinally polarized target. The twist-2 and twist-3 contributions are calculated by the , which is related to the color electric and magnetic polarizabilities of nucleon, plays a pivotal role in phenomenological studies of the high-twist contributions. f p−n 2 is sensitive to Q 2 and its Q 2 -evolution satisfies [77,78] where γ 0 /8β 0 = 32/81 with n f = 3. The magnitude of f p−n

2
(1) shall be fit by comparing with the data. Moreover, it has been argued that the O(Q −4 )-term µ p−n 6 may also have sizable contribution, so we take µ 6 /Q 4 -term into consideration to have a better fit of the data.

C. The strong coupling constant αs
The α s -running behavior in perturbative region is governed by the RGE (15). Its solution can be written as an expansion over the inverse powers of the logarithm L = ln µ 2 r /Λ 2 ; and up to four-loop level, we have [80] where Λ is the scheme-dependent asymptotic scale, which could be fixed by matching the measured value of α s at a reference scale such as M Z or m τ to its predicted value under a specific scheme.
In infrared region, when the scale is close to Λ or even smaller, α s becomes large whose magnitude cannot be well described by the RGE. To make the QCD prediction more reliable, we shall adopt four low-energy models for the α s to do our calculation.
The first low-energy model is based on the analytical perturbation theory (APT) [81,82], and we call it as the APT model. In APT model, its strong coupling constant α APT s is described by applying the perturbation theory directly to the spectral function, which takes the following form, where µ is the energy scale, y = µ 2 /Λ 2 with where φ(z) satisfies 1/φ(z) Its freezing value is close to α APT s (10 −10 )/π ≈ 0.43. The second low-energy model is an alteration of Eq. (27), we call it as the WEB model [83], which is suggested to suppress the nonperturbative power corrections of the APT model, and it takes the following form where these phenomenological parameters b = 1/4 and p = c = 4. The obtained corresponding approximate freezing value ∼ α WEB (10 −10 )/π ≈ 0.21. The third low-energy model is based on the "massive analytic pQCD theory" (MPT) [84,85], which takes the phenomenological glue-ball mass m gl = √ ξΛ as the infrared regulator, and we call it as the MPT model. It takes the following form whose freezing value at the origin satisfies a cr = π/(β 0 ln ξ).
Under the Landau gauge, we have a cr | ξ=10±2 = 0.61∓0.05, which leads to the freezing point α MPT s (0)/π = 0.19 −0.01 +0.02 . The fourth low-energy model is based on the continuum theory [86] and we call it as the CON model, where the exchanging gluons with effective dynamical mass m g is adopted and the non-perturbative dynamics of gluons is governed by the corresponding Schwinger-Dyson equation. It takes the following from [86,87], which leads to the freezing point α CON III. NUMERICAL RESULTS To do the numerical analysis, we take the nucleon axial charge ratio g A = 1.2724 ± 0.0023 [88]. The asymptotic QCD scale Λ can be fixed by using the α s -value at the reference point such as α MS s (m τ ) = 0.325±0.016 [88], which gives Λ MS | n f =3 = 0.346 +0.028 −0.029 GeV by using the four-loop RGE. Using the relation (28), we obtain Λ APT | n f =3 = 0.244 +0.033 −0.031 GeV. In Fig. 1, we present the typical running 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 -running behavior derived from the RGE under MS-scheme is given as a comparison. Fig. 1 shows the importance of the using of low-energy models in the region of small energyscale. Using the criteria suggested in Ref. [93] for the analytic matching of α s in perturbative and nonperturbative regimes, we obtain the transition scales (Q 0 ) for various low-energy models, which are ∼ 1.77 GeV, ∼ 1.78 GeV, ∼ 1.78 GeV, ∼ 1.19 GeV for APT, WEB, MPT and CON models, respectively. As a subtle point, because the transition scales Q 0 for the cases of WEB and MPT are slightly bigger than m τ , and for self-consistency, we use the low-energy α WEB/MPT s (m τ ) = 0.325 ± 0.016 to fix Λ, which is 0.206 ± 0.022 GeV or 0.294 +0.033 −0.032 GeV, respecitvely.
For later convenience, in the following discussions, we simply use α MS s to stand for the case of using MS-scheme α s in all Q 2 -region, α APT s to stand for the case of using APT model in low-energy region (Q < Q 0 , as mentioned above, Q 0 is different for different low-energy model) and MS-scheme α s in large Q 2 -region, α WEB The perturbative contributions to the leading-twist part E ns (Q 2 ) has been known up to N 3 LO. Under conventional scale-setting approach, the pQCD series is scale dependent, and by setting µ r = Q, we obtain On the other hand, the pQCD series becomes scale invariant by applying the PMC, and we obtain for any choice of renormalization scale, where Q ⋆ is of perturbative nature, which can be determined up to NNLL accuracy ln Q * 2 Q 2 = −1.08 − 1.87a(Q) − 24.06a 2 (Q). (34) One may observe that the perturbative coefficients in PMC series (33) are much smaller than those of conventional series (32), especially for those of high-orders, which are due to the elimination of divergent renormalon terms as n!β n 0 a n s . This indicates that a much more convergent perturbative series can be achieved by applying the PMC. At the same time, the PMC scale Q ⋆ also shows a fast convergent at high Q-range, e.g. the relative absolute values of the LL, the NLL and the NNLL terms are 1: 0.064 : 0.030 for Q = 100 GeV. Thus the residual scale dependence due to unknown even higher-order terms can be greatly suppressed.
Using the convergent PMC perturbative series, one can obtain a reliable prediction of unknown O(a 5 )-term by using the PAA, e.g. by using Eq.(22), we obtain We present the predicted leading-twist part of the spin structure function Γ p−n 1 (Q 2 ) under four low-energy models in Fig. 2, where the results under conventional and PMC scale-setting approaches are presented. The experimental data are from SLAC [4][5][6][7]10], DESY [20-24], CREN [25][26][27] and JLab [32,34,35]. The PMC predictions are independent to any choice of µ r , and the shaded band shows the conventional renormalization scale uncertainty by varying µ r ∈ [Q/2, 2Q]. Under conventional scale-setting approach, the spin structure function Γ p−n 1 (Q 2 ) shows large scale dependence, especially in low-energy region. In low-energy region, the results by using the IR-fixed couplings are much more reliable. And since couplings behaves differently in low-energy region, the spin structure function Γ p−n 1 (Q 2 ) behaves quite differently for Q → 0. When the energy scale is large enough, such as Q > 1.5 − 2.0 GeV, the perturbative leading-twist terms could explain the experimental data well. Fig. 2 also shows that in low-scale region, the leading-twist terms alone cannot explain the data and one must take the high-twist terms into consideration. By comparing with the data, this fact inversely provides us a good platform to achieve reliable predictions on the magnitudes of high-twist contributions.

B. Analysis of high-twist contributions under various low-energy models
Following the discussions of Sec.II.B, we need to fit two parameters, f p−n 2 (1 GeV 2 ) and µ 6 , so as to determine the high-twist contributions. We adopt the most recent data listed in Refs. [34,35] to do the fitting, whose momentum transfer lies in the range of 0.054 GeV 2 ≤ Q 2 ≤ 4.739 GeV 2 . We adopt the APT, WEB, MPT, and the CON couplings in doing the fitting. The quality of fit is measured by the parameter of χ 2 /d.o.f , e.g.

2
(1 GeV 2 ) and µ 6 are presented in Table. I. The right-most column shows the smallest χ 2 /d.o.f for the predictions before and after applying the PMC under four α s models. The magnitudes of those two parameters are small, which agree with the usual consideration that at large Q 2region, the high-twist terms are power suppressed and are negligible. However in low Q 2 -region, they will have sizable contributions; especially f p−n 2 (1 GeV 2 ) is important for a reliable theoretical prediction on Γ p−n 1,the. (Q 2 ) in low Q 2 -region. Table. I shows that the fitted parameters under conventional scale-setting approach have strong scale dependence, whose quality of fit χ 2 /d.o.f varies from tens to hundreds, and the optimal fit are achieved for the case of µ r ∼ Q. This, together with a better pQCD convergence due to the elimination of divergent log-terms ln µ 2 r /Q 2 , in some sense explain why µ r = Q is usually taken as the preferable renormalization scale for conventional scale-setting approach. On the other hand, the fitted parameters for the PMC scale-setting approach is independent for any choice of renormalization scale, thus a more reliable and accurate prediction is achieved.
At present, the twist-4 coefficient f p−n
sponse of the color magnetic and electric fields to the spin of the nucleon [94,95]. Using the PMC predictions for the hard-part of the leading-twist contributions, we obtain where the errors are squared average of those from ∆d p−n 2 = ±0.0036 and ∆f p−n 2 for the four low-energy α s models (e.g. Table. I).
We present the prediction of Γ p−n 1 (Q 2 ) with both leading-twist and high-twist contributions in Fig. 4. Comparing with Fig. 2, Fig. 4 shows that a more reasonable prediction can be achieved by including hightwist contributions. Under conventional scale-setting approach, the large scale dependence for the leading-twist prediction of Γ p−n 1,Conv. (Q 2 ) can be greatly suppressed by including high-twist terms due to the cancellation of scale dependence among different twist-terms. Under PMC scale-setting approach, the scale-invariant Γ p−n 1,PMC (Q 2 ) under APT, MPT and CON α s models are close in shape, which as shown by Table. I also have close quality of fit χ 2 /d.o.f ; while the PMC prediction under WEB model is slightly different from those of other α s models.
As a final remark, to improve the quality of fit, as suggested by Ref. [37], we use the JLab data points with Q 2 > 0.268 GeV 2 to do fit. By using the scale-invariant PMC pQCD series, the quality of fit χ 2 /d.o.f improves to be ∼ 34 for APT model, ∼ 52 for WEB model, ∼ 34 for MPT model and ∼ 38 for CON model, respectively, which correspond to the p-value around 95% − 99% [88].
C. An analysis of high-twist contributions with massive high-twist expression As shown by Fig. 4, the predictions drops down quickly in very small Q 2 -region, and the quality of fit is greatly affected by the data within this Q 2 -region, indicating the twist-expansion could be failed in very small Q 2 -region. It has been suggested that by using the "massive" hightwist expansion to do the data fitting, cf. [37,38,[89][90][91][92], one may obtain a better explanation of the data in very low Q 2 region. As an attempt, we take the following "massive" high-twist expansion to do the fit [37] where the parameter m represents a dynamical effective gluon mass, whose square satisfies Here we have set the initial scale of the squared mass as 1 GeV, and we shall take the parameters M 2 = 0.5 GeV 2 and p = 0.1 to do the calculation, which are within the suggested range of Ref. [92]. At present, to fit the magnitude of the "massive" high-twist terms, the parameters f p−n

2
(1 GeV 2 ) and m 2 (1 GeV 2 ) are used to fit with the data [34,35]. When doing the fitting with the experiments data within the range of 0.054GeV 2 ≤ Q 2 ≤ 4.739GeV 2 , we adopt four α s models. The results for the two parameters f p−n 2 (1 GeV 2 ), m 2 (1 GeV 2 ) and their corresponding quality of fit χ 2 /d.o.f are presented in Table. II. Those two parameters are obtained by considering the systematic error σ j,sys at each data point Q 2 j into the fitting; We put the details of fitting in the end of Appendix B. Comparing the smallest χ 2 /d.o.f listed in Table. I and Table. II, one may observes that the "massive" BSR shows a better behavior with smaller quality of fit χ 2 /d.o.f . The conventional predictions for twist-4 f p−n 2 (1 GeV 2 ) apparently depends on the choice of µ r . The quality of fit χ 2 /d.o.f for conventional predictions with α s under WEB model varies from tens to hundreds, while similar χ 2 /d.o.f for conventional predictions with α s under APT, MPT and CON models are along with different fit parameters f p−n 2 (1 GeV 2 ) and m 2 (1 GeV 2 ), respectively. If using the PMC scale-independent series and the "massive" high-twist term, we can obtain the corresponding color polarizability χ p−n E and χ p−n B : where the errors are squared average of those from ∆d p−n 2 = ±0.0036 and ∆f p−n 2 for the four low-energy α s models (e.g. Table. II).
To compare with Fig. 4, Fig. 5 shows that by using the "massive" high-twist term with the fitted parameters f p−n 2 (1 GeV 2 ) and m 2 (1 GeV 2 ), a better prediction in agreement with the experiments data for Q 2 below 0.5 GeV 2 can be achieved, which results as a smaller Table. II. Different from the PMC predictions, the scale-dependence for conventional predictions is enhanced in small Q 2 region. Then, without renormalization scale dependence, the PMC predictions for the twist-4 contribution are more reliable; more explicitly, we observe that the quality of fit χ 2 /d.o.f can be improved as ∼ 28 for APT model, ∼ 45 for WEB model, ∼ 29 for MPT model and ∼ 36 for CON model, respectively, all of which correspond to a p-value ≥ 99%.

IV. SUMMARY
In the paper, we have applied the PMC single-scale approach to deal with the perturbative series of the leadingtwist part of Γ p−n 1 (Q 2 ) up to N 3 LO level. The pQCD series for both Γ p−n 1 (Q 2 ) and the PMC scale Q * are convergent in large Q 2 -region. We have also provided a prediction on the uncalculated N 4 LO by using the more convergent and scheme-and-scale invariant PMC confor-  Basing on the PMC predictions on the perturbative part, we then provide a novel determination of the hightwist contributions by using the JLab data, whose momentum transfer lies in the range of 0.054 GeV 2 ≤ Q 2 ≤ 4.739 GeV 2 . In large Q 2 -region, the high-twist contributions to Γ p−n 1 (Q 2 ) are power suppressed and negligible, which are however sizable in low and intermediate Q 2region; Fig. 2 shows that in low Q 2 -region, the leadingtwist terms alone cannot explain the JLab data. The high-twist term is necessary and it can fix this prob- lem with two fit parameters as Fig. 5 shows. Taking the high-twist contributions up to twist-6 accuracy, we have fixed the twist-4 coefficient f p−n 2 and the twist-6 coefficient µ 6 by using four typical α s -models, which give f p−n 2 | APT = −0.120 ± 0.013 and µ 6 | APT = 0.003 ± 0.000, f p−n 2 | WEB = −0.081 ± 0.013 and µ 6 | WEB = 0.001 ± 0.000, f p−n 2 | MPT = −0.128 ± 0.013 and µ 6 | MPT = 0.003 ± 0.000, f p−n 2 | CON = −0.139 ± 0.013 and µ 6 | CON = 0.002 ± 0.000, respectively. Here the errors are squared averages of those from the statistical and systematic errors of the measured data. As an attempt, by taking the "massive" high-twist expansion such as Eq. (45) to do the fit, we have shown that a better explanation of the data in very low Q 2 range can be achieved. (1 GeV 2 ) and µ6 According to Ref. [37], it is straightforward to deduce the squares of the standard deviation of f p−n 2 (1 GeV 2 ) and µ 6 , based on the minimization of χ 2 /d.o.f . Comparing the experimental data Γ p−n 1,exp (Q 2 ) and theoretical prediction Γ p−n 1,the (Q 2 ), we describe this difference at a specific point Q 2 j by using the following symbol: The quality parameter χ 2 is rewritten as where z j ≡ 1/Q 2 j and w j ≡ 1/σ 2 j,stat from the squared statistical uncertainties of experimental values. The valuesμ 4 andμ 6 can be obtained by the condition: the simultaneous minimization of χ 2 (µ 4 , µ 6 ). Here the reduced values ofμ 4 andμ 6 are defined aŝ with the unnormalized "average" Simplifying f p−n 2 and µ 6 as constants, there are the following approximations: the statistical uncertainties at different points are considered as uncorrelated; the systematical uncertainties at different point are considered as uncorrelated between different experiments but correlated under the same experiment.
After using Taylor expansion of χ 2 (µ 4 , µ 6 ) around the point (μ 4 ,μ 6 ) up to the terms quadratic in the deviations, the approximate relations are [37] χ 2 (μ 4 + σ(μ stat 4 ),μ 6 ) = χ 2 min + ) from Eq.(23). Moreover, the calculation of systematical uncertainties of the parameters µ 4 and µ 6 at different point should consider the weighted mean values of different experiments. From this, we firstly express the quantities yz 2 and yz in form ofμ 4 andμ 6 yz =μ 4 z 2 +μ 6 z 3 , yz 2 =μ 4 z 3 +μ 6 z 4 , Considering the two different experimental group from Refs. [34,35], theμ 4 andμ 6 redefined bŷ where i = 1, 2 and weighted factorsα i ,β i ,k i andh i satisfying that The D all is the unnormalized averages over two experiments. Then, we estimate the systematical uncertainty by averaging the deviations where N = 4, 6, symbols "UP" and "DO" refer to the valuesμ (i) N extracted from the experimental data plus or minus the uncertainty σ j,sys at momentum Q j , respectively. Finally, the systematical uncertainty ∆μ sys