Optimized QCD determination of the polarized Bjorken sum rule

We discuss and improve the determination of polarized Bjorken sum rule $\Gamma_1$ in two different ways: phenomenological determination of $\Gamma_1^\text{exp}$ by means of integration of the spin structure function $g_1$ within the truncated moment approach; and the optimization of the truncated perturbation series for the coefficient function of the corresponding DIS process by means of the renormalization group. Both obtained values $\Gamma_1^\text{exp}$ and $\Gamma_1^\text{th}$ become closer to each other when processing COMPASS data. Meanwhile, this optimization approach is universal and can be applied to any of the DIS sum rules.

The theoretical ground for the deep inelastic scattering (DIS) is the "factorization theorem" for hard exclusive processes at large transferred momentum q, −q 2 = Q 2 ≫ P 2 = m 2 h , where m 2 h is the hadron mass. In this framework the DIS of leptons on nucleons can be derived from basic assumptions of quantum field theory and the scaling violation for the DIS is described by the renormalization group (RG) equation providing a direct test of QCD. In terms of parton distribution functions (pdfs) f p/h (x, µ 2 ) on the Bjorken x variable, 1 > x = Q 2 /(2P q) > 0, the idea of formalism describes the scaling violation with pdfs being objects of nonperturbative nature, while their evolution with scale µ 2 up to the scale Q 2 is governed by the well-known DGLAP equations [1][2][3][4]. This is a beautiful manifestation of the short-and long-distance factorization in hard processes. Experimental verification of the DIS sum rules faces however the difficulty that in any realistic experiment one cannot reach arbitrarily small values of the Bjorken x, x x 0 ≡ Q 2 min /(2(P q) max > 0). This is a serious obstacle in the determination of the sum rules that involve the Mellin moments of pdfs, i.e. integrals of pdfs weighted with x n over the whole range (0,1) of x. We propose a method based on the truncated (cut) Mellin moments (CMM) approach which is realized in the range (x 0 , 1) and can mitigate this x 0 problem. The CMM of the parton densities were applied to QCD analysis almost twenty years ago [5][6][7][8] and then developed and successfully compared with the experimental DIS spin data [9][10][11][12][13][14][15][16][17]. The CMM approach provides not only a natural framework of DIS analysis in the restricted kinematic region of x x 0 but also allows one to effectively study a very low x regime of the sum rules. In [16,17] we elaborated the generalized Bjorken sum rule (gBSR) and demonstrated how, with the help of this construction, one can determine the Bjorken sum rule (BSR) Γ 1 = 1 0 g 1 (x)dx from experimental data in the available restricted x region. This paper is divided into two parts: the analysis of BSR data within the framework of the optimized CMM approach, which we then combine with optimization of the perturbative expansion for the coefficient function C Bjp (α s ) based on the solution of RG. In the next section, we present the generalized Bjorken sum rule and its advantages in determining the Bjorken sum rule value from experimental data for a constrained range of x. In Sec. III, we estimate Γ 1 from the recent COMPASS data [18,19] using the two-point and multi-point versions of the gBSR approach. This latter allows us to incorporate different uncertainties for each kinematic bin of the experiment that can be interpreted as optimization of integration. Another kind of optimization applied to the perturbation series for the RG invariant quantity BSR is considered in the next sections. The criteria for the optimized analysis of the C Bjp (α s ) within the perturbation QCD and the RG are discussed and fixed in Sec. IV. Following these criteria we find numerically the admissible domains for the corresponding new normalization scales µ 2 (µ 2 = Q 2 ). In Sec. VI, the results of optimization are presented and discussed. Finally, we optimize the perturbation expansion for the BSR at the COMPASS kinematics and compare one with the result obtained here with the help of the generalized Bjorken sum rule. The important technical issues are presented in two Appendices.

