Global analysis of polarized DIS & SIDIS data with improved small-x helicity evolution

,

The proton spin puzzle has been one of the most intriguing and profound mysteries in our understanding of the proton structure for over three decades (for reviews see Refs.[1][2][3][4][5][6][7][8][9]).The main challenge is to determine, both qualitatively and quantitatively, how the proton spin is distributed among the spins and orbital angular momenta (OAM) of its quark and gluon constituents.The question is usually formulated in terms of spin sum rules, such as the Jaffe-Manohar sum rule [10] (see also the Ji sum rule [11]), that decompose the proton spin of 1/2 (in units of ℏ) into the sum of the quark (S q ) and gluon (S G ) spins and the OAM carried by the quarks (L q ) and gluons (L G ), Each of the contributions in Eq. ( 1) can, in turn, be written as the integral of a partonic function over the longitudinal momentum fraction x carried by the parton.For example, with similar expressions for the OAM contributions [12][13][14][15][16], where ∆Σ(x, Q 2 ) is the flavor singlet combination of the quark helicity parton distribution functions (hPDFs) ∆q(x, Q 2 ) (quark flavor q) and ∆G(x, Q 2 ) is the gluon hPDF [10].The goal of current research in the field of proton spin physics is to determine ∆Σ(x, Q 2 ), ∆G(x, Q 2 ), L q (x, Q 2 ), and L G (x, Q 2 ) across a broad range of x and Q 2 in order to quantify how much of the proton spin is carried by the partons in different kinematic regions.
The standard way to address the proton spin puzzle is by extracting the hPDFs ∆q(x, Q 2 ) and ∆G(x, Q 2 ) from experimental data using collinear factorization along with the spin-dependent Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [17][18][19] to relate observables at different Q 2 .There has been a number of very successful extractions of hPDFs over the years within this approach [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34].Nevertheless, the DGLAP-based methodology has a drawback: since the DGLAP equations evolve PDFs in Q 2 , they cannot truly predict the x dependence of PDFs.The x dependence is greatly affected by the functional form of the PDF parametrization at the initial momentum scale Q 2 0 , which gives the initial conditions for the DGLAP evolution.The parameters are then determined by optimizing agreement between the theoretical calculations to the experimental measurements.In this way, the experimental data, in the x range where it is available, make up for the inability of DGLAP evolution to predict the x dependence of PDFs.Conversely, in the x region which has not yet been probed experimentally, DGLAP-based predictions typically acquire a broad uncertainty band due to extrapolation errors.This is particularly true in the small-x region.Since no experiment, present or future, can perform measurements down to x = 0, further theoretical input is needed to constrain the hPDFs at low x.The benefit of small-x helicity evolution is it makes a genuine prediction for the hPDFs at small x given some initial conditions at a higher x 0 .Due to the integrals in Eq. (2), precise control over the behavior of hPDFs at small x is mandatory to resolving the proton spin puzzle.

B. Proton spin at small x
The first resummation of hPDFs at small x was performed in the pioneering work by Bartels, Ermolaev and Ryskin (BER) [35,36], who employed the infrared evolution equations (IREE) formalism from Refs.[37][38][39][40][41].The BER IREE resummed double logarithms of x, i.e., powers of the parameter α s ln 2 (1/x) (with α s the strong coupling constant), which is referred to as the double-logarithmic approximation (DLA).The leading small-x asymptotics for the flavor singlet combination of quark hPDFs and the gluon hPDF can be written as with α h the helicity intercept.BER found α h = 3.66 αsNc 2π in the pure-gluon case and α h = 3.45 αsNc 2π for N f = 4 (the numbers 3.66 and 3.45 were calculated numerically, the latter for N c = 3, with N c /N f being the number of quark colors/flavors).These intercepts are numerically large, with α h > 1 for realistic coupling α s = 0.2 − 0.3, making the integrals (2) divergent as x → 0. One may hope that the higher-order corrections in α s , when calculated, would lower the intercept α h below 1, making the integrals (2) convergent.In addition, at very small x, parton saturation corrections (see Refs. [42][43][44][45][46][47][48][49] for reviews) are likely to significantly modify the asymptotics (3) by slowing down (or completely stopping) the growth of hPDFs with decreasing x (see, e.g., [50] for the impact of saturation effects on the unpolarized flavor nonsinglet evolution).Phenomenological applications of the BER IREE approach were developed in Refs.[51][52][53][54][55][56].Recently, the BER approach has been applied to the OAM distributions as well [57].
The equations developed in Refs.[58,60,63,64,66,71] were also derived in the DLA.Similar to the unpolarized evolution equations [75][76][77][78][79][80][81][82][83][84], the helicity evolution equations [58,60,63,64,71] only take on a closed form in the large-N c [99] and large-N c &N f [100] limits.In that case they become the evolution equations for the so-called "polarized dipole amplitudes," which are dipole scattering amplitudes with an insertion of one gluon or two quark operators at the sub-eikonal level into the light-cone Wilson lines [63,64,71,92].The earlier version of this evolution, constructed in Refs.[58,60,63] (which we will refer to as KPS) led to an intercept of α h = 4 in the large-N c limit [61,62], significantly smaller than the intercept of α h = 3.66 αsNc 2π found by BER in the same limit.The KPS evolution has recently been augmented [71] by inclusion of the operators which couple what can be interpreted as the OAM of the gluon probe (in A − = 0 light-cone gauge of the projectile) to the spin of the proton. 1he revised evolution equations, which we will refer to as the KPS-CTT equations [58,64,71], have been solved at large N c both numerically [71] and analytically [101].While the former reference found the numerical value of the intercept to be α h = 3.66 αsNc 2π , appearing to agree with BER, the analytic solution [101] found that the BER and KPS-CTT intercepts at large N c disagree in the third decimal point.Very recently, a numerical solution of the large-N c &N f version of the KPS-CTT evolution [102] established a disagreement with BER (in the same limit) at the 2-3% level, with the discrepancy increasing with N f .While the observed differences between the two sets of results appear to demand further theoretical investigation, they are sufficiently small to allow one to proceed with rigorous phenomenological applications of the KPS-CTT evolution equations [58,60,63,64,71].
The first phenomenological application of the polarized dipole amplitude formalism, more precisely its KPS version, was performed by a subset of the present authors in Ref. [103].In that work a successful "proof-of-principle" fit of the world polarized DIS data for x < 0.1 and Q 2 > m 2 c (with m c the charm quark mass) based solely on small-x helicity evolution was performed.Since the analysis of Ref. [103] was limited to DIS data, only the g 1 structure functions of the proton and neutron were extracted instead of the individual flavor hPDFs.The impact of DIS data from the EIC on our ability to predict the g 1 structure function at small-x was also estimated.In addition, in order to demonstrate that it is possible to extract the combinations ∆q + (x, Q 2 ) ≡ ∆q(x, Q 2 ) + ∆q(x, Q 2 ) for q = u, d, s using small-x helicity evolution, parity-violating DIS EIC pseudodata was utilized.We refer to ∆q + (x, Q 2 ) as the C-even hPDFs, whereas the flavor nonsinglet C-odd hPDFs are similarly defined as ∆q

C. Subject of this work
In the present paper we perform, for the first time, a phenomenological analysis based on the KPS-CTT version of small-x helicity evolution with several other significant new features beyond the work of Ref. [103].Instead of the large-N c limit of evolution employed in Ref. [103], we base our analysis on the large-N c &N f limit.In addition to the polarized DIS data, we also include in our analysis polarized SIDIS data.Since the SIDIS data is sensitive to the individual quark and anti-quark helicity PDFs, ∆q(x, Q 2 ) and ∆q(x, Q 2 ), it is not sufficient to just use the flavor singlet helicity evolution from Ref. [71], which only yields the ∆q + (x, Q 2 ) combination (in addition to the gluon hPDF ∆G(x, Q 2 )).One also needs the flavor nonsinglet quark hPDFs ∆q − (x, Q 2 ).Those are constructed using the large-N c small-x helicity evolution equation for the flavor nonsinglet case from Ref. [60].Finally, to make the calculation more realistic and avoid the integrals (2) diverging at x → 0, we include running coupling corrections into the kernel of the evolution equations (both flavor singlet and nonsinglet).We make the coupling run with the daughter dipole size, which ends up effectively reducing the intercept α h for ∆q + and ∆G below 1. (The intercept of the flavor nonsinglet hPDFs is smaller than 1 even at fixed coupling in the realistic α s = 0.2 − 0.3 range; still, for consistency, we apply running coupling corrections to the flavor nonsinglet helicity evolution as well.)The analysis of SIDIS data also requires input for fragmentation functions, which are not specific to the small-x evolution at hand; therefore, we employ the existing JAM fragmentation functions for pions, kaons, and unidentified hadrons from Ref. [34].
The paper is structured as follows.We begin in Sec.II by outlining the polarized dipole amplitude formalism developed in Refs.[58,60,63,64,71] and explicitly writing out the flavor-singlet KPS-CTT large N c &N f DLA smallx helicity evolution equations with running coupling corrections, along with the flavor nonsinglet helicity evolution equation derived in Ref. [60].We also present the details of our numerical methodology in solving these evolution equations.We describe the calculation of observables (double-longitudinal spin asymmetries) in DIS and SIDIS, particularly detailing the calculation of the polarized SIDIS cross section at small x.We explain our analysis of the world polarized DIS and SIDIS low-x data and describe the implementation of the KPS-CTT evolution within the JAM Bayesian Monte Carlo framework.The results of our analysis are presented in Sec.III, which include plots of data versus theory, the hPDFs, and the g 1 structure function as well as an estimate of how much of the proton spin is carried by the net spin of partons at small x.We also conduct an EIC impact study on the aforementioned quantities.Conclusions and an outlook are given in Sec.IV.

