Impact of recent COMPASS data on polarized parton distributions and structure functions

We perform a new extraction of polarized parton distribution functions (PPDFs) from the spin structure function experimental data in the fixed-flavor number scheme (FFNS). In this analysis, we include recent proton and deuteron spin structure functions obtained by the \texttt{COMPASS} collaboration. We examine the impact of the \texttt{new COMPASS} proton and deuteron data on the polarized parton densities and compare with results from our previous study (KATAO PPDFs), which used the Jacobi polynomial approach. We find the extracted PPDFs of the proton, neutron, and deuteron structure functions are in very good agreement with the experimental data. The results for extracted PPDFs are also compared with available theoretical models from the literature.


I. INTRODUCTION
One of the principal goals of Quantum Chromodynamics (QCD) has been the detailed investigation of the spin structure of the nucleon and nuclei, as well as the determination of the partonic composition of their spin projections. The extraction of polarized or spin-dependent parton distribution functions has been recognized as a longstanding issue of physical interest [1,2], and theoretical studies on the spin structure of the nucleon have been discussed extensively in several reviews [3][4][5][6][7].
Determinations of polarized parton distribution functions (PPDFs) with an estimate of their uncertainties have been presented in multiple studies . The variation among these PPDFs sets can be due to a number of factors including the choice of experimental data sets, the form of the parameterization and uncertainty calculation, as well as the details of the QCD analysis such as the treatment of heavy quarks or higher-twist corrections.
The results from various calculations can lead to a wide range of expectations for the polarized observables; hence, it is illuminating to compare the results of different methodologies to the experimental measurements. In our previous analysis, we performed the detailed pQCD analysis of PPDFs using the orthogonal Bernstein and Jacobi polynomial methods at next-to-leading-order (NLO) [44][45][46]. Other theoretical studies implementing a QCD analysis on the spin structure of the nucleon using orthogonal polynomials have been reported in Refs. [47][48][49][50]. Thus, one goal of our investigation is to revisit this topic using a Mellin moment approach instead of the orthogonal polynomial approach.
Recently, the COMPASS collaboration [68,69] extracted new DIS measurements of the polarized proton and deuteron structure functions for the region of 0.0035 < x < 0.575, 1.03 < Q 2 < 96.1 GeV 2 and 0.0045 < x < 0.569, 1.03 < Q 2 < 74.1 GeV 2 , respectively. Thus, we will combine the data sets used in Ref. [44] with the COMPASS16 and COMPASS17 data sets to extract improved polarized structure functions and PPDFs.
The plan of this paper is as follows. In Sec. II, we review the theoretical framework and basic formalism of the polarized structure function analysis based on the inverse Mellin technique. In Sec. III we outline the parameterization of PPDFs and the selection of the data sets. In Sec. IV, we present the structure functions, PPDFs, and moments obtained in our fit, and compare these both to our earlier KATAO fit (using orthogonal polynomials) as well as other results from the literature; this also includes an evaluation of the impact of the new COMPASS data sets. Finally, in Sec. V, we provide a summary and concluding remarks.

II. THEORETICAL FORMALISM
The QCD formalism allows us to express the spin dependent nucleon structure function g 1 (x, Q 2 ) in terms of a Mellin convolution of the polarized non-singlet δq N S i , the polarized singlet δΣ, and the polarized gluon δg distributions with the corresponding Wilson coefficient functions δC N S q , δC S q , and δC g . The polarized structure function is then given by [4] where e j denotes the charge of the jth quark flavor, n f is the number of light flavors, x is the Bjorken variable, Q 2 = −q 2 is the four-momentum transfer, and the symbol ⊗ denotes the Mellin convolution. The coefficient functions δC i which we use in the present analysis are calculated in the MS-scheme at next-to-leading order [70][71][72][73]; in particular, we make use of the Pegasus routines [74]. The spin dependent flavor non-singlet distribution δq N S j evolves independently, while the spin dependent singlet δΣ and gluon distributions δg are coupled in the QCD evolution.
In the above equation, the polarized non-singlet and singlet PPDFs are expressed by the individual spin dependent quark flavor contributions as where δq j is the polarized quark distribution function of the jth light flavor. In our fits, we will take the strong coupling constant α s (Q 2 0 ) at initial scale Q 2 0 as a free parameter to be fit. The evolution of the strong coupling constant α s (Q 2 ) can be obtained from the QCD renormalization group equation and is determined by the β-function, β(Q 2 ): Here we have expanded the β-function in powers of α s out to NLO, and the first two coefficients can be computed in the MS-scheme to be β 0 = 11− 2 3 n f and β 1 = 102− 38 3 n f . Thus, given the value of α s (Q 2 0 ) at the initial scale Q 2 0 , we can numerically solve the differential equation in Eq. (4) for any Q 2 scale [74]. For the present analysis, we will work in the FFNS with n f = 3 light partonic flavors {u, d, s}.
For our fit, we will use the spin dependent proton, neutron, and deuteron structure functions. The spin dependent deuteron structure function xg d 1 (x, Q 2 ) can be represented in terms of the proton and neutron structure functions, xg p 1 (x, Q 2 ) and xg n 1 (x, Q 2 ) using the relation where ω D = 0.05 ± 0.01 is the D-state wave probability for the deuteron [75]. For comparison with the data, we will need to compute the PPDFs and structure functions at a variety of Q 2 scales. The evolution in Q 2 is performed using the wellknown DGLAP collection of integro-differential evolution equations [76,77] which can be solved analytically after a conversion from x-space to Mellin N -moment space.
The N 'th Mellin moments of the spin dependent parton densities δf (x) are defined to be The Mellin transform will decompose the convolution of parton densities δf (x) of Eq. (1) into a product of Mellin moments: To invert the Mellin transform, the argument N is analytically continued to the complex plane. Note that the basic method of solving the spin dependent non-singlet, singlet, and gluon evolution equations in Mellin space is reported in the literature in detail [71,72,78,79]. The solution of the flavor non-singlet, singlet and gluon evolution equations at NLO are given by where a s ≡ a s (Q 2 ), a 0 = a s (Q 2 0 )/4π, δP NS and δP (1) NS denote the LO and NLO non-singlet splitting functions.
Here, the matrices U 1 and L are evolution matrices, for some details see Ref. [79].
The basic framework of this method is described in the literature [30,37,80]. The resulting δf (x) for all PPDFs depends on the initial value of α s (Q 2 0 ) and unknown parameters of the spin dependent parton distributions; we will now discuss our parameterization form.

III. INPUT PARAMETERIZATION AND DATA SETS
To study the impact of the recent COMPASS16 [68] and COMPASS17 [69] data on the spin dependent parton distribution functions, we will start by comparing to our previous KATAO [44] results; hence, our initial parameterization and χ 2 minimization will be based on this work.