II. GENERALIZED BJORKEN POLARIZED SUM RULE
The well-known Bjorken polarized sum rule [20,21] enables one to study the internal spin structure of the nucleon. Referring to the first moment of the nonsinglet spin dependent structure function g 1 , BSR is a powerful tool in understanding the proton helicity decomposition in QCD. Here we present a novel method for determining the BSR at experimental constraints. To this aim, we construct the generalized Bjorken sum rule being one of the results of the CMM approach [14,16]. The generalized truncated moment f (z, n, ω), obtained as a Mellin convolution of the pdf f with any normalized function ω(x), obeys the DGLAP evolution equation with the rescaled kernel [16]: In the case of the non-singlet polarized structure function g 1 and for n = 0, one obtains the generalized structure function g 1; ω , with the same DGLAP evolution kernel as for g 1 , namely P (y). In this way, we define the truncated generalized Bjorken sum rules (gBSR), Γ 1; ω (x 0 ), which is equal to the ordinary BSR in the limit x 0 → 0: We showed in [16] that the truncated generalized BSR Γ 1; ω (x 0 ) approaches the limit Γ 1 (0) more quickly than the truncated ordinary BSR, Γ 1 (x 0 ), This property of the gBSR allows very effective determination of Γ 1 (0) with the use of a few first orders of Taylor expansion: Below we present the main expressions for the successful determination of Γ 1 (0) from the experimental data in the two-and multi-point gBSR approaches. For more details, see [16].

A. Two-point approach
One can estimate the value of Γ 1 (0) from the smooth extrapolation of the truncated moments Γ 1; ω (x 0 ) in x 0 , where Γ 1; ω (x 0 ), Eqs. (3) and (4), is based on the simple sign-changing normalized function ω(x) depending on three parameters z 1 , z 2 , A, Here the ω-model parameters are 0 < x 0 < z 1 < z 2 and 0 < A for the sign change. It is convenient to rewrite the approach in terms of the experimental parameters x min and r, where x min = x 0 /z 2 denotes the smallest x available in the experiment and r = z 1 /z 2 , x min < r < 1 is the ratio of two experimental points from the set 0 < x min ≡ x 1 < x 2 < · · · < x max < 1. The idea of the gBSR leads to a "shuffle" of the initial function g 1 in x variable and the truncated generalized BSR Γ 1; ω (x min , r). Indeed, substituting Eq. (8) for the ω-model into Eq.(4) for Γ 1; ω and integrating, one obtains [16] that saturates the limit Γ 1 (0) more quickly than the ordinary Γ 1 (x min ), In other words, the use of the gBSR "mimics" the extension to lower values of x in the experimental kinematic regime. We have checked that the BSR limit Γ 1 (0) can be determined very effectively with the use of the first order of Taylor expansion independently of the small x behavior of g 1 [16]. For a number of cases even the zero order of Taylor expansion can give a satisfactory result. As we will see, the analysis of COMPASS data is that case. Figure 1 illustrates the features of Γ 1; ω , where we plot Γ 1; ω (x min , r) Eq. (9), as a function of x min for "constant" and "quasi-linear" behavior corresponding to the special values of the model parameter A, A = A 00 or A = A 0 , respectively. The A 00 (x, r) ensures simultaneous vanishing of the first and second derivatives, Γ ′ 1;ω (x) = Γ ′′ 1;ω (x) = 0 while the A 0 (x, r) ensures separately the second condition, Γ ′′ 1;ω (x) = 0. These special values of the parameter A take the following forms [16]: We show also the standard truncated BSR Γ 1 (x min ), Eq. (10), (lower thick solid curve).
In the case of A 00 (x, r), the requirement Γ ′ 1;ω (x) = Γ ′′ 1;ω (x) = 0 relates the model ratio r and the argument x via the parameters β and γ of the auxiliary function g f it 1 fitted to the experimental data, namely x is a non-negative solution of the quadratic equation: (βγ − 1)x 2 + (β + 1)(r + 1)x + (β/γ − 1) r = 0 (14) and, inversely, Then, the total BSR limit Γ 1 (0) can be determined in the zero or the first order approximations, respectively: with A 00 (x min , r) given by Eqs. (11), (15) and A 0 (x min , r) given by Eq. (12). In our test plot, we choose g 1 ∼ x −0.4 (1 − x) 3 (1 + 5 x) reflecting theoretical predictions on the small x behavior of the nonsinglet structure function g 1 [22]. We fix A 00 at x = 0.04 what, according to Eqs. (11) and (15), implies r = 0.8 and A 00 = 5.2. Using the same r, we obtain for chosen value of x = 0.02 A 0 = 5.8. One can see that saturation of  The presented above idea of the generalized BSR can be modified to analyze the experimental data incorporating different uncertainties for each measurement. Usually, the results for g 1 are extracted from the data for kinematic bins {Q 2 i , x i } with uncertainties ∆(g(x i )). In order to compute moments of g 1 and verify the BSR, bins must be evolved to a common scale Q 2 . Then, the experimental bins {x i } n 1 can be used to construct multi-point weight function ω in Eq. (1), where This multi-point ({z i , A i } n 1 ) model with z 1 < z 2 < . . . < z n is a generalization of the previous two-point ({z 1 , z 2 , A}) one, Eq. (8), and leads to a new presentation for Γ 1; ω , In order to use the generalized BSR to determine Γ 1 (0), one only needs to construct the set of {A i }. Choosing the most appropriate weights w i and hence A i , one is able to tune the analysis to the real experimental constraints. In terms of the experimental parameters the multi-point gBSR, Eq. (20), has the form where A natural way to implement the experimental uncertainties of g 1 into the model of Γ 1; ω is to use weights to increase the contributions of the data with smaller uncertainties, for instance, taking ones inversely proportional to the relative uncertainties ∆(g 1 )/g 1 at each x i . We choose the weightsw i ,w i ≡ w n−i , where U i is the relative statistical uncertainty of g 1 (x i ), reflecting the increase of the experimental uncertainties of g 1 with decreasing x.
Then we arrive at the first order approximation for Γ 1 (0) in the multi-point case which is a generalization of the two-point one, Eqs. (17), (12):

III. ANALYSIS OF COMPASS DATA
In this section, we determine the BSR limit, Γ 1 (0), from the recent COMPASS data [18,19,23], where x min = 0.0036, using the two-point and multi-point versions of the gBSR approach described in the previous section. We follow the method utilizing Eqs. (11), (15), (16) and Eqs. (12), (17) for the two-point zero and first order approaches, respectively. For the multi-point analysis, incorporating the experimental uncertainties in each measurement of the spin function g 1 , we use the first order approximation, Eqs. (22)- (27). The contributions to Γ 1 (0) in all our approaches are calculated from the experimental data evolved to Q 2 0 = 3 GeV 2 . In our analysis we used two auxiliary, global for all x, fit functions. The first one was used to evolve g 1 (x i , Q 2 i ) to the common value Q 2 0 . The second auxiliary function, fitted to the data after evolution, was used only to calculate the derivatives g ′ 1 in the equations for A 0 , to find r in A 00 and to obtain small rest contributions to the BSR, . The latter arises from the relation between x and r in the zero order approximation, Eq. (15), when we arrive at x i /r not being the experimental point x k . Our fit function for Q 2 0 , where we did not make any assumption on the normalization, has the form Instead of this global fit, one can use a series of polynomials g i 1 (x) locally fitted for each range (x i ; x i+1 ) additionally enhancing the impact of the experimental measurements themselves on the final estimations of Γ 1 (0). Our results are presented in Fig. 2 and Table I parameter A: two of them are A 00 (x, r), Eqs. (11), (15), and one is A 0 (x, r), Eq. (12), together with the truncated ordinary BSR, Eq. (10). The arguments x i in A 00 and A 0 are chosen from COMPASS x points, shown in column 2 of Table 1: x 1 = 0.0036, x 2 = 0.00459, etc. The set of experimental data points x 1 , x 2 , ..., x n determines the ratio parameters r. For the two-point IAPX plot, the r is calculated for x 1 and the closest next x i , r = x 1 /x 2 , while for the two-point 0APX plots, the ratios r are given by Eq. (15) at two values of x min : x 6 and x 7 . One can see from Fig. 2 that, similarly to our test plot in Fig. 1, Γ 1; ω (x min ) obtained from COMPASS data approaches the BSR limit Γ 1 (0) visibly quicker than the initial BSR Γ 1 (x min ). The quasi-constant behavior of the gBSR begins far from the small x limit contrary to the standard BSR. This is a strong motivation for using in our approach apart from the first order also zero order approximation to find the BSR value. In Table I we present our estimations for the Bjorken sum  (0) for Q 2 0 = 3 GeV 2 in the generalized BSR approach from the COMPASS data [18,19,23]. The experimental {xi} n 1 set is shown in column 2. We present the two-point zero order (0APX) results, Eqs. at Q 2 0 from reanalyzed COMPASS data is displayed in the last row.  [19] rule Γ 1 (0) for Q 2 0 in the two-and multi-point approaches as well. For both cases we use the first order approximation Γ IAPX 1 , which is a very effective method, particularly when experimental x min 0.2 [16]. The quasi-constant behavior of the gBSR for COMPASS data, see Fig. 2, should ensure the applicability not only of the first order approximation (IAPX) but also the zero order one (0APX). Therefore, within the two-point case we present additionally the zero order analysis. This is shown in column 4, where Γ 1 (0) is obtained from Eqs. (11), (15), (16) for all admissible experimental x points related via ratios r, Eq. (15), (column 3). In turn, using the first order approximation we proceed step by step from x 1 = 0.0036, true lowest x for COMPASS, to larger values of x min (column 2) simulating the experimental restriction on the x-range. For the two-point case, Eqs. (12), (17), we choose two closest x points, x min and the next to it, obtaining the corresponding ratio r = x i /x i+1 (column 5) and the result for Γ 1 (0) (column 6). The multi-point estimations of Γ 1 (0), Eqs. The multi-point estimations and the two-point IAPX approach give very similar values of Γ 1 (0) up to x min ≈ 0.1. The zero order two-point results, which can be found only for the first seven experimental x points, are also in very good agreement with the rest of our estimations of Γ 1 (0) for this x region. All our results agree well with the value provided by COMPASS from reanalyzed data [19] (last row of Table I). One can also see that for x min ≈ 0.2, for which g 1 suffer from larger experimental uncertainties, Eq. (25), the multi-point approach incorporating these uncertainties via suitable weights provides a more stable estimation of Γ 1 (0) than the two-point ones. We can conclude that this version of our approach is universal, reliably fixes Γ exp-opt 1 at Γ exp-opt and looks most promising in determining the DIS sum rules.

