Helicity at Small $x$: Oscillations Generated by Bringing Back the Quarks

We construct a numerical solution of the recently-derived large-$N_c \&N_f$ small-$x$ helicity evolution equations with the aim to establish the small-$x$ asymptotics of the quark helicity distribution beyond the large-$N_c$ limit explored previously in the same framework. (Here $N_c$ and $N_f$ are the numbers of quark colors and flavors.) While the large-$N_c$ helicity evolution involves gluons only, the large-$N_c \&N_f$ evolution includes contributions from quarks as well. We find that adding quarks to the evolution makes quark helicity distribution oscillate as a function of $x$. Our numerical results in the large-$N_c \&N_f$ limit lead to the $x$-dependence of the flavor-singlet quark helicity distribution which is well-approximated by \begin{align} \Delta \Sigma (x, Q^2)\bigg|_{\mbox{large-}N_c \&N_f} \sim \left( \frac{1}{x} \right)^{\alpha_h^q} \, \cos \left[ \omega_q \, \ln \left( \frac{1}{x} \right) + \varphi_q \right]. \end{align} The power $\alpha_h^q$ exhibits a weak $N_f$-dependence, and, for all $N_f$ values considered, remains very close to $\alpha_h^q (N_f=0) = (4/\sqrt{3}) \sqrt{\alpha_s N_c/(2 \pi)}$ obtained earlier in the large-$N_c$ limit. The novel oscillation frequency $\omega_q$ and phase shift $\varphi_q$ depend more strongly on the number of flavors $N_f$ (with $\omega_q =0$ in the pure-glue large-$N_c$ limit). The typical period of oscillations for $\Delta \Sigma$ is rather long, spanning many units of rapidity. We speculate whether the oscillations we find are related to the sign variation with $x$ seen in the strange quark helicity distribution extracted from the data.