II. METHODOLOGY
A. Flavor singlet evolution at small x The small-x helicity formalism in the light-cone operator treatment (LCOT) framework along with the large-N c &N f small-x evolution equations for helicity were revised in Ref. [71].In the new formalism, the (DIS) g 1 structure function is given by where e q is the quark electric charge as a fraction of the magnitude of the electron's charge.The C-even quark hPDFs in the DLA take the form [64,71] ∆q The gluon hPDF in the DLA is [63] ∆G(x, Note that the quark and gluon hPDFs ∆q + and ∆G are expressed in terms of the impact-parameter-integrated polarized dipole amplitudes Q q and G 2 , whose operator definitions can be found in Refs.[58,64,71] and Ref. [63], respectively.The dipole amplitudes depend on the transverse size of the dipole x 10 = |x 1 − x 0 |, where the "polarized" (sub-eikonally interacting) line is located at x 1 and the unpolarized (standard) Wilson line is at x 0 in the transverse plane.The amplitudes also depend on the center-of-mass energy squared s of the projectile-proton scattering.The dimensionless longitudinal momentum fraction z can be thought of as the momentum fraction of the softest of the two lines in the dipole.(However, this definition is somewhat imprecise, and it is more accurate to think of zs as the effective energy of the dipole-proton scattering [58,60,70].)The momentum scale Λ denotes our infrared (IR) cutoff and is the scale characterizing the proton.No dipole can be larger than 1/Λ, that is, the transverse size x 10 < 1/Λ.
At small x, Eq. ( 4) was derived in Refs.[58,60,61].However, the contribution of G 2 to ∆q + in Eq. ( 5) was recognized only recently [71].Given that G 2 is closely related to the gluon hPDF ∆G, as follows from Eq. ( 6), Eqs. ( 4) and (5) show that in our LCOT approach the contribution of ∆G to g 1 comes in through ∆q + [71,102] (see more on this below).We have also expanded the definition of the amplitude Q q to include dependence on the quark flavor q = u, d, s, such that we have three different amplitudes Q u , Q d and Q s for the light flavors, which is necessary since the quark spinor field operators are flavor dependent.The operator definition for the three flavors is the same, but the flavor dependence can enter through the initial condition of the dipole amplitude evolution.
While Eq. ( 4) appears to correspond to the leading-order (LO) expression in the collinear factorization approach to polarized DIS (see, e.g., Eq. (4.5) in Ref. [104]), in the LCOT framework it contains more information than that.In collinear factorization at the next-to-leading order (NLO) and beyond, the expression for the g 1 structure function also involves the contribution of ∆G.More precisely, one can write [18,19,[105][106][107][108][109][110][111][112][113][114] with the coefficient functions ∆c q (z) and ∆c G (z) calculated order-by-order in perturbation theory.In the MS scheme, the small-x large-N c &N f coefficient functions are [105] (see also [114] for the three-loop contribution, which we do not show explicitly here) Note that after the z-integration in Eq. ( 7), the contribution from the order-α s terms in Eqs. ( 8) becomes of the order α s ln 2 (1/x), the contribution from the order-α 2 s terms in Eqs. ( 8) becomes of the order [α s ln 2 (1/x)] 2 , etc. Consequently, in the collinear factorization power counting, the contributions from ∆c q (z) and ∆c G (z) in Eq. ( 7) are NLO and beyond, allowing one to truncate the expansion at a given order in α s determined by the accuracy of the calculation.In our DLA small-x power counting, the leading small-x parts of ∆c q (z) and ∆c G (z) are already included to all orders in the powers of α s ln 2 (1/x).This is precisely what Eq. ( 4) accomplishes [102].While it appears to be just the LO part of Eq. ( 7), the fact that ∆q + in it is evolved with the DLA small-x helicity evolution [58,60,63,64,71], resuming powers of both α s ln 2 (1/x) and α s ln(1/x) ln Q 2 /Q 2 0 , implies that Eq. ( 4) contains both the DLA DGLAP evolution of ∆q + , which mixes it with ∆G (by resumming the powers of α s ln(1/x) ln Q 2 /Q 2 0 ), and the leading small-x parts of the coefficient functions ∆c q (z) and ∆c G (z), resummed to all orders in α s ln 2 (1/x), bringing in the ∆G and additional ∆q + contributions into g 1 , as expected from Eq. (7) (see [102] for a more detailed discussion).The fact that all these contributions are contained in Eq. ( 4), which looks much simpler than Eq. ( 7), appears to suggest that we are working in the "polarized DIS scheme" [102] for our hPDFs (cf.[115] for the standard DIS scheme), where ∆G does not contribute to g 1 directly, unlike the more widely used MS scheme from Eq. (7).Other small-x calculations, such as the NLO BFKL evolution [116,117] (in the small-x power counting), result in the spin-independent GG anomalous dimension in the DIS scheme [108].This appears to be similar to our calculation giving a polarized DIS scheme result, with the difference between the anomalous dimensions in different schemes being proportional to N f [102,108].
The polarized dipole amplitudes Q q and G 2 , which enter Eqs. ( 4), ( 5) and ( 6), are found by solving the small-x evolution equations.The DLA large-N c &N f revised evolution equations at fixed coupling are given by Eqs.(155) in Ref. [71] (see also Refs.[58,64]).Its existing numerical solution [102] (with fixed coupling) leads to a large intercept α h for the flavor singlet hPDFs and for ∆q + (see Eq. ( 3) with the intercept values in the text following that equation), making the integrals in Eq. (2) divergent as x → 0. As we discussed above, this divergence may be regulated by higher-order corrections and/or by the onset of saturation, which is likely to slow down the growth of hPDFs as x → 0. As the unpolarized small-x evolution [72][73][74][75][76][77][78][79][80][81][82][83][84] is single-logarithmic, resumming powers of α s ln(1/x), a consistent inclusion of saturation effects is beyond the double-logarithmic approximation employed here.While, strictly-speaking, phenomenology based on small-x evolution in the DLA should work with the high intercepts found in Ref. [102], it appears unphysical to perform an analysis of experimental data with a formalism that would yield an infinite amount of spin at small x.While we cannot include the single-logarithmic (resumming powers of α s ln(1/x)) corrections to the revised DLA evolution equations (155) from Ref. [71], since they have not been fully calculated yet (see Ref. [70] for the single-logarithmic corrections to the earlier KPS evolution), we can include running-coupling corrections into the DLA evolution.A similar approximation was employed in the BER framework [53,55] and for the spin-independent eikonal small-x evolution [118,119], resulting in successful phenomenology.
In the DLA equations (155) from Ref. [71], the scale of the coupling could be given by either the "parent" (x 10 ) or the "daughter" (x 21 or x 32 ) dipole.The running coupling corrections to the (un-revised) KPS evolution, calculated in Ref. [70] (along with other single-logarithmic corrections), indicate that at DLA the coupling runs with the daughter dipole size.For the neighbor dipole amplitudes Γ, Γ, and Γ 2 , introduced in Refs.[58,60,63,64,66,71] and also entering helicity evolution equations, the coupling runs with the dipole size x 32 , which determines the next emission's lifetime and is integrated over in the kernel [70].Therefore, we proceeded by running the coupling with the daughter dipole size (or, more precisely, with the dipole size that we integrate over in the kernel) in all the terms of the KPS-CTT evolution.(See Refs.[120][121][122][123][124] for calculations and analyses of the running coupling corrections in the unpolarized small-x evolution case.)The resulting running-coupling version of the large-N c &N f helicity evolution equations (155) from [71] reads The running coupling in Eqs. ( 9) is given by the standard one-loop expression, with Λ QCD the QCD confinement scale.We have also modified Eqs.(9) compared to Eqs. (155) in Ref. [71] in two additional ways: first, we are now treating the momentum scale Λ as the infrared cutoff (assuming that Λ > Λ QCD ); second, since the amplitude Q q is now flavor dependent, we replaced the N f factors from Ref. [71] by flavor sums q .Eqs. ( 9) also include the dipole amplitude G, which is defined in Ref. [71]: as one can see from Eqs. (4), ( 5) and ( 6), the g 1 structure function and hPDFs do not depend on this dipole amplitude: this will affect our analysis below.Following Refs.[58,60,63,64,66,71] we have introduced the impact-parameter integrated "neighbor dipole amplitudes" Γ q (x 2 10 , x 2 32 , zs), Γ(x 2 10 , x 2 32 , zs) and Γ 2 (x 2 10 , x 2 32 , zs) for the amplitudes Q q , G and G 2 , respectively, with physical dipole transverse size x 10 and lifetime ∼ x 2 32 z.This lifetime for the neighbor dipole amplitudes depends on the transverse size of another (adjacent) dipole, giving rise to the "neighbor" amplitude name.
The inhomogeneous terms (initial conditions) in Eqs. ( 9) can be calculated at the Born level for a longitudinally polarized massless quark target instead of the proton.This gives [58,60,63,71] where ) is the Casimir operator in the fundamental representation of SU(N c ).These expressions will motivate our choice of the initial conditions for our phenomenological analysis.(While strictly-speaking we should have included running coupling corrections into the expressions (11) as well, the fixed-coupling form has a sufficient variety of dependence on the relevant variables zs and x 10 to motivate a fairly broad class of initial conditions we will implement below.) Flavor nonsinglet evolution at small x As one can see from Eq. (4) in the previous subsection, measurements of the g 1 structure function in DIS off a nucleon are only sensitive to a specific linear combination of ∆q + (x, Q 2 ).Such DIS measurements were the topic of our previous study [103].However, the polarized SIDIS process, as we will see below, provides information on the individual flavor hPDFs ∆q(x, Q 2 ), or, equivalently, on both ∆q + (x, Q 2 ) and ∆q − (x, Q 2 ) ≡ ∆q(x, Q 2 ) − ∆q(x, Q 2 ).The above evolution equations (9) only allow us to calculate ∆q + (x, Q 2 ).To perform the polarized SIDIS data analysis we need to supplement them with the small-x helicity evolution in the flavor nonsinglet channel.
A closed evolution equation at small x yielding ∆q − (x, Q 2 ) in the LCOT framework can be obtained in the large-N c limit, which is equivalent to the large-N c &N f limit for the flavor nonsinglet helicity evolution in DLA.(In the DLA, the flavor nonsinglet evolution is N f -independent, since virtual quark bubbles do not contribute.Thus, the large-N c and large-N c &N f limits are identical for flavor nonsinglet evolution.)Employing Eq. (54b) of [60] we write in the DLA We see that ∆q − (x, Q 2 ) only depends on one (impact-parameter integrated) polarized dipole amplitude, G NS q (x 2 10 , zs), for each flavor q = u, d, s.The definition of this dipole amplitude can be found in Eqs.(55) of Ref. [60].Just as in the flavor singlet case, the nonsinglet dipole amplitude can be determined by solving the small-x evolution equation, which reads [60] To be consistent with the flavor-singlet evolution, we have also inserted a running coupling into Eq.( 13), modifying it slightly compared to the fixed-coupling flavor nonsinglet evolution equation derived in Ref. [60].The inhomogeneous term in Eq. ( 13) can also be calculated at Born level for a quark target [60]: This expression will again motivate our choice of the flavor nonsinglet initial conditions in phenomenology.
C. Numerical implementation of the flavor singlet and nonsinglet evolution Similar to our previous works [61,67,71,102], small-x helicity evolution equations simplify if one performs the following change of variables, Here z (n) = z, z ′ , z ′′ , . .., while η (n) = η, η ′ , η ′′ , . ... Note that this form, in contrast to the earlier works, removes the factor √ α s from the definition of the variables η and s ij , so that the one-loop running of the coupling can be implemented via (cf.Eq. ( 10)) Since we assume that Λ > Λ QCD , we have s 0 > 0. As all our dipole sizes are smaller than 1/Λ, we see that s 21 > 0, thus avoiding the Landau pole at s 21 = −s 0 < 0 in the coupling.(In general, having an IR cutoff for the dipole sizes, x ij < 1/Λ, implies that all s ij > 0.) Before discretizing our evolution equations, we need to impose the starting value of x for our evolution (cf. Ref. [103]).For z = 1 and x 10 = 1/Q, we have the "rapidity" variable y ≡ η − s 10 = Nc 2π ln 1 x .Hence, if our evolution starts at some value of x labeled by x 0 , then the x < x 0 condition implies that η − s 10 > Nc 2π ln 1 x0 ≡ y 0 .Regarding the value of x 0 , it was observed in Ref. [103], using the older (KPS) version of our helicity evolution, that good-χ2 fits of the polarized DIS data can be obtained with x 0 = 0.1 (and even for a slightly higher values of x 0 ).This is in contrast to the x 0 = 0.01 starting point of the evolution [75][76][77][78][79][80][81][82][83][84] for phenomenological analyses of the unpolarized observables (see, e.g., Refs.[118,119]).As discussed in Sec.III A below, it was speculated in Ref. [103] that such a discrepancy could be attributed to the helicity evolution resumming the double-logarithmic parameter α s ln 2 (1/x) while the unpolarized evolution [77-84, 125, 126] resums single logarithms α s ln(1/x).This way, the resummation parameter for helicity evolution is larger at small x, making the helicity evolution start at larger x values.We thus put x 0 = 0.1 in all our analyses below. 2  The full process of discretizing our flavor singlet and nonsinglet evolution equations with running coupling is detailed in Appendix A. In the end, the discretized version of Eqs.(9) written in terms of the variables (15) reads where the numerical step sizes are chosen such that ∆η = ∆s 10 = ∆s 21 ≡ ∆, and the indices are defined by {η, s 10 , s 21 } → {j, i, k} • ∆.Eqs.(17) allow us to compute the numerical solution for the flavor singlet evolution equations (9).Note that it is only necessary to loop over the ranges dictated by our physical assumptions, 0 ≤ i ≤ k ≤ j ≤ j max and i < j.Furthermore, it is useful to notice that the neighbor dipole amplitudes reduce to their dipole-amplitude counterparts when k = i, that is, We can continue this convention and write the quark and gluon hPDFs from Eqs. ( 5) and ( 6) in the new variables, and where the only difference compared to ∆G from Eq. ( 6) is the running coupling.
The last pieces to consider are the inhomogeneous terms.According to the Born-level initial conditions (11), they can be re-written using our new logarithmic variables as Since Eqs.(21) are linear in η and s 10 , we follow Ref. [103] and employ the linear-expansion ansatz, i.e., Thus, for the three light flavors we consider, q = u, d, s, the full set of initial conditions for the flavor singlet evolution depends on 15 parameters a u , b u , c u , a d , . . ., c 2 which we will fit to the data.Moreover, because the evolution equations we are solving are linear, their solution can be written as a linear combination of 15 "basis" dipole amplitudes, each of which is constructed by performing the iterative calculation outlined above while setting one parameter (from all the a's, b's and c's) in Eqs. ( 22) to be 1 and all the other parameters to 0. Furthermore, since all hPDFs and the g 1 structure function depend linearly on the polarized dipole amplitudes, they are also linear combinations of their corresponding basis functions as well.
For example, ∆u + (x) can be expressed as a linear combination of the 15 "basis hPDFs" shown in Fig. 1.Since ∆u + (x) depends directly on the linear combination Q u + 2G 2 (see Eq. ( 5)), one may expect that Q u and G 2 have the largest contributions to ∆u + (x) at moderate x.This is indeed the case, with the top and bottom panels in Fig. 1 having the largest-magnitude contributions to ∆u + (x).Some of the other amplitudes contribute more significantly at lower x's, as their magnitudes begin to influence those of Q u and/or G 2 through evolution.At the smallest values of x in Fig. 1, the largest contributor is G 2 , followed by G, while the contributions from Q d and Q s remain small for all values of x.
A consequence of this observation, that we will return to later, is that the sign of the g 1 structure function is influenced mainly by the sign of G 2 (or, equivalently, the sign of ∆G) and the sign of G.A challenge for phenomenology presents itself: G is slow to grow and hence less sensitive to available data near x = x 0 , but it has a potentially large effect on the small-x asymptotics.Unless we have sufficient data from an observable that is directly sensitive to G, constraining that amplitude will be difficult.
Similar to the singlet evolution, the discretization of the nonsinglet evolution equation ( 13) reads (again, see Appendix A for details) The corresponding flavor nonsinglet quark hPDF is given by with the integrals also discretized and evaluated numerically.Interested readers are directed to Appendix C for a discussion about convergence testing the numerical solutions of the flavor (non-)singlet evolution equations and the discretized versions of the hPDFs.The Born-level approximation ( 14) is linear in the logarithmic variables (15), so we make a linear expansion ansatz for the inhomogeneous term in the flavor nonsinglet evolution, for each of the three light flavors, q = u, d, s.This means that flavor nonsinglet hPDFs can be reconstructed as a linear combination of 9 flavor nonsinglet basis functions, generated by putting one of the 9 parameters (a NS u , b NS u , . . ., c NS s ) to 1, while setting all others equal to 0. Combining this with the 15 parameters from Eqs. (22) describing the inhomogeneous terms for the flavor singlet dipole amplitudes, we have 24 parameters (and associated basis functions) for the eight amplitudes , which we will fit to describe the world polarized DIS and SIDIS experimental data at low x.

D. SIDIS cross section at small x
We will now derive a formula for the SIDIS structure function g h 1 (x, z) at small x.Using the notation of Ref. [71], we start with the DIS structure function g 1 (x) and write it as where σ ⃗ γ * +⃗ p→X (λ, Σ) is the total virtual photon-proton cross section for the proton with helicity Σ and for the transversely polarized virtual photon with polarization λ, and α em is the fine structure constant.The virtual photonproton cross section is always inelastic at this order in α em , as the virtual photon has to decay into a quark-anti-quark pair, with the quark and anti-quark fragmenting into hadrons in the final state.Consider producing a hadron with a fixed value of z ≡ P • P h /P • q, where P and q are the 4-momenta of the proton and virtual photon, respectively, while P h is the momentum of the detected hadron, as shown in Fig. 2. At high energy/small x we can work in the frame where the proton has a large P + momentum component, while the virtual photon has a large q − momentum component.Then z ≈ P − h /q − is the fraction of the virtual photon's minus momentum carried by the produced hadron.All other components of the hadron's momentum are integrated over.
The SIDIS process at small x.An incoming virtual photon with momentum q decays into a quark-antiquark pair which interacts with the target proton carrying momentum P .The quark and antiquark then fragment into hadrons, and one of these hadrons is detected with momentum P h .
We then write, by analogy to Eq. ( 26), in the collinear approximation [127-129] where k ⊥ and P h⊥ are the transverse momentum vectors for the quark and produced hadron in Fig. 2, while D h/q 1 (z, Q 2 ) is the collinear fragmentation function.The sum q,q goes over the produced quarks and antiquarks.While only quark fragmentation is depicted in Fig. 2, an antiquark could instead fragment there, by reverting the particle number flow direction on the quark line in the diagram.
In arriving at Eq. ( 27) we have employed the aligned jet configuration, dominant in DLA [58,71], in which k − ≈ q − , such that the produced hadron carries the fraction P − h /k − ≈ P − h /q − = z of the quark's momentum.Consequently, we assume that z is not very small, such that the hadron is produced in the forward (virtual photon) direction/current fragmentation region and arises from the fragmentation of the forward-moving quark with 4-momentum k in Fig. 2, and not from the fragmentation of the antiquark, which is separated from the quark by a large rapidity interval.This is similar to the hybrid factorization approach to particle production [130][131][132].(The fragmentation of the antiquark in Fig. 2 would contribute to small-z hadron production, and is neglected here since we are interested in order-one values of z.)In addition, the scale in the argument of the fragmentation function could be chosen to be k 2 ⊥ .However, in our small-x kinematics, the typical value of Integrating Eq. ( 27) over k ⊥ and P h⊥ we obtain Comparing this to Eqs. ( 26) and ( 4), we arrive at reproducing the result in Eq. ( 2) of Ref. [30] (see also Refs.[127,133,134]), derived in the collinear factorization framework.(As we mentioned above, since quarks and antiquarks have different fragmentation functions, the righthand-side of Eq. ( 29) cannot be expressed solely in terms of the ∆q + linear combinations of hPDFs, and the ∆q − functions will enter as well.)We conclude that the expression (29) for the polarized SIDIS structure function is the same in the collinear and small-x formalisms for large z.However, we emphasize that a similar discussion as that surrounding Eqs. ( 4) and ( 7) applies to Eq. ( 29) regarding its interpretation in the LCOT framework as implicitly including higher-order α s corrections. 3

E. Global analysis
Our goal is to describe the world data on the longitudinal double-spin asymmetries in DIS and SIDIS at low x using small-x helicity evolution.We start with the longitudinal DIS asymmetry, A ∥ (see, e.g., Refs.[29,135]), where the arrow ↑ (↓) denotes the lepton spin along (opposite to) the beam direction, and the arrow ⇑ denotes the target polarization along the beam axis.The kinematic variables are given by where y = ν/E is fractional energy transfer of the lepton in the target rest frame, γ 2 = 4M 2 x 2 /Q 2 , and R = σ L /σ T is the ratio of the longitudinal to transverse virtual photoproduction cross sections.When 4M 2 x 2 ≪ Q 2 (γ 2 ≪ 1), we have η ≪ 1 and the virtual photon-target asymmetries are implying Similarly, in polarized SIDIS for the production of a hadron h, the asymmetry A h 1 can be expressed as (see, e.g., Refs.[23,30]) In principle there is another observable in the DIS/SIDIS family that could help constrain hPDFs: parity-violating DIS.This process is sensitive to the g γZ 1 structure function which is approximately proportional to ∆Σ [136,137].Unfortunately there is little to no data for g γZ 1 in the small-x (x < 0.1) region (see, e.g., Ref. [138]), not allowing us to employ this observable in our analysis.
Between the two scattering processes, we have ten unique observables: two in DIS (proton or deuteron/ 3 He target) and eight in SIDIS (proton or deuteron/ 3 He target with charged pion or kaon final states) from which in principle we can constrain the eight polarized dipole amplitudes (five associated with the C-even and flavor singlet hPDFs , and three with the flavor nonsinglet hPDFs (G NS u , G NS d and G NS s )).In our formalism, the g 1 and g h 1 structure functions are calculated in terms of hPDFs using Eqs.( 4), (29), respectively.(Note that ∆q = (∆q + +∆q − )/2 and ∆q = (∆q + − ∆q − )/2.)This is the bridge connecting small-x helicity evolution to the experimental data.Fitting the hPDFs to A ∥ , A 1 and A h 1 at moderate x ≲ 0.1 allows us to determine the initial conditions of the polarized dipole amplitudes (22), (25).We then evolve the polarized dipole amplitudes toward lower values of x using Eqs.( 9) and ( 13) to obtain hPDFs in that region, and compare with existing data as well as make predictions at smaller x.We mention that the structure functions F 1 and F h 1 involve the unpolarized PDF q(x, Q 2 ) and, for the latter, the unpolarized fragmentation function (FF) D h/q 1 (z, Q 2 ).We compute F 1 and F h 1 up to next-to-leading order using 3 Strictly speaking, for consistency the fragmentation functions D h/q 1 (z, Q 2 ) should also be taken in the polarized DIS scheme, but since the only presently available fragmentation functions are given in the MS scheme, we make use of the existing extractions.
collinear factorization and DGLAP evolution, based on the JAM analysis in Ref. [34].(To be consistent, strictly speaking one should include small-x evolution also for F 1 and F h 1 .However, for us the results of Ref. [34] serve as a faithful proxy of the experimental data for these structure functions.A more comprehensive analysis that also utilizes small-x evolution for F 1 and F h 1 is left for future work.)Let us present a short discussion about our ability to constrain G 2 and G, which are two important polarized dipole amplitudes driving the small-x evolution of the hPDFs.The polarized dipole amplitude G 2 is directly related to the gluon hPDF, per Eq.(20).However, the observables we consider here do not directly couple to the gluon hPDF.Instead, as we saw above, they couple only to quark hPDFs.The dipole amplitude G 2 enters the quark hPDFs ∆q + along with the dipole amplitude Q q .Moreover, they always enter in the same linear combination, Q q + 2 G 2 for q = u, d, s (see Eq. ( 19)).We see that while G 2 and Q q couple directly to the spin-dependent structure functions for DIS and SIDIS, we do not have an observable (or a linear combination of observables) in this analysis which separately couples only to G 2 or only to Q q .
What may help us to separate G 2 and Q q is the fact that these dipole amplitudes have a different pre-asymptotic form.While it is established numerically that at asymptotically small x, both polarized dipole amplitudes G 2 and Q q are proportional to the same power of x with the same intercept [102] and are, therefore, probably hard to distinguish, in the pre-asymptotic region where the asymptotic form has not yet been reached, their contributions to the quark hPDFs may be quite different.This can be studied by comparing the Q u and G 2 basis functions for ∆u + in Fig. 1, shown in the top and bottom panels of that figure, respectively.If these functions were identical, they could be freely interchanged against each other while still producing the same structure functions: in such a case it would be impossible to separate G 2 and Q u from the data.Since the contributions of different amplitudes to quark hPDFs differ from each other, as follows from Fig. 1, these basis contributions cannot be adjusted at one value of x while maintaining the same value for the observables at all other x.Therefore, we may be able to separate G 2 and Q u using the polarized DIS and SIDIS data.However, since the Q u and G 2 basis functions have similar shapes, per Fig. 1, it might be the case that the uncertainties in the resulting extractions of Q u and G 2 will be large.
The polarized dipole amplitude G, on the other hand, does not couple to any of the polarized DIS or SIDIS observables we consider here.Rather, it mixes with other polarized dipole amplitudes only through evolution (see Eqs. ( 9)).This is why the G basis function of ∆u + (second from the bottom panel in Fig. 1) appears to be vanishingly small above x > x 0 .The consequence of this is that in the region of x where the polarized DIS and SIDIS data exist, 5 × 10 −3 < x < 0.1, the G amplitude is very small and is, therefore, much less constrained by the data than the Q q and G 2 dipole amplitudes.At small x, however, the G amplitude is quite large, second only to G 2 (see Fig. 1).As we will see below, G, unconstrained by the existing polarized DIS and SIDIS data, will dominate over the other polarized dipole amplitudes at small x, adversely affecting our ability to make precise predictions at even smaller x.Nevertheless, it is possible that G might be constrained with slightly more leverage in x.We will discuss this in Sec.III D when we explore the impact of the future EIC data on our uncertainties.
In our global analysis we use the JAM Bayesian Monte Carlo framework (see, e.g., [29,139,140]) to randomly sample (roughly 500 times) the space of 24 parameters a, b, c from Eqs. ( 22) and ( 25), namely a u , b u , c u , a d , . . ., c NS s .For each combination of these parameters, we solve our evolution equations ( 9) and ( 13) to determine the polarized dipole amplitudes The actual numerical solution is facilitated by the basis functions introduced above.)Next, using Eqs.(19) and (24), we calculate the quark hPDFs at small x, which, via Eqs.( 4) and (29), can be used to determine the structure functions g 1 and g h 1 that enter the numerator of the asymmetries A ∥ , A 1 (Eqs.(32), (33)) and A h 1 (Eq.( 34)), respectively.The χ 2 -minimization procedure allows us to construct the posterior distributions of the parameters, and the corresponding solutions of our evolution equations then allow us to infer the quark and gluon hPDFs (the latter via Eq.( 20)).We confirmed that the posterior distributions of the parameters are distributed more narrowly than the initial flat sampling and are approximately Gaussian, indicating a convergence in their values.These extracted quark and gluon hPDFs, and the quantities that can be computed from them, are the main results of our work, which we present below.