IV. RENORMALIZATION GROUP ANALYSIS OF QCD PT SERIES FOR BSR
The estimates of Γ exp 1 = Γ 1 (0) discussed before and presented in Table I for COMPASS measurements can be compared with the QCD predictions Γ th 1 (Q 2 ) for the BSR. The latter includes radiative corrections obtained in the MS -scheme in O(α n s ), n = 1, 2, 3 and 4 approximation in [24][25][26] and [27], respectively, and the nonperturbative higher twist effects (HT), Here in the RHS of Eq.(30) C Bjp (a s ) is the leading twist nonsinglet coefficient function of the polarized BSR, |g A /g V | = 1.29 ± 0.05 stat ± 0.1 syst [19], while µ p−n 4 is the scale of the first power correction to the high twists whose effects become essential in the small/moderate Q 2 region. Below we shall investigate the QCD radiative corrections to C Bjp (a s ) based on the renormalization group transform.
A. The problem of PT optimization for coefficient function C Bjp (as) The perturbation expansion for C Bjp (a s ), the known coefficients of which are presented in Appendix A, reads: One can see that the convergence in (32a) becomes expectedly worse with order growth. Below we are going to perform optimization of expansion in (31a) choosing an appropriate new normalization scale µ → µ ′ following the renormalization group (RG). The value of the truncated series for BSR in (31) starts "breathing" under the mentioned variation of µ around the norm scale µ = Q; this is the inevitable effect of the truncation that we shall use for optimization. The corresponding approach goes back to the generalization of the Brodsky-Lepage-Mackenzie (BLM) [28] method suggested in [29,30] for the RG invariant quantities. The approach is based on the {β}-expansion for the PT coefficients [29] 1 , see here Appendix A for details, and allows managing their values in great detail. An alternative approach to the PT optimization named PMC is elaborated and applied to BSR in [32], we will return and mention its results in Sec.VI. However, for a kind of practical tasks it looks like that one does not need to use all such detailed information. Here we shall try to avoid the details of {β}-expansion, instead we shall apply first rather brute and direct the approach to make smaller the contribution of radiative corrections in (31a) for BSR. In other words, we shall reform 4 successive orders of radiative corrections in the parentheses in the RHS of Eq.(32a) to make its sum in Eq.(32b) minimal following the RG transform. In the next subsection, we will remind the reader of appropriate elements of the corresponding formalism that can be applied to any RG invariant quantities, see [29,30,33].