I. INTRODUCTION
Understanding the partonic (quark and gluon) structure of the proton is an essential part of our understanding of Quantum Chromodynamics (QCD). One of the big open questions in the studies of proton structure is the proton spin puzzle. Proton is a composite spin-1/2 particle made of quarks and gluons. Consequently, the spin of the proton should be a sum of all the spins and orbital angular momenta (OAM) carried by the quarks and gluons comprising the proton. This statement is formalized in terms of helicity sum rules, one due to Jaffe and Manohar [8] and another one due to Ji [9]. The former reads S q + L q + S G + L G = 1 2 (2) where S q and S G are the total spin carried by the quarks and gluons, respectively, whereas L q and L G are the OAM carried by the quarks and gluons in the proton. One can write S q and S G , which are functions of the momentum scale Q 2 , as integrals of helicity distribution functions over the Bjorken x variable, where is the flavor-singlet quark helicity distribution function. Here ∆f and ∆f denote the helicity distribution functions of quarks and antiquarks, respectively. In Eq. (3), ∆G is the gluon helicity distribution. The quark and gluon OAM, L q and L G , can also be written as x-integrals of the OAM distributions [10][11][12][13][14]. Note that the latter cannot be expressed as expectation values of twist-2 operators in the proton state, unlike the helicity parton distribution functions (hPDFs). The proton spin puzzle started with the discrepancy between the prediction of the (perhaps, naïve) constituent quark model stating that all the proton spin is carried by the constituent quarks (such that S q = 1/2) and the experimental measurements by the European Muon Collaboration (EMC) [15,16] that reported a much lower number for S q . More recently, the experiments at the Relativistic Heavy Ion Collider (RHIC) measured a non-zero gluon spin S G [17,18]. The current experimental values for quark and gluon spin are S q (Q 2 = 10 GeV 2 ) ≈ 0.15 ÷ 0.20, integrated over 0.001 < x < 1, and S G (Q 2 = 10 GeV 2 ) ≈ 0.13 ÷ 0.26, integrated over 0.05 < x < 1 (see [19][20][21][22][23] for recent reviews). One can see that these reported S q and S G do not add up to 1/2 and the proton spin puzzle is still open. (Note also that none of the terms in Eq. (2) is positive-definite.) We conclude that the missing part of the proton's spin must come either from the smaller-x regions in the integrals of (3) for the spin terms S q and S G than the x-ranges reported above or from the OAM terms L q and L G .
No experiment, present or future, can measure helicity PDFs all the way down to x = 0, since this would require infinite energy. In addition, at very small (but non-zero) values of x higher-twist corrections become comparable to the leading-order contribution for many observables, making it hard (if not impossible) to extract twist-2 hPDFs from the data. (The way such corrections may come in is described using the parton saturation framework, see [24][25][26][27][28][29] for reviews.) Therefore, to assess the amount of parton spin and OAM coming from the small-x region, one has to develop a robust theoretical formalism, which, if able to describe the existing experimental data and predict the values of the future longitudinal spin measurements, can be extrapolated down to x = 0 with a (hopefully) good level of confidence.
Helicity evolution equations were derived in [1,40,42] using the DLA. Similar to the BK evolution, the equations close in the large-N c limit, with N c the number of quark colors. In addition, and different from the unpolarized case, the helicity evolution equations [1,40,42] also close in the large-N c &N f limit, where N f is the number of quark flavors. Owing to the complexity of the large-N c &N f equations [1,42], only the large-N c equations have been solved to date, resulting in the following power-law small-x asymptotics for the quark and gluon helicity distributions in DLA [2,3,41] ∆Σ(x, Q 2 ) large-Nc (When writing equations like (5) we suppress the potentially non-trivial Q 2 dependence of the involved observables, along with the possible sub-dominant x-dependence in the prefactor, and concentrate on the leading x dependence of the quantities.) The small-x asymptotics of quark and gluon OAM distributions at large N c were found in [38] using the same formalism. The goal of this work is to numerically solve the large-N c &N f helicity evolution equations derived in [1,42] and use the solution to assess possible corrections to the large-N c result (5) for the asymptotic behavior of ∆Σ(x, Q 2 ) at small x. The results can be (and are) combined with the values of ∆Σ extracted from the experiment [4][5][6][7] to extrapolate ∆Σ to smaller values of x and estimate the total quark contribution to the proton spin S q . In addition, our results can be employed to find the small-x asymptotics of ∆G and the OAM distributions for quarks and gluons in the large-N c &N f limit, that is, beyond the large-N c expressions for ∆G in Eq. (5) and for the OAM small-x asymptotics derived in [38].
The paper is structured as follows. In Section II, we rewrite and simplify the large-N c &N f helicity evolution equations from [1,42]. The equations are written in terms of the quark and gluon polarized dipole amplitudes, which are defined below as correlators of the polarized and regular light-cone Wilson lines. Our numerical algorithm is presented in Sec. III, where we rewrite the equations in terms of the variables convenient for our numerical approach, discretize the equations and cast them in the form suitable for implementing our algorithm. The resulting numerical solutions of the large-N c &N f helicity evolution equations for N f = 2, 3, 6 are presented in Sec. IV, see Figs. 2 and 3 there. The surprising new result, best visible in Fig. 3, is that the quark and gluon polarized dipole amplitudes oscillate as functions of energy (or of ln(1/x)), changing the sign back and forth. These oscillations of dipole amplitudes, in turn, result in oscillations in ∆Σ(x, Q 2 ) as a function of ln(1/x), as demonstrated in Eq. (32) or, equivalently, Eq. (1) above, and depicted in Fig. 4. Such oscillating behavior was absent in the large-N c (and N f = 0) analyses carried out in the previous works [2,3,38,41], which only saw a power-of-1/x growth (5) of helicity distributions at small x. The oscillations of ∆Σ(x, Q 2 ) with ln(1/x) are the main qualitative result of this paper.
The parameters α q h , ω q , and ϕ q from Eq. (1), describing ∆Σ(x, Q 2 ) at small x and at large-N c &N f , are extracted from our numerical solution in Appendix B and are summarized in Eq. (33). In Sec. IV we also present a fit (29) for the N f -dependence of the oscillation frequency ω q which increases with N f , meaning that that the oscillation period gets shorter and, therefore, the oscillations become more pronounced, in the regime where more quark flavors become dynamically relevant. In Sec. V we follow the method from [2] to construct a preliminary estimate of the impact of the asymptotic form for ∆Σ(x, Q 2 ) given by Eq. (1) on the amount of spin carried by the quarks in the proton. The results are given in Figs. 5 and 6, similar to [2] indicating a potential for a large correction. We conclude in Sec. VI and speculate whether the oscillations we find in our solution are related to the oscillation of the strange quark hPDFs extracted from the experimental data in [4][5][6][7].

