Bjorken sum rule with hyperasymptotic precision

We obtain an improved determination of the normalization of the leading infrared renormalon of the Bjorken sum rule: $Z_B (n_f=3)= -0.407\pm 0.119 $. Estimates of higher order terms of the perturbative series are given. We compute the Bjorken sum rule with hyperasymptotic precision by including the leading terminant, associated with the first infrared renormalon. We fit the experimental data to the operator product expansion theoretical prediction with $\hat f_{3}^{\rm PV}$ as the free parameter. We obtain a good agreement with the experiment with $\hat f_{3}^{\rm PV}\times 10^3=32^{+187}_{-196}\;{\rm GeV}^{2}$ for $Q^2 \geq 1$ GeV$^2$.


Introduction
The Bjorken sum rule [1] is one of the cleanest observables, from the theory point of view, on which to apply the operator product expansion (OPE) in its nonperturbative version [2]. At low orders in the 1/Q 2 expansion, it has the following form: where the definitions are the following: 1 The first few terms of the perturbative expansion of C B (Q) are known. The O(α s ) correction has been computed in [4], the O(α 2 s ) correction in [5], the O(α 3 s ) correction [6], and the O(α 4 s ) correction in [7][8][9].
We also consider the correction due to the charm quark (with finite mass) to the perturbative series of C B . It was first computed in Ref. [10], and later corrected in Ref. [11]. Arguments from renormalons in the context of heavy quark physics suggest that the charm decouples at those scales for the high order coefficients of the perturbative expansion [12]. Therefore, this is the situation we will consider in this paper. In any case, the effect of these corrections is tiny.
Power corrections were calculated in [13,14]. The LO renormalization group running of the twist-four operators has been computed in Ref. [15]. The nonperturbative matrix elements are defined in the following way: |g A |s σ = 2 p, s|J 5,3 σ |p, s , where J 5,a σ (x) = ψγ σ γ 5 t a ψ(x) is the nonsinglet axial current, t a is a generator of the flavor group, and J 5 is the singlet axial current. |g A | is the absolute value of the constant of the neutron beta decay, |g A /g V | = F + D = 1.2754 ± 0.0013 [16].
The dimension two condensate f 3 can be related with the expectation value (sandwiched between the proton state) of a local operator. f 3 is scale dependent and it is defined at Q 2 0 , i.e., f 3 is the reduced matrix element of R 3 2σ , renormalized at Q 2 0 , which is defined for the general flavor indices, with t i being the flavor matrices, as R i 2σ = gψG σν γ ν t i ψ, p, s|R i 2σ |p, s = f i s σ (i = 0, · · · , 8), andG µν = 1 2 ε µναβ G αβ is the dual field strength. Each term of the OPE expression written in Eq. (1) is ill defined. Therefore, as such, the OPE is a formal expression, and, without further qualifications, it is of little use. The origin of this problem comes because, even though the Bjorken sum rule is an observable, and, therefore, is well defined, the splitting between the different terms of the OPE is ambiguous. When working in dimensional regularization using minimal-like subtraction schemes, this reflects in the fact that the perturbative series of C B (α s ) is asymptotically divergent. Therefore, its sum does not converge to a number, and a method has to be used to regularize the perturbative sum. A suitable one is to first construct the Borel transform of the perturbative sum and, afterwards, to do the inverse of the Borel transform (also named Borel sum or integral). Such inverse is still ill defined, due to singularities in the positive axis of the Borel plane. The location and character of such singularities are determined by the OPE [17,18], but not the overall normalization. This information is enough to determine the divergence pattern of the perturbative series up to an overall normalization (for the case of the Bjorken sum rule, it has been computed with logarithmic accuracy in [3]). Such divergent pattern of the perturbative expansion associated with the OPE is usually referred to as renormalons [19].
To fully fix the leading term of the OPE, we have to specify how we handle the singularities in the Borel plane when doing the Borel sum. Following the discussion in [20], we use the principal value (PV) prescription for the Borel integral (the median resummation above and below the real axis). The reason is that the outcome is expected to be real and scale/scheme independent (see the discussion in Refs. [20,21]).
Defining C B (α s ) also defines f 3 . Therefore, we rewrite Eq. (1) as where (often we will work with the variable u = β 0 t 4π instead of t) andf For convenience, we have absorbed the prefactor [α s (Q 2 0 )] It is usually stated that, after introducing such prefactor,f 3 is a renormalization group invariant (with the precision the running of the nonperturbative condensate is known). We would like to emphasize that this is not necessarily so. In orderf 3 to be renormalization scale and scheme independent, it is necessary that the regularization of the perturbative sum of C B is made in such a way that it is explicitly scale and scheme independent (something that the PV delivers but not other regularizations of the perturbative sum), so that the difference between the Bjorken sum rule and g A 6 C B (α s (Q)) is also scale and scheme independent. Assuming the validity of the OPE in its nonperturbative version (which we take for granted),f PV 3 can be written asf PV 3 =f PV 3,MS Λ 2 MS , wheref PV 3,MS is a dimensionless constant. Let us emphasize that such equality holds irrespectively of the scheme used for the strong coupling constant. Therefore, we can also generically writê where X stands for the renormalization scheme of the strong coupling constant and we define b = β 1 /(2β 2 0 ) (note that the definition of b that we use in this paper is different from the one used in Ref. [3]). Therefore, providedf PV 3,X is obtained in one scheme, one can easily transform it to a different scheme using the very same conversion factor one uses to transform Λ QCD from one scheme to another. In the present work X = MS.
Whereas the above procedure yields unambiguous and convenient definitions of the different terms of the OPE, this does not mean that we have the means to compute them. In practice, we only have a relatively small set of the first order terms of the perturbative expansion, and the experimental (or lattice) data. Nowadays, it is not possible to compute the nonperturbative corrections from first principles via analytic methods. At present, it is only possible to compute them in some cases, numerically, from lattice simulations. Nevertheless, such computation is unavoidably plagued by perturbative corrections. Actually, those are the dominant contribution to the observable. A paradigmatic case is the computation of the gluon condensate in the lattice, where its value is orders of magnitude bigger than the actual size of the nonperturbative gluon condensate [22,23]. Overall, to fit the nonperturbative condensate to lattice or experimental data, one first needs to compute the perturbative series with exponential [in −1/α(Q)] accuracy (or with power in 1/Q 2 accuracy). Whereas it is not possible to obtain the exact expression of C PV B , it can be computed approximately, and, more importantly, with a well-defined method to quantify the error in a parametric way [20,24,25]. This method adapts the hyperasymptotic expansion used in ordinary differential equations [26] (see also [27]) to the case of quantum field theories with marginal operators, and it has successfully been applied to a variety of observables [23,25,28]. The hyperasymptotic approximation has not yet been confronted directly to experimental data. The Bjorken sum rule yields a fantastic opportunity to check the OPE in its nonperturbative setup directly into experiment, without resorting to lattice. Therefore, it is our plan to apply such a method to the Bjorken sum rule. Nowadays, the perturbative series of C B (α s (Q) is known to high enough orders to start showing its asymptotic nature for the range of energies for which such sum rule has been measured. Therefore, this opens the venue to determine the Borjen sum rule with exponential accuracy. This is very interesting, since it will allow us to determine the leading nonperturbative correction with hyperasymptotic precision.
Overall, we have now turned the problem into evaluating C PV B with the highest possible accuracy, i.e., including the leading terminant. We address this goal in Sec. 2. One necessary ingredient in the evaluation of the terminant is the determination of the normalization of the renormalon. A previous determination can be found in [3]. In the present paper, we give an improved determination. On the one hand, we can use the new term of the perturbative expansion that is now known. Another improvement is a new way to determine the normalization of the renormalon, which has proven to be more stable [29]. As this method had not been applied to the GLS [30] and Ellis-Jaffe [31] sum rule, we profit this paper to do the same analysis for those sum rules. For the case of the GLS sum rule, the perturbative expansion is also known to an order higher [9], but not for the Ellis-Jaffe sum rule. The latter has also been criticized as being less fundamental than the Bjorken and GLS sum rules (see [32]).
The paper is organized as follows. In Sec. II, C PV B will be computed with hyperasymptotic precision, and an updated determination of the normalization of the leading infrared renormalon, as well as new estimates of the higher order terms of the pertubative series, are given. Updated determinations of the normalization of the leading renormalons of the Ellis-Jaffe and GLS sum rules are also given. In Sec. III, the comparison with the experimental data will be done allowing us the extraction off PV 3 . Finally, the conclusions are presented in Sec. IV.
2 Hyperasymptotic approximation to C PV B

Renormalons
We first need to know the renormalon structure of the perturbative series. We take the results relevant to our case from the analysis made in Ref. [3].
The ultraviolet renormalon structure of the moments of the deep inelastic scattering structure functions has been computed in Ref. [33]. For the case of the Bjorken sum rule, the ultraviolet renormalon formally dominates for n f > 2 and n → ∞. Nevertheless, at low orders in perturbation theory, the infrared renormalon appears to be dominant. This can be seen from the fact that the sign of the known terms of the perturbative series is equal, whereas if the ultraviolet renormalon were to be dominant, then we would find a sign alternating series. In any case, we will perform the conformal mapping [34]: to guarantee that the u = −1 renormalon does not play a significant role. This transformation maps the first infrared renormalon to w = 1/3 and all other singularities to the unit circle |w| = 1. In the conformal mapping, the expansion parameter is w = 1/3. In practice, the effect of doing the conformal mapping is small, which, again, points to the fact that the effect of ultraviolet renormalons is small in comparison with the effect of the infrared renormalon located at u = 1. A similar conclusion was obtained in Ref. [35] using Pade approximants. We will only give numbers for the conformal mapping case for theoretical reasons. Nevertheless, as we have already mentioned, they will be quite similar to the computation without conformal mapping. The Borel transform near the closest infrared renormalon singularity has the following structure (we have changed N X → Z X compared with the notation in Ref. [3] to avoid confusion with N P ): where S reg (u) is an analytic function at u = a. a = 1 and [18]: b B dictates the strength of the singularity. It is interesting to study its dependence on n f . In the Bjorken sum rule, if n f ∈ (0, 6) ⇒ 1 + 2b + b B ∈ (1, 2) so, formally, one could just keep the first two terms of the series in Eq. (11), since the next term would go to zero for u → 1.
If the Wilson coefficients multiplying the higher twist operators were known exactly, then one could also fix the coefficients d B r . Unfortunately, only their leading log running is known. Nevertheless, by performing the matching at a generic scale ν, one can resum the terms of the type (1 − u) n ln n (Q 2 /ν 2 ) and obtain the logarithmically dominant contribution to d B r ∼ ln r (Q 2 /ν 2 ). This was done in [3], and we have little to add in this respect in this paper. The outcome for the asymptotic expression of the coefficients was the following: The above expression contains subleading terms in the 1/n expansion. In the strict 1/n expansion, it simplifies to: An alternative that we will consider in this paper is to work with the quantitỹ the coefficients of which have a more compact expression for its asymptotic expansion: As a final comment, the whole previous discussion also applies to the GLS and Ellis-Jaffe sum rules changing any subscript B by the GLS and EJ subscripts respectively. In Note also that for the case n f = 0, the perturbative expansions of the three sum rules coincide, as differences between singlet (S) and nonsinglet (N S) have to do with light fermion effects, which disappear in the n f = 0 approximation.

Determination of the normalization constant
In this subsection, we will obtain the normalization of the leading infrared renormalon for the Borjken, GLS, and Ellis-Jaffe sum rules: Z B , Z GLS , and Z EJ . For short, we will use the notation Z X , where X = {B, GLS, EJ}. A previous determination of these quantities was obtained in Ref. [3]. In such reference, the ultraviolet renormalon was neglected, and the normalization of the leading infrared renormalon was obtained using the following equality (actually, using the analogous expression after doing the conformal mapping of Eq. (10)): where It has been observed in Refs. [12,29] that a more stable result is obtained by determining the normalization directly from the ratio of the exact and asymptotic expression. Therefore, this is the path we will follow in this paper instead. The results will not be very different.
To eliminate the anomalous dimension of the leading infrared renormalon, we consider the perturbative series ofC B (for consistency, the prefactor is expanded with one-loop accuracy), defined in Eq. (15), and take the asymptotic expression given in Eq. (16). To both expressions we do the conformal mapping and consider the ratio of the corresponding coefficients at x = 1 with NNNLO precision, where x = ν/Q. These we will take as our central values.
Actually, for the asymptotic expression, we already take the expected asymptotic expression after conformal mapping, i.e., with the leading singularity located at w = 1/3 [this is nothing but multiplying Eq. (16) by 3 n ]. We have checked that the difference is indeed small, and well inside of what would be our errors. We have also checked that not doing the conformal mapping does not change our determination significantly, which is consistent with the interpretation that the ultraviolet renormalon is not important. Another check of this interpretation comes from taking the n f → ∞ limit of our results. These should converge to the large β 0 results obtained in [36]. Indeed, this is what approximately happens: the result converges to the result expected from the large β 0 analysis for the infrared renormalon. Actually, this also happens even if we do not do the conformal mapping. The error in our determination of the normalization of the infrared renormalon is due to the incomplete knowledge of the perturbative series. This reflects in that the result will depend on the scale, the order we truncate, and the explicit expression we use to determine the normalization. We use these three methods as indicators of the error. In particular, we do the following: (i) Scale variation between x ∈ (1/2, 2).
(iii) We consider the difference between using (the conformal transform of) Eqs. (16) and (17) for the asymptotic coefficient (as they generate different 1/n corrections).
The three methods yield different ways to measure the fact that n is not infinity. We take the largest of the three as our estimate of the error. The scale variation will yield the largest error that we then symmetrize to its maximum value. Our best values for N X , and the associated error, can be found in Table 1. They have been computed with NNNLO accuracy for Z B and Z GLS (and with NNLO for Z EJ , except for n f = 0, where all three sum rules are equal), after conformal mapping, at x = 1. The scale dependence of the results, as well as the convergence, are shown in Fig. 1 for some selected values of n f . In view of the figures, we believe that our error is conservative. To reinforce this conclusion, we also consider different possibilities for the expression of the asymptotic coefficient after the conformal mapping: either to do the conformal mapping to the original asymptotic expression, or to directly write what should be the new expression  Table 1: Values of the infrared renormalon residue of the Bjorken, Ellis-Jaffe and GLS leading-twist Wilson coefficient C X . For n f = 0 the three are equal. Note that the determination for the Ellis-Jaffe sum rule is obtained at one order less than the other normalizations (except for the case n f = 0).
for the asymptotic expression (changing the location of the singularity). The differences are well inside the error. The error we quote in this paper is bigger than the error quoted in Ref. [3]. The reason is that we are more conservative here. If we had used the same method, we had obtained a larger error in that reference. Therefore, as we already said, we consider our error estimate to be conservative. On the other hand, the central values are quite close, except for large n f , where we now observe a faster trend to get smaller values of the normalization as we increase n f .
In Fig. 1, we can see that the scale dependence becomes smother as we go to higher orders (for n f = 0 or 3). The convergence depends on the number of flavors. It is optimal for n f = 3, which actually happens to be the most interesting case from the physical point of view, and it deteriorates for large n f . Also, for large n f , the value of the normalization is consistent with zero. This fits with the picture that the renormalon is less important when the number of flavors grows and one can reach to the point where the infrared renormalon disappears. For the Ellis-Jaffe perturbative series with n f = 0 (which is known to one order less), the same discussion applies but the flip of sign of the normalization happens at n f = 4, and for n f = 6, it starts to blow up.
The values of Z B and Z GLS are consistent with each other within errors. This is consistent with the interpretation that the light-by-light term does not contribute to the renormalon, as it was done in Ref. [34].
We are then able to give some estimates for the coefficients of the perturbative series. We provide them in Tables 2, 3, and 4. We should stress that our numbers incorporate the right asymptotic behavior, which is not the case for large-β 0 estimates. For the Bjorken and GLS sum rules, they were calculated in [36].
Previous estimates of the coefficients of the perturbative series can be found in the literature for the Bjorken (GLS) sum rules [34,35,37,38]. Nevertheless, those estimates do not profit from the knowledge of the O(α 4 ) coefficient, which was not known back then.    Table 1.     GLS for ν = Q and for different number of flavors. We do not display the column with n f = 0 since the numbers are equal to the Bjorken case. We use the expression from Eq. (13).

B
As we have already mentioned, the divergent behavior of the perturbative series is regulated using the PV prescription for the Borel sum [see Eq. (7)]. The exact expression of the Borel transform is unknown. Therefore, one has to use approximations. In this context, it is of paramount importance to have a parametric control on the error with exponential accuracy. Consequently, we apply the hyperasymptotic approach developed in [20,24,25]. The Borel sum can then be split into a partial sum truncated at the minimal term, plus the leading terminant, and plus a leftover of the perturbative expansion (the perturbative expansion is not known with high enough accuracy to go beyond that):  Table 4: Renormalon-based estimates of the perturbative coefficients C (s) EJ for ν = Q and for different number of flavors. We do not display the column with n f = 0 since the numbers are equal to the Bjorken case. We use the expression from Eq. (13). where In the weak coupling limit, this expression can be approximated by where is Eq. (13). N max = 3. By default, we will take N P = N max . This will fix c. For illustration, we show the relation between N P and c in Fig. 2. In principle, we would like c to be small or of order 1. For Q around 1 GeV, N P = 3 is okay. For larger Q one should ideally take a larger value for N P . As N P = 3 is the maximum value known nowadays, it is not possible to take a largest value. In our analysis, we will always take N P = 3. We have checked that if we take N P = 2, and the corresponding c, the variation is minimal.
Note the prefactor multiplying Ω in Eq. (20). This characteristic is novel compared to earlier analyses using the hyperasymptotic expansion, and resums the large ln Q logarithms associated with the anomalous dimension. Note also that, in practice, Ω has a hidden dependence in Q through c. In the weak coupling approximation, this hidden dependence is encoded in K B . We show the dependence of K on Q 2 in Fig. 3. This makes that Ω effectively does not behave (leaving aside the anomalous dimension) as α(µ)Λ 2 MS . For numerics, we will use Eq. (23). The difference with the weak coupling expressions in Eq. (24) is sizable for Q 2 < 1 GeV. We will discuss this issue further later.
We will combine the statistical and systematic errors of each experiment (in case they are given separately by the experiment) in quadrature. Nevertheless, for the more recent (and more precise) datasets from JLAB [45][46][47][48], some correlations are expected. To account for those, we follow the same procedure as in Refs. [47,48] and combine part of the systematic error with the statistical one in quadrature, such that the reduced χ 2 is fixed to be one. The remaining systematic error is taken to be completely correlated.
In the left-hand side of Eq. (1), target-mass effects have been included using the Nachtmann variable [49], see also [50,51]. They read where N = p, n (remember that M B 1 = M p 1 − M n 1 , and accordingly Γ B 1 = Γ p 1 − Γ n 1 and so on), ξ = 2x/ 1 + 1 + 4m 2 N x 2 /Q 2 is the Nachtmann scaling variable, m N is the nucleon mass. The quantity M 1 is the first Nachtmann moment of g 1 that absorbs all the target mass corrections, ∼ (m 2 N /Q 2 ) n , and where a N 2 is the target mass correction given by the x 2 -weighted moment of the polarized leading-twist g 1 structure function, and d N 2 is a twist-three matrix element given by For these quantities we use the values given in [47] (a recent lattice determination can be found in [52]): d p−n 2 = 0.0080(36) , a p−n 2 = 0.031 (10) .
Note that the elastic contribution has to be included in the experimental numbers in order the sum rule to be fully inclusive, i.e., The empirical parametrization of the elastic form factors is taken from the recent paper [53]. For large momentum this contribution is completely negligible. The difference starts to be sizable for Q 2 < 2 GeV 2 . For Q 2 < 1 GeV 2 , it becomes very important.
In comparing theory and experiment, we will letf PV 3 to be a free parameter, which will then be fitted to the experimental data. The value of α s (M z ) = 0.0179(9) is taken from the PDG [54]. After running down to scales of the order of 1 GeV using [55], it yields We assess the error associated with the parametrization of the elastic form factors by performing the same fits using the older parametrization [56]. We then take the difference as the associated error. It will be small in comparison with the other experimental errors.
In our final fits, we will neglect the error associated to the coefficients g A , a 2 and d 2 , as their impact is negligible in comparison with other uncertainties.
We now consider the theory part of the problem. We first illustrate the problem of convergence of the perturbative series by drawing the perturbative series in the MS scheme at different orders in α s in Fig. 4. The perturbative series has a relative good convergence. However, this convergence deteriorates when we approach to low energies. We also see how the perturbative theoretical result diverges from the experimental numbers at low energies. This is a reflection that, at these scales, we need a description of the experiment with exponential accuracy. In other words, power corrections have to be incorporated, and the perturbative sum needs to be evaluated with exponential accuracy. We do so by using Eq. (6) with C PV B evaluated using the expression in Eq. (20), which has exponential precision. We can see that the agreement with the experimental data improves. This is specially so around the 1 GeV region. We can also see a qualitative change in the figure for scales smaller than 1 GeV (see the black continuous line). Nevertheless, this is an artifact of using Eq. (23) for the terminant. If we instead use the weak coupling expressions for it one has in Eq. (24), the behavior is different (see Fig. 5). This illustrates that we can not trust our predictions below 1 GeV. If we eliminate the terminant and only keep the dimension two condensate (blue dashed line in Fig. 4), the agreement deteriorates. The black line has been obtained from a fit lettinĝ f PV 3 to be free using data fulfilling that Q 2 ≥ 1 GeV 2 . As just mentioned above, we have used Eq. (6) with C PV B evaluated with the exponential precision using the expression in Eq. (20) with N P = 3. For Ω, we use Eq. (23). The value obtained is listed in Table 5.
We now turn to the discussion of how robust this result is, and to determine the error. In Table 5, we disclose in detail the error budget.
The experimental errors have already been discussed above. We have separated them into the statistical, systematic, and the one associated with different parametrizations of elastic terms. By far, the dominant error is the systematic error associated with the JLAB data. We also include the error due to α s (see Eq. (32) in Table 5), which is also much smaller than the systematic error. We then combine them in quadrature to give our final experimental error.
The theoretical error is generated by our incomplete knowledge of C PV B and the incomplete knowledge of the Wilson coefficient of the dimension-two term of the OPE. Theoretically, there is also the error associated with higher order terms of the OPE. These are of O(1/Q 4 ). This is parametrically smaller than the other errors associated with our incomplete knowl-   (1) We allow for a variation of N B according to the error given in Table 1.
(3) Approximating Ω by its weak coupling expansion at leading order. In other words, working with the second equality in Eq. (24).
We now discuss these methods to estimate the error in more detail.
(1) The error of N B is due to our incomplete knowledge of the perturbative series, and it can be, to some extent, be mixed with the O(α 3/2 s e − 4π β 0 αs ) corrections of the terminant.
(2) The scale dependence is a reflection of the incomplete knowledge of the terminant and of the Wilson coefficient of the dimension two term of the OPE. If we keep the running of the coupling constant consistently with the precision of the computation, there would not be scale dependence. We implement the scale dependence by considering the running of the coupling with the maximal known accuracy.
(3) Measures the size of the O(α 3/2 s e − 4π β 0 αs ) correction to the terminant. An alternative way to measure these corrections is to approximate Ω using the first equality of Eq. (24). This indeed produces larger errors but still much smaller than the error associated with N B , or ν. Another alternative consists of changing N P from 3 to 2. This changes c. If the terminant were known with infinite precision, the dependence in c would disappear, in the sense that the dependence on c of the terminant would cancel with the dependence on c of the finite sum. Therefore, the variation of c assesses such dependence. The effect happens to be very small.
Let us now come back to the issue of the O(1/Q 4 ) corrections. As we have already mentioned, they are parametrically subleading than the errors we have considered. Therefore, we will neglect them. 2 Still, fits with different energy ranges may give indirect hints of their possible importance. Related to this discussion, a potentially worrisome issue is the size of the elastic terms. For large momentum this contribution is completely negligible. The difference starts to be sizable for Q 2 < 2 GeV 2 . For Q 2 < 1 GeV 2 becomes very important. In this respect, it is important to mention that the short distance behavior of the elastic form factors is expected to produce corrections to the Bjorken sum rule with powers bigger than 1/Q 2 [58] (for the present discussion, see also [53,56]). Therefore, one may argue that they are a measure of power corrections of order 1/Q 4 or higher. Nevertheless, such prediction is not in such a robust status as fully inclusive sum rules, as they somehow rely on local quark-hadron duality (one expects them to fail in the large N c limit). Therefore, we will not consider this as a source of error, although it would deserve further studies.
For the final estimate of the theory error, we consider items 1-3 as different estimators of the incomplete knowledge of the OPE. Even though, there is some overlapping between them (they measure similar things), we will consider them as independent, and combine them in quadrature to give our final theory error off PV 3 . The sensitivity to the theoretical and experimental error may vary depending on the range of energies used for the experimental data. Therefore, we use different datasets for the fits to see the dependence on the dataset. The different datasets we consider are the following: set I: Q 2 ≥ 1 GeV 2 and set II: Q 2 ≥ 2 GeV 2 . The set I has 31 points, and the set II has 19 points (and are less precise). We refrain from considering smaller values of Q, since the dependence on the specific parametrization of Ω becomes important. As a general trend, when we increase the number of points, the statistical and systematic errors off PV 3 decrease. On the other hand, the error associated with the parametrization of the elastic form factors moderately increases but to a lesser extent. On the theory side, we observe the following dependence on the dataset: As expected, the fit is more dependent on N B , as we go to smaller energies (set I). We find a large dependence on ν with the fit to the set I if we lower ν since we take ν 2 = Q 2 /2. On the other hand, we find very little dependence if we take higher values of ν. For set II the error is more symmetric and, overall, smaller. Finally, the error associated with Ω is very tiny, and basically zero when using set II for the fit.
For our final numbers we use the results obtained with the dataset I. This dataset yields smaller errors. For convenience, we display the results obtained for the different datasets: Set II :f PV 3 × 10 3 = 123 +417 −417 (exp) +17 −24 (th) = 123 +417 −418 GeV 2 .
For the final error, we combine the experimental and theoretical error in quadrature. The experimental error is the dominant one. As we have already mentioned, we take the numbers obtained from the fit to the dataset I as our final numbers. Our final number forf PV 3 reads: We also give a number forf PV 3,MS [defined in Eq. (9)]. Its error is overwhelmingly dominated by the error off PV 3 , being the error associated to Λ  To illustrate the quality of the fit and the impact of the error, we plot our final results in different ways. In Fig. 6, we plot our final predicted curve for Γ B 1 , including the error associated withf PV 3 . In Fig. 7, we directly show our prediction off PV 3 (and its associated error) and how it compares with the experimental data. We observe that the fit (the experimental data) approximately follows a constant value forf PV 3 , as expected from the use of a summation scheme of the perturbative series that is scheme and scale independent.
The results obtained in this paper can be compared with the more traditional renormalon subtracted (RS) scheme [59], as used in Ref. [3]. We obtain f RS (1; GeV) = 0.152 GeV 2 , f RS (0.8; GeV) = 0.105 GeV 2 using the dataset I and f RS (1; GeV) = 0.260 GeV 2 , f RS (0.8; GeV) = 0.239 GeV 2 using the dataset II. To translate these determinations of f 3 from the RS scheme to the PV scheme, one has to use the following expression (N P = N max and Q = ν): (37) Note that the left-hand side is independent of Q and ν f . If we take ν f = Q = 1 GeV in the right-hand side, then the fits with the RS scheme give very similar numbers to the ones we have obtained in this paper. Using set I, one obtainsf PV 3 = −0.003 GeV 2 , and with set II, one obtainsf PV 3 = 0.141 GeV 2 . This shows the stability of the fit to using a different scheme for handling the renormalon. Within errors, the result is stable to variation of ν f and Q, as shown in Fig. 8, for dataset I.
We can also compare with the result obtained in [3] using the RS scheme. We only aim to make a qualitative comparison. In that paper, it was obtained f 3,RS (ν f ) = −0.124 +0.137 −0.142 by a global fit to the Bjorken, Ellis-Jaffe, and the GLS sum rules to existing data back then with Q 2 ≥ 0.66 GeV 2 . A pure fit to the very same data to only the Bjorken sum rule (and treating slightly different the strong coupling) yields -0.102 GeV 2 . If we use the new experimental data, the new value of N B , and the new known coefficient of the perturbative expansion, then the result does not change much, one obtains -0.091 GeV 2 . Note that all the experimental systematic error has been combined with the statistical one in quadrature, since the reduced χ 2 is 1.6, bigger than 1. Actually, this is the reason the result is similar to the one obtained in [3], as it makes that the new data weights similarly to the old one. If one keeps part of  the systematic error as independent, then the new data weighs more, and shifts the central value (but then the rest of the error goes into the systematics, not affecting the comparison within errors). This effect is seen in the fits with data constrained to be Q 2 ≥ 1 GeV 2 or Q 2 ≥ 2 GeV 2 . Taking into account the error of the determinations and of the conversion, the difference is of the order of one sigma, which we consider reasonable. Moreover, in this paper, we prefer to play a more conservative attitude, in view of the different behavior of the terminant (depending on the approximation used for it) for Q 2 ≤ 1, and restrict to fits with Q 2 ≥ 1 GeV 2 .

Conclusions
In this paper, an updated determination of the normalization of the leading infrared renormalon of the Bjorken, GLS and Ellis-Jaffe sum rule has been obtained. Compared with the evaluation of Ref. [3], this evaluation profits from new terms of the perturbative expansion that are known now, as well as from a more optimal way to pin down the renormalon, as highlighted in Ref. [29]. The results are summarized in Table 1 Using these results, estimates of higher order terms of the perturbative series are also given (see Tables 2, 3, and 4). For the first time ever, the hyperasymptotic approximation has been applied and validated directly in experimental data (up to possible concerns on the extrapolations to x → 0 of the dispersive integrals of the sum rules, where there is no direct experimental data). We computed the Bjorken sum rule with hyperasymptotic precision by including the leading terminant, associated with the first infrared renormalon. One advantage of the present method is that it minimizes the dependence on the normalization of the renormalon. Another major advantage is that it allows us to obtain the "real" size of the nonperturbative corrections, in the sense that the nonperturbative correction scales as powers of Λ QCD .
After fitting the experimental data to the OPE theoretical expression lettingf PV 3,MS as a free parameter, we obtain a reasonable agreement with experiment for Q 2 ≥ 1 GeV 2 with We emphasize that this result is independent of the renormalization scale and scheme used for the strong coupling. The error off PV 3,MS is overwhelmingly dominated by the error off PV 3 , and the latter is overwhelmingly dominated by the experimental error. There is some correlation in the experimental error of the most precise data from the JLAB experiments [45][46][47][48]. We follow the methodology of [47,48]. This correlation yields the most important source of error of the present analysis from the experimental side. Any improvement in this respect will immediately lead to a reduction of the errors of the numbers obtained in this paper.
As we can see in Figs. 4 and 7, the fit complies with expectations. The difference between the leading twist contribution and the experimental result is reasonably fitted by a quadratic correction with a Q-independent value off PV 3 . It approximately follows a straight line (see Fig. 7), as expected. Its value is small and consistent with zero within one sigma significance.
The introduction of the terminant makes the result more stable. The inclusion of the terminant introduces a qualitative change in the perturbative behavior around the 1 GeV region, making it much closer to the experimental figure. One may think that the terminant is basically equivalent to the nonperturbative condensate (with the same anomalous dimension), albeit α(µ) suppressed. Note, however, that, effectively, Ω has an implicit Q dependence, since it depends on c, which is fixed such that N P = 3.
The theoretical error off PV 3 is dominated by the scale variation and N B and ν, and it is bigger when fitting to the dataset I than to the dataset II. In any case, it is still much smaller than the experimental error. All these theoretical errors would benefit from higher order computations. To know the anomalous dimension of the condensate to a higher order would also be of help.
We let to future research the application of the hyperasymptotic approximation to other sum rules. gramme under Grant Agreement No. 824093. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya.