Quark spins and Anomalous Ward Identity

We calculate the intrinsic quark spin contribution to the total proton spin using overlap valence quarks on three ensembles of $2+1$-flavor RBC/UKQCD domain-wall configurations with different lattice spacings. The lowest pion mass of the ensembles is around 171 MeV which is close to the physical point. With overlap fermions and topological charge derived from the overlap operator, we verify the anomalous Ward identity between nucleon states with momentum transfer. Both the connected and disconnected insertions of the axial-vector current are calculated. For the disconnected-insertion part, the cluster-decomposition error reduction (CDER) technique is utilized for the lattice with the largest volume and the error can be reduced by $10\%\sim40\%$. Nonperturbative renormalization is carried out and the final results are all reported in the $\overline{{\rm MS}}$ scheme at 2 GeV. We determine the total quark spin contribution to the nucleon spin to be $\Delta\Sigma=0.401(25)(37)$, which is consistent with the recent global fitting result of experimental data. The isovector axial coupling we obtain in this study is $g_A^3=1.256(16)(30)$, which agrees well with the experimental value of 1.2723(23).


I. INTRODUCTION
The decomposition of the proton spin into its quark and glue constituents has long been a puzzle ever since the first deep inelastic scattering (DIS) experiment around three decades ago [1,2] revealed that not all the proton spin originates from the quark intrinsic spin as depicted in the naive quark model, leading to the so-called "proton spin crisis." Now we understand that the proton spin, consisting of quark spin, quark orbital angular momentum, glue spin and glue orbital angular momentum, is the result of complicated QCD dynamics which cannot be described by the quark model. However, the precise proportion of the total proton spin carried by these components remains unclear. On the experimental side, since the integration of the spin-dependent parton distributions over the momentum fraction x gives the fraction of the proton spin which is carried by the corresponding flavor, that is, where µ is the MS scale, the global fit to the experimental data of DIS or Drell-Yan processes for extracting the parton distributions will provide us knowledge about the quark spin contribution to the proton spin. Three recent experimental results from D. de Florian et al. [3], the NNPDF collaboration [4] and the COMPASS collaboration [5] determine the total quark intrinsic spin contribution ∆Σ to be 0.366 +0.015 −0.018 , 0.25 (10) and [0.26, 0.36], respectively. On the lattice side, a recent calculation [6] is carried out with physical pion mass, but with only one single lattice ensemble of N f = 2 clover-improved twisted mass fermions. More careful lattice studies with ensembles of different lattice spacings and different lattice volumes are imperative to push the results to the physical limit and to control the corresponding systematic uncertainties.
In this work, we use overlap fermions on three domain-wall ensembles to calculate the quark spin contribution to the nucleon spin. Since for each quark flavor the intrinsic spin is actually half of the corresponding axial coupling of nucleon, we need to calculate the axial coupling for the flavor-diagonal case. Thus, both the connected-insertions and disconnectedinsertions of the correlation functions need to be included. The anomalous Ward identity is carefully checked to see if any normalization due to lattice artifacts needs to be applied to the axial-vector current to make the identity hold. We actually find that the same normalization constant for the local axial-vector current as used in the isovector case to satisfy the chiral Ward identity also satisfies the anomalous Ward identity. This is not true in general for non-chiral fermions. For the disconnected-insertion part, the cluster-decomposition error reduction (CDER) technique [7] is utilized for the lattice with the largest volume to reduce the statistical error. For the connected-insertion part, an improved axial-vector current is employed such that the finite lattice spacing effect can be reduced. All of our results are matched to the MS scheme at 2 GeV using nonperturbative renormalization. We propose a new renormalization pattern where we separate the connected-insertion part and the disconnected-insertion part from the beginning which is more natural for the lattice calculation and offers more information than the conventional flavor irreducible representation approach.
This paper is organized as follows. The formalism of quark spin and anomalous Ward identity are discussed in Sec. II. In Sec. III we describe all the numerical details of our simulation. Then in Sec. IV, we check the axial Ward identity to address the normalization issue. The bare results of the disconnected contribution are shown in Sec. V. The detailed results of the connected contribution come in Sec. VI. We discuss the renormalization in Sec. VII and make global fits to get the final results in Sec. VIII. A short summary is given in Sec. IX.

II. FORMALISM OF QUARK SPIN AND ANOMALOUS WARD IDENTITY
The quark spin contribution to the nucleon spin is associated from the nucleon matrix element of the flavor-singlet axial-vector current, where A 0 µ is the flavor-singlet axial-vector current The flavor u, d and s contributions to g 0 A are denoted as ∆u, ∆d and ∆s in Eq. (2), so that g 0 A = ∆u + ∆d + ∆s.
A special property of the flavor-singlet axial-vector current is that it satisfies the anomalous Ward identity (AWI) where the Adler-Bell-Jackiw anomaly appears from the Jacobian factor of the fermion determinant due to the U (1) chiral transformation [8] where the pseudoscalar density P f and the topological charge density operator q representing the anomaly are where G a µν is the gauge field strength tensor andG a µν = µνρσ G a ρσ . Note, notations are in Euclidean space and the coupling constant g is absorbed in the definition of the gauge potential A a µ . As far as renormalization is concerned, it is shown that A 0 µ has a two-loop renormalization [9,10] and the topological charge has a one-loop mixture with ∂ µ A 0 µ [10] so that the renormalized AWI in the dimensional regularization scheme becomes with the anomalous dimension γ = −(α s /π) 2 3 8 C F . m R and P R are renormalized quark mass and pseudoscalar density. We see that, the α 2 s renormalization term on the left is the same as that on the right from mixing. Thus, mP and ∂ µ A 0 µ + 2iN f q are renormalization group invariant (the latter to second order at least) and the form of AWI is the same with or without renormalization.
On the lattice, the AWI is preserved by the overlap fermion which is chiral and satisfies the Ginsparg-Wilson relation [11]. The m f P f is renormalization group invariant since Z m Z P = 1 for the chiral fermion and the local version of the topological charge q(x) derived from the overlap operator is equal to 1 16π 2 tr c G µνGµν (x) in the continuum [12][13][14][15], i.e.
where D ov is the overlap operator. In the overlap case, the chiral axial-vector current is derived [16] and one can directly proceed to carry out the renormalization of the chiral axial-vector current perturbatively or non-perturbatively. However, this chiral axial-vector involves a non-local kernel K µ = −i δDov(Uµe iαµ(x) ) δαµ(x) | α=0 and is somewhat involved to implement numerically. We shall use the local current in the present study. As such, it invokes a normalization constant Z 0 A which warrants that the unrenormalized AWI in its 'semiclassical' form (Eq. (5)) is satisfied on the lattice and is itself scale independent. Therefore the normalization and renormalization takes two steps. First, one needs to find the normalization Z 0 A for the local axial-vector current which satisfies the unrenormalized AWI where A 0 µ = f =u,d,s ψ f iγ µ γ 5ψf and P f = ψ f iγ 5ψf are the local axial-vector current we use on the lattice andψ = (1 − 1 2 D ov )ψ is for giving rise to the effective quark propagator which conforms to the form in the continuum.

After the normalization constant Z 0
A is determined, one then takes on the renormalization procedure. We shall discuss the determination of Z 0 A in Section IV after we give the numerical details of the calculation and will carry out the renormalization in Section VII.
Before we check the AWI on the lattice, we shall first give some numerical details of the lattice calculation.

III. NUMERICAL DETAILS
We use overlap fermions [17] as valence quarks to perform our calculation. Since the overlap action preserves chiral symmetry at finite lattice spacing via the Ginsparg-Wilson relation [18], there is no additive renormalization for the quark mass. The effective quark propagator of the massive overlap fermion is the inverse of operator D c +m [19,20] where D c satisfying {D c , γ 5 } = 0 is exactly chiral and can be defined from the original overlap operator D ov as D c = ρDov 1−Dov/2 . The overlap operator can be expressed as D ov = 1 + γ 5 (γ 5 D w (ρ)) where is the matrix sign function and D w is the Wilson kernel with κ = 0.2 (corresponding to parameter ρ = 1.5). As discussed above, another great feature of the overlap operator is that the local version of the topological charge of the gauge field can be defined as q(x) = [12][13][14][15]. The Atiyah-Singer index theorem [21] is satisfied which relates the total topological charge to the index of zero modes of the overlap operator so no multiplicative renormalization is needed for this definition of q. These two features help us to feasibly check the AWI which we can use as a normalization condition in the disconnectedinsertion case. We use multiple partially-quenched valence quark masses to cover a wide range of pion mass using the multi-mass algorithm. More details regarding the calculation of the overlap operator and eigenmodes deflation in the inversion of the fermion matrix can be found in [22].
The three lattice ensembles we use for the calculation are 2+1-flavor domain-wall fermion (DWF) ensembles generated by the RBC/UKQCD collaboration [23,24]. They are labeled as 24I, 32I and 32ID and the detailed parameters of the ensembles can be found in Table I. We have three different lattice spacings and the lowest pion mass at 171 MeV is close to the physical one. To calculate the quark spin or, in practice, to calculate the axial coupling, we need to construct 3-point correlation functions where χ is the nucleon interpolation field, G denotes the source grid and A µ =ψiγ µ γ 5ψ is the local axial-vector current withψ = (1 − 1 2 D ov )ψ for giving rise to the effective quark propagator (D c + m) −1 . The correlation function can have two kinds of current insertions, i.e., the connected-insertion (CI) and the disconnected-insertion (DI), corresponding to two ways of Wick contractions. They are depicted in Figure 1.
For the CI calculation, we use the stochastic sandwich method (SSM) with low-mode substitution (LMS) [25] to better control the statistical uncertainty. We use Z 3 noise grid sources with Gaussian (24I and 32I) or block smearing (32ID) [26] coherently at t src = 0 and t src = 32 in one inversion. The sinks are block smeared and located at different positions with different separations in time from the source. Setups regarding the valence simulation of the CI case are listed in Table II. Technical details regarding the LMS of random Z 3 grid source with mixed momenta and the SSM with LMS for constructing 3-point functions can be found in references [25][26][27]. Due to the fact that multi-mass inversion algorithm is uniquely applicable to the overlap fermion with eigenvector deflation, we calculate 5 ∼ 6 valence masses each for the three lattices. Table II. The details of the overlap simulation in the valence sector for the CI case, including the name of the lattice, the grid type of source G src (the notations such as 12-12-12 denote the intervals of the grid in the three spatial directions; see reference [25] for more details), the number of source grids N src , the positions of sources t src , the grid type of sink G sink , the number of the noises for the sink grids N sink , the source-sink separations (t sink − t src ) and the bare valence quark masses m v q a. For disconnected-insertion calculations, we use the low-mode average (LMA) technique to calculate the quark loops which improves the signal-to-noise ratio particularly for the pseudoscalar and scalar currents. The low-mode part is calculated exactly while the highmode part is estimated with 8 sets of Z 4 -noise on a 4-4-4-2 space-time grid with even-odd dilution and additional time shift. The same Z 3 -noise grid source with smearing as in the CI case is used in the production of the nucleon propagators. We make multiple measurements by shifting the source time-slice to improve statistics; the spatial position of the center of the grid is randomly chosen for each source time-slice to reduce autocorrelation. References [27][28][29] contain more details regarding the DI calculation. When constructing quark loops, we include more valence quark masses to cover the strange region. The bare strange quark mass is determined by the global-fit value at 2 GeV in the MS scale calculated in our previous study [30] and the nonperturbative mass renormalization constant calculated in [31].
To obtain the axial coupling, we construct a ratio of the 3-point correlation function to the nucleon 2-point function where f k is a kinematic factor which is related to the Lorentz index of the current, Γ p is the polarized projector of the nucleon spin, Γ e is the non-polarized projector and C 2 (t f ) = x χ(t f , x)χ(0, G) . The matrix element g A can then be obtained asymptotically g A = R(t f τ, τ 0). However at finite t f and τ , the excited states will contribute to the ratio and we need to extract g A by fitting the ratio to more complicated function forms. A commonly used form with two-state fit reads which assumes only the first excited state has effects and δm is the energy difference between the ground state and the first excited state. In practice, higher excited-states' contribution can alter the value of δm, making it a free parameter, accounting for an effective energy difference.

IV. ANOMALOUS WARD IDENTITY ON THE LATTICE
To verify the AWI in Eq. (9), we note that there is no flavor-mixing in this unrenormalized form. Thus, one can check it for individual flavors and, furthermore, since the lattice calculations of matrix elements are separated in the CI and DI cases as shown in Figure 1, one can separately check the chiral Ward identity for the connected matrix elements In this case, the matrix elements are for the u or d valence quark with quark mass m q . Here, the normalization constant Z A (CI) is due to the fact that we use the local current in this calculation.
Similarly, the AWI for the matrix elements in the DI case is In principle, the normalization constants Z A (CI) and Z A (DI) can be different, especially when non-chiral fermions are used in the lattice calculation and also when the topological charge q is not calculated from the overlap operator D ov as in Eq. (8). We shall check them in the following.

A. Disconnected insertion (DI) case
We shall check the disconnected insertion (DI) case first which has been investigated in our previous study [28]. The anomalous Ward identity (AWI) in the DI case in Eq. (14) relating the nucleon matrix element of the divergence of the axial-vector current A µ to those of the product of the quark mass m q and the pseudo-scalar current P and also to the topological charge term q is an important check for lattice spin calculations involving the flavor-diagonal matrix elements (MEs) of the axial-vector current. This is especially true for the strange quark as it only contributes in the DI. Only properly extracted MEs plus correct lattice normalization will make this identity hold. Our previous work [28] utilized the AWI via the form factors which is actually an extended form of the Goldberger-Treiman relation for the flavor-singlet case at finite momentum transfer q 2 and is expressed as where g A and h A are the axial and induced pseudoscalar form factors respectively from the nucleon matrix element of the axial-vector current and g P (q 2 ) and g Q (q 2 ) are the pseudoscalar and anomaly form factors defined in Eqs. (16) and (17) have the same form separately for the CI and DI, while Eq. (18) for the topological form factor only appears in the DI. Equation (15) can be derived by inserting the currents between nucleon states with momenta p and p and applying the divergence to the nucleon states where E and E are the energies of the source and sink nucleons and q i is the momentum In the earlier study [28] we found that a normalization factor of κ A ∼ 1.4 on the axial-vector side is needed in order to satisfy the identity which is much larger than the isovector normalization constant 1 Z 3 A (24I) = 1.111(6) computed by using the chiral Ward identity in the pion 2-point function case [31]. Since on the right hand side of Eq. (15) we have Z m Z P = 1 and there is no multiplicative renormalization of the topological charge defined by the overlap operator, and, as shown above, the renormalized 2 AWI is the same as the un-renormalized one at two-loop level [28], the factor κ A was believed to be a necessary normalization factor in the DI case for compensating the violation of the AWI induced by lattice artifacts and was used to normalize the DI axial-vector MEs. In this study, we shall have a critical reexamination of this issue. We also make a similar check for the light quarks case of 24I and 32I and for the new 32ID lattice by calculating the following ratio Note that we keep the dependence of the sink time t f and the current time τ for all the form factors and, therefore, the excited-state effects are not handled until we fit the final ratio.
We also check the AWI more carefully at the ME level, in other words, treating ∂ µ A µ as an operator insertion between the nucleon states p and p . The lattice version of the AWI reads 1 We choose to call it normalization constant rather than renormalization constant since it is a finite renormalization which has no scale dependence and deviates from unity only because of finite lattice spacing effects. 2 When we say renormalization, we mean there is a non-zero anomalous dimension and therefore is scale dependent.
where A i (τ, x) is an abbreviated form of the ME p |A i (τ, x)|p ,î is the unit vector along the ith direction and the continuum partial derivative is replaced by the backward-difference on the lattice. This equation cannot be checked directly since the momentum projection is always done before we can take the spatial difference is a good approximation if q i is small enough. So we have a simplified form which can be checked numerically. A second ratio R 2 is thus defined as Furthermore, for the temporal part, by inserting complete sets of intermediate states and using the time evolution operator, the time dependence of the ME can be formulated as up to exponential suppression of the excited-states contamination, where ∆E = E − E is the energy difference between the sink and source nucleon states, such that leading to if the excited-states contamination can be ignored or completely removed by fit. This is actually the counterpart of Eq. (19) and we now have the third ratio to check the AWI   (6) where the value from a constant fit is 1.074 (24). The momentum transfer dependence of R 1 and R 2 are plotted in Figure 3. Each point on the plot of R 1 comes from a two-state fit while the points on the plot of R 2 come from constant fits. In the R 1 plot, except for the first two | q| 2 (the point of the first | q| 2 is not shown in the figure since the corresponding two-state fit fails), the values are basically flat within errors and the fitted value using the middle three points is 1.806(46). If we believe that the ratio R 1 is a proper check of the AWI, this value should be used as a normalization factor. In the R 2 case, all the points lie on a constant line within errors and a constant fit excluding the 4th point gives 1.090 (18), quite consistent with Z 3 A (24I) = 1.111 (6). The problem now is to understand the discrepancy and to determine which one is correct.
Since R 1 and R 3 in Figure 2 are quite similar, we shall only compare R 3 and R 2 . It is easy to see that the only difference between R 2 and R 3 is to use Points of different t f are shifted slightly to enhance the legibility.
∆E A 4 (τ, q) . We have proven that they are exactly the same if there are no excited-state effects. So it is useful to see what A 4 (τ, q) looks like. The left panel of Figure 4 shows the behavior of A 4 (τ, | q| = 2π L ) as a function of τ ; no obvious plateau is discernible even we go to relatively larger t f , which means that the excited-states contamination are large and, in this case, even a two-state fit cannot extract the ME reliably. To be more specific, all the components of the AWI are plotted in the right panel of Figure 4. In this | q| = 2π L case ∆Ea ∼ 0.03, so the values of ∆E A 4 (τ, q) are very close to 0; however the ∂ 4 A 4 values are of order 0.1 which pinpoints the problem. One can ask why the ∂ 4 A 4 case is so different since there should also be some excited-states contamination. The answer is that the AWI is actually a relation of the current operators and it should hold regardless of whether the currents are inserted between two nucleon ground states or excited states. The only problem is that when we use A 4 (τ, q)−A 4 (τ − 1, q) = ∆E A 4 (τ, q) , we are assuming that the ME is the ground-state ME and the ∆E is the energy difference between two ground states which is, apparently, not the case. We can thus conclude that if the conditions τ 0 and t f τ are satisfied, the three ratios will be the same; for finite τ and t f , the ratio of   (6) shows that the values of R 2 are consistent with this normalization factor except for the boundary points. R 1 shows large discrepancy and approaches the horizontal line with large t f .
Points of different t f are shifted slightly to enhance the legibility. quark masses. In other words, we have Z 0

