On the small-$x$ behavior of the orbital angular momentum distributions in QCD

We present the numerical solution of the leading order QCD evolution equation for the orbital angular momentum distributions of quarks and gluons and discuss its implications for the nucleon spin sum rule. We observe that at small-$x$, the gluon helicity and orbital angular momentum distributions are roughly of the same magnitude but with opposite signs, indicating a significant cancellation between them. A similar cancellation occurs also in the quark sector. We explain analytically the reason for this cancellation.


I. INTRODUCTION
Over the past decades, tremendous effort has been poured into determining the partonic helicity contributions to the nucleon spin. It has been known for quite some time that quarks' helicity ∆Σ accounts for only a quarter of the nucleon spin. A recent NLO global QCD analysis has found a nonzero contribution from the helicity of gluons ∆G [1]. When combined, these two contributions come closer to, but still fall short of the expected value of 1 2 . One might expect that the remaining discrepancy could be resolved by a precise future determination of the gluon helicity distribution ∆G(x) in the small-x region where the current theoretical uncertainties are overwhelmingly large.
However, a priori there is no reason to expect that the nucleon spin entirely originates from partons' helicity. As is clear from the Jaffe-Manohar sum rule [2], the resolution of the spin puzzle cannot be complete without a full understanding of the orbital angular momentum of quarks L q and gluons L g [3,4]. Unfortunately, at the moment very little is known about the actual value of L q,g , and the community still has a long way to go in extracting them from experiments. The recent proposals of observables for L q,g [5][6][7] as well as the first lattice QCD computation of L q [8] are particularly encouraging in this direction.
In this paper, we investigate the QCD evolution of the orbital angular momentum. The four terms in (1) actually depend on the renormalization scale Q 2 [9]. Moreover, they can be written as the integral over Bjorken-x of the corresponding partonic distributions. For ∆Σ and ∆G, these are the usual polarized parton distributions ∆Σ(x) and ∆G(x). A less known fact is that the x-distributions for L q,g can also be defined [3,[10][11][12][13] L q,g (Q 2 ) = 1 0 dxL q,g (x, Q 2 ). ( The Q 2 -evolution of L q,g (x, Q 2 ) has been previously studied in [14] by solving the renormalization group equation for the moments L q,g (j, Q 2 ) = We shall be particularly interested in the small-x behavior of L q,g (x) which was not a focus of interest in [14], but is phenomenologically important in view of the recent controversy over the large uncertainties in ∆G. Indeed, we find that there is a significant cancellation between the helicity and orbital angular momentum distributions at small-x both in the quark and gluon sectors. We explain analytically that such a cancellation is a robust feature of the evolution equation.
In this work, we restrict ourselves to the Q 2 -evolution equation, and do not discuss the evolution equation in x. The latter requires an intricate resummation of double logarithms α s ln 2 1/x which has recently enjoyed renewed interest for the helicity distributions [15][16][17][18][19], but not yet for the orbital angular momentum distributions. While we push our numerical simulation down to very small values of x, one has to keep in mind that at some point the DGLAP equation breaks down, and should be smoothly taken over by the small-x equation. Where and how exactly this transition occurs is presently not understood.

II. EVOLUTION EQUATION
Let us express the four terms in (1) as the first moment in Bjorken-x of the corresponding partonic distribution functions where f is the flavor index. The polarized quark ∆q f (x), antiquark ∆q f (x) and gluon ∆G(x) helicity distributions are standard, whereas the orbital angular momentum distributions for quarks L f (x), antiquarksL f (x) and gluons L g (x) are perhaps less familiar. Their operator definitions in the light-cone gauge (A + = 0) were first introduced in [11] (see also [10]). The fully gauge invariant definitions of L q,g (x) in QCD have been obtained in [12] 1 (see, also, [3,13]) where it has been shown that L q (x) and L g (x) are not usual twist-two distributions, but have both the twist-two ('Wandzura-Wilczek') and genuine twist-three components, so like the g 2 (x) structure function for the transversely polarized nucleon. In this paper we only consider the singlet distributions 1 In Ref. [12], Lq,g(x) have been defined such that These are related to the present convention as (Note that Lg(x) is an even function of x.) The four distributions ∆Σ(x, Q 2 ), ∆G(x.Q 2 ). L q,g (x, Q 2 ) satisfy the renormalization group equation in Q 2 . For ∆Σ(x) and ∆G(x), this is the well-known DGLAP equation which reads, to leading order, where 3 , n f being the number of flavors and β 0 = 11 − 2n f 3 . The corresponding equation for L q,g (x, Q 2 ) has been implicitly derived in [10] to one-loop and explicitly written down in [5]. Because L q,g (x) have a twist-two component, they mix with ∆q(x) and ∆G(x) under renormalization as whereP Integrating both sides of (12) over x from 0 to 1, we obtain in agreement with [9]. One can check that 2 that is, the total angular momentum 1 2 is conserved.

III. NUMERICAL RESULTS
In this section we show the result of our numerical simulation of the coupled equations (7) and (12) directly in the x-space. We set n f = 3 throughout and use the one-loop running coupling constant α s (Q 2 ) = 4π β 0 ln Q 2 /Λ 2 with Λ = 0.25 GeV. The initial scale is chosen to be Q 2 0 = 1 GeV 2 . As for the initial condition, in principle one should draw guidance from the existing global analysis fits [1,[20][21][22][23][24]. However, this is not practical for the present purpose due to several reasons. Firstly, since our primary interest is in the small-x region, the details of parameterization in the large-x region are presumably not very important. Besides, currently the uncertainties of ∆G(x) in the small-x region are enormous. Moreover, there has been no global analysis study of the orbital angular momentum distribution. (There are, however, some model calculations [26][27][28] and an estimate in the Wandzura-Wilczek approximation [29].) We therefore restrict ourselves to the following very simple models: The nucleon spin is equally distributed to the four terms in (1) at the initial scale. The distributions are nonsingular as x → 0. For simplicity, we assume they are constant in x.
• Helicity dominance model Initially the helicity contributions alone saturate the sum rule. We try where A q and A g are fixed by the conditions ∆Σ(Q 2 0 ) = 1 4 and ∆G(Q 2 0 ) = 3 8 . We vary the parameter a g to explore different possibilities. 3 Fig. 1 shows ∆Σ(x), ∆G(x), L q (x), L g (x) in the democratic model as a function of x (left) and rapidity Y ≡ ln 1/x (right) at Q 2 = 10 GeV 2 . We find that at small-x, ∆Σ(x) and L g (x) turn negative and significantly cancel L q (x) and ∆G(x), respectively. In Fig. 2, we plot the quantity as a function of x min at Q 2 = 10 GeV 2 (left) and Q 2 = 100 GeV 2 (right). For comparison, we also plot the helicity part alone Depending on the value of Q 2 , the helicity contribution undershoots (small-Q 2 ) or overshoots (large-Q 2 ) the total spin 1/2. In either case, in this model the helicity contribution becomes the dominant part of the total spin, although initially it has the same magnitude as the orbital angular momentum contribution. Fig. 3 shows the result for the helicity dominance model with a g = −0.6 at Q 2 = 10 GeV 2 (left) and Q 2 = 100 GeV 2 (right). Although ∆Σ(x) is initially positive and has a power-law divergence, after the evolution again it turns negative and partly cancels L q (x). L g (x) is initially zero, but it quickly develops a strong singularity and significantly cancels ∆G(x). Fig. 4 is the plot of (25) in this model. We see that even though the orbital angular momentum vanishes at Q 2 = 1 GeV 2 ,  already at Q 2 = 10 GeV 2 it plays a crucial role to fulfill the spin sum rule. We also notice that very little orbital angular momentum is generated in the large-x region, x > 0.1. The result for a g = −0.3 is qualitatively similar. ∆Σ(x) turns negative even though initially it is as strongly divergent as ∆G(x).
Comparing the two models, we see that the growth of orbital angular momentum at small-x is mainly governed by the parameter a g , rather than the initial value L g (Q 2 0 ). If the initial helicity distribution is singular a g < 0, a significant amount of orbital angular momentum is generated at small-x. It is thus important to better constrain the value a g in future global analyses. Finally, in Fig. 5 we plot the quark and gluon contributions to the nucleon spin separately in the two models at Q 2 = 10 GeV 2 . Due to the cancellation, the two curves flatten quickly in the democratic model. In the helicity dominance model, a gradual rise of the gluon angular momentum is observed down to x ∼ 10 −3 , implying that the cancellation is not complete.

IV. ANALYTICAL INSIGHTS
The significant cancellation between ∆G(x) and L g (x) and also between ∆Σ(x) and L q (x) we found numerically in the previous section is phenomenologically important and calls for a theoretical explanation. Let us try to understand by analytical means how such a cancellation can arise from the structure of the evolution equation.
Since |L g (x)|, |∆G(x)| |∆Σ(x)|, |L q (x)| at small-x, to first approximation we may ignore ∆Σ(x) and L q (x) altogether. The equation then reads Let us first consider the double logarithmic approximation (DLA) which is familiar in the context of the unpolarized distributions but can be readily generalized to the helicity distributions [30][31][32].
In deriving (29), it has been assumed that the DLA saddle point j = j 0 is to the right of all poles in the complex j-plane. This may not be the case if the initial condition (24) has a singularity a g ≡ −c < 0. When j 0 > c, the DLA saddle point rules and the initial power-law is washed out. On the other hand, when c > j 0 , an extra term appears, and this will dominate over (29). The boundary j 0 = c forms a two-dimensional surface in the (Y, ln Q 2 , c) space as illustrated in Fig. 6. A recent global analysis [25] found 4 a rather large value 1 c 0.5 for ∆G(x) at Q 2 = 4 GeV 2 and x < 10 −3 . This suggests that the power-law regime (32) indeed exists and is accessible in the present and future experiments. We now show that in the power-law regime, the ratio L g (x)/∆G(x) is directly related to the 'Regge intercept' c. For this purpose we look for a solution whose leading singularity is of the form In order for the first moments to be finite, we assume 1 > c. (The boundary value c = 1 is actually interesting, see below.) Substituting (34) into (27), we can perform the z-integral analytically. 5 In practice, the lower limit of the integral can be sent to zero since the neglected terms are subleading in x. We thus obtain where is the harmonic number.) α, β, δ are the anomalous dimensions [10] analytically continued to noninteger values j → c. 6 Eq. (35) can be readily solved as 5 To be more precise, since (34) is valid only for x x0 with some x0 < 1, one has to divide the z-integral into different regions x/x 0 dz. It is easy to see that only the second integral, where the form (34) can be used, contributes to the leading singularity. 6 The entries in (28) are the first two terms of α, β, δ in the limit c → 0.
where C, C are the integration constants. Since δ > 0 and δ > α for 1 ≥ c ≥ 0, the second term in A is formally subleading (though it may not be numerically subleading in practice). This gives for sufficiently large Q 2 . We see that L g (x) and ∆G(x) have opposite signs at small-x, and that L g (x) is larger in magnitude than ∆G(x). The DLA corresponds to the limit c → 0 where the anomalous dimensions diverge and the distributions depend logarithmically on x (instead of a power-law).
Similarly, the small-x behavior of ∆Σ(x) is governed by the exponent c. Keeping only ∆G(x) from (34) on the right hand side of (7), one finds This immediately gives The ratio is plotted in Fig. 7 (lower curve). ∆Σ(x) has an opposite sign to ∆G(x) and its magnitude is strongly suppressed as c → 1. This is consistent with the numerical results in the previous section. Finally, for the quark orbital angular momentum, we find so that This is also plotted in Fig. 7 (upper curve). The limit c → 1 is particularly interesting. In this limit, We see that the boundary value c = 1 is permissible from the sum rule point of view. While ∆G(x) and L g (x) are both logarithmically divergent, the divergent parts cancel exactly. Interestingly, the first relation in (45) formally agrees with the argument in [5] based on an operator analysis without any reference to the small-x behavior of ∆G(x), L g (x). 7 It has been also observed in an explicit model calculation in [28]. In practice, it is difficult to numerically confirm the ratios and exponents obtained above. This is because the approach to the asymptotic regime is slow, especially due to the running of the coupling. For realistic values of Q 2 and x, the subleading corrections which depend on the initial condition are still not negligible. For example, in Fig. 3 we find |∆G(x)| > |L g (x)|, although asymptotically the inequality should be reversed, see (40). (Note that in this model L g (x) = 0 initially, so the second term of A in (39) can be comparable to the first term.) Also, whether the Y -dependence is exponential e #Y or DLA-like e # √ Y , is difficult to determine because of the limited Y -range. We however checked that as c = −a g is increased, a DLA-like fit becomes less and less favorable. 8 The picture emerging from our analysis is that, for any value of c between 0 and 1, among the four terms in the Jaffe-Manohar decomposition in the x-space, there is a significant cancellation between the first two terms (quark sector) and also between the last two terms (gluon sector). Note that when c = 1 exactly, in practice one is computing the anomalous dimensions relevant to the first moments (3), see (21). Therefore, a similar cancellation occurs among the integrated quantities 1 2 ∆Σ + L q and ∆G + L g as was observed in [9,33]. We have shown that such a cancellation already occurs in the density space at small-x, and the deviation from exact cancellation is controlled by the Regge intercept c.

V. CONCLUSIONS
In this paper we numerically solved the QCD evolution equation for the orbital angular momentum distributions. Compared with the previous work [14], our work is focused on the small-x region where an interesting cancellation occurs between L g (x) and ∆G(x) and also between ∆Σ(x) and L q (x). For ∆G(x), such a cancellation has been previously suggested in [5,28] from different arguments. As we demonstrated analytically, in the present approach this naturally follows from the structure of the evolution equation.
Our finding has an important implications for phenomenology. On one hand, the precise value of ∆G is of intrinsic interest in QCD, and it is certainly imperative to reduce the uncertainties of ∆G(x) in the small-x region in future experiments such as at the planned Electron-Ion Collider (EIC). On the other hand, this is not sufficient to solve the nucleon spin puzzle because a good fraction of the would-be spin from ∆G(x) at small-x is canceled by the orbital angular momentum 7 In [5], there was a mistake in the normalization of Lg(x) from Section III on. Because of this, the correct relation Lg(x) ≈ −∆G(x) was incorrectly presented as Lg(x) ≈ −2∆G(x) and this coincided with the DLA result in the appendix of [5] (where the normalization is correct), creating an apparent 'consistency. in the same x-region. This suggests that one has to look into the orbital angular momentum in the large-x region [6,7]. After all, this is a very natural and obvious future direction of research.
As already mentioned in the introduction, the DGLAP-type evolution equation considered in this paper eventually breaks down and should be superseded by the small-x evolution equation which resums double logarithmic contributions (α s ln 2 1/x) n . (Not to be confused with the DLA which resums powers of α s ln 1/x ln Q 2 .) For the helicity distributions, it is known that such a resummation dynamically generates a power-law behavior [15][16][17][18][19]. Furthermore, there may be a regime where nonlinear evolution equations come into play, as is the case for the unpolarized distributions. Unfortunately, at the moment very little is known about the small-x resummation for the orbital angular momentum distributions. Admittedly, we may have pushed our numerical solution to too small values of x for which the present approach is not justified and an alternative approach is needed. Still, one can naturally expect that the rapid growth of the distributions, either due to the DLA or the Regge behavior augmented by the QCD evolution, is smoothly connected with the power-law generated by the small-x resummation. Clearly this issue deserves further study. It is straightforward to generalize (28) by including the quark distributions ∆Σ(x) and L q (x).
The following linear combinations diagonalize the evolution equation Asymptotically, S 4 (x) dominates and we have the relation This may be compared with the c → 0 limit of (42), ∆Σ(x) ≈ −0.25∆G(x) (for n f = 3). The difference is because (42) is obtained by first approximating |∆Σ(x)| |∆G(x)|.