III. RESULTS
In this section we present the results of our numerical analysis.We will concentrate on the proton g 1 structure function, and the quark and gluon hPDFs (along with quantities, such as net spin, that can be computed from them).

Data set (A1)
Target Npts χ 2 /Npts SLAC (E142) [141] 3 He 1 0. Our analysis (JAMsmallx) of the world polarized DIS and SIDIS data at low x utilizes measurements from SLAC [141][142][143][144][145], EMC [146], SMC [147][148][149], COMPASS [150][151][152], and HERMES [153,154] for DIS, and SMC [155], COMPASS [156,157], and HERMES [158,159] for SIDIS.The data of interest falls in the Bjorken-x range of 5 × 10 −3 < x < 0.1 ≡ x 0 , and the Q 2 range is 1.69 GeV 2 < Q 2 < 10.4 GeV 2 .Since x ≈ Q 2 /s, the minimum cut on Q 2 determines the minimum accessible x in the data set (for a given experimental center-of-mass energy), and conversely the maximum cut on x determines the maximum Q 2 .The upper limit on x (denoted by x 0 ) was chosen based on our previous (DIS-only) work [103], as (almost) the highest value of x which gave a good-χ 2 fit.This x 0 is the point where we start the small-x helicity evolution.The fact that our small-x approach was able to describe data up to such a high value of x could be due to the fact that, unlike the unpolarized Balitsky-Fadin-Kuraev-Lipatov (BFKL) [125,126], Balitsky-Kovchegov (BK) [75][76][77][78] and Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) [79][80][81][82][83][84] small-x evolution, which resums powers of α s ln(1/x) at the leading order, our helicity evolution has a different (larger) resummation parameter, α s ln 2 (1/x).For α s ≈ 0.25, our resumation parameter becomes of order 1 for x ≈ 0.1, potentially justifying our use of x 0 = 0.1 as the starting point for our evolution.Note that the value of our resummation parameter α s ln 2 (1/x) at x = x 0 = 0.1 is comparable to (and even slightly larger than) the value of the resummation parameter α s ln(1/x) for the unpolarized small-x evolution at x = 0.01, which is where the latter evolution is usually initiated in phenomenological analyses [118,119].The lower limit of Q 2 is set by the charm quark mass, m 2 c = 1.69 GeV 2 .This is also the cut placed by the JAM FF set we use [34], which has independent functions for π + , K + , h + (π − , K − , h − are found through charge conjugation) that we evolve through the DGLAP equations.By analogy to [103], we choose our IR cutoff to be Λ = 1 GeV.Also, in the Q 2 range specified above, the strong coupling in Eq. ( 16) is taken with N f = 3 (and N c = 3).
The range of the outgoing hadron momentum fraction z in polarized SIDIS is 0.2 < z < 1.0, and we do not place any explicit cut on this variable.In practice, the data (after all the appropriate cuts) generally has values of 0.4 < z < 0.6 ; some data sets integrate z ∈ [0.2, 1] while others cover z ∈ [0.2, 0.85].After all the cuts we are left with 122 polarized DIS data points and 104 polarized SIDIS data points, for a total N pts = 226.The overall χ 2 /N pts of our fit, based on the central theory curves, is 1.03.(We have also performed fits with cutoffs of x 0 = 0.08 and x 0 = 0.05, which produced no significant change in χ 2 /N pts .)The breakdown of the data by experiment, along with our χ 2 /N pts for those individual data sets, is shown in Table I for DIS and Table II for SIDIS.The plots of the experimental data versus our JAMsmallx theory are shown in Fig. 3 for polarized DIS and Fig. 4 for polarized SIDIS.Overall, our results demonstrate very good agreement with the existing world data.

B. Proton g1 structure function
We now examine our result for the g 1 structure function of the proton to analyze the predictive capability of our formalism.Our calculation of g p 1 for all replicas is given in Fig. 5.This is the result of 500 individual fits of the experimental data where the (quark and gluon) hPDFs were extracted and then (the quark ones) used to compute We color code each replica by its asymptotic sign at small x in order to clarify the structure of the plot as well as to help establish correlations with the hPDFs below.While g p 1 is well constrained in the region where there is experimental data (5×10 −3 < x < 10 −1 ), it is largely unconstrained at smaller x.The major difficulty in constraining g p 1 is caused by the insensitivity of the data to the G 2 and G amplitudes described above.
That being said, the asymptotic solution of the large-N c &N f evolution equations [102] guarantees that the small x behavior of g p 1 must be exponential in ln(1/x).This implies that it has to pick a sign (positive or negative) when x → 0. Our results indicate (see Fig. 7) that, given the existing experimental data constraining our formalism, the asymptotic sign is likely to be picked by x = 3.5×10 −4 with 10% uncertainty, with the uncertainty decreasing to 5% at approximately x = 2.5 × 10 −5 .Currently, 70% of the replicas are asymptotically positive and 30% are asymptotically negative.These percentages are stable as the number of replicas increases.The primary source of uncertainty is how low in x one must go to determine the sign, as some replicas that appear positive may undergo a sign change at smaller x.Interestingly, our observation of a preference for g p 1 to be positive at small x agrees with the recent papers analyzing (unpolarized and polarized) DIS structure functions using the anti-de Sitter space/Conformal Field Theory (AdS/CFT) correspondence [160][161][162] that make an even stronger statement that g p 1 clearly grows positive at small x.This behavior also has implications for the net parton spin expected at small x, as we discuss in Sec.III C. FIG. 5.The small-x calculation of the g1 structure function of the proton.The black curve is the mean of all the replicas with the green band giving the 1σ uncertainty.Red (blue) curves are solutions that are asymptotically positive (negative).

Sign of g p 1 and quantifying numerical ambiguity
From Fig. 5 alone, one can make the qualitative observation that indeed each replica of g p 1 grows exponentially with ln(1/x) as we suggested earlier, and the color indicates the asymptotic sign of g p 1 for that given replica.We mentioned in the previous section that the exponential behavior of helicity functions in our theory makes it difficult for a given replica to maintain a near-zero value, and thus it must eventually choose to (rapidly) increase in magnitude towards positive or negative values.Given the numerical nature of our global analysis, we cannot compute each fitted replica down to x = 0 (corresponding to ln x → −∞), so the color-coding and sign assignment is determined by the slope of a replica at the lowest-computed value of x: if the slope increases (decreases) as x goes to zero, then it is considered "asymptotically" positive (negative).To balance our time and computational resources, the results discussed in this section use replica data computed down to x asymp = 10 −7.5 .One may realize potential issues with this system: a given replica may have multiple different "asymptotic" signs depending on the lowest computed value of x.
Any given replica is defined by its specific combination of basis functions, and since our Bayesian analysis samples parameters (Eq.( 21)) that may be either positive or negative, competition between basis functions can result in nodes.Replicas with two nodes in g p 1 (x), such as the one illustrated in Fig. 6, can occur for linear combinations of similar basis functions with opposite signs, as in the top/bottom panels of Fig. 1.These changes in sign can occur at various values of x depending on the initial conditions, making the prediction of the asymptotic sign dependent on what x value is used to make the prediction.
Careful readers may have already noticed this from Fig. 5, where there are a few red-coded replicas that appear to be growing negative (and a blue-coded replica that appears to be growing positive) at x = 10 −5 .This is due to each of these replicas having a delayed critical point ( = 0) that occurs at x < 10 −5 , where a different basis function takes over the growth and the replica changes the sign of its slope.These critical points also are connected to the issue of ambiguity, where at a specific value of x we may be able to measure that a replica is growing positive (or negative) but has a magnitude that is actively negative (or positive), leaving its asymptotic sign unconfirmed.Luckily, investigations of these incidents show that they occur in a statistically small portion of replicas from the perspective of our considerably small x asymp .
Since our goal is predictability at small x, we decided to quantify the amount of ambiguity by its probability density in x.That is, for each replica we count the smallest-x instance of ambiguity, and take note of where in x it occurred.For example, Fig. 6 shows a replica that begins positive (true for all replicas) and evolution drives it more positive until it reaches a critical point, after which the replica then grows negative.After the critical point (in the gray region), the replica will be considered ambiguous until it crosses g p 1 (x 1 ) = 0, and then it is considered asymptotically negative (in the blue region).Only when the sign of g p 1 and the sign of its first derivative (as x decreases) agree can the replica be considered asymptotically positive or negative.If we wanted to predict the asymptotic sign of the FIG. 6.An example replica of g p 1 (x, Q 2 = 10 GeV 2 ) that demonstrates how the asymptotic sign is dependent on x pred .If x pred resides in the red (blue) region then the replica will be considered asymptotically positive (negative) according to the sign of the first derivative (for decreasing x), and its agreement with the sign of the magnitude.If x pred resides in either gray region then the asymptotic sign is ambiguous due to a contradiction between the sign of the slope and the sign of the magnitude.replica based on an observation at x = x pred that resides in this (blue) region, then we would predict that this replica is "asymptotically negative" as x → 0. However, this same replica has a small-x critical point (around x = 10 −4 ) that causes the sign of its slope to change; the replica observed in the (gray) region (on the left) between the critical point and g p 1 (x 2 ) = 0 would be considered ambiguous again.After crossing zero a second time, a prediction made at x pred < x 2 would therefore designate the replica to be "asymptotically positive."The smallest-x instance of ambiguity is thus counted in a bin at x 2 .In this way, each replica is counted exactly once, and replicas that oscillate multiple times about the g p 1 = 0 axis only have their most delayed ambiguity counted.We can define the number of replicas that have their smallest-x instance of ambiguity in a particular bin of x as C A (x) (the counts of ambiguities) and make a histogram.The ambiguity count C A (x) is normalized such that it sums to the total number of replicas N ambig containing at least one ambiguity: Because some replicas are always unambiguous across the entire range of x, the ambiguity count is less than the total number of replicas: N ambig ≤ N tot .Now suppose we want to predict the asymptotic behavior of g p 1 at small x based on the behavior of the function at some value x pred .Knowledge of the ambiguity count C A (x) allows us to estimate the accuracy of this prediction by estimating the probability that an unobserved ambiguity remains at x asymp < x < x pred .This probability is given by a summation as in Eq. ( 35), but over the truncated range in x: From the normalization condition (35), we see that Eq. ( 36) implies the truncated moment is normalized at x pred = x 0 to the total fraction of replicas containing at least one ambiguity: From the left panel of Fig. 7 we see that the number of smallest-x ambiguities decreases greatly as x approaches zero.
The right panel shows we must go down to approximately x = 3.5 × 10 −4 , 2.5 × 10 −5 , and 6 × 10 −7 to capture the asymptotic sign with 10%, 5%, and 1% uncertainty, respectively.This is strong justification that x asymp = 10 −7.5 is reasonably low enough to capture the asymptotic sign of our replicas with low uncertainty.Due to Eq. ( 37) we also know how many replicas are completely unambiguous; since we impose our evolution to begin at x 0 = 0.1, the running integral at that point quantifies the total ratio of replicas that have at least one ambiguity.According to the right panel of Fig. 7, approximately 50% of replicas choose their asymptotic sign immediately as evolution begins.
Note that the data constrains the initial condition for g p 1 to be positive, so all completely unambiguous replicas are asymptotically positive.
Furthermore, splitting the replicas by their asymptotic sign (not shown in Fig. 7) allows us to also investigate how early (or late) the different solutions are chosen relative to each other.We gather that ambiguously negative replicas tend to choose their sign earlier than their positive counterparts, with the caveat that the majority of asymptotically positive replicas do not have any ambiguities at all.Approximately 75% of asymptotically positive replicas are completely unambiguous, and the remaining 25% are determined by x ≈ 2 × 10 −5 with 5% uncertainty.Though fewer in number, a still significant portion of replicas are asymptotically negative, 95% of which are confirmed by x ≈ 4.3×10 −4 .This suggests that using a lower x pred will affect the positive-identified and negative-identified solutions differently.In particular, a lower x pred is likely to identify a greater number of asymptotically positive solutions by correcting replicas that would have been misidentified as asymptotically negative at a higher x pred .This asymmetric impact on positive-identified versus negative-identified solutions can be traced back to constraints from the data at large x, which strongly prefers g p 1 > 0. The fact that this positive preference persists down to small x suggests that the polarized dipole(s) which dominate the small-x asymptotics are partially (but not fully) constrained by the large-x data.This will be discussed in detail in Sec.III B 3. We performed a similar analysis of the smallest-x critical points of each replica (rather than the ambiguities).On average, the smallest-x critical point occurs 4% earlier in ln(1/x) than its smallest-x zero.Since the ambiguous region of a replica is precisely the region in x between its critical point and zero, this small 4% difference indicates that any remaining ambiguities are quickly resolved at small x.Thus, we conclude that, from the perspective of Fig. 7, if we had data down to x ≈ 10 −5 we could determine the asymptotic sign of g p 1 with high certainty (> 95%).

Asymptotic behavior of g p 1
Collectively utilizing the information in Figs. 5 and 7 paints a curious picture: there are many more g p 1 replicas that adopt their asymptotic forms early than there are replicas that change their signs at small x.This results in some clustering behavior, e.g., in the left panel of Fig. 7 there is a cluster of replicas around x = 5 × 10 −3 , implying that these replicas share similar critical points and rates of growth.As mentioned previously, the majority of replicas have no ambiguities and adopt their asymptotic growth rather quickly, effectively clustering their critical points at x = x 0 ≡ 0.1 (not explicitly shown).This behavior supports the idea that early adoption of asymptotic growth is preferred, whereas replicas with late critical points are fewer in nature.Consequently, we expect that there should be a form of bimodality in g p 1 between the rapidly growing positive solutions versus the rapidly growing negative solutions.This is a novel result, which we quantitatively analyze below.
While Fig. 5 may appear to show the anticipated bimodality (red versus blue curves), upon closer inspection the values of g p 1 are normally distributed, both at small x (x = 10 −3 ) and very small x (x = 10 −7.45 ), as depicted in Fig. 8.To uncover the bimodal behavior it is necessary to construct a new observable related to the curvature of g p 1 which is sensitive to how quickly our evolution equations drive the g p 1 replicas toward the asymptotic limit.The emphasis, therefore, is not so much on g p 1 as on the exponent of its power-law behavior at small-x, i.e., g p 1 (x) ∼ x −α h .
The generalized x-dependent exponent α h (x) can be extracted through the logarithmic derivative of g p 1 : where g p (0) 1 = const.Examining the distribution of α h (x) across replicas can provide complementary information to the distribution of g p 1 (x) itself.Notably, the exponent provides a meaningful way to scale the solutions: if they have the same α h (x), they have the same curvature, whether the magnitude of g p 1 (x) is large or small.To further capture the signed behavior of g p 1 (x) and distinguish between solutions trending positive or negative at small x, we can generalize the logarithmic derivative (38) to reflect the sign of g p 1 itself: Both the effective exponent α h (x) (38) and its signed generalization (39) are shown in Fig. 9 at varying values of x (from the same global fit that produced Fig. 5).(We remark that if a g p 1 replica has a delayed critical point it will result in a delayed zero that may cause an artificially large ratio if g p ′ 1 (x) ≫ g p 1 (x) ≈ 0. In order to avoid these statistical outliers, any replica with a ratio value outside of 5σ from the average are omitted from the results in Fig. 9.) The distribution in the right panel at x = 10 −2 (blue histogram) is skew-normal, which is expected since we are definitively outside of the asymptotic regime.However, at x = 10 −3 (yellow histogram) we already see the formation of two separated peaks, one positive and one negative.As x continues to decrease down to x = 10 −5 (green histogram), the two peaks become more refined as the evolution equations predict specific curvature related to the intercept α h (see Eq. ( 39)).Without the sign dependence, as displayed in the left panel of Fig. 9, as x → x asymp , a single peak emerges that approaches the expected asymptotic value for α h .The decreasing uncertainties are a consequence of our small-x evolution, where the predictive power constrains the value of α h (x).
From the perspective of the right panel of Fig. 9 it appears that data sensitive to this curvature at x as large as x = 10 −3 may be enough to identify which bimodal peak g p 1 belongs to.Unambiguously identifying this curvature will provide us the asymptotic sign of g p 1 as well as the asymptotic sign of all the (flavor singlet and C-even) hPDFs, as will be discussed below.The fact that such a conclusion could be made at x ≈ 10 −3 by analyzing the curvature of g p 1 (x), compared to x ≈ 10 −5 by studying g p 1 (x) itself (see the discussion around Fig. 7), makes the idea of curvature a useful quantity to consider once future low-x data is available from the EIC.

Origins of asymptotic behavior
To understand what differentiates the positively and negatively growing solutions for g p 1 displayed in Fig. 5, we examine the polarized dipole amplitude parameters themselves, defined in Eqs.(22).We note that the experimental data is only sensitive to the polarized dipole amplitudes as a whole, and not to any specific basis function.For example, combining Eqs. ( 4), (5), and (32) shows that A 1 is constructed from the dipole amplitudes Q q and G 2 , and any combination of parameters that reconstructs the experimental data with good χ 2 is equally valid.An appropriate change of variables can reorganize the basis hPDFs to increase the sensitivity to their overall sign.We can then   (38) showing that as x decreases, the intercept α h (x) becomes more constrained as a consequence of the small-x evolution equations.(Right) Keeping information on the sign dependence by using Eq. ( 39) produces bimodal peaks at ±α h (x).At large x there is no asymptotic behavior, and for smaller values of x two refined peaks emerge.
classify which of these parameters are most correlated with the asymptotic sign of g p 1 .We find enhanced sensitivity to the asymptotic sign of g p 1 from the linear combinations a ′ ≡ (a + b)/2 and b ′ ≡ (a − b)/2.Then the dipole initial condition G (0) = a η + b s 10 + c can be written as These new basis functions are displayed in Fig. 10.Compared to Fig. 1, the alternative parameters a ′ , b ′ change the shape of the basis hPDFs.In particular, we note that this greatly increases the separation between the a ′ 2 = 1 and b ′ 2 = 1 basis functions at large x where the data provides constraints.When we bin the replicas into asymptotically positive/negative g p 1 at small x, we find that the parameter with the largest difference between the solutions is ã′ .The asymptotically-positive solutions preferred a negative parameter ã′ = −1.56 ± 2.32, while the asymptotically-negative solutions preferred the positive ã′ = 1.42 ± 2.34.No other systematic differences in parameters were observed.
We can understand from the basis hPDFs shown in Fig. 10 why asymptotically-positive/negative g p 1 correlates, respectively, with negative/positive values of ã′ , and why ã′ shows the greatest discrimination power.First, we note that the basis hPDFs themselves are negative-definite functions of x for positive values of the initial parameters a ′ , b ′ , c, which is simply a consequence of the explicit minus sign in Eq. ( 5).Second, we note that the hPDFs arising from both the G (0) (with parameters ã′ , b′ , c) and G Comparing the effect G and G2 has on the overall sign of g p 1 (x) at small-x.Top row: the priors are restricted so that (left) G ≤ 0 and (right) G ≥ 0. Bottom row: the priors are restricted so that (left) G2 < G = 0, and (right) G2 > G = 0.All other parameters initially are randomly sampled just as they were in the fit shown in Fig. 5.We see that controlling the sign of G strongly influences the sign of g p 1 , and that the sign of G2 will also influence the sign of g p 1 .
large at small x; the a ′ 2 = 1 basis function also being sizeable at large x, whereas the a ′ basis function only contributes meaningfully at small-x.The large-x behaviour means that the parameter a ′ 2 , while important for determining the small-x asymptotics, is constrained by higher-x experimental data, and it specifically prefers negative values: a ′ 2 = −0.98 ± 1.00.The origin of the different asymptotic behaviors seen in Fig. 5 therefore appears to be due to the dipole G, which makes no contribution to the basis hPDFs at larger x, and, thus, the sign of ã′ evades experimental constraints.
To test this hypothesis, we ran fits where all of the G initial condition parameters (ã, b, c) were restricted to either be negative-definite or positive-definite, with all other parameters unchanged.All g p 1 replicas in the negative-definite G fit were asymptotically positive.The positive-definite G fit was slightly less selective but still generated a 73% majority preferring asymptotically negative g p 1 replicas (recall the original fit in Fig. 5 had a 70% positive preference).The results, shown in the top row of Fig. 11, clearly demonstrate that the sign of the G dipole amplitude determines the small-x asymptotics of g p 1 , as anticipated by the basis functions in Figs. 1 and 10.The reason G leads to a g p 1 that is poorly constrained at small x can be seen directly from Eqs. ( 4)-( 6), ( 9) and Eqs. ( 12), ( 13), (29): G does not contribute directly to any hPDF.Whereas all the other (non-neighbor) polarized dipole amplitudes directly enter a DIS/SIDIS observable, the effects of G are only felt indirectly through its impact on the evolution of the other amplitudes.As a result, hPDFs mediated by G only become large at very small x (see the top panel of Fig. 10), where there are no constraints from data.
While G is the driving factor in determining the small-x asymptotics of g p 1 , G 2 also plays a role.In fact, if G was removed, G 2 would be the most important amplitude in controlling the small-x asymptotics of g p 1 .We see this explicitly when setting the initial conditions for G all to zero (ã = b = c = 0) and repeating the previous analysis of now restricting the G 2 initial condition parameters to be always positive or always negative.The result, shown in the bottom panel of Fig. 11, confirms that, although constrained by large-x data, G 2 plays the second-most important role after G in determining the small-x asymptotics of g p 1 .The negative-definite G 2 fit was 100% selective of asymptotically 12. Color coding the hPDF replicas according to the asympotic sign of g p 1 shows that there is a novel correlation: at small x, quark hPDFs (left) have the same sign as g p 1 (only ∆u + is shown) while the gluon hPDF (right) has the opposite sign as g p 1 .
positive g p 1 replicas, while the positive-definite G 2 fit was 96% selective of asymptotically negative g p 1 replicas.Fig. 11 then compactly summarizes the origin of the asymptotic behavior seen in Fig. 5.The origin of the huge uncertainty band at small x is due to the inability to constrain the sign of G from large-x data, and the overall preference of the central curve in Fig. 5 favoring positive solutions is due to the fact that there is an experimental constraint which prefers G 2 < 0, leading to g p 1 > 0. Knowing now that the dipole amplitude G controls the small-x asymptotics of g p 1 gives us powerful insight into the hPDF correlations which characterize the fits.Comparing Eqs. ( 4), ( 5), and ( 6) we can draw the conclusion that at asymptotically small x these quantities are simply related by where the last step in each line represents the fact that the evolution of Q q and G 2 is driven by G (see Eqs. ( 9)).At small x, the two hPDFs ∆q + and ∆G are both driven by the same polarized dipole amplitude G, but have opposite signs.Since g p 1 is proportional to ∆q + (weighted by quark electric charge squared and summed over flavors), it follows that if the quark hPDFs for all flavors have the same sign, then, at small x, g p 1 will have the same sign as the quark hPDFs and opposite sign as the gluon hPDF.These anticipated (anti)correlations among the hPDFs are shown in Fig. 12, where we plot only ∆u + and ∆G for brevity.Note that the color coding used for the replicas in Fig. 12 indicates the ultimate asymptotic sign of g p 1 , not the hPDF itself.That is, an hPDF replica is colored red (blue) if the corresponding g p 1 replica is asymptotically positive (negative).The fact that the asymptotic signs of ∆q + and ∆G are, respectively, correlated and anti-correlated to the sign of g p 1 at small x is a robust, novel prediction of the small-x helicity evolution framework. 4, 5Thus, in order to better predict the asymptotic sign of g p 1 , ∆q + and ∆G, we need to better constrain the polarized dipole amplitude G.One option is data from the future EIC, discussed in the Sec.III D. We also outline several additional ways in Sec.III E.
C. Extracted helicity PDFs and calculation of net parton spin and axial-vector charges at small x Our results for the hPDFs are shown in Fig. 13.Since our small-x analysis is only valid for x < x 0 = 0.1, we restrict the plots to that region.As with the g p 1 structure function shown in Fig. 5, the hPDFs themselves also exhibit broad uncertainty bands at small x. 6 The uncertainty bands for all four hPDFs span zero below x ≲ 10 −3 , indicating that the hPDFs in that region may be positive, negative, or consistent with zero.By far the largest uncertainty is seen in ∆G, which, unlike ∆q + , is not directly sensitive to inclusive DIS constraints on g p 1 (Eq.( 4)).As shown in Figs.11  and 12, the large uncertainty in ∆G is due to the lack of sufficient constraints on the dipole amplitudes G and G 2 that dominate both ∆q + and ∆G at small x.This conclusion is further supported by the left panel of Fig. 13, where ∆u + , ∆d + and ∆s + exhibit approximately the same error band below x ≈ 10 −4 .At larger x, where the hPDF behavior is driven more by the Q q dipole amplitudes, we can observe flavor separation between the three quarks.The uncertainty of the ∆s + distribution then becomes much larger than that for ∆u + and ∆d + , most likely due to the limited SIDIS kaon data.The similar error bands at small x for ∆u + , ∆d + and ∆s + are in contrast to markedly distinct error bands for ∆u − , ∆d − and ∆s − , shown in the right panel of Fig. 13, which exhibit significant flavor separation even down to small x.Recall that the flavor nonsinglet hPDFs are driven by a different polarized dipole amplitude, G NS (see Eq. ( 12)), which is sensitive to flavor separation through the SIDIS data.As a result of the different evolution, the x∆q − distributions converge quickly to zero at small x, unlike the x∆q + distributions, due to the smaller intercept at small x (see also Appendix B).The similarity of the error bands for ∆u + , ∆d + and ∆s + appears to be driven by the error band of the polarized dipole amplitude G 2 , which affects all quark flavors in the same way, per Eq. ( 19).Consequently, additional input which can better constrain G and/or G 2 may well reduce this uncertainty by forcing the hPDFs to choose a definite sign at small x.We discuss possible strategies to achieve this in Sec.III E.
One feature of note in our hPDFs from Fig. 13 is that ∆s + and ∆G are much larger in magnitude than the same hPDFs obtained in the JAM framework using the DGLAP-based approach [30,33,34].In particular, our extracted ∆s + distribution is below zero at about the 1σ level at x ≈ 10 −2 .This is to be compared with Fig. 6 of Ref. [33], which exhibits a ∆s + consistent with zero across the entire considered range 5 × 10 −3 ≤ x ≤ 0.9.Note that the global analyses conducted in Refs.[30,33,34] are quite different than the one we present here, e.g., they use DGLAP evolution within collinear factorization, include data across the full range of x, and in some cases impose SU(2) and SU(3) flavor symmetries.Nevertheless, it is a valuable cross-check to see whether zero strangeness polarization is consistent with our results as well.To that end, we have separately re-fit the data, setting the strangeness polarization identically to zero: ∆s + (x, Q 2 ) = ∆s − (x, Q 2 ) = 0.The overall quality χ 2 /N pts = 1.04 of the zero-strangeness fit is slightly worse than the quality χ 2 /N pts = 1.03 of the default fit, with the asymmetries A h 1 from tagged kaon SIDIS being the most affected by the change.For that subset of the data, the quality-of-fit degraded from χ 2 /N pts = 0.81 in the default fit to χ 2 /N pts = 1.05 in the zero-strangeness fit.This marginal degradation of the fit quality is consistent with the 1σ departure of ∆s + from zero preferred by the default fit in Fig. 13, with the tagged kaon data only accounting for 26/226 data points in total.Therefore, we conclude that small ∆s + is indeed consistent with our formalism, and that there is a real (but weak) preference from the data for nonzero ∆s + at x ∼ 0.01 within our small-x framework.
Next, we address the contribution to the proton spin and axial-vector charges from small x.The flavor singlet quark helicity distribution is given by for the light flavors considered in this work.Using the hPDFs in Fig. 13, we show x∆Σ(x, Q 2 ) in Fig. 14.Again, the uncertainty band at small x based on current experimental data is rather wide, spanning zero so that the sign of ∆Σ is uncertain.
From ∆Σ(x, Q 2 ) and ∆G(x, Q 2 ) we can determine how much net parton spin (see Eq. ( 2)) resides at small x by computing truncated moments of the distributions.We can similarly determine the small-x contributions to the triplet g A and octet a 8 axial-vector charges from truncated moments of the appropriate linear combinations of quark hPDFs.Focusing on the x region 10 −5 ≤ x ≤ 10 −1 of our analysis, we consider the following truncated moments: Here we consider two representations of the truncated moments: either as a function of the upper limit x max with fixed lower limit 10 −5 , or as a function of the lower limit x min with fixed upper limit 0.1.That is, in the notation of Eq. ( 43), we have (x 1 , x 2 ) = (10 −5 , x max ) for [x max ] and (x 1 , x 2 ) = (x min , 0.1) for [x min ].We have also dropped the Q 2 dependence of the truncated moments on the left-hand side of Eq. ( 43) for brevity.
Both [x max ] and [x min ] representations of the truncated moments are plotted in Fig. 15.From the truncated moment of the total parton helicity 1  2 ∆Σ + ∆G [x max(min) ] , we conclude that, despite the sizable uncertainties, the amount of the proton spin coming from the net spin of small-x partons could be quite large.The outer bounds of these truncated moments also allow for the possibility that the net quark and gluon spin contained within the small-x region may be even more significant than what has been computed at large x.We observe that, despite the wide error bands in ∆G(x, Q 2 ) and ∆Σ(x, Q 2 ) separately, the error in the truncated moment 1  2 ∆Σ + ∆G is narrower than if the two were uncorrelated.Because of the replica-by-replica anticorrelation between ∆q + (x, Q 2 ) and ∆G(x, Q 2 ) seen in Fig. 12, there is a systematic cancellation between them, resulting in a truncated moment 1  2 ∆Σ + ∆G which skews net negative and is more tightly constrained than either ∆Σ(x, Q 2 ) or ∆G(x, Q 2 ) alone.In addition, the nonzero slope of 1  2 ∆Σ + ∆G [xmax] as one approaches x max = 10 −5 indicates that this truncated moment has not fully saturated at that point in x.In contrast, the small-x contribution to g A and a 8 appears to saturate around x = 10 −4 , giving a finite, non-negligible contribution from small-x partons.
Taken at face value, our formalism strikingly predicts a negative contribution to the proton spin from the net spin of small-x partons even when accounting for the 1σ error band.In this scenario favored by our default fit, a significant positive contribution from orbital angular momentum would be needed to satisfy the Jaffe-Manohar sum rule (1).Interestingly, similar observations have been made in using AdS/CFT to analyze g p 1 [160][161][162][163].We also predict that approximately 15-21% of the known value of g A and 12-77% of the known value of a 8 are generated from partons with 10 −5 ≤ x ≤ 10 −1 , where the values of the moments over the full range x ∈ [0, 1] are known from neutron and hyperon β-decays [24]: g A = 1.269(3) and a 8 = 0.586 (31).
However, we caution the reader that our small-x analysis is strongly dependent on the large-x initial conditions to our evolution, and that the error bands shown throughout this work are strictly statistical in nature.These are an accurate representation of the uncertainty coming from the experimental data and from the Monte Carlo sampling procedure, but in particular they do not reflect the systematic bias that comes from omitting large-x data that cannot be captured in this formalism.Combining our small-x evolution equations with external input from large x can therefore possibly result in large, systematic changes to the extracted hPDFs beyond the 1σ statistical error bands.This suggests that an appropriate matching procedure onto hPDFs extracted from a large-x, DGLAP-based analysis like JAM [30,33,34] will be crucial to determining the proton spin budget.Moreover, since JAM found both viable positive ∆G(x, Q 2 ) and negative ∆G(x, Q 2 ) solutions [33,34], the predictions for the small-x truncated moments may even depend on which large-x solution is chosen for the matching.Indeed, as we show in Fig. 18 below, matching to the positive gluon hPDF solution could lead to a substantially different outcome for ∆G(x, Q 2 ), deviating beyond the 1σ error band over a significant range of x.Clearly a rigorous implementation of such a matching will be an important aspect of future analyses; a first attempt is detailed in Section III E below.Having emphasized this vital caveat, we summarize our results for the small-x truncated moments of ( 12 ∆Σ + ∆G)(x, Q 2 ), g A (x, Q 2 ) and a 8 (x, Q 2 ) over the small-x window x ∈ [10 −5 , 0.1] for Q 2 = 10 GeV 2 : a few options to constrain it.The first such constraint is positivity, which is the statement that the number densities for positive and negative helicity partons are positive.In particular, for gluons this leads to where G(x, Q 2 ) is the unpolarized gluon PDF.(We will set aside issues as to whether Eq. ( 45) is strictly satisfied under (MS) renormalization [166][167][168].)We impose this constraint by checking the value of ∆G(x, Q 2 ) in the region x < x 0 = 0.1, and punishing the χ 2 of the fit if the positivity constraint is violated.Unfortunately, by the time our evolution begins, our baseline fit for ∆G(x, Q 2 ) and the JAM DGLAP-based G(x, Q 2 ) [33,34] are of comparable size.The latter grows much faster at small x than our extraction for ∆G(x, Q 2 ), causing the positivity constraint to have a negligible effect.This is perhaps not surprising, given that at small x the unpolarized gluon distribution G(x, Q 2 ) is eikonal, while ∆G(x, Q 2 ) is sub-eikonal, and, hence, suppressed by a power of x.
Another constraint on ∆G(x, Q 2 ) that we explored was a preliminary matching onto the (large-x) JAM DGLAPbased extraction of ∆G(x, Q 2 ) in Refs.[33,34], in particular the SU(3)+positivity scenario.The result is shown in Fig. 18; the red box is bounded by 10 −1.3 < x < 10 −1 and 0.05 < ∆G(x, Q 2 ) < 0.2.The motivation is that any complete description of ∆G(x, Q 2 ) should agree with DGLAP extractions in this region.The matching is performed in a simple way, by choosing an intermediate region in x and forcing our fit of ∆G(x, Q 2 ) to qualitatively agree with the JAM DGLAP-based extraction.This is done in a similar way to the positivity constraint described above, whereby we punish the χ 2 whenever ∆G(x, Q 2 ) strays outside of the matching region (red rectangle in Fig. 18).This constraint causes our extracted ∆G(x, Q 2 ) to take on mostly positive values at small x, seemingly changing sign from our original extraction.However, note that while the baseline extraction uncertainty band grew negative for large x, there were still a significant number of replicas (with good χ 2 ) that grew positive at large x and overlapped with the red region.Forcing ∆G(x, Q 2 ) to pass through that area then preferentially selects those replicas.Consequently, the whole uncertainty band for ∆G(x, Q 2 ) remains shifted upward even in the small-x region.Given that g p 1 (x, Q 2 ) ∝ −∆G(x, Q 2 ) (see Eq. ( 41)), the matching constraint leads to a quantitative change to the distribution of g p 1 replicas: they are now 40% positive and 60% negative.As we emphasized previously, input on hPDFs from large x can have a significant effect on predictions made at small x, motivating future work into a more rigorous matching to DGLAP-based hPDF fits.
Furthermore, the issue with constraining G could be alleviated by a more rigorous way of handling the starting point of evolution x 0 .In this work, we chose x 0 = 0.1 and then used experimental data to fit initial conditions for the polarized dipole amplitudes in order to obtain the correct starting values for all of the extracted hPDFs.Only after these starting values have been determined do we then evolve the distributions in a region dominated by our double logarithmic resummation.In reality, evolution in x begins at x = 1, but is sub-leading, with the dominant contribution at large x given by DGLAP-driven large-x dynamics.The method of matched asymptotic FIG. 17.Relative uncertainty for both this work (red) and a JAM DGLAP-based extraction [165] (blue) for EIC impact studies using the high g p 1 scenario.Dotted lines denote extrapolating beyond lowest x for which pseudodata was generated.For this work, pseudodata was generated down to x = 10 −4 .For the JAM DGLAP-based fit, pseudodata was generated down to x = 2 × 10 −4 [165].
expansions [169,170] suggests that we start the evolution at x 0 = 1, include the DGLAP contributions, but subtract the double counting of logarithms that are present in both resummations.By starting evolution earlier, G might become more sensitive to the data.As discussed at the end of Sec.II E, the challenge in constraining G stems from the fact that it has a small magnitude in the region where there are measurements (see Fig. 1).The magnitude of the G contribution to ∆u + is so small at larger x partly because G enters only through evolution, and evolution is delayed until x 0 = 0.1.If x 0 = 1, G will start growing sooner, and it might then have a large enough contribution to be sensitive to the experimental data.
Moreover, perhaps the most direct way to constrain ∆G is to include in the analysis an observable directly sensitive to it.(Recall that in the polarized DIS and SIDIS processes considered here the contribution from the gluon hPDF is suppressed by a factor of α s .)Two possibilities, which have been used in DGLAP-based extractions [23,26,27,33,34], are jet and hadron production in polarized proton-proton collisions.The numerator of the double-longitudinal spin asymmetry A LL in ⃗ p + ⃗ p collisions takes the following form where ∆f is the parton hPDF for either the quarks or gluon, a(b) is the parton coming from proton A(B), and σ ab is the partonic cross section of parton a interacting with parton b.For hadron production, Eq. ( 46) needs also to be convoluted with the D 1 FF.More work is needed to derive an analogue of Eq. ( 46) in the KPS-CTT small-x evolution framework, and initial developments can be found in Ref. [98].Lastly, in the future, it will also be interesting to attempt to constrain the large-x behavior of the hPDFs by direct matching onto nonperturbative calculations from lattice QCD.Such matching in the vicinity of x ∼ 0.1 is actually feasible for the double-logarithmic helicity evolution, unlike for the case of single-logarithmic unpolarized small-x evolution which would require reliable lattice data down to much smaller x.In addition, recently a new approach to determining the initial conditions for small-x evolution by starting at the level of the proton wave function has been developed in Ref. [171].While that work was done in the context of unpolarized small-x evolution, it is possible that it could be extended to the polarized case, helping us constrain the initial conditions for helicity evolution at hand.

IV. CONCLUSIONS
In this paper we have presented the first phenomenological implementation of the KPS-CTT theoretical framework [58,64,71] for the evolution of hPDFs.This work represents a significant improvement over our previous study [103] by utilizing the revised evolution equations instead of the original KPS equations.On top of that, we have adopted the large-N c &N f limit, which enables a more realistic description of the physics, now including quarks FIG.18.The result of matching onto the ∆G(x) extraction from DGLAP [33,34] at intermediate x.The green band is our baseline fit.The blue band is the result of matching.The light red square is the region where we enforce matching.
in addition to gluons.Another key advancement of this research is an expansion of our analysis beyond just polarized DIS data by also incorporating polarized semi-inclusive DIS measurements.This allowed us to extract both the C-even and C-odd quark hPDFs ∆q + and ∆q − , along with the gluon hPDF ∆G.To extract ∆q − we had to, for the first time, implement the numeric solution for the KPS evolution of the nonsinglet hPDFs.Moreover, we have included running coupling corrections in the evolution of ∆q + , ∆q − , and ∆G, which is another feature of the analysis that makes our approach more rigorous.
Through the application of the JAM Bayesian Monte Carlo framework, we have successfully described all available polarized DIS and SIDIS data below the threshold x 0 = 0.1, achieving a very good fit with χ 2 /N pts = 1.03.However, when attempting to extend our predictions to lower values of x, the uncertainty associated with our results was found to be substantial.This large uncertainty arises from the inherent insensitivity of the data to the polarized dipole amplitudes G 2 and G.To address this challenge, we discussed several potential future improvements, among which investigating jet or hadron production in longitudinally polarized proton-proton collisions emerges as a promising medium-term solution.However, more theoretical developments are desirable in the short term, where one must identify the observables which can be expressed in terms of the polarized dipole amplitudes G 2 and G.
Another issue which needs to be clarified in the medium term is the impact of the axial anomaly on the g 1 structure function and hPDFs at small x.The role of the axial anomaly in the polarized structure functions, originally pointed out in Refs.[10,172,173], has been recently revisited in Refs.[174][175][176][177].The effect appears to be distinct from the DLA of BER and KPS-CTT evolution.Developing the corresponding phenomenology is left for future work.
Based on current experimental data, we find that there could be significant negative net spin, as well as non-negligible contributions to the triplet and octet axial-vector charges, coming from small-x partons.However, there are large uncertainties in our estimates, including unaccounted-for systematics in matching onto large-x DGLAP-based fits, which will be important to implement in future work.Nevertheless, in such a scenario (negative net parton spin), significant OAM would be needed to satisfy the (Jaffe-Manohar) spin sum rule.The inclusion of EIC data in the long term would greatly enhance our understanding of hPDFs, as our impact study showed, and enable more precise statements about the distribution of (spin and orbital) angular momentum within the proton.