B. CI case
We also check the chiral Ward identity in the CI case. In fact, the violation of the chiral Ward identity in terms of form factors at small momentum transfers in the CI case is also observed and reported in [32], where their formula is equivalent to checking the ratio R 1 . In the CI case, the definitions of R 1 and R 2 are the same as those in the DI case but without the topological charge term. The results of R 1 and R 2 at different t f of the 24I lattice are plotted in Figure 5 as The results of the CI case are similar and the conclusion is the same. The ratio R 2 shows well established chiral Ward identity while R 1 shows violation. The difference of R 1 and R 2 reflects how we treat the A 4 term. Even if the form factors g A (q 2 ) and h A (q 2 ) are calculated using A i only, the ratio R 1 is still problematic since Eq. (19) is used in the derivation and it assumes that the ME of A 4 gives the same form factors without excited-states contamination.
A cure to this problem is to go to large enough t f where the excited-states contamination can be ignored. Unfortunately, this requires much larger statistics. We will test this in the future.

V. DISCONNECTED-INSERTION CONTRIBUTION
As is mentioned above, we use the two-state fit to extract the MEs. Examples of the fitting on the 32ID lattice at the unitary point can be found in Figure 6; both the light and strange quark results are included. We use three source-sink separations t f = 6a, 7a and 8a, which correspond to 0.86 fm, 1.00 fm and 1.15 fm respectively, to carry out the fit. The  (Table III).
The cluster-decomposition error reduction (CDER) technique [7] is used in order to better control the statistical uncertainties for the 32ID lattice where the CDER technique may improve the signal more significantly since the size of this lattice L ∼ 4.6 fm is relatively large, while we do not use this technique for the 24I and 32I lattice due to their small sizes (L ∼ 2.7 fm and L ∼ 2.6 fm respectively). In order to use the CDER technique, the 3-point functions can be rewritten as where we put a cutoff R to the distance between the quark loop and the sink of the nucleon propagator and we can vary R to obtain different 3-point functions. As demonstrated in [7], the signal will saturate after R is larger than the corresponding correlation length but the noise will keep growing due to the the fact that the variances of the two disconnected operators in their vacuum expectation values are independent of each other. Therefore an optimal cutoff R can be found if the lattice size is larger than the correlation length between the operators whereupon the signal-to-noise ratio is improved. However, as is shown in Figure 7 where the g A (DI) from Eq. (29) is plotted as a function of R, no clear plateau shows until at very large R, especially for the light quark case, which is probably because the correlation length is not much smaller than the lattice size. So we cannot find an optimal R in this case. Nevertheless a correlated fit using the following asymptotic form [7] with k and M being free parameters helps in extracting C 3 (∞) properly. The blue bands in the figure show the results of the correlated fit while the green bands show the results of an uncorrelated fit in contrast. The reason why we need this comparison is because the data points of different R are strongly correlated and an uncorrelated fit will underestimate the error very much. To keep the correlation of different R, we cannot do single two-state fits respectively for each R. Instead, a simultaneous two-state fit combining all the R to keep the whole correlation matrix is carried out. The error of the correlated fit, which is not much smaller than the error of the data points with large cutoff, is believed to be a reasonable estimation. In this way, the final statistical uncertainties of the MEs on the 32ID lattice can be reduced by 10% ∼ 40% for different quark masses.

