The Spin Budget of the Proton at NNLO and Beyond

We revisit the scale evolution of the quark and gluon spin contributions to the proton spin, $\frac{1}{2}\Delta \Sigma$ and $\Delta G$, using the three-loop results for the spin-dependent evolution kernels available in the literature. We argue that the evolution of the quark spin contribution may actually be extended to four-loop order, and that to all orders a single anomalous dimension governs the evolution of both $\Delta \Sigma$ and $\Delta G$. We present analytical solutions of the evolution equations for $\Delta \Sigma$ and $\Delta G$ and investigate their scale dependence both to large and down to lower"hadronic"scales. We find that the solutions remain perturbatively stable even to low scales, where they come closer to simple quark model expectations. We discuss a curious scenario for the proton spin, in which even the gluon spin contribution is essentially scale independent and has a finite asymptotic value as the scale becomes large. We finally also show that perturbative three-loop evolution leads to a larger spin contribution of strange anti-quarks than of strange quarks.


Introduction
The decomposition of the proton spin in terms of the contributions by quarks and anti-quarks, gluons, and orbital motion is a key focus of modern nuclear and particle physics. As has become well-known [1,2,3], in a gauge theory the decomposition is not unique. The two physically most relevant spin sum rules for the proton are the Ji decomposition [2], which ascribes the proton spin to gauge-invariant contributions by quark spins and orbital angular momenta, and total gluon angular momentum, and the Jaffe-Manohar decomposition [1], in which there are four separate pieces corresponding to quark and gluon spin and orbital contributions, respectively. The two sum rules have in common only the quark spin piece, and there is no relation among the other pieces. In particular, in the gauge-invariant definition of Ref. [2] only the total gluon angular momentum is well defined and cannot be split in a physically meaningful way into helicity and orbital angular momentum contributions.
The Jaffe-Manohar sum rule corresponds to the canonical decomposition of the proton's angular momentum. It may be regarded as a "partonic" spin sum rule, since both the quark and gluon spin pieces are related to parton distributions measurable in inelastic ℓp or pp scattering processes. The sum rule reads where 1 2 ∆Σ and ∆G are the quark and gluon spin contributions and L q and L g the orbital ones. ∆Σ and ∆G may be obtained from the first moments of the helicity parton distributions ∆q(x, Q 2 ), ∆q(x, Q 2 ) (where q = u, d, s, . . .) and ∆g(x, Q 2 ) of the proton: where in the first line the sum runs over all active quark flavors whose number we denote by N f . In the following we will mostly use the simplified notation and likewise for the anti-quarks. We note that although the parton distribution ∆g(x, Q 2 ) and hence its first moment are gauge-invariant, the identification with the gluon spin contribution is only valid in the light-cone gauge. The same is true for the orbital angular momentum pieces.
As indicated in Eq. (1), the contributions to the proton spin are all scale dependent, although the dependence cancels in their sum. The dependence on Q 2 is given by spin-dependent QCD evolution equations. The kernels relevant for the evolution of the first moments ∆q(Q 2 ), ∆q(Q 2 ), ∆G(Q 2 ) have been derived to lowest order (LO) in Refs. [4,5], to next-to-leading order (NLO) in Refs. [6,7,8,9], and recently to next-to-next-to-leading order (NNLO) in Refs. [10,11,12]. The kernels for the separate evolution of L q and L g in the Jaffe-Manohar decomposition are known only to LO [13,14,15], although the evolution of their sum is known from Eq. (1) to the same order as that of 1 2 ∆Σ + ∆G, that is, to NNLO.
In this letter, we will discuss some features of the higher-order scale evolution of the first moments, starting from the NNLO results of [10,11,12] for the evolution kernels. We will first study the singlet evolution, where we will extend previous arguments by Altarelli and Lampe [16] to show that the evolution of ∆Σ may actually be determined even to next-to-next-to-next-to-leading order (N 3 LO). The evolution equations for ∆Σ and ∆G may straightforwardly be decoupled and solved in closed form. We present the solutions in analytical form and show numerical results for their evolution at various perturbative orders. We note that studies along these lines were first presented to LO in Refs. [17,18,19,20]. The paper [21] considered the NNLO evolution of ∆Σ. In our paper we go beyond the previous work by extending the results for the evolution of ∆Σ to N 3 LO and that of ∆G to NNLO. Apart from the intrinsic value of this, we believe that our results could also have interesting applications in comparisons to models and in studies of nucleon spin structure in lattice QCD [22]: Although nowadays renormalization on the lattice is typically performed at nonperturbative level, comparison to high-order perturbative evolution should be valuable as a benchmark.
As is well-known [16,23,24] the gluon spin contribution ∆G(Q 2 ) in general evolves as the inverse of the strong coupling constant and thus rises logarithmically with growing scale, either to large positive or negative values, depending on the input ∆G(Q 2 0 ), ∆Σ(Q 2 0 ) at some scale Q 0 . However, in between there is a unique solution for which the gluon spin contribution remains almost flat in Q 2 and tends to a finite asymptotic value. Such a "static" solution in fact occurs [25] in the early NLO DSSV analysis [26,27]. We show that static solutions may be found at every order in perturbation theory and determine the asymptotic values for the gluon spin contribution at LO, NLO, and NNLO.
We will finally also study the evolution in the flavor non-singlet sector. Higher-order evolution is known to generate interesting patterns of flavor-or charge-symmetry breaking in the nucleon sea. It was shown a long time ago [28,29,30] that NLO evolution leads to an asymmetryū =d both in the unpolarized and the helicity parton distributions. At NNLO, a new type of valence splitting function emerges [31,12,32], which gives rise to a difference in the strange and antistrange parton distributions in the nucleon [33], just from the fact that the nucleon carries net up and down valence distributions. In Ref. [33] estimates for the spin-averaged s(x, Q 2 ) −s(x, Q 2 ) were given that showed that the asymmetry resulting from evolution is not as small as might be expected from a three-loop effect. Of course, non-perturbative physics may well be the dominant source of the strangeness asymmetry in the nucleon [34]. In the present paper we will extend the perturbative study in [33] to the spin-dependent case. An interesting difference with respect to the spin-averaged asymmetry is that the first moment 1 0 dx ∆s(x, Q 2 ) − ∆s(x, Q 2 ) does not have to vanish, whereas 1 0 dx s(x, Q 2 ) −s(x, Q 2 ) = 0 due to the fact that the nucleon does not carry net strangeness. As a result, strange quarks and anti-quarks may make different contributions to the proton spin. Indeed, as will be a result of this paper, such a net strangeness helicity asymmetry arises from NNLO evolution.