(A3)
The expansions for other (neighbor) dipole amplitudes are similar.Note that the transverse sizes in neighbor dipoles are always ordered such that x 32 < x 21 < x 10 , which implies that s 32 > s 21 > s 10 .Neglecting order-∆ 2 terms for small step sizes ∆ ≪ 1, Eqs. (A1) can be written as Q q (s 10 , η + ∆) = Q q (s 10 , η) + Q (0) q (s 10 , η + ∆) − Q (0) q (s 10 , η) the amplitudes, a sign change in the s 10 contributions due to the positive starting point and negative growth, and the asymptotic behavior at small x.The last property is also useful for checking the implementation of our hPDF calculation, since the dipole amplitudes and hPDFs should have the same asymptotics.
We show in Fig. 19 high-resolution (small step-size) numerical solutions of the polarized dipole amplitudes, as functions of η for a fixed s 10 , compared to their analytic counterparts.The general shape and growth of the flavor nonsinglet amplitudes (see the left panels in Fig. 19) shows a good agreement between the numerical and analytic solutions with a reasonably small step-size of ∆η = ∆s 10 = ∆ = 0.03.One can see that the analytic solution grows in magnitude slightly faster than the numeric solution.The logarithm of the absolute value of the dipole amplitudes, plotted in the right panels of Fig. 19, reveals further quantitative agreement, where we see that the numerical intercept α h converges to within 1.4% of the analytic solution.The logarithmic scale also allows us to compare the two solutions' large-x (low η) behaviors using the location of the sign change (the cusp) in the b N S contribution (the middle right panel of Fig. 19).The lower the fixed s 10 value, the lower the sign change.We see in Fig. 19 that when s 10 = const = 0.3 the sign changes coincide just above η = 2.5, implying that our numerical solution is equally valid as x → x 0 .Furthermore, we can delay the sign-change by increasing s 10 for these plots, and that will allow us to to determine the necessary resolution for retaining agreement as x becomes small.This test is given by the left-hand panel of Fig. 20, which informs us that a resolution of ∆η = ∆s 10 = ∆ < 0.06 will retain analytic agreement at the dipole amplitude level.We routinely use ∆ ≤ 0.025 for our numerics and global analysis.
The polarized dipole amplitude-level agreement gives us confidence to compare how each solution impacts our observables ∆q − .We employ the plots on the right-hand panel of Fig. 20 to extract the intercept of the ln |∆u − | basis functions and confirm that the hPDFs asymptotics given by the analytic and numerical dipole amplitudes match within 1%, and are consistent with the intercept that was computed at the dipole amplitude level.This completes the cross-check of our numerical solution for the flavor nonsinglet evolution equations.

Appendix C: Convergence testing of numerical solutions
The discretization defined in Appendix A is very useful for solving complicated integral equations which are very difficult if not impossible to solve analytically.The numerical solution is rather straightforward to derive, but it has the same faults as any discrete function, namely the fact that the accuracy of a numerical solution is dependent on the resolution, i.e., the step size.In our case, we have two different variables to work with (η, s 10 ), which results in a two-dimensional grid (G[i, j]) for our numerical solution to compute.To simplify the discretization, we defined the step sizes for η and s 10 to be the same, ∆η = ∆s 10 ≡ ∆.The requirement we impose on our numerical solution to confirm its validity is that as the step-size decreases, the computed values should converge to a single output.
We have tested each of our flavor singlet basis functions (Fig. 1) as well as the flavor nonsinglet basis functions (not shown).However, the results can be summarized by their subsequent implementation in calculating the hPDFs ∆q + (x) and ∆q − (x).The left-hand panel of Fig. 21 shows x∆u + (x) for a "test state" of initial conditions.We define a test state simply as any replica that has been confirmed to fit data with χ 2 Npts ≈ 1.This hPDF was plotted multiple times for varying step-sizes, and it is clear that as the step-size decreases the solutions converge to a single output.
The same convergence test was conducted on x∆q − (x) and is displayed in the right panel of Fig. 21.In this case there is also an analytic solution, as discussed in Appendix B. We find not only a convergence of the numerical solution to a single output as ∆ becomes smaller, but also that the converged output is exactly that of the analytic solution.We note here that Fig. 21 is a demonstration of the convergence.The results discussed in Sec.III were computed using much higher resolutions, ∆ ≈ 0.02.

behavior of g p 1 20 3 . 1 28
singlet evolution at small x 4 B. Flavor nonsinglet evolution at small x 7 C. Numerical implementation of the flavor singlet and nonsinglet evolution 8 D. SIDIS cross section at small x 11 Origins of asymptotic behavior 21 C. Extracted helicity PDFs and calculation of net parton spin and axial-vector charges at small x 24 D. Impact of EIC data on g p

FIG. 1 .
FIG. 1.The u-quark hPDF, x∆u + (x), constructed solely out of each basis function in the range x ∈ [10 −5 , 1].The legend in each panel shows which basis function was used for which curve.For example, the blue curve in the top panel corresponds to x∆u + (x) constructed from the initial conditions Q

FIG. 3 .
FIG.3.Comparison of the experimental data and the fit based on our small-x theory for the double-spin asymmetries A1 and A ∥ in polarized DIS on a proton (red), deuteron (blue) and 3 He (green) target.

FIG. 7 .
FIG. 7. (Left) Histogram that counts the number of replicas with a smallest-x ambiguity at a given value of x. (Right) The running sum of the ambiguity histogram, telling us what percentage of replicas have an ambiguity below a given value of x.

FIG. 9 .
FIG. 9. (Left) Histograms utilizing Eq.(38) showing that as x decreases, the intercept α h (x) becomes more constrained as a consequence of the small-x evolution equations.(Right) Keeping information on the sign dependence by using Eq.(39) produces bimodal peaks at ±α h (x).At large x there is no asymptotic behavior, and for smaller values of x two refined peaks emerge.

ds
FIG. 20. (Left) A plot of (the logarithm of) the s10 contribution to G NS u (parameterized by b NS u ) as a function of η.Each color represents a different fixed value of s10.The location of the sign-change in the amplitude, indicated by the cusp, appears to vary with s10.Smaller step sizes lead to convergence of the sign-change between the analytic and numeric solutions, and ∆η = ∆s10 = ∆ < 0.06 retains small-x agreement.(Right) A plot of (the logarithm of) each ∆u − basis functions (parameterized by a NS u , b NS u , and c NS u ) as a function of log (x).Each plot depicts the asymptotic agreement between the numeric and analytic solutions, as well as a measure of the intercept α h .

TABLE II .
Summary of the polarized SIDIS data on A h 1 included in the fit, along with the χ 2 /Npts for each data set.