B. General basis of optimization
Re-expansion of the running couplingā s (t) = a s (∆, a ′ s ) and its powers in terms of the logs, t = ln(µ 2 /Λ 2 qcd ), t − t ′ = ∆ = ln µ 2 /µ ′2 and the new coupling a ′ reads which is the way to write the corresponding RG solution for a(t) through the operator exp (−∆β(a)∂ a ) [. . .] | a=a ′ (see [29] and refs therein). The shift of the logarithmic scale ∆ in its turn can be expanded in perturbation series in powers of a ′ β 0 where the argument of the new coupling a ′ depends on t ′ = t − ∆. It is sufficient to take this renormalization scale for a ′ , which corresponds to the solution in the previous step of iteration, rather than to solve the exact equation a s (t − ∆(a ′ s )) ≡ a s (t ′ ) = a ′ s . Re-expansion a s in terms of a ′ s and ∆ i leads to rearrangement of the perturbation series for the RGI quantity C Bjp = i a i s c i → i (a ′ s ) i c ′ i , the elements in the r.h.s. can be expressed as c ′ i = B ij c j where B ij is a triangular matrix. In this notation C Bjp transforms to The elements B ij appear as a composition of transforms in Eq.(33) and Eq.(34), In the square brackets below we write the elements of the triangle matrix B explicitly: Recall that the standard BLM [28] is based on the decomposition  (36).
In this way, we have constructed a device to improve perturbation expansion by means of Eqs. (34,36); the prize we should pay is the demand to control simultaneously both the expansions for ∆ in (34) and for c i in (36). We shall discuss the admissible domains of {∆ 0 , ∆ 1 , . . .} that satisfy the above conditions in the next subsection.

V. THE ADMISSIBLE DOMAINS OF {∆} PARAMETERS
Here we consider the application of the general scheme to the actual observable quantity C Bjp in order to obtain the truncated PT series (35) with small corrections at the appropriate convergence. At the end we will fit the appropriate values for the parameters {∆ 0 , ∆ 1 , ∆ 2 }. Let us introduce the notation A ′ = β 0 a ′ s for the scaled renormalization group solution for the brevity and simplicity. To satisfy the reliability of the PT expansion, we set natural inequalities for its successive terms shown in the items (i -iii) below: (i) for ∆(t) in (34) it reads (ii) for PT expansion in (35) we set a similar condition with respect to c ′ We consider these Eqs. (37,38)  (iii) To fix the PT domain of applicability, we put for the logarithmic variable t ′ = t − ∆(t ′ ) the appropriate lower bound at µ 2 0 ≃ 1 GeV 2 that corresponds to t µ0 ≃ 2.3 We shall scan t in the practically interesting interval 2.  Fig. 3. The constraints in Eqs. (38,39) are much more restrictive for the parameters in the right half plane at ∆ 0 > 0; therefore, the corresponding domains are significantly smaller than on the left one. For all that the "BLM value" {∆ 0 = 2, ∆ i = 0}, see Eq.(A.3b) in App. A, belongs to the admissible domain. Obviously, the larger is t, the larger (and lighter in color in the figure) is the corresponding admissible domain. Let us mention first the black triangle that corresponds to the conditions c ′ 2 = c ′ 3 = 0; it belongs to the region of applicability and was predicted in [30], The conditions in Eqs.(40) correspond to the new norm scale µ ′2 = µ 2 exp −∆(a ′ s ) = −1.56 + 0.396β 0 a s (µ ′2 ) > 0.21µ 2 , see Fig.1(Right) in [30]. If one imposes the fourth term in the condition (38) for {∆ 0 , ∆ 1 } parametrization, then the domains {∆} in the right half plane become slightly vuggy and get a cut along the ∆ 0 direction. This effect can be seen from the corresponding cross section at ∆ 2 = 0 of 3D admissible domains presented in Fig.4. 3D . Two examples of admissible 3D-domains {∆} 2 0 in the order O(a 4 s ) are calculated at t = 5 (µ 2 ≈ 15.0 GeV 2 ) and at t = 8 (µ 2 ≈ 301.0 GeV 2 ) and shown in Fig.4; see the left and right panels respectively. Again, the larger is t, the larger is the corresponding admissible domain. Pay attention to the difference of the axes scales in the left and right panels. One can notice that the point {∆ 0 , ∆ 1 , ∆ 2 } corresponding to the conditions c ′ 2 = c ′ 3 = c ′ 4 = 0, black ball there, is not contained in the admissible region. The reason is that this condition at t 9 contradicts inequality (37) for perturbation expansion of ∆. We shall discuss the task of optimization within this frame for C Bjp in the order O(a 4 s ) in the next section.