VI. CONNECTED-INSERTION CONTRIBUTION
As for the CI case, we use the improved axial-vector current following our previous work on the 24I and 32I lattice [26] to reduce the discretization errors on the 32ID lattice as well.
For the 24I and 32I lattice, we reanalyze the data with u quark and d quark separately. Twostate fits are also applied to the MEs of A i =ψiγ i γ 5ψ , and the fitting setups are also listed in a similar table (Table IV). We plot the fitting results of the 32ID lattice at the unitary point in Figure 8 as an example. To implement the improvement, we also need to fit for the MEs of three more currents: For these currents, the signal-to-noise ratio is not as good as that for the A i case and no obvious excited-state contribution can be observed; we are only able to make a constant fit combining several separations. A example of D 4 is plotted in Figure 9, note that we drop three points on each of the source and sink sides.
The spatial and temporal components of the improved axial-vector current are defined as A im i = A i + gD i and A im 4 = A 4 + gD 4 , where the factor g is determined by assuming the final g A calculated from the two components of the improved current are identical [26].
Although the results of the currents with derivative are noisy and the constant fit may not be a perfect choice, it is enough for this calculation since the improvement itself is only around 3% or less. Plots of the improvement are shown in Figure 10. For the d quark case,

VII. RENORMALIZATION
The renormalization of the axial-vector current is indispensable for comparing our result with experiment and phenomenology. The scale-independent isovector normalization con-stant Z A (CI) can be calculated by imposing the chiral Ward identity in CI as in Eq. (13) or between the vacuum and a pion state [31]. There is no difference between the u and d quark in this case, as can be seen in the RI/MOM non-perturbative procedure. Hence, Z A (CI) = Z 3 A , the isovector normalization constant. Since we adopt the mass-independent renormalization scheme, it is also the same as the octet renormalization Z 8 as conventionally used in the literature. After checking the AWI in the DI, we concluded in Sec. IV that axial-vector current with the normalization of Z A (CI) satisfies the AWI, thus there is no additional normalization factor for the AWI and Z A (DI) = Z A (CI) = Z A is the only normalization constant as far as tree-level AWI is concerned. Through the chiral Ward identity in the CI, we can determine Z A to a high precision. Since we have calculated MEs of both CI and DI, the disconnected part of vertex functions also needs to be computed and the corresponding renormalization can be obtained by the lattice nonperturbative approach in the RI/MOM scheme [33]. This part contains a scale-dependent DI piece and also mixing effects and is referred to as renormalization, to be distinguished from the normalization, discussed so far, in upholding the AWI at the tree level.