A. Parameterization of the polarized parton densities
For the parameterization of the spin dependent parton densities in x space at our initial scale Q 2 0 = 4 GeV 2 , xδq j (x, The free parameters are {η j , α j , β j , γ j }, and we use the common notation δq j = {δu v , δd v , δq, δg} for the partonic flavors up-valence, down-valence, sea, and gluon. In this yields a total of 16 degrees of freedom, but we will introduce some constraints to reduce the number of free parameters in order to achieve a stable and reliable minimum. The A j and η j parameters are not independent. Since the first moment of polarized parton densities plays an important role, the normalization constants A j are selected such that η j are the first moments of spin dependent of parton densities δq j (x, Q 2 0 ); specifically η j = 1 0 dxδq j (x, Q 2 0 ). Thus, the normalization factors A j can be computed to be: where B(m, n) is the Euler beta function. We will presume a SU (3) flavor symmetry such that δq ≡ δu = δd = δs = δs. As we mentioned before, by including only inclusive DIS data in the QCD fit, it is not possible to separate polarized quarks from polarized anti-quarks. In fact, inclusive polarized DIS data con-strain the total polarized quarks and anti-quarks combinations. 1 Thus, we will focus on the PPDF combinations (δq + δq) as displayed in Fig. 5.
Using the above results, we can analytically compute the Mellin-N space transform of the polarized parton densities at the input scale of Q 2 0 : The first moments of the polarized valence distribution, δu v and δd v , can be fixed by utilizing the parameters F and D as measured in neutron and hyperon β-decays [81,82]. In fact q 3 and q 8 are the non-singlet combinations of the polarized parton densities: δq 8 = (δu + δu) + (δd + δd) − 2(δs + δs) . (13) captionComparison of the parameter values and their statistical errors at the input scale Q 2 0 = 4 GeV 2 in the different cases: KATAO (Jacobi polynomial method) [44], Base (without COMPASS16 and COMPASS17), Fit A (with COMPASS16), Fit B (with COMPASS16 and COMPASS17) obtained from the best fit to the data.
The first moments of the above distributions are found to be:ˆ1 Using F =0.464±0.008 and D= 0.806±0.008 from the literature [31,83], we find the first moments of δu v and δd v to be η uv = +0.928 ± 0.014 and η dv = −0.342 ± 0.018; in our QCD fit we will fix {η uv , η dv } to these central values.
The first moments of δq and δg do not have prior constraints, and these will be determined in the fit by the free parameters ηq and η g . The above value for the octet axial charge assumes a good SU (3) symmetry. It was noted in Refs. [84,85] that this symmetry can be broken by about 20% which would then yield F ∼ 0.43 and D ∼ 0.84, and thus η uv ∼ +0.865 and η dv ∼ −0.405. We have also run our fit with these modified values and observed that the variation due to these changes is small and well within our PPDF uncertainties.
The factor of (1 + γ j x) in Eq. (9) provides flexibility of the parameterization in the intermediate x region. This flexibility is beneficial for fitting the the polarized valence distributions δu v , δd v . In contrast, we find that the parameters γq and γ g have a very mild impact on the fit and it is sufficient to set them to zero and remove these degrees of freedom. (We note that the QCD analysis of polarized SIDIS data [80] is sensitive to the γq and γ g parameters. ) We have now reduced the number of free parameters from 16 to 12. Preliminary fits indicate that some of the parameters such as {γ uv , γ dv , βq, γ g } are very weakly constrained by the present data set and have very large uncertainties. In fact, the precision of the data which we used is not high enough to constrain these mentioned parameters sufficiently. We found that, altering them within these uncertainties does not obtain a significant change of χ 2 . Therefore we will also fix the values of these parameters, and we now have a remaining 8 free parameters for the PPDFs in addition to the QCD coupling constant α s (Q 2 0 ) to fit from the data.
In this analysis we will evolve the PPDFs from the initial scale Q 2 0 = 4 GeV 2 up to arbitrary scales to compare our theoretical predictions with the data across the full kinematic range. We construct a global χ 2 function using the experimental measurements g Exp 1 , the experimental uncertainty (statistical and systematic added in quadrature) ∆g Exp 1 , and theoretical prediction g T heory 1 . Our χ 2 is constructed as follows: where the i-index sums over all experimental data sets, and in each experimental data set the j-index sums over all data points. We introduce a weight w i which allows us to apply separate weights to different experimental data sets; for the present analysis we choose all weights to be unity, w i = 1. These data sets include statistical and systematic errors which we combine in quadrature. There is also a normalization for each experiment K i and an associated uncertainty ∆K i . The normalization shifts K i are fitted at the start of our procedure, and then fixed. We present these values in Table I, and find that all the K i shifts are less than 1% except for a single value; for the SLAC/E155 experiment we find K i = 1.024.
As outlined in Sec. III A, we have a total of 9 unknown free parameters: 8 parameters describing the PPDFs at Q 2 0 , and also α s (Q 2 0 ) as another free parameter. We will use the CERN library MINUIT package [86] to minimize χ 2 by varying the free parameters and obtain a best fit.    [62], E155 [58,59], E143 [57], and HERMES06 [51].
We are now ready to extract the polarized parton densities.

IV. RESULTS OF THE QCD ANALYSIS
In this section, we will demonstrate how inclusion of the new COMPASS proton g p 1 data [68] and deuteron g d We will divide our analysis into three steps. As a first step, we perform a fit with all the data of Table I with the exception of the COMPASS16 [68] and COMPASS17 [69] experimental data; this totals 379 data points, and we identify this as our "Base" fit. We then include the COMPASS16 proton data, and this is our "Fit A" which contains 430 data. Finally, we include the COMPASS17 deuteron data, and this is our "Fit B" with the full 473 data points. As Fit B contains the complete data set, we will use this for comparisons in Figs. 2, 3, 4, and 5 where it is identified as "This Fit." In Table III, the final values of the fit parameters for the different data sets are summarized. We find that χ 2 /d.o.f is less than unity in all cases indicating a good quality of fit. Additionally, our fits compare well with our previous KATAO analysis where we find We will begin with the comparison of the xg N 1 structure functions as this is the primary input to our fit. In Figs. 2, 3, and 4, we display the comparison of our theoretical predictions with the structure function data for xg p 1 , xg d 1 and xg n 1 , respectively. The figures are given as a function of Q 2 at different values of x and are compared to all of the experimental data that we used in the present analysis. The theoretical predictions are in good agreement with the experimental measurements across the fill x-range. In the following sections, we will investi-
We derive the uncertainties of the polarized parton distributions for the different polarized observables using the covariance matrix elements of the QCD fit.
Examining Fig. 5 we find that the spread of results for the x(δu + δū) distribution is comparatively narrow indicating this flavor component is well constrained. The results of "Fit B" are comparable to our previous analysis using the Jacobi polynomial expansion method (KATAO) [44], as well as many of the other results from the literature. Our results are slightly larger than those of BB in the larger x region (x ∼ 0.2). The x(δd + δd) distribution is also comparatively narrow suggesting this too is well constrained. Again, our results of "Fit B" are generally comparable to the other results from the liter-ature, with "Fit B" yielding a slightly larger x(δd + δd) than BB in the region x ∼ 0.1. For the x(δs + δs) distributions (or 2δq in our notation), we find a broader spread of both our results ("Fit B" and KATAO) and the other fits from the literature suggesting this component is less constrained. Specifically, "Fit B" roughly coincides with many of the other predictions, but the DSSV and LSS10 results yield a changes sign as a function of x and LSS14 yields a larger result. Of all the components we examine, clearly the gluon distribution xδg has the widest spread of predictions and the greatest uncertainty. "Fit B" is similar to the KATAO results, but yields a smaller result in the region x ∼ 0.3; compared to the other curves, these results generally give a smaller xδg than the other predictions. In particular, in the region x ∼ 0.1 AAC give the largest result and DSSV gives a negative results. Clearly, the xδg distribution leaves much room for improvement and it will be interesting to see which predictions are favored by future data sets. Presumably, the choice of data sets (such as SIDIS) may contribute to these differences. are showen with error bands compared to results obtained by KATAO [44], BB [37], AAC [31], DSSV [35], GRSV [27] and LSS [13,14]. For clarity, we only present our Fit B results, labeled here as "This Fit."

Comparison of {Base,Fit A, Fit B} on PPDFs
Since it is the new COMPASS data on xg p 1 and xg d 1 that represent the important new additions to our data set, we want to focus on the variations among our fits: {Base, Fit A, Fit B}. To investigate the specific impact of COMPASS16 and COMPASS17 data sets, we compare our results for our individual fits: "Base" (without including COMPASS16 and COMPASS17), "Fit A" (including COMPASS16) and "Fit B" (including COMPASS16 and COMPASS17). These results are shown in Fig. 6 where we have displayed both the absolute value of the PPDFs and also the ratio compared to our base fit.
As suggested by the results of Fig. 5, in Fig. 6 we find that xδu v (x) and xδd v (x) appear to be strongly constrained with little variation among the separate fits. Specifically, the variation is on the order of a percent except for the region at large x where the PPDFs vanish and there are no data constraints.
In contrast, xδq(x) and xδg(x) do display some differences amount the fits due to the addition of the COM-PASS data; the variations of "Fit A" and "Fit B" of Fig. 6 are quite similar, and these differ from the "Base" fit. The xδq(x) function displays some variation in the small x region 10 −1 while the variation of xδg(x) function is generally at larger x 10 −1 ; again, the very large x region should be discounted as before. To examine how the fits change with the inclusion of the COMPASS data, we examine the partial χ 2 contributions to COMPASS16 and COMPASS17 data set for each of our fits: {Base, Fit A, Fit B}. If we compute χ 2 for the COMPASS16 data set using the "Base" fit (which does not include this data), we find a total χ 2 value of 34.67 for the 51 COMPASS16 data points, and when we include this data in the fit ("Fit A") it improves slightly to 33.48. Correspondingly, if we fit the COMPASS17 data set using the "Fit A" (which does not include this data), we find a total χ 2 value of 27.43 for the 43 COMPASS17 data points, and in the fit ("Fit B") this is quite similar at 27.22. Thus, both the COMPASS16 and COMPASS17 data set are in reasonable agreement to the initial "Base" fit. The changes among the {Base, Fit A, Fit B} sets is most evident in the ratio plots of Fig. 6.
Finally, in Figs. 7 and 8, we directly compare our "Fit B" with the proton and deuteron polarized structure functions from COMPASS16 [68] and COMPASS17 [69] experimental data in a composite plot; as the individual data range over Q 2 , we display our predictions with selected values of Q 2 to illustrate the evolution effects. This allows us to see the comparison of data and theory in a compact, albeit approximate, manner.

C. Moments and Sum Rules
We now turn to integrated moments and sum rules. Note that the calculation of the moments integrates over the full range x = [0, 1], so this requires some extrapolation outside the x range where the structure functions have been measured.

PPDF Moments
We start by computing the PPDF moments, as these will be the necessary ingredients for the other moments and sum rules that follow.
In Table II, we compare the results of the first moments of the polarized parton densities for our fits with results from the literature at NLO in the MS-scheme at Q 2 = 4 GeV 2 . Comparing our "Base" fit with "Fit A" and "Fit B" we see the moments are generally stable with the exception of the ∆g which varies by ∼ 30%. Including the other PPDF moments from the literature, we see the results for {∆u v , ∆d v } are quite stable (∼ 1%) while {∆u, ∆d} show a bit more variation (∼ 10%), and finally {∆q, ∆g} a larger spread (> 100%). We will now look at the influence of the above PPDF moments on the experimentally measurable structure functions.

Structure Function Moments
We next examine the first moment of the xg N 1 (N = p, d, n) structure functions defined to be: In Table III, we compare the results for Γ N 1 (Q 2 ) of Fit B with the COMPASS measurements. We observe that the fit agrees with the COMPASS results within ∼ 1σ of the experimental uncertainty.
Next, in Table IV, we compare our first moment results with those from the literature. The theoretical results for Γ p 1 are uniform within ±2%, while the range on Γ d 1 increases to ±5%, and the range on Γ n 1 further increases to ±15%.
3. Bjorken Sum Rule, xg N S 1 (x, Q 2 ) and Γ N S 1 (Q 2 ) Following Ref. [51], in the scaling (Bjorken) limit we have We can isolate the a 3 term by taking the difference between the proton and neutron terms, and we will identify this as the non-singlet (NS) contribution. Thus, Γ N S 1 (Q 2 ) enters the polarized Bjorken sum rule [100] and is related to the ratio of the axial and vector coupling constants (g A,V ). Here, C N S 1 (Q 2 ) is the non-singlet coefficient function.
In a similar manner we define g N S 1 (x, Q 2 ) as the difference between the proton and neutron structure functions: In Fig. 10 we compare our results for xg N S 1 (x, Q 2 ) with the HERMES experimental data [51] for selected bins of Q 2 . We find minimal variation among our different theoretical fits (including the previous KATAO fit), and these curves compare well with the experimental results.
From Eq. (20) we can also relate Γ N S 1 (Q 2 ) to the previously computed proton and neutron first moments as: These results are presented in Table IV and with the COMPASS results. The result of our "Fit B" is comparable to COMPASS16, and below (but within uncertainties) to COMPASS17.
We can also calculate the structure function g N 2 (x, Q 2 ) via the Wandzura-Wilczek relation [101,102]: Figure 11 shows the polarized structure function xg p 2 and xg d 2 as a function of x for different cases of Base, Fit A, Fit B and our previous KATAO results [44] in  Figure 10. NLO non-singlet polarized structure function xg N S 1 (x, Q 2 ) as function of x in comparison with the results of KATAO [44] and HERMES experimental data [51].  [103], HERMES [104], and SMC [105] experimental data at Q 2 = 5 , 6 GeV 2 . As the data actually span over a range of Q 2 , in Fig. 12 we display the Q 2 evolution of the polarized structure function xg 2 (x, Q 2 ) for the proton and deuteron as function  Fig. 11 we see that our "Base" and "Fit A" coincide throughout the x range suggesting a minimal impact from the COMPASS16 data on this observable; conversely, our "Fit B" does differ, especially in the larger x region, suggesting a stronger influence of the COMPASS17 data on xg d 2 (x, Q 2 ).  Figure 11. NLO polarized structure function xg p 2 (x, Q 2 ) and xg d 2 (x, Q 2 ) as a function of x for Q 2 = 5, 6 GeV 2 compared to E143 [57], E155 [103], HERMES [104], and SMC [105] experimental data. We present also the results of different Base (dashed-dotted-dotted), Fit A (dashed) and Fit B (solid line) QCD fits which are compared with our previous KATAO (dashed-dotted) results.

D. The Proton Spin
It is important for us to understand the decomposition of the proton spin in terms of the separate contributions from the quarks, gluon, and the orbital angular momentum components. The spin of the proton can be computed from the first moment of the polarized parton densities together with the quark and gluon orbital momentum (L q , L g ) is as following [106]   Here L z (Q 2 ) = L q (Q 2 ) + L g (Q 2 ) is the total orbital angular momentum of all the quarks and gluons, ∆g(Q 2 ) = 1 0 dx δg(x, Q 2 ) is the first moment of the polarized gluon distribution, and ∆Σ(Q 2 ) =´1 0 dx δΣ(x, Q 2 ) with δΣ ≡ δu v + δd v + 6δq is the first moment of the polarized singlet distribution. In Eq. (22), we note that the spin sum ( 1 /2) is actually independent of Q 2 even though each individual term is dependent on Q 2 . In Table V we compute { 1 /2∆Σ, ∆g} using "Fit B" at Q 2 = 4 GeV 2 , and then infer the value of L z (Q 2 ) assuming Eq. (22). As we observed in Table II the Figure 13. NLO polarized singlet parton density xδΣ(x) at Q 2 0 = 4 GeV 2 for this fit ("Fig B") as a function of x, compared with results from the literature including KATAO [44], BB [37], AAC [31], DSSV [35], GRSV [27] and LSS [13,14]. literature were displayed in Fig. 5, and there is quite a bit of variation. In contrast, Figure 13 shows our NLO singlet polarized parton density xδΣ(x) (≡ δu v + δd v + 6δq) compared with other results from the literature. The results of this fit ("Fit B") with the previous analysis KATAO [44] are quite similar as suggested by Table II. Generally, the singlet polarized distributions are negative for x 0.04 for most of the models, but there is some slight variation in the range x 0.06 to x 0.02. Overall, the variation of xδΣ(x), as compared to xδg(Q 2 ), is reduced; this is notable as xδΣ(x) is a combination of both valence and sea PPDFs.
In contrast to our previous polarized analysis (KATAO), we did not use the Jacobi polynomial expansion method. Our results for the PPDFs are comparable to other extractions, and generally it appears that xδu v and xδd v are comparatively well determined in contrast to xδq and xδg which display a larger variation across the x range.
We also computed various structure functions and moments for the proton, neutron, and deuteron, and these also compare well with the both the COMPASS data, as well as other determinations from the literature. Again the results from this fit are comparable to the previous KATAO [44] results using orthogonal polynomials; it is reassuring to see that the results are generally independent of the underlying calculational methodology.
The strong coupling constant α s (Q 2 0 ) was extracted from the fits, and the uncertainty is slightly decreased compared to the KATAO analysis. This α s (Q 2 0 ) can be evolved up to α s (M 2 Z ) by assuming an evolution order (LO, NLO, ...) and heavy quark mass thresholds; we find values that are low compared to the world average, but within uncertainties.
From this analysis, it appears that both the various theoretical analyses using a variety of techniques and (x-space, N -space, orthogonal polynomials) are generally converging to yield a homogeneous set of predictions which are in good agreement with the diverse sets of experimental measurements. While there is still room for further improvements, such studies provide a strong validation of the underlying QCD theoretical framework.
A standard LHAPDF library file of our polarized PDFs {xδu v (x, Q 2 ), xδd v (x, Q 2 ), xδq(x, Q 2 ), xδg(x, Q 2 )} and their uncertainties can be obtained via e-mail from the authors upon request.