II. HELICITY EVOLUTION EQUATIONS AT LARGE-Nc&N f
At small x, in the DLA, the flavor-singlet quark helicity PDF can be written as [1,40,42] where s ≈ Q 2 /x is the center-of-mass energy squared, Λ is a transverse momentum scale characterizing the target before small-x evolution, and the impact parameter-integrated polarized quark dipole amplitude is defined by where the polarized quark dipole amplitude is The dipole consists of a quark and an anti-quark located at transverse positions x 1 and x 0 depending on the term in Eq. (9). In our notation the transverse vectors are denoted by In addition, x ij denotes the magnitude of the vector x ij . The light-cone variables are defined by x ± = (x 0 ± x 3 )/ √ 2. Our longitudinally polarized proton moves in the x + direction. The quantity z is the fraction of some original (x −direction-moving) probe's momentum carried by the softest (anti-)quark in the dipole. The angle brackets in Eq. (9) denote the averaging in the proton wave function in the small-x/saturation sense [24][25][26][27][28][29], albeit now taking into account that the proton is longitudinally polarized.
The polarized dipole amplitude (9) is defined in terms of the light-cone fundamental Wilson lines and the so-called polarized Wilson lines V pol [1,[40][41][42], consisting of regular light-cone Wilson lines along with the sub-eikonal helicity-dependent local operator insertion(s), (Note that the definition (11) is different by an extra factor of s in overall normalization as compared to that used in earlier works [41,42].) We also use an abbreviated notation where As usual, T in Eq. (9) denotes time ordering of the operators. In equations (10) and (11), A + is the eikonal gluon field of the shock wave (often referred to as the target, or the target proton), F 12 is the helicity-dependent sub-eikonal fieldstrength tensor component of the target gluon field (both in the fundamental representation), while ψ andψ are the (sub-eikonal) quark and anti-quark fields of the shock wave. In addition, p + 1 is the large light-cone momentum of the parton in the shock wave generating those sub-eikonal fields. In Eq. (11) we have also employed the adjoint light-cone Wilson line The matrices t a and t b are the fundamental generators of SU(N c ) with the indices a, b = 1, . . . , N 2 c − 1. The large-N c &N f helicity evolution equations we are about to study mix the polarized quark dipole amplitude (9) with the polarized gluon dipole amplitude [41,42] defined in terms of the adjoint polarized Wilson lines [42] (U pol FIG. 1: Diagrammatic representation of the quark (Q) and gluon (G) polarized dipole amplitudes. The shaded rectangle is the shock wave, while the square represents insertions of sub-eikonal operator(s).
(Note again an additional s in the normalization of Eq. (14) as compared to [42]; in addition, G 10 (z) in Eq. (13) was labeled G adj in [42].) The impact parameter-integrated polarized gluon dipole amplitude is defined similar to Eq. (8) The polarized quark and gluon dipole amplitudes are depicted diagrammatically in Fig. 1. The polarized proton, along with all the partons to be produced by the small-x evolution, are depicted by the rectangular shock wave. The square on one of the lines depicts insertions of sub-eikonal operators present in Eqs. (11) and (14).
As follows from Eq. (7), we only need to find the amplitude Q(x 2 10 , z) in order to construct the quark helicity PDF. However, the evolution equations (16) mix all four involved amplitudes, and have to be solved for all the amplitudes in order to obtain Q(x 2 10 , z). In this paper we are chiefly interested in the asymptotic behaviors of the dipole amplitudes at high energies. In [2] it was shown that the high-energy asymptotics of helicity amplitudes at large-N c is largely independent of the initial conditions/inhomogeneous terms in the evolution, which only affect the overall prefactor in the asymptotic expression. (In addition, it is well-known that the small-x asymptotics of the Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution [56,57] depends on the initial conditions only in its prefactor as well.) Inspired by these examples of independence of the initial conditions, we assume it to be the case in the large-N c &N f approximation as well, and, for convenience, we take the inhomogeneous terms in Eqs. (16) to be instead of using the proper Born-level initial condition [2].
In the regime of our interest, 1 z s allowing us to rewrite the evolution equations (16) as Our aim is to solve this system of integral equations numerically. Its solution for Q(x 2 10 , z), via Eq. (7), would give us the quark helicity distribution at small Bjorken x and in the regime where both the gluon and quark contributions to evolution are relevant (due to the large-N c &N f limit employed in deriving Eqs. (20)).

III. NUMERICAL SOLUTION: DISCRETIZATION AND ALGORITHM
To numerically evaluate the integrals in Eqs. (20), we first make the following changes of variables where z (n) can be z, z or z which relates to η, η or η respectively. In addition, kl = 10, 21 or 32. In terms of these new variables, Eqs. (20) can be re-written as Note that the integrals in equations (22) include the regions with s 21 < 0, s 32 < 0. This means that Λ is not an infrared (IR) cutoff, but a perturbative transverse momentum scale characterizing the proton at the start of our evolution. Next, we discretize the resulting equations (22). In particular, let i (n) and j (n) be the discretized version of s kl and η (n) , with step sizes ∆s and ∆η, respectively, such that η (n) = j (n) ∆η ≡ η j (n) along with s 10 = i ∆s ≡ s i , s 21 = i ∆s ≡ s i , and s 32 = i ∆s ≡ s i . The discretized amplitudes are defined as (22) are self-contained over the following region in (s 10 , η)- where η max is some arbitrary positive upper value of the η-range. This is the region where we will solve them numerically. Discretized Eqs. (22) are In principle, Eqs. (23) can already be solved numerically, similar to [2,41]. Instead we will simplify these equations which will allow us to implement a faster algorithm in the numerical solution. We replace j by j − 1 in each of the four equations in (23). For instance, in equation (23a) involving G ij , this replacement gives Comparing Eq. (23a) to Eq. (24), we obtain Since zs ≥ 1/x 2 10 , we have η ≥ s 10 , which leads to j ≥ i. For i = j we have Q j(j−1) = 1 in the last term of Eq. (25): this is determined by the initial conditions. Including this additional contribution from a single point in the i, j grid does not significantly affect the numerical solution.
Writing each equation in (23) in this recursive form allows for a numerical evaluation with one fewer layer of loops, resulting in much shorter computation time for smaller step sizes, ∆η and ∆s, and for a larger η-range, defined by η ∈ [0, η max ]. In order to write down the recursive equations similar to Eq (25) for Q ij , Γ ikj and Γ ikj , notice that where we have used the fact that i < k < j, a consequence of the regime in which Γ (and Γ) are defined and employed, 1/z s x 2 21 x 2 10 , or, equivalently, η s 21 s 10 . Equations (23) and (26) imply that (again for i < k < j) We solve Eqs. (25) and (27) numerically in steps along the η-axis. To obtain the values for 0 ≤ η ≤ η max , we start from η = ∆η, which is equivalent to j = 1, at which we use Eqs. (25) and (27) to determine G i1 and Q i1 for 1 − ηmax ∆η ≤ i ≤ 1, together with Γ ik1 and Γ ik1 for 1 − ηmax ∆η ≤ i ≤ k ≤ 1, assuming that G i0 = Q i0 = Γ ik0 = Γ ik0 = 1 are determined by the inhomogeneous term in Eqs. (23). Afterward, we repeat the same steps by applying Eqs. (25) and (27) for j = 2, then j = 3, and so on, until j max = ηmax ∆η . At each j, we compute the dipole amplitudes G ij and Q ij for i in the range j − ηmax ∆η ≤ i ≤ j and also find the amplitudes Γ ikj and Γ ikj for all pairs of i and k satisfying j − ηmax ∆η ≤ i ≤ k ≤ j. This process determines the values of Q(s 10 , η) in the η ∈ [0, η max ], η − η max ≤ s 10 ≤ η region, which we will use below to determine the high energy asymptotics of the polarized quark dipole amplitude. We always take ∆s = ∆η for simplicity.

IV. NUMERICAL SOLUTION: RESULTS
In this Section we present the results of our numerical solution of the large-N c &N f evolution equations (22) for the quark and gluon polarized dipole amplitudes Q(s 10 , η) and G(s 10 , η), along with the corresponding quark helicity PDF ∆Σ(x, Q 2 ). The plots of ln |G(s 10 , η)| and ln |Q(s 10 , η)| in the (s 10 , η)-plane are shown in Fig. 2 for N f = 2, 3, 6. Throughout this paper, we take the number of colors to be N c = 3.
The plots in Fig. 2 demonstrate an approximately linear rise of ln |G(s 10 , η)| and ln |Q(s 10 , η)| with η, similar to the large-N c case studied previously in [2,3]. The only difference is that now, in Fig. 2, the rise of those functions is not monotonic, and appears to be periodically interrupted by lines of sharp local minima.
To illustrate the origin of this non-monotonicity, we plot sgn[G(0, η)] ln |G(0, η)| and sgn[Q(0, η)] ln |Q(0, η)| as functions of η in Fig. 3 for N f = 3. From these plots we see that G(0, η) and Q(0, η) oscillate with η. The oscillations explain the non-monotonic behavior we saw in Fig. 2. These oscillations are the main qualitative difference of the small-x asymptotics for the quark helicity distribution in the large-N c &N f limit, as compared to the large-N c case. It appears that introducing the quarks back into the helicity evolution equations generates oscillations. While the absolute values of Q(s 10 , η) and G(s 10 , η) still grow exponentially with η, both dipole amplitudes also oscillate. Moreover, from Fig. 2 one can see that the oscillation period appears to get smaller (and the oscillation frequency appears to get larger) with increasing N f , which is also consistent with the fact that these oscillations are absent in the gluon-only large-N c equations solved in [2,3].
While an analytic solution of Eqs. (22) is beyond the scope of this work, we can try to find analytic formulas approximating our numerical results in Fig. 3, at least in the large-η asymptotics. Combining the oscillations with the exponential growth of the maxima of |Q(0, η)| and |G(0, η)| with η we propose the following asymptotic forms for the polarized dipole amplitudes: The oscillation frequencies are denoted by ω G and ω Q , while the initial phases are denoted by ϕ G and ϕ Q . As one can already see from Fig. 2, both the frequencies and the initial phases depend on N f . This is confirmed by the detailed analysis of our numerical solution in Appendix A. Furthermore, the amplitudes of oscillations in Q(0, η) and G(0, η) grow exponentially with η, while the exponents α G and α Q (the intercepts in Regge terminology) also appear to depend on N f . The results of the analysis carried out in Appendix A, where we fit our numerical solution with the ansatz (28) and extract the corresponding frequencies, phases, and intercepts, are summarized in Table I. From Table I we see that, for each N f , the two frequencies are equal within the numerical accuracy, ω G = ω Q . The frequencies increase with N f : while the exact analytic dependence of ω G = ω Q on N f is yet to be determined, we construct a Padé approximant to write    The resulting intercept α, frequency ω, and initial phase ϕ, for the gluon and quark polarized dipole amplitudes G(0, η) and Q(0, η) at N f = 2, 3, 6 and N c = 3.
The intercepts, α G and α Q , of the polarized dipole amplitudes' exponential growth given in Table I, are also equal to each other for each N f with the precision of our numerical solution, α Q = α G . For all N f studied they remain close to α q h (N f = 0) = 4 √ 3 ≈ 2.309, which is the intercept for G(0, η) in the large-N c pure-glue limit with N f = 0 (in units of α s N c /(2π)), derived analytically in [3]. It appears, though, that α Q = α G is decreasing slowly with N f .
The initial phase ϕ in Table I is always between 0 and π 2 for G(0, η) and between − π 2 and 0 for Q(0, η). The values of ϕ Q and ϕ G do not display a clear relation with each other. Their functional dependence on N f is non-monotonic, and its form is also not obvious from Table I. In fact, the initial phase depends greatly on the choice of the initial condition (the inhomogeneous term) in Eqs. (16). For instance, if one performs a similar computation using the Born-level-inspired dipole amplitudes as initial conditions (cf. [2]), while still taking the inhomogeneous terms for G and Q to be equal for simplicity,  Table I, which were obtained for the initial condition (18). The difference between ϕ G = −1.146 and ϕ Q = 1.530 and the phases in Table I is too large to be attributed to discretization errors. We thus conclude that the phases ϕ Q and ϕ G are indeed very much dependent on the initial conditions to the evolution, and their values listed in Table I are To determine the small-x asymptotics of ∆Σ we need to evaluate the integral in Eq. (31). Note that for Q ≥ Λ the integration region in Eq. (31) lies within the area η ∈ [0, η max ], η − η max ≤ s 10 ≤ η where we found the numerical solution for Q(s 10 , η) if we choose η max = αs Nc 2π ln Q 2 xΛ 2 . Since Q(s 10 , η) is known numerically, we perform a numerical integration to obtain the values of ∆Σ(x, Q 2 ) as a function of Bjorken x at fixed Q 2 = 10 GeV 2 for the case of 3 quark flavors, N f = 3, while choosing, for simplicity, Λ 2 = Q 2 , such that η max = αs Nc 2π ln 1 x . The resulting ∆Σ(x, Q 2 ), similar to the polarized dipole amplitudes Q and G, is an oscillating function of ln(1/x), with the oscillation amplitude growing exponentially with ln(1/x). In Fig. 4 we plot sgn[∆Σ(x, Q 2 )] ln |∆Σ(x, Q 2 )| as a function of x, demonstrating the oscillations explicitly. Inspired by the success of the ansatz (28) for the dipole amplitude, we propose the following ansatz for the small-x asymptotics of ∆Σ(x, Q 2 ), The parameters α q h , ω q and ϕ q are extracted from our numerical results in Appendix B by applying the parameter fitting process outlined in Appendix A. We use α s (10 GeV 2 ) ≈ 0.25. This gives (again, for N f = 3) 2π . This shows that the intercept and frequency of the dipole amplitude Q determine the intercept and frequency of ∆Σ: both are uniquely determined by the evolution. (Note that one cannot simply substitute Eq. (28b) into Eq. (31) to obtain this result analytically, since the former is valid only at s 10 = 0, while the latter has an integral over a range of non-zero values of s 10 .) However, as remarked previously for the dipole amplitudes, the initial phase, ϕ q , is, in fact, a by-product of our choice of initial conditions, G (0) and Q (0) . In practical phenomenological applications, G (0) and Q (0) , and, hence, the initial phase ϕ q , have to be determined from the data.

V. QUARK HELICITY: AN ESTIMATE
In this Section we follow the strategy employed in [2] to estimate the possible impact of our new functional form for the small-x asymptotics of ∆Σ(x, Q 2 ) (32) on the amount of the proton spin carried by the small-x quarks. Our analysis below should be understood as a rough estimate, with the anticipation of a more detailed phenomenology to be done in the future.
As the asymptotic form (32) only holds for small x, we use that expression to extrapolate the quark helicity distribution from the DSSV14 result in [6], starting at a particular connecting point, x 0 1, into the small-x region, x min ≤ x ≤ x 0 . In particular, we write the small-x quark helicity distribution for x < x 0 as with α q h and ω q given in Eq. (33). For x > x 0 we take the quark hPDF to be given by the DSSV14 parameterization [6], ∆Σ DSSV (x, Q 2 ). We, therefore, need to match our ∆Σ(x, Q 2 = 10 GeV 2 ) from Eq.
Unlike [2], where the ansatz for ∆Σ at small x contained only the power of 1/x, we now have the function in Eq. (34) with two unknown parameters, the overall normalization factor K and the phase ϕ q . We assume that a more complete phenomenological approach would be able to uniquely determine ϕ q from the large-x (x > x 0 ) data. We further assume that the value of ϕ q obtained from a more complete approach would generate a smooth matching of ∆Σ at x = x 0 , ensuring continuity of both ∆Σ(x, Q 2 ) and its derivative ∂∆Σ(x, Q 2 )/∂x at x = x 0 . This assumption is not very reliable for our purposes, since the functional form (34) is asymptotic, and is, therefore, valid only for x x 0 : using it to ensure the continuity of ∆Σ(x, Q 2 ) and ∂∆Σ(x, Q 2 )/∂x at x = x 0 is somewhat questionable. Strictly-speaking we have to admit that ϕ q is an almost arbitrary parameter, whose value appears to be very important for assessing the amount of quark helicity at small x. While the more detailed phenomenology should better constraint the allowed ranges of K and ϕ q , we will fix them here by simply requiring the continuity of ∆Σ(x, Q 2 ) and ∂∆Σ(x, Q 2 )/∂x at x = x 0 between our asymptotics (34) and ∆Σ DSSV (x, Q 2 ) from [6].
We perform this computation for x 0 = 0.01 and x 0 = 0.001, doing a separate matching at each x 0 . The resulting x∆Σ(x, Q 2 = 10 GeV 2 ) is plotted versus x in Fig. 5, which depicts the original DSSV14 [6] curve, extrapolated by a power-law in x to very small x, along with the two curves resulting from matching our asymptotics (34)  To see more clearly the implication of this result on the quark helicity, S q (Q 2 = 10 GeV 2 ), we consider the integral similar to Eq. (3) but with the lower limit at some small finite x min , such that 0 < x min < x 0 (cf. [2]) From Eqs. (3) and (35) we see that S q (Q 2 ) = (1/2) ∆Σ [xmin] (Q 2 ) in the x min → 0 limit. In Fig. 6, we plot ∆Σ [xmin] (Q 2 = 10 GeV 2 ) versus x min for the DSSV14 parameterization [6] along with the two curves resulting from matching our (34) to DSSV14 at x 0 = 0.01 and x 0 = 0.001. In addition, for comparison, we show two extra curves in Fig. 6 (dashed lines, labeled KPS16), resulting from using the power-law ansatz for small-x ∆Σ obtained from the large-N c pure-glue evolution [2,3], and adjusting its overall normalization to ensure the continuity of ∆Σ(x, Q 2 ) with DSSV14 at x 0 = 0.01 and x 0 = 0.001, as was done in [2]. The thin (red) line is DSSV14 [6], while the thick (green) and medium-thick (blue) solid lines show the results of matching our asymptotics (34) to DSSV14 at x 0 = 0.01 and x 0 = 0.001, respectively. The two dashed lines represent the small-x extrapolation done in [2] using a power-law expression for ∆Σ(x, Q 2 ) at small x resulting from the large-N c evolution, also matched onto DSSV14 at x 0 = 0.01 (thick dashed, green) and x 0 = 0.001 (medium-thick dashed, blue). Figure 6 illustrates the potential of small-x evolution to significantly affect the amount of the quark spin content in the proton, which was already observed in [2]. We conclude from Fig. 6 that the small-x extrapolation of ∆Σ(x, Q 2 ) at Q 2 = 10 GeV 2 varies significantly with the matching point, x 0 . For the larger x 0 -value, x 0 = 0.01, we also see a strong variation between using large-N c (KPS16) and large-N c &N f small-x evolution for helicity. It is also interesting to see in Fig. 6 that for x 0 = 0.001 both the KPS16 curve and the curve based on Eq. (34) are rather close to the DSSV14 curve and to each other. Again, let us stress that our use of the asymptotic expression (34), valid for x x 0 , to perform the matching of both ∆Σ(x, Q 2 ) and ∂∆Σ(x, Q 2 )/∂x at x = x 0 , which is probably outside of its region of applicability, gives us only a rough estimate of quark hPDF at small x. We expect that future more detailed studies would place this matching under a more solid theoretical control.

VI. CONCLUSIONS AND DISCUSSION
In this work, we have numerically computed the asymptotic high-energy behavior of the polarized quark dipole amplitude resulting from the double-logarithmic small-x helicity evolution [1,42] in the limit of large number of quark colors and flavors. The obtained amplitude Q is plotted in Figs. 2 and 3 as a function of η and s 10 defined in Eq. (21). The amplitude Q in the large-N c &N f limit displays an oscillatory pattern as a function of the center-of-mass energy and of the dipole transverse size, on top of the exponential growth seen before in the large-N c limit with N f = 0 [2,3], when all the evolution was gluon-driven. We conclude that these oscillations appear after including quarks back into the small-x helicity evolution of [1,42].
Our numerical results for Q are well-approximated by Eq. (28b). The frequency ω Q and initial phase ϕ Q of the oscillations depend on the number of flavors, with this dependence summarized in Table I and Eq. (29). The intercept α Q , also given in Table I This is the main result of this work. It is illustrated in Fig. 4. The intercept, α q h , and the frequency, ω q , given by Eq. (33), are within the margin of error from the corresponding parameters, α Q and ω Q , for the dipole amplitude, Q. The initial phase, ϕ q , on the other hand, differs greatly from ϕ Q , due to the relation (31) between the two quantities. In general, the phase ϕ Q , and, therefore, the phase ϕ q , exhibit a very strong dependence on the initial conditions (inhomogeneous terms) for our large-N c &N f evolution. This is in contrast to α q h ≈ α Q αsNc 2π and ω q ≈ ω Q αsNc 2π which are independent of the initial conditions and are universal properties of the evolution. This ambiguity in the initial phase unfortunately trickles down to the ambiguity in our prediction for the quark helicity inside a proton. Fixing the phase should be done by determining the initial conditions for the evolution. Resolving this issue would require a more detailed phenomenological work, which is beyond the scope of this paper. Note again that the oscillation frequency ω q vanishes in the N f = 0 limit, such that the oscillation is a property of having quarks in the evolution. The N f = 0, N c → ∞ limit itself warrants a little further discussion. In [1][2][3][40][41][42] the large-N c limit was understood as N f = 0 gluons-only evolution. In this regime the polarized Wilson lines are given only by the first (∼ F 12 ) terms in Eqs. (11) and (14), such that Q(x 2 10 , z) ≈ G(x 2 10 , z)/4 [42]. The small-x asymptotics of the quark and gluon amplitudes Q and G were, therefore, given by the same power of 1/x. However, strictly-speaking one needs to clarify what is implied by the quark dipole amplitude at N f = 0. The quark dipole amplitude obtained as a part of the gluon dipole amplitude in the large-N c limit, as employed in [1][2][3][40][41][42] and given by Q(x 2 10 , z) ≈ G(x 2 10 , z)/4, is not exactly the right object to determine the quark helicity PDF ∆Σ in Eq. (7). To properly define ∆Σ as the quark helicity distribution in the large-N c limit, one needs to have a non-zero N f , at least N f = 1. The N f = 1, N c → ∞, α s N c = const limit of equations (16) is obtained by dropping all the N f -terms, that is, formally by putting N f = 0 in the equations. We report here that the solution of the resulting equations performed as part of this work leads to the intercept α Q = 2.39. If we compare this intercept to those listed in Table I we see that this result appears to support our earlier conclusion about mild N f -dependence of the intercept. Note that the N f = 1, N c → ∞, α s N c = const limit is further complicated by the fact that in Q(x 2 10 , z) the Born-level interaction with the quark target, which should be used for calculating Q (0) (x 2 10 , z), consists of a sum of two terms with different N c -scaling [40] (see Eq. (30) above): the t-channel quarks exchange is N c -enhanced compared to the t-channel gluon exchange contribution. Careful implementation of this observation shows that imposing the N f = 1, N c → ∞, α s N c = const limit on equations (16) by putting N f = 0 and obtaining the intercept α Q = 2.39, is correct only for a pure-glue shock wave, that is, for probing a polarized glueball state instead of the proton. For the scattering on the actual proton, or on a single polarized quark, the limit is more subtle, and includes one iteration of the N fterms in the solution of Eqs. (16). This regime has not been explored in this work due to the higher phenomenological relevance of the large-N c &N f approximation considered here.
Finally, let us comment on the potential phenomenological implications of our main qualitative result: at small-x the flavor-singlet quark helicity distribution ∆Σ(x, Q 2 ) oscillates in ln(1/x). Similar oscillation has been found in the strange quark helicity distribution ∆s extracted from the experimental data by the PDF collaborations [4][5][6][7]. This oscillation is the driving force behind the sign change of ∆Σ DSSV (x, Q 2 ) in Fig. 5 above. If the ∆s oscillation in x is confirmed by the future data extractions, it appears reasonable to ask a question whether it is related to our oscillating result (36) for ∆Σ(x, Q 2 ). While we do not separately consider individual quark flavors, the frequency of ∆Σ oscillations increases with N f , and should be more pronounced if more flavors are included in the helicity evolution. This may be related to the oscillation observed in ∆s, and not in ∆u or ∆d. On the other hand the period of these oscillations in ln(1/x) is Using ω q = (0.469) αsNc 2π from Eq. (33) with α s = 0.3 and N c = 3, we obtain T (N f = 3) ≈ 35, which is a very large number for rapidity or for ln(1/x). However, the first sign flip one encounters depends on the initial phase of the oscillation, which is hard to determine. If we start from the maximum of the cosine function in Eq. (36) at some x 0 , then the first sign flip would happen at x/x 0 ≈ e −T /4 ≈ 10 −4 , which is a more phenomenologically-reasonable number, but is still rather low to be relevant for the upcoming Electron-Ion Collider (EIC) experimental program [19,23,58]. On yet another hand, the period of oscillations found above may be significantly affected by the higherorder corrections in α s , even by the running of the coupling. In addition, we are encouraged by the similarity of the shapes between our curves and the DSSV14 line in Fig. 5, indicating that some of the physics behind the DSSV14 line might be accurately described by our small-x evolution. We, therefore, leave the final verdict on the issue of the phenomenological relevance of the ∆Σ oscillations found in this work for the future investigations.

VII. ACKNOWLEDGMENT
The authors would like to thank Mr. Daniel Adamiak for providing his code to help us construct Figs. 5 and 6, which, in turn, was based on the code provided by Prof. Daniel Pitonyak, to whom we are also grateful. We also thank Prof. Pitonyak for a discussion of numerical simulations for helicity at small x. YK would like to thank Mr. Mohammed Karaki for his work on this project in its very early stages.
This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award Number de-sc0004286.
where we also adopt a convention in which ω > 0. We use the largest-η maximum available in our simulation to extract ω using Eq. (A3). (Indeed the extracted value of ω can be cross-checked by comparing π/ω to the spacing between the positions of the local maxima along the η-axis in the numerical solution.) The phase ϕ can then be determined from the second derivative maximum condition ω η * + ϕ = π n, where η * is the numerically-extracted position of the same largest-η maximum of the second derivative (A2) and n is an integer. The value of n is adjusted so that ϕ ∈ (−π, π], where the choice between ϕ ∈ (0, π] and ϕ ∈ (−π, 0] is done by making sure that the corresponding f (η * ) given by Eq. (A1) is positive or negative (assuming that K > 0) in agreement with the numerical value of the function f (η) at η = η * . Finally, a linear regression on ln f (η) cos (ω η + ϕ) = α η + ln K (A4) allows us to extract α from the slope of this function. For the numerical values of G(0, η) and Q(0, η) found in the range η ∈ [0, η max ], we only use η ∈ [0.75η max , η max ] to extract the intercepts α G and α Q using Eq. (A4), in addition avoiding the values of η close to the nearest cosine zero, η n = 1 ω π 2 − ϕ + π n , by at least 5% of the cosine's period, (0.05)T = π 10ω . This is done in order to obtain the intercept as close as possible to the asymptotic value, minimizing the errors due to numerical artifacts and oscillation.
Following the above prescription, in Fig. 7 we plot d 2 dη 2 ln |G(0, η)| and d 2 dη 2 ln |Q(0, η)| as functions of η for N f = 3. The corresponding plots for other values of N f also display the same qualitative behavior. For large η, the shapes of the graphs approach that of the function in Eq. (A2), displaying periodic local maxima below the η-axis. This provides another justification for the proposed asymptotic form, (28), for G(0, η) and Q(0, η).
The method outlined above is employed to determine the asymptotic forms for G(0, η) and Q(0, η) at N f = 2, 3, 6. In particular, for each N f , ∆η and η max , we find the values of ω G , ω Q , ϕ G and ϕ Q from the largest maximum (η * ) of the graphs in Fig. 7, which, in turn, correspond to the function in (A2), in order to get as close as possible to the asymptotic behavior at large η. The frequencies are found by using Eq. (A3), while the phases are extracted using ω η * + ϕ = π n, with the integer n adjusted as described above. The parameters α G and α Q can be deduced from the slope of the function in Eq. (A4).
As expected for any numerical and asymptotic solution, the resulting intercepts, α G and α Q , frequencies, ω G and ω Q , and initial phases, ϕ G and ϕ Q , all vary slightly with the step size ∆η and the maximum η max of the computation range η ∈ [0, η max ]. Since the exact continuum and asymptotic solution corresponds to the limit where ∆η → 0 and η max → +∞, we perform the computation with various ∆η's and η max 's, obtaining these parameters for each computation. Then, these values of the parameters are fitted with second-order polynomials in ∆η and 1/η max (cf. [2]). For each parameter, the value of its best-fit quadratic surface at ∆η = 0, 1 ηmax = 0 is taken to be our numerical estimate for the parameter. This technique is employed to obtain the estimates of α G , α A , ω G , ω A , ϕ G and ϕ A for N f = 2, 3, 6 in the limit ∆η → 0, 1 ηmax → 0. Fig. 8 displays by the dots the values of these parameters at N f = 3 for each pair of ∆η and η max that we performed the computation for, together with the best-fit quadratic surfaces. For completeness, let us list the equations describing the best-fit quadratic surfaces (for N f = 3): The qualitative features of these plots and quadratic fit functions are similar for N f = 2, 6, but we omit them for brevity. The resulting values of the parameters extracted with the quadratic fit are given in Table I of the main text. For each of the parameters, the quadratic fit is also tested with Akaike information criterion, giving the R 2 of at least 0.999. The numerical errors shown in Table I are derived by taking the difference between the ∆η → 0 and η max → +∞ extrapolations of the quadratic and linear fits in ∆η and 1/η max to the data points in Fig. 8. One may wonder why α G and α Q in Fig. 8 approach ∆η = 0 and 1 ηmax = 0 from below, while in [2] the intercept approached the same limit from above. To better understand the difference, we re-ran the large-N c evolution simulations done in [2] with the initial conditions different from those used in [2]. We used G (0) = 1, while in [2] Born-level initial conditions (30) were employed. Using the trivial initial condition G (0) = 1 resulted in the intercept approaching the ∆η = 0 and 1 ηmax = 0 limit from below for the large-N c evolution. We thus conclude that, while the asymptotic and continuum value of the intercept appears to be independent of the initial conditions, the approach to this intercept in finite-step-size and finite-η-range numerical simulations appears to depend on the initial conditions.  8: Plots of parameters α G , α Q , ω G , ω Q , ϕ G and ϕ Q as functions of ∆η and 1/η max for N f = 3 and N c = 3. The dots represent our numerical evaluation, while the solid surfaces depict the best fits (quadratic in ∆η and 1/η max ) used for extrapolating to the continuum asymptotic values at ∆η = 0 and 1 ηmax = 0.
As a final cross check for the asymptotic form (28), we plot e −α G η G(0, η) and e −α Q η Q(0, η) versus η in Fig. 9. We see that the functions display clear sinusoidal pattern for η 10, demonstrating sinusoidal oscillation in the large-η asymptotics, as expected from the ansatz (28).