Evolution equations
We start by considering the generic evolution equation for the first moment of a spin-dependent parton distribution a, b ≡ u,ū, d,d, s,s, . . . , G: where ∆P ab describes the splitting b → a (it is the first moment of the usual x-dependent splitting function). The ∆P ab are perturbative in the strong coupling α S ; their perturbative series starts at with a s ≡ α S /(4π). The running coupling obeys the renormalization group equation where with C F = 4/3 and C A = 3. Keeping just the first term in each of Eqs. (5) and (6) yields the leading order (LO) evolution of the parton distributions. Taking into account also the second, or the second and third, terms corresponds to next-to-leading order (NLO) and next-to-next-toleading order (NNLO) evolution, respectively.
The evolution equations may be simplified by introducing non-singlet and singlet combinations of the quark and antiquark distributions; see e.g. Ref. [29]. Following the notation of Ref. [31] and using charge conjugation invariance and flavor symmetry of QCD, we first write the evolution kernels ∆P ab as The splitting functions ∆P S qq and ∆P S qq thus describe splittings in which the flavor of the quark changes. Starting from NNLO, ∆P S qq and ∆P S qq differ [29,33,12].
We now introduce three types of flavor non-singlet combinations of parton densities: which turn out to diagonalize the evolution equations in the non-singlet sector. (Up to NLO it would be sufficient to consider only two non-singlet combinations; owing to ∆P S qq = ∆P S qq , it becomes necessary to introduce a third combination at NNLO and beyond.) Each of the three combinations evolves in a simple closed form: where the corresponding evolution kernels are The decoupled non-singlet equations are trivial to solve; we will present the solutions later.
In the singlet sector defined by Eq.
(2) we have coupled evolution equations for ∆Σ and ∆G: where and with the first moments of the splitting functions involving gluons, ∆P qG , ∆P Gq , ∆P GG . As we shall discuss below, thanks to the simplicity of the evolution kernels in the spin-dependent case the singlet equation may also be solved analytically in a simple way.
The evolution of the helicity parton distributions is in itself closed and not affected by contributions from orbital angular momentum. On the other hand, L q and L g in Eq.
(1) are both scale dependent and have their own evolution equations. As it turns out, their evolution is not closed but is partly driven by ∆Σ(Q 2 ) and ∆G(Q 2 ). This has to be the case since the left-hand-side of Eq. (1) needs to remain independent of the scale. Presently, the evolution of L q and L g is known only to lowest order [13,14,15] * . Beyond LO, we therefore cannot separate the evolution of the two orbital components. However, we can still consider the evolution of the total orbital angular momentum L ≡ L q + L g by simply taking the derivative of Eq. (1): This relation will serve as an important cross-check for future calculations of the separate evolution of L q and L g at higher orders. We note that, like for the helicity quark and antiquark distributions, there will be a separate angular momentum piece L q i for each flavor, and full evolution equations will require introduction of non-singlet and singlet combinations. L q as appearing in the spin sum rule is the singlet. * We note that the evolution for the total quark and gluon angular momenta in the Ji decomposition may be derived by profiting from the relation between the total angular momentum operators and the quark and gluon energy momentum tensors [2] and is actually known up to NNLO accuracy [35,36]. Unfortunately, since there is no direct connection between the Ji and Jaffe-Manohar spin decompositions (except for the quark spin piece), it is not possible to use these results to obtain the higher-order evolutions of L q and L g .