A. Formalisms
The axial-vector coupling has conventionally been classified as the isovector g 3 A = ∆u − ∆d, the octet g 8 A = ∆u + ∆d − 2∆s through the diagonal SU (3) chiral transformation and the singlet g 0 A = ∆u + ∆d + ∆s through the U(1) transformation and their renormalization follows. One can obtain the renormalized ∆u, ∆d and ∆s in term of their unrenormalized counterparts through these flavor-irreducible representations and the details are given in Ref. [34]. On the other hand, the lattice calculations are carried out in terms of flavors and MEs in the CI and DI. It is natural to use them as the basis in renormalization. As we shall see, this has the advantage of preserving the the CI piece which is scale independent and can be compared in different lattice calculations. Moreover, it is physical and can be extracted from the global fitting of the polarized PDF.
In the RI/MOM renormalization scheme, the renormalized quantities are related to the unrenormalized ones through the vertex and the field renormalization. The most general form from the lattice classification is the following where ∆f (f = u, d, s) is the bare axial-vector current matrix element for a particular flavor f and ∆f RI is the corresponding renormalized one in the RI scheme. Σ C or Σ D in the matrix is defined by the following trace indicating the renormalization condition.
where Z q is the quark field renormalization constant, Λ C/D (p) is the connected or discon- the DI MEs receive contributions from both CI and DI. These equations are defined in the quark massless limit so that the RI-MOM is a mass-independent renormalization scheme. In practice, we do calculations at finite quark masses and then extrapolate to the chiral limit.
In principle, Z q can be determined by considering the derivative of the quark propagator with respect to the discretized momenta. However, Z q so determined is known to have large discretization error. We shall use Z A from the chiral Ward identity as an input; therefore, we have Σ C = 1 Z A and Z q is determined via Eq. (32) instead as employed in Ref. [31]. The renormalization constants come from the inverse of the matrix in Eq. (31). The renormalized quark spins in the RI scheme are where In the present calculation, N f = 3.
To compare with experiments, we need to match the results from the above RI scheme to that of MS at 2 GeV. As we see from Sec. II, it entails a two-loop perturbative calculation of the axial-vector current in the RI and MS schemes respectively. On the other hand, it is shown in Eq. (7), this is the same renormalization constant from the one-loop mixing of the topological charge. We carry out the simpler one-loop mixing calculation of the topological charge for the matching factor from the RI scheme at momentum p to the MS scheme at scale µ based on Package-X [35,36] and this matching ratio can be represented as a matrix R m which needs to be where f m = αs Thus, after R m is multiplied to the renormalization matrix in Eq. (34), the renormalized quark spin in the MS scheme is where the notations of ∆u N (CI) and ∆d N (CI) mean they have normalization only and In practice, Z D,MS A are to be evolved to a given scale such as 2 GeV for each p 2 a 2 in RI and extrapolated to p 2 a 2 = 0. This involves an evolution For the axial-vector current, the anomalous dimensions are It is shown in [34] that, at two-loop order, the evolution of the flavor-singlet renormalization factor is given by where Z A + 3Z D,MS A (µ) is the renormalization constant for the flavor-singlet case which we will show later, and the relevant constants are β 0 = 1 The the evolution of α s at two-loop level is given in [37], where W −1 is the lower branch of the Lambert function and Λ is set to be the PDG value The final results of the renormalized u/d quark spin can be decomposed into the CI part and the DI part for each flavor where the connected insertion part is scale independent. We should caution that this is true for the axial-vector case due to the chiral Ward identity. This is not true in general, such as for the case of the scalar and the energy-momentum tensor matrix elements where the CI parts are also scale dependent. On the other hand, the disconnected insertion parts depend on the MS scale of µ where Σ = ∆u + ∆d + ∆s = ∆u(CI) + ∆d(CI) + (∆u + ∆u + ∆s)(DI).
This decomposition of the quark spin in terms of flavor and CI and DI is common for all renormalization schemes, and is not limited to the RI or MS scheme. When the CI and DI components are added together from Eq. (37) to get the total matrix elements, one arrives at a simpler expression where f = u, d, s. In terms of the flavor irreducible representations, they are Eq. (47) can be derived by starting the renormalization from the combined CI and DI matrix elements so that Eq. (37) becomes Similarly, Eq. (47) for the renormalized quark spin for each flavor can be derived from the the basis of flavors irreducible representations g 3 A , g 8 A and g 0 This has been worked out in [34] with the same results.
It is not surprising that one arrives at the same renormalized results in Eq. (47) [38,39].
It is found that there is a connected-sea parton which is in the connected-insertion of the current-current correlator in addition to the disconnected-sea partons in the corresponding disconnected-insertion. The former is responsible for the Gottfried sum-rule violation [38].
These two sea partons have not been separated in the global fittings so far. However, it is demonstrated, in one example, that by combining the strange parton distribution from the semi-inclusive DIS experiment of HERMES, the global fitting result ofū(x) +d(x) and the lattice calculation of the ratio of x s x u(DI) , one can separate the connected-sea from the disconnected-sea distribution of the u and d partons [40]. It is shown that in the operatorproduct-expansion, it is the moments of the combined connected-sea and valence parton distributions that correspond to the local matrix elements of the CI in the lattice calculation [39]. The parton evolution equations with separate connected-and disconected-sea parton is formulated [41]. Provided that future global fitting take this separation into account when fitting experiments at different Q 2 , one can obtain the moments of the valence and connected-sea to extract ∆u(CI) and ∆d(CI) and other moments from the unpolarized and polarized partons and compare directly with the lattice calculation of moments.

B. Numerical results of the renormalization
The results of Z A on the 24I and 32I lattice have been obtained in our previous study [31] to be 1.111(6) and 1.086(2) at the massless limit for both valence and sea quarks. The The corresponding chiral extrapolation is shown in the right panel.
Z A on the 32ID lattice is calculated in this study using the same strategy:   22a, 38a and 15a for the 24I, 32I and 32ID lattices respectively. The improvement of the signal-to-noise ratio is ∼ 50% or less. The criterion of choosing the cutoff is based on the χ 2 of the linear fit with respect to a 2 p 2 which is described in detail in Ref. [42].

VIII. GLOBAL FITTING AND RESULTS
Having the bare MEs and the renormalization constants we obtained in the above sections, we can now carry out the global fitting to push our results to the physical pion point, the continuum limit and the infinite volume limit. The functional form used is where m 2 π,v means the valence pion mass, m 2 π,s means the sea pion mass and L is the size of the lattice. We have two m 2 π terms in the fitting since we are using partially-quenched valence quark masses. We use the form like m 2 π,v − m 2 π,p where m 2 π,p is the physical pion mass in order to let c 0 = g phy A to be the value in the physical limit. However, not all the coefficients in the fitting function have statistical significance during the fit, meaning that the lattice data has no constraint on the corresponding term, or in other words, the effect of the corresponding term is weak enough to be ignored with the current statistical uncertainty. The fitting results of the coefficients of the DI case are listed in Table V term. The χ 2 /d.o.f. stays basically the same or even gets smaller which is another proof that the c 3 term is negligible. The reason why we choose to exclude this term is because this term will enlarge the final uncertainty unreasonably due to the overfit of the data. All the other four terms are kept in our final fit of the DI case. For the CI case, the situation is similar. However, we need to further exclude the c 1 term to have reasonable control of the uncertainty potentially because we use the improved axial-vector current and the finite lattice spacing effects are small. The fitting results of the coefficients of the CI case are listed in Table VI.    The final results of global fitting are shown in Figure 15 and Figure 16 [4] are at Q 2 = 10 GeV 2 and the integration range over the momentum fraction is from 10 −3 to 1. The COMPASS results [5] are at scale Q 2 = 3 GeV 2 . All the lattice results are calculated in the MS scheme at µ = 2 GeV. Since the evolution of ∆Σ   if the divergence of the axial-vector current is inserted as an operator between nucleon states. This is an important check indicating that the lattice artifacts are under control.
For the disconnected-insertion part, the CDER technique is used for the 32ID lattice when constructing 3-point functions and the statistical error can be reduced by 10% ∼ 40%.
The DI contributions to the light and strange quark ∆l(DI) and ∆s(DI) are determined to be −0.073(13)(15) and -0.035 (8)(7), respectively. For the connected-insertion part, we use the improved axial-vector current aiming to reduce the finite lattice spacing effects. The results of the CI contribution to u and d quarks ∆u(CI) and ∆d(CI) are 0.919(13)(28) and −0.337(10)(10) respectively. As we mentioned in Sec. VII, they are scale independent due to the chiral Ward identity and can be compared to other lattice calculations. They can be extracted from deep inelastic scattering, provided the connected-sea and disconnected-sea partons are separated in the global fit [41]. Nonperturbative renormalization is carried out so the reported results are all in the MS scheme at 2 GeV scale. The numerical results are collected in Table VII; the total intrinsic quark spin contribution is ∆Σ = 0.401 (25)(37), which is consistent with the recent global fitting results of experimental data [3][4][5]. The isovector g 3 A = 1.256 (16)(30) with 3% combined statistical and systematic error is within one sigma of that of the experimental value at 1.2723 (23).
When checking the axial Ward identity, we find that the effects of the excited states are crucial to understand the violation of the extended Goldberger-Treiman relation and even a two-term fit cannot always extract the MEs unbiasedly, so our estimations of the systematic uncertainties are relatively large. Our results can be further improved by carrying out the same calculation at the physical point directly and by using larger source-sink separations to reduce the excited-states contamination.