VI. THE RESULTS OF RENORMALIZATION GROUP OPTIMIZATION FOR BSR
The problem of optimization of QCD radiative corrections in the parameter space {∆} was put in [30], see Sec.VI there. We resolve the problem here by means of direct numerical calculation of the minimum of radiative corrections to BSR, i.e. , the minimum of the function |f Rad |, obtained within the admissible domain {∆} at every t. The approach to PT optimization when the bare minimum of |f Rad (t; {∆}| is restricted by the set of inequality constraints in (37,38,39) is universal and does not depend on knowledge of details of the {β}-expansion for the quantity. Indeed, we will not use in the numerical analysis below any information about the intrinsic structure of the PT coefficients c i .

A. The numerical results of renormalization group optimization of C Bjp
The quantity |f Rad (t; {∆}| is equal to α s π (1 + O(α s )) for C Bjp and was interpreted in [32] as auxiliary "effective α g1 s /π", which we will use for convenience. Different optimization results below will be compared with the initial result in Eq.(32a), Let us start with optimization in the {∆ 0 , ∆ 1 }-space, the optimal points in Fig. 3 (blue points there) are located in the left part mostly on the boundary of admissible domains {∆} 1 0 . For t τ = ln(m 2 τ /Λ 2 qcd ) and t τ − ∆ = t ′ τ we have As it follows from the comparison of Eq.(43) and Eq.(44c), the effective α g1 s /π changes from 0.171 to 0.149, respectively. Therefore, the advantage of the optimization looks substantial in this case, while the value of ∆ in (44a) is moderate. To hold a correspondence with the original BLM result (∆ BLM 0 = 2 > 0) and qualitatively close hereto the PMC results in [32], we shall also check a minimum of |f Rad | at the additional condition ∆ 0 > 0. The corresponding points are shown in Fig. 3 (red in color) within the admissible domains in the right half of the plane, do not provide significant advantages in comparison with the corresponding 2D results in Eqs.(44c, 45c). On the other hand, the PT convergence is spoiled for this case, see the underlined terms there, and the final results in Eqs.(46c, 47c) do not seem reliable. One can only admit that the effective α g1 s /π 0.144. Therefore, we admit and will use further the single 2D result in Eq.(44c) which is sufficiently good.
The disadvantage of this brute force approach is that the "blind analysis" based only on inequality constraints (37,38,39) can lead to the qualitatively unsatisfactory result. The minimum of the radiative corrections in the cases (45b, 46b, 47b) appear on the boundary of constraint (38), see the underlined terms there, where the PT convergence is deteriorating. Another lesson from these results is that the "BLM/PMC preferable" scales at ∆ 0 > 0 appear not close to the bare minimum of radiative corrections.
These conclusions following from the results of numerical analysis require either to return to the accurate consideration of the {β}-expansion of the c i coefficients, or to propose more restricted criteria for the optimum. In this connection smaller values for α g1 s /π obtained in [32] within PMC correspond rather to the condition ∆ 0 > 0, although the direct comparison is difficult due to the different PT schemes and presentations. The behavior and value of this α g1 s /π and the corresponding normalization scales in [32] qualitatively agree with the condition, while the bare minimum is not the explicit goal of the used PMC procedure. Let us mention here that we do not agree with the construction of {β}-expansion used in this PMC, see our criticism in [30,34].

B. Optimized Bjorken Sum Rule for COMPASS measurements
Here we compare the results for the BSR obtained in both presented approaches, see the beginning of Sec. IV. In the result of COMPASS data analysis in Sec. III based on the optimized generalized BSR incorporating experimental uncertainties on the spin function g 1 we obtained Γ exp-opt 1 = 0.191 for Q 2 0 = 3 GeV 2 , while the original COMPASS estimations [19] provide Γ c-ss 1 = 0.192 ± 0.007 stat ± 0.015 syst . Bellow we shall miss the high twist correction µ p−n 4 /Q 2 in the theoretical RHS part of Eq. (30), in the discussion of COMPASS result. The reason is that in processing COMPASS data for Γ c-ss 1 , this effect has not been included in an explicit form due to high Q 2 for the most of these data. In the previous section, we discussed the optimized results for the QCD radiative corrections at the world reference scale m 2 τ . Below, we provide similar results starting with C Bjp 1, a s (Q 2 ) at the COMPASS reference scale Then, using the already checked 2D {∆ 0 , ∆ 1 } optimization in Eq.(44), we find the optimized value of C Bjp , that is visibly larger than the non-optimized result in Eq.(48). This value leads to the estimate for Γ th 1 in (30), in comparison with the result Γ th 1 (3GeV 2 ) = 0.177 from Eq.(48) that is taken at a s (Q 2 0 ) ≈ 0.0268. Comparing our results for experimental determination of Γ exp-opt 1 = 0.191 and theoretical prediction Γ th-opt 1 = 0.182, we can conclude that the difference is a bit reduced due to both kinds of optimization. The result for the bare minimum, Γ th-bare 1 ≈ 0.183, does not improve situation distinctly.

VII. CONCLUSIONS
Here we have discussed and improved the determination of the polarized Bjorken sum rule Γ 1 in two different aspects: (i) phenomenological determination of Γ 1 = 1 0 g 1 (x)dx by means of integration of the spin structure function g 1 within the truncated moment approach [16] improved here; (ii) the optimization of the truncated perturbation series for the coefficient function C Bjp (α s ) of the corresponding DIS process by means of the renormalization group. Using the truncated Mellin moments approach, we constructed the generalized Bjorken sum rule that allows one to determine the BSR value from the experimental data in a restricted kinematic range of the parton fraction x. We also showed how to incorporate different uncertainties for each kinematic bin from experiment. We applied our analysis to recent COMPASS data on g 1 and demonstrated better convergence in low limit of integration and more stable estimation of Γ exp-opt 1 ≈ 0.191. We provided certain criteria for the optimized analysis of C Bjp (α s ) within the perturbation QCD and the renormalization group. Based on these criteria we found the admissible domains for the corresponding new normalization scales µ 2 for α s . Within these domains we found numerically the minimum of the QCD radiative corrections to C Bjp (α s ), which leads to optimum values of theoretical predictions for Γ th-opt 1 ≈ 0.182 in the order O(α 4 s ). As a result, both values Γ exp 1 and Γ th 1 found in this way become closer to each other when processing the data of COMPASS experiment. with the SU c (N )-group fundamental fermion invariants The exact solution of Eq. (B.4) can be expressed in terms of the Lambert function W (z), [35] (see also [36,37]) defined by z = W (z) e W (z) .

(B.5)
This solution has the form , where z(t) = (1/b 1 ) exp (−1 + iπ − t/b 1 ) and the branches of the multivalued function W are denoted by W k , k = 0, ±1, . . .. A review of the properties of this special function can be found in [35]. The known second-iteration solution of Eq. (B.4) that provides us with sufficient accuracy with the following result 3. The approximate solution of the renormalization-group equation in the four-loop of QCD [38], where the βfunction is given by Eq. (B.3), assumes the asymptotic expansion where l = ln(t).