First moments of the splitting functions
We now collect the various ∆P ab as available from the literature. At lowest order we have [4,5] The second-order results in the MS scheme may be found in Refs. [7,8,9].
with ζ i ≡ ζ(i) the respective value of Riemann's zeta function. Finally, at NNLO we have from Refs. [6,11,12]: In the above equations, C F = 4/3, C A = 3, N f is the number of flavors, and d abc d abc /C A = 5/18.
There are systematic patterns among the above results which may be understood from general arguments. First of all, ∆P + has to vanish to all orders in the strong coupling. As follows from Eqs. (9,10), ∆P + governs the evolution of combinations such as ∆u + ∆ū − (∆d + ∆d), which correspond to matrix elements of flavor non-singlet axial currents. Such matrix elements are not renormalized and are hence scale independent [37] to all orders, consistent with the Bjorken sum rule [38]. † As is well known (see Eqs. (12,13)), the vanishing of ∆P + immediately implies that the evolution of the singlet ∆Σ is to all orders driven by the "pure-singlet" anomalous dimension N f ∆P S qq + ∆P S qq . The explicit results shown in the above equations suggest that or, generalized to all orders, Furthermore, we deduce from Eqs. (16,17) ∆P qG = 0 , The all-order results just given may in fact be understood by an argument given in Ref. [16]. The quark singlet combination ∆Σ corresponds to the proton matrix element of the flavor-singlet axial current, where S is the proton's polarization vector. Because of the axial anomaly, the singlet axial current † Evidently, in a perturbative calculation, one could in principle choose a factorization scheme in which ∆P + becomes non-zero at NLO or beyond. However, such a scheme would be unphysical [30].
is not conserved: In the first line, F µν is the gluonic field strength tensor andF µν its dual. In the second line we have used that Tr[FF ] may be written as the divergence of the "anomalous current" that we denote by K. From Eq. (22) we conclude that j µ 5 − 2N f a s K µ is conserved: The relation ∂ µ j µ 5 = 2N f a s Tr[FF ] holds to all orders in perturbation theory. As a result, Eq. (23) holds to all orders as well. As was discussed in [16], in perturbation theory we may relate matrix elements of K µ to the gluon spin contribution: S µ ∆G = − P, S | K µ | P, S . Although K depends on the choice of gauge, its forward proton matrix element is gauge invariant, except for topologically nontrivial gauge transformations that change the winding number. (The latter feature makes the identification of ∆G with the matrix element of K impossible beyond perturbation theory [1].) From the conservation law in (23) we may thus conclude Inserting the general evolution equations for ∆Σ and ∆G in (12), as well as the renormalization group equation for a s (Q 2 ) in (6), we find, on the other hand, The right-hand-side vanishes when the all-order relations given in Eqs. (19) and (20) are satisfied. Equivalent results are found when studying the renormalization of the axial anomaly in dimensional regularization in [39]. One may object that Eqs. (20) do not follow from (25) in a strict mathematical sense; however, there is little (if any) freedom physically to obtain results other than (20) from the last term in (25). In particular, the C A parts in P GG can only be canceled by those in the β-function. The explicit verification to three loops by the results of [11,12] is of course a strong argument for the all-order validity of (20). In addition, the vanishing of ∆P qG in any physical scheme is a consequence of helicity conservation.
We note that as seen in Refs. [11,12,30] relations like ∆P + = 0 and (19) and (20) may not emerge automatically in an actual higher-loop calculation of the splitting functions, where dimensional regularization and a prescription for γ 5 and the Levi-Cività tensor have to be adopted. They may then be reinstated by a factorization scheme transformation, so that the correct physical splitting functions are obtained.
It is now clear that in the MS scheme a single anomalous dimension, ∆P ΣΣ , resulting from the axial anomaly, governs the evolution of the quark and gluon spin contributions and (via Eq. (14)) of the total orbital angular momentum. Inserting our findings into Eq. (12), we obtain where we have dropped the ubiquitous argument Q 2 . We may further simplify this equation by defining [16] ∆Γ From (26) we then have The lower right entry of the evolution matrix now vanishes since in the product a s (Q 2 )∆G(Q 2 ) the evolution of the strong coupling exactly cancels the ∆P GG part of the evolution of ∆G. Clearly, Eq. (28) is straightforward to solve, and we will return to the equation shortly.
Thanks to Eq. (19) we may now determine the four-loop (N 3 LO) contribution to the anomalous dimension ∆P ΣΣ from the three-loop value ∆P (2) Gq computed in Ref. [12]: 4 Higher-order solutions in the singlet sector We now proceed to solve the singlet evolution equation (28). Changing d ln Q 2 to da s via Eq. (6), we have for the evolution of ∆Σ, up to N 3 LO: ΣΣ + a 2 s ∆P where in the second line we have used that ∆P This equation is readily solved analytically. The solution gives the first moment of the singlet at scale Q in terms of its boundary value at the "input" scale Q 0 : where a Q ≡ a s (Q 2 ) and a 0 ≡ a s (Q 2 0 ). For completeness, we have included the LO term, which predicts a constant ∆Σ.  Figure 1 shows the quark singlet evolution factor on the right-hand-side of Eq. (32), assuming a fixed number N f = 3 in the anomalous dimensions and the beta function, and using the full NNLO evolution of the coupling constant ‡ . We choose a relatively low input scale Q 0 = 1, with a value α s (Q 0 ) = 0.404 § . One can see that the NLO evolution affects the quark spin content of the proton by up to 7% while NNLO evolution adds an extra ∼ 1 − 2% effect. The numerical impact of the four-loop term ∆P ΣΣ reaches only O(0.2%) at the highest scale. ‡ Alternatively, one could use at each order a coupling constant defined by truncating the QCD β-function to that order. This approach would mostly affect the LO results, since the coupling constant at this order is larger. It would also slightly limit the range of applicability of the LO calculation in the backward evolution since the nonperturbative regime would be reached already at higher scales. Nevertheless, the main results of this paper would remain basically unchanged since the values of the coupling constant are quite similar at NLO, NNLO and beyond.
§ That value corresponds to the conventional α s (M Z ) = 0.1181 for N f = 5. In this paper we always use the NNLO expression for the coupling constant, independently of the order considered, as a way of isolating the effect of the higher order splitting functions in the corresponding evolution. Ultimately, as discussed in Ref. [21,40], one may want to compare helicity parton distribution functions extracted from experiment or computed on the lattice [22] with calculations performed in QCD-inspired models of nucleon structure. The latter typically are formulated at rather low momentum scales of order of a few hundred MeV. Given the high order of perturbation theory now available for evolution, it is therefore interesting to evolve the singlet spin contributions not only to large perturbative scales, but also "backward" towards the limit of validity of perturbation theory [21]. In Fig. 2 we show the evolution of ∆Σ at LO, NLO, NNLO and N 3 LO down to Q ∼ 0.35 GeV, starting from the initial scale Q = 2 GeV. Since at low scales the approximate analytical expressions for the running of the coupling constant α s start to deviate from the exact result, we rely on the accurate numerical solution of Eq. (6) to NNLO accuracy. As can be observed, and as is expected, the higher order terms affect the evolution of the singlet in a significant way, much more strongly than what we found for the evolution to larger scales. On the other hand, a striking feature is that the evolution remains relatively stable even down to such low scales as considered here: At the lower end of Figure 2 the N 3 LO contribution enhances the singlet by a modest 8% compared to the previously known NNLO result, despite the fact that at Q = 0.35 GeV the coupling constant becomes α s ∼ 1.3, precariously past the boundaries of the perturbative domain. In addition, all higher orders (NLO, NNLO, N 3 LO) go in the same direction. We note that the upturn of ∆Σ toward small scales -in the direction of large quark and anti-quark spin contributions to the proton spin -was already observed to NLO and NNLO in Refs. [40] and [21], respectively. We also remark that results on high-loop evolution may be useful for lattice-QCD studies of nucleon structure, possibly allowing cross-checks of the nonperturbative renormalization carried out on the lattice.
The solution of the evolution equation for the gluon spin contribution now follows directly. From the lower row in Eq. (28) we have by simple integration and using again d ln where again a Q ≡ a s (Q 2 ) and a 0 ≡ a s (Q 2 0 ), and where in the second line we have used Eq. (19) to replace ∆P ΣΣ by ∆P Gq , which is more natural in the case of the gluon distribution.
An immediate observation is that the integral on the right-hand-side of (33) starts at order a s (Q 2 ) and a s (Q 2 0 ). Therefore, we arrive at the well-known result [17,18,23] that the leading term in ∆Γ is a constant in Q 2 , so that the first moment of the gluon spin contribution evolves as the inverse of the strong coupling. Inserting the solution for ∆Σ(Q 2 ) from Eq. (32) into (33) and carrying out the integration, we find the full NNLO analytical solution where with We note that in contrast to ∆Σ, we can only give the NNLO evolution of ∆G here. This is due to the fact that ∆Γ is shifted by one power of a s relative to ∆G. In order to obtain the N 3 LO solution for ∆G one would need the four-loop splitting kernel ∆P Gq which is presently still unavailable. Figure 3 shows the NNLO evolution of the gluon spin contribution to the proton spin, starting from the values ∆G = 0.102 and ∆Σ = 0.254 at Q 0 = 1 as realized in the global analysis [46]. We also show the evolution of 1 2 ∆Σ and the evolution of the total orbital angular momentum L resulting from (14). Notice that both ∆G and L have a divergent behaviour at large scales, resulting in a rather unphysical cancellation of two very large contributions to fulfill the spin sum rule. Figure 3: Evolution of the quark and gluon spin contributions 1 2 ∆Σ and ∆G at NNLO, starting from the inital scale Q = 1 GeV. We also show the evolution of L following from Eq. (14).
As we discussed above for the singlet contribution, it is also interesting to analyze the behavior of the gluonic spin contribution at lower scales. In Fig. 4 we show the backward evolution of ∆G at LO (dashes), NLO (dots) and NNLO (solid line) for three different scenarios, corresponding to setting ∆G(Q 0 ) = +1, 0.1, −1 at the initial scale Q 0 = 1 GeV ¶ . For each scenario, we observe a striking convergence of the fixed order results down to very low scales, always towards small gluonic contributions. Even though the "F " term in Eq. (34) contains corrections proportional to positive powers of α s that could spoil the convergence of the expansion in the non-perturbative region, the evolution of the gluonic contribution is completely dominated by the leading 1/a s term in Eq. (34), as can be observed in Figure 4 where we also present this term separately for each scenario. Our findings set a strong constraint on the proton spin content carried by gluons at hadronic scales. Within the rather extreme scenarios analyzed here (for which the gluon contribution accounts for as much as twice the spin of the proton at Q 0 = 1 GeV!), we obtain the requirement |∆G(Q ∼ 0.35 GeV)| 0.3. Indeed, the few available model estimates of ∆G suggest values of the order 0.2 − 0.3 [41,42,43,44,45] at a low hadronic scale.

"Static" value of ∆G
As we have discussed, ∆G(Q 2 ) in general evolves as 1/a s (Q 2 ) for large scales. As inspection of Eq. (34) shows, depending on the input values of ∆G(Q 2 0 ) and ∆Σ(Q 2 0 ) the evolution can be towards large positive or negative values. This implies that there is a specific input, a "critical point", for which ∆G(Q 2 ) actually remains almost constant [26] and tends to a finite asymptotic value as Q 2 → ∞. This "static" value of ∆G is expected to change from order to order in perturbation theory. To determine it at a given order, we only need to tune the input such that the ¶ ∆G(Q 0 ) = 0.1 corresponds to the result of the DSSV fit in [46].
To LO, using (36), this condition becomes where we have used N f = 3 flavors and again ∆Σ(Q 2 0 = 1 GeV 2 ) = 0.254. The gluon spin contribution then remains constant at the value ∆G LO stat (Q 2 0 ).
At NLO, the necessary input value for the static solution becomes The NLO "static" solution is no longer completely constant in Q 2 . However, by construction it does converge asymptotically to a finite value, given by We note that a value of similar size was in fact found in the early NLO DSSV analysis [25,26,27].
Finally, at NNLO, the corresponding values are with an asymptotic value given by Numerical results for the "static" solutions for ∆G are shown in Fig. 5. Figure 5: Evolution of the "static" gluon solutions, starting from the inital scale Q 0 = 1 GeV.
At NNLO the sum of the contributions by quarks and gluons starts at ∆Σ(Q 2 0 )/2+G NNLO stat (Q 2 0 ) = 0.0025 with an asymptotic result of ∆Σ(∞)/2 + G NNLO stat (∞) = 0.01335. In that particular scenario the total orbital angular momentum almost accounts for the entire proton spin, L q + L g ≃ 1/2, and is almost constant in Q 2 (from Q 2 0 = 1 to ∞ it varies by less than 3%).
As mentioned above, in our study of the "static" ∆G we have for simplicity chosen a fixed number of flavors, N f = 3. This will not be entirely adequate when considering the limit of large Q 2 , and a matching to N f = 4 and N f = 5 should be performed at the charm and bottom mass scales, respectively. For the inputs ∆G stat (Q 2 0 ) given explicitly above, each matching would slightly upset the cancelation of the 1/a s term in the solution for ∆G, so that the resulting gluon spin contribution would not be entirely "static" anymore. However, this is expected to be a small effect. We have checked that using a fixed number of N f throughout changes the asymptotic value of the static ∆G by less than 10%. In any case, it is clear that a static solution for ∆G exists even if one performs a full matching to N f = 4 and N f = 5: One could always start with an input at the bottom mass, Q 2 0 = m 2 b , that creates a static solution for all higher Q 2 with N f = 5. This solution could then be evolved backward to any lower Q 2 one desires, even to Q 2 = 1 GeV 2 where N f = 3. This result at Q 2 = 1 GeV 2 would then be the input to be used to obtain a static solution with full matching.
We believe these solutions, especially because of the fact that they have a well behaved asymptotic limit at large scales, deserve further attention since they arise as strong boundaries on non-perturbative physics from almost purely perturbative considerations.

Non-singlet evolution of the valence quark spin contribution
We finally turn to the evolution in the non-singlet sector. As discussed in the Introduction, we focus here on the strangeness "valence" spin contribution (∆s−∆s)(Q 2 ) generated by three-loop evolution.
Each of the non-singlet evolution equations in Eq. (10) has the solution where ∆q (A) (Q 2 0 ) is the corresponding input non-singlet combination and the evolution operator U (A) is given by .
We may readily use (43) with A = − and A = V to obtain the solution for a valence quark contribution ∆q − ∆q, resulting in [33] (∆q − ∆q) ( where ∆q (V ) (Q 2 0 ) = q (∆q − ∆q ) (Q 2 0 ) is the total spin-dependent valence distribution in the nucleon as defined in Eq. (9), at the initial scale. The first term on the right represents the homogenous component of the evolution of the valence distribution, which starts at NLO; its explicit expression is identical to the one in Eq. (32) with the change ∆P ΣΣ → ∆P (−) . As follows from Eq. (11), ∆P (V ) − ∆P (−) = N f ∆P S qq − ∆P S qq , which (see Eqs. (15)- (17)) becomes nonzero starting from NNLO. Therefore, the second term on the right of (45) will in general be non-vanishing as well, as long as ∆q (V ) (Q 2 0 ) = 0, which of course is the case. We conclude that NNLO evolution generates an asymmetry ∆s = ∆s in the contribution by strange and anti-strange quarks to the proton spin, even if ∆s = ∆s at the initial scale Q 0 . This is at variance with the spin-averaged case where the first moment of s −s is protected by the fact that there can be no net valence strangeness in the proton and remains zero to all orders.
To NNLO accuracy, the evolution factor in the second term in Eq.(45) reduces to where, as before a Q = a s (Q 2 ) and a 0 = a s (Q 2 0 ). Therefore, assuming ∆s(Q 2 0 ) = ∆s(Q 2 0 ) in order to estimate the purely perturbative effect, we have, to NNLO where in the second line we have inserted the explicit value of ∆P (2)S qq − ∆P (2)S qq from Eq. (17). The last factor on the right is of course just the total valence spin contribution at the initial scale.
(48) Figure 6 shows ∆s V ≡ ∆s − ∆s as a function of Q. The difference reaches −0.001 at Q = M Z . As expected for a three-loop effect, it is small. On the other hand, for the latest extractions of (∆s + ∆s) [47] the relative asymmetry |∆s − ∆s|/|∆s + ∆s| would be of order 1%. Evidently, non-perturbative contributions [34] may well be the dominant source of the polarized strangeness asymmetry. However, the effect we describe here would certainly need to be taken into account in a full analysis. We emphasize that the perturbative asymmetry is robustly predicted to be negative, so that ∆s > ∆s.

Conclusions
We have presented a set of studies of the evolution of the quark and gluon spin contributions to the proton spin at higher orders in perturbation theory, motivated by the recent calculations of the helicity splitting functions at full NNLO [10,11,12]. We have argued that the evolution of ∆Σ is known even to four loops, which may prove valuable for lattice studies, as well as for comparisons to models residing at lower "hadronic" scales. The anomalous dimension relevant for the evolution of ∆Σ and related to the axial anomaly also turns out to generate the evolution Figure 6: Evolution of the perturbatively generated "valence" strange contribution ∆s V ≡ ∆s − ∆s, starting from the inital scale Q 0 = 1 GeV.
of ∆G. The same must then be true for the total orbital angular momentum L q + L g in the Jaffe-Manohar sum rule, although the separate evolutions of L q and L g are presently only known to LO.
We have obtained analytical higher-order solutions for ∆Σ and ∆G and presented numerical results for their evolution. These show a stable upturn of ∆Σ toward low scales, bringing it actually closer to quark model expectations that favor a large quark spin contribution to the proton spin. The gluon spin ∆G, when evolved backwards, shows a remarkable focus towards low values, again in line with quark model assumptions, setting a strong constraint on the gluon contribution at hadronic scales. We have also shown that at every order of the perturbative evolution, there is a unique solution for which ∆G tends to a finite asymptotic value as the scale becomes large. We have estimated the values of ∆G in such a scenario.
We have finally also examined the size of the new effect arising from three-loop evolution in the flavor non-singlet sector, the generation of an asymmetry in the strange and anti-strange contributions to the proton spin. We have found that perturbative evolution predicts ∆s − ∆s to be negative, with a magnitude of order 1% of the total ∆s + ∆s.