Lattice Calculation of Nucleon Isovector Axial Charge with Improved Currents

Lattice Calculation of Nucleon Isovector Axial Charge with Improved Currents

Jian, Yi-Bo Yang, Keh-Fei,
Andrei, Terrence Draper and Raza Sabbir Sufian (QCD Collaboration)
Department of Physics and Astronomy, University of Kentucky, Lexington, KY 40506, USA Department of Physics, The George Washington University, Washington, DC 20052, USA

We employ dimension-4 operators to improve the local vector and axial-vector currents and calculate the nucleon isovector axial coupling with overlap valence on -flavor Domain Wall Fermion sea. Using the equality of from the spatial and temporal components of the axial-vector current as a normalization condition, we find that is increased by a few percent towards the experimental value. The excited-state contamination has been taken into account with three time separations between the source and sink. The improved axial charges and are obtained on a lattice at pion mass of 330 MeV and a lattice at pion mass 300 MeV are increased by 3.4% and 1.7% from their unimproved values, respectively. We have also used clover fermions on the same DWF configurations and find the same behavior for the local axial charge as that of the overlap fermions.

11.15.Ha, 12.38.Gc, 14.20.Dh

I introduction

Lattice QCD, as a non-perturbative method for solving QCD problems, has achieved enormous success in the field of hadron physics. However, there exist several quantities whose lattice results still do not agree with experimental values. One outstanding example is the nucleon isovector axial charge . The value of this charge, Olive et al. (2014), has been well determined by the neutron -decay experiments (see, for example, Mendenhall et al. (2013); Mund et al. (2013)). However, many lattice simulations (for recent results, please refer to Alexandrou et al. (2013); Bhattacharya et al. (2014); Owen et al. (2013); Junnarkar et al. (2015); Syritsyn (2015); Ohta (2015); Yang et al. (2016a); Yoon et al. (2016)) yield results as much as lower than the experimental value. We note that some of lattice calculations do have results consistent with experiments Horsley et al. (2014); Bali et al. (2015); Abdel-Rehim et al. (2015), due to, for example, the use of large source-sink time separations, which is encouraging, but we think it is essential to have consistent results from different fermion actions at the physical pion mass with systematics such as the continuum and large volume extrapolations under control so that there is consensus from the community. So it is still important to pursue lattice studies on and nucleon structure in general both theoretically and technically to test the whole lattice methodology critically.

Furthermore, decomposing the nucleon spin into quark and gluon constituents, i.e. quark spin, quark orbital angular momentum, gluon spin and gluon angular momentum, has long been an important issue. The axial charge gives the intrinsic quark spin in the nucleon of flavor . So to solve the spin problem via lattice QCD, an accurate from lattice must be demonstrated Liu (2016). However, for the strange quark, there have been a number of lattice calculations Bali et al. (2012); Engelhardt (2012); Babich et al. (2012); Abdel-Rehim et al. (2014); Alexandrou et al. (2016) with their ( is the nucleon spin polarization direction) in the range from to , whose absolute values are several times smaller than those of the global fits of DIS experiments, e.g., de Florian et al. (2009), Nocera et al. (2014) and Adolph et al. (2016), leading to a large discrepancy. So a lot of effort is still needed to close the gap and the accurate and precise calculation of can serve as a benchmark for all the lattice calculations.

The deviations between lattice results and the experimental values of can originate from lattice artifacts and lattice systematic errors, such as the lack of chiral symmetry, finite volume effects, finite lattice spacing effects and excited-state contamination, if we believe that QCD is the correct description of the strong interaction. One issue that escaped scrutiny is the fact that practically all the recent lattice calculations use only local currents instead of improved currents or conserved currents for the chiral fermions such as domain wall fermion or overlap fermion. With an improved fermion action, one needs to adopt the correspondingly improved current to remove the error. And even with the conserved current, one needs improvement as well to obtain the improved conserved current to make the matrix elements calculated to be , especially for the off-forward case Martinelli et al. (1991). As we shall see from the present study, we find that the two values obtained by inserting point current , or are obviously different on each of the two lattice ensembles with and fm respectively: the one coming from , denoted by , is smaller than the one coming from , which is denoted as . Since must be calculated in a moving frame, and its signal-to-noise ratio (SNR) is worse compared to , most previous lattice studies focus on only, so this may be the first observation of this kind of deviation. This deviation, manifesting the asymmetry between spatial and temporal components on the lattice, is due to the finite lattice spacing artifact and could be resolved by using the conserved or improved conserved lattice axial current. And hopefully, using the improved current may also improve the final result at finite lattice spacing, leading to a more accurate result to be compared with the experimental value.

One can build a conserved vector current easily for Wilson-like fermion actions, while constructing a conserved axial vector current is only viable for chiral fermions, such as the overlap fermion in our case Hasenfratz et al. (2002). Our future goal is to employ the conserved axial vector current (or its improved version) to calculate of the nucleon. In the meantime, as an exploratory study, we shall use simple improved currents in this work, to see whether this kind of improvement can result in degenerate and with common coefficients on different valence quark masses, and lead to a value closer to experiments as well.

Since, in the isovector channel, the disconnected insertions are canceled between the two degenerate light quark flavors, we focus on the the connected insertion (CI) only in this work. As a cautious test of our codes and results, we also carry out a test calculation using clover fermions as the valence.

This paper is organized as follows. In Sec. (II), we present the numerical details, including lattice setup, smeared-to-smeared 3-point function construction, and a new low-mode substitution scheme for the nucleon propagator. Sec. (III) introduces the improved currents used in this work. Sec. (IV) contains the numerical results. Then, we show the results of the test calculation on the clover fermions in Sec. (V). A summary is given in Sec. (VI).

Ii Numerical details

ii.1 Lattice Setup

In this work, we use overlap fermions for the valence quark on -flavor Domain-Wall fermion (DWF) sea Aoki et al. (2011) to carry out the calculation. The effective quark propagator of the massive overlap fermion is defined as the inverse of the operator  Chiu (1999); Liu (2005), where is exactly chiral, i.e. Chiu and Zenkin (1999), and can be expressed in terms of the overlap Dirac operator as


where is the matrix sign function and is the Wilson Dirac operator with =0.2 (corresponding to parameter ). The RBC/UKQCD DWF gauge configurations used are from the (24I) and (32I) ensembles Aoki et al. (2011). The parameters of the ensembles are listed in Table (1). We use 5 different quark masses with corresponding ranging from MeV to MeV on each of the two ensembles to study the pion mass dependence of and . We compute 4 different source-sink separations, , , and on 24I and three separations , , on 32I, to handle the excited-state contamination. The largest separations are fm for 24I and fm for 32I, respectively.

Symbol a (fm) (MeV)
24I  0.1105(3) 0.005 330 203
32I  0.0828(3) 0.004 300 309
Table 1: The parameters for the RBC/UKQCD configurations: label, spatial/temporal size, lattice spacing Blum et al. (2016), the degenerate light sea quark mass, the corresponding pion mass and the number of configurations used in this work.

To better control the SNR, we use two 12-12-12 (16-16-16 for 32I) noise grid sources with Gaussian smearing at and and four 2-2-2 (1-1-1 for 32I) grid sources at which are 8, 10, 11 and 12 (12, 14 and 15 for 32I) time-slices away from the source positions with block smearing. The notations such as 12-12-12 denote the intervals of the grid in the 3 spatial directions (please see reference Yang et al. (2016a) for more details). The block smearing function can be defined as


where is the smearing size, is a path of gauge links between spatial coordinate and , for instance,


and is the component of the coordinate , is a product of gauge links in the direction,


We average 6 types of paths to avoid bias on the order of the gauge links in the path. Numerically, we have an algorithm to speed up the smearing process with such a definition, making the cost of the smearing proportional to , rather than Yang et al. (2016b). The size of is tuned very close to the size of the Gaussian smearing at the source, so the source and sink are approximatively symmetric. Different from the block case, the cost of Gaussian smearing is proportional to the iteration times Gong et al. (2013), which is around 10 times slower than the block smearing for the smearing size used in this paper, so we choose to employ the block smearing at the sink. More technical details regarding the calculation of the overlap operator, eigenmode deflation in the inversion of the fermion matrix, low-mode substitution (LMS) of random grid source with mixed momenta, and the stochastic sandwich method with LMS for constructing 3-point functions can be found in references Li et al. (2010); Gong et al. (2013); Yang et al. (2016a). In the next subsection, we will discuss some improved numerical techniques of LMS.

ii.2 Stochastic sandwich method with smeared sink

In our previous paper Yang et al. (2016a), we use point-sink when we calculate 3-point functions by means of the stochastic sandwich method with LMS, but at the source position, we employ a Gaussian smeared source, leading to an asymmetry between the source and sink. In that case, the excited-state contamination will be larger on the sink side and that asymmetry also causes difficulties when we try to do the fit.

Figure 1: The quark diagram of a nucleon correlator from position to with a connected insertion at . The product of the quark propagators in the shadowed region constitutes the current inserted propagator . The propagator from the current position to the sink is decomposed into its low- and high-mode contributions ( and respectively).

The 3-point function of the nucleon can be expressed formally as Yang et al. (2016a),


where is a effective propagator coming from source , going through the inserted current , and ending at the sink (the shadowed region in Fig. 1); is a functional of the other two propagators and (the black and red lines in Fig. 1) without current insertion. When performing LMS on the quark propagator between the sink and the current, we decompose the sink part of into high-mode and low-mode parts,


where is the exact low-mode propagator propagating from to , is the noise-estimated high-mode propagator propagating from to , is the noise vector helping to pick up the starting point of and here is the index of noise vector. If we focus on the sink spatial position only, omit all other coordinate indices, and denote , we can rewrite Eq. (5) and Eq. (6) in a more compact form,


The simple summation on illustrates that it is a point-sink. Since the low-mode propagator is an all-to-all propagator, we can pick up exactly any sink point without noise estimation. This is why LMS can increase the SNR greatly especially for the low-mode dominated cases.

To implement the smeared-sink, we actually want


where denotes the sink grid, is the smearing function from point to , either Gaussian or block. In our practical calculations, we firstly apply smearing on the sink of


Eq. (8) then becomes


Note that the function is symmetric, i.e., . So we can exchange the two summations and rewrite the above equation as


Here is the point in the sink grid. So for each , what we actually do is: first “anti-smearing” on , i.e., , and then complete the final nucleon contraction.

It is well known that the smeared-to-smeared 2-point functions have poor SNR compared to that of the smeared-to-point correlators. However for the smeared-to-smeared 3-point function to 2-point function ratios, we find that the SNR is almost the same as the smeared-to-point ones, which should be attributed to the cancelations of the fluctuations when taking the ratio. The fitting of the ratios are more stable and reliable as expected now for the symmetric source and sink.

ii.3 New LMS contraction scheme

As shown above, the stochastic sandwich method of 3-point function construction with LMS will eventually turn into several 2-point function contractions of effective propagators. In our previous implementation, the nucleon correlation function with LMS is expressed as (copying Eq. (5) of reference Yang et al. (2016a) here for convenience),


where the functional means a normal contraction operation, denotes the high-mode part of the noise grid propagator, is the low-mode propagator, is the noise phase at the position which belongs to grid , and is the full noise grid propagator which is a combination of the high-mode part and the low-mode part. The above equation is actually an expansion by dividing into these two parts: . In this old scheme, we have 4 functional in line 2 and line 3 of Eq. (12). Since each one entails a normal contraction, we need to do 4 times of nucleon contraction for these two lines. While in the last two lines of Eq. (12), functional are within the summation of , and each summation consists of times of contraction operation, where is the number of points in the grid , so there are actually times of nucleon contractions contained in these two lines. Therefore, totally we need times of normal contraction operations for a complete LMS contraction.

A better choice is


Since the order of the grid summation and the functional can be exchanged optionally in the low-high-high case, i.e., , this expression is exactly equivalent to Eq. (12),


However, Eq. (13) needs only times of normal contraction operation: times for the first term and time for the second one. This new LMS contraction scheme reduces considerably the complication of coding, and takes only 1/4 of the computer time as compared to that in Eq. (12).

Iii improved currents

Lattice QCD is a system based on a hypercube lattice and a lattice action is a discretized version of the continuum QCD action. Therefore, a simple local current, which is conserved in the continuum limit, is no longer conserved when the lattice spacing is finite. For the clover fermions, even through the action is improved to , one still needs to have improvement in addition to the conserved current to remove the lattice artifact Martinelli et al. (1991). In this case, the best option to calculate a matrix element on the lattice is to use the so-called improved conserved current.

For the vector case of Wilson-like fermions, the conserved current is the point-split one Doi et al. (2009); Burger et al. (2014), which can be expanded in lattice spacing to have a local current and a derivative term in the next order in :


where . To the lowest order, this is actually an modification to the point current. For the clover fermions which is improved to owing to the use of the clover term , one needs to add an extra term to improve, either the expanded point-split current or the conserved one, to at the same time if one wants the final matrix element calculated to be Martinelli et al. (1991); Heatlie et al. (1991). So in the clover case, both and should be used in addition to the point current to implement a full improvement. For the overlap case, though it is not easy to derive the expansion of the conserved current Hasenfratz et al. (2002), we notice that a chiral lattice Dirac operator can always have a expansion as Niedermayer (1999)


where a clover term with is automatically generated. So in this sense, we can assume that the improvement terms and also apply to the overlap case with point current. Therefore, the general form of the improved vector current in our case can be expressed as


Whereas, for the axial vector case, only chiral fermions, like overlap Hasenfratz et al. (2002) and Domain Wall fermions Boyle (2014), can have a conserved current and both are harder to calculate numerically. For the time being, we shall assume that the common improvement terms and of axial current for general improved actions Luscher et al. (1996) apply to the overlap case as well with free coefficients. Therefore, the general form of improved axial current can be expressed as


and in the above two equations are the improved vector and axial vector currents; and are the normalization (finite renormalization) constant, which gives rise to the continuum-like effective quark propagator ; , , and are the coefficients of the improvement terms. We omit the lattice spacing in the above formulas for simplicity. However, the term does not contribute to the forward matrix element when calculating and the term does not contribute to the unpolarized matrix element for . And for chiral fermion in our case, we have , so the final improved currents used in this work are


The charges corresponding to the improved currents and are marked as and ; similarly, we also have the notations as and . Here and are averaged over values of . In the following, when using a latin letter, e.g., , as the Dirac index, it ranges from to , while greek letters range from to . Specifically, to calculate the connected 3-point functions of the improved currents, we need to carry out computation for each of the following currents


We use superscript or to denote the point currents and the dimension-4 currents with derivative, respectively. We mark the correspoding charges of these currents as , , , , , , and for further convenience. Our goal is to calculate the improved isovector axial charge and to eliminate the deviation between the spatial and temporal parts, that is to say, we will demand, after our improvement, and as our normalization conditions. Using these two equations, we can solve for the coefficients and as


After we have these two coefficients, we need the normalization constant to calculate the final improved charges. We have determined the axial normalization constant before through the chiral Ward identity for the pion Liu et al. (2014). It that case, none of the derivative terms contributes since the pion is at zero momentum, and thus we determine only . In this work, we will use the proton vector charge to check whether it is unity with which is set to equal to the from the pion, such that we can know whether the same applies to the proton case.

To calculate the above charges, we need to calculate the forward nucleon matrix element . This can be obtained via the 3-point function to 2-point function ratios ,


where is the smeared proton interpolating field, and is the same except for using instead. In the vector case, is the non-polarized projector of the nucleon spin; in the axial vector case, is the polarized projector. When is large enough, there should exist a plateau, which is denoted as . is a product of the desired matrix element and a kinematic factor . To extract the matrix elements, we need to compute these factors first. For example if , is some gamma matrix,


and and are the nucleon mass and energy. The matrix element can then be expressed as .

current factor
Table 2: Kinematic factors in the 3-point function to 2-point function ratios. The last two columns show the results in two special cases, which are used in this work.

All the factors are listed in Table (2). For the axial vector case, we choose the polarization index the same as the Dirac index for currents and . The index of the dependence for currents and comes from the polarized projector. It is found that all improvement currents have the same structures as the local ones, and the kinematic factors of the improvement currents just have one additional “” multiplied to the factors of the corresponding local currents. For and , we need to carry out the calculations in a moving frame because the factors are proportional to the nucleon momentum. For and , we can do the calculation in the rest frame of the nucleon. This is the reason why they are separated into two parts. For the axial case, it is the other way around.

Iv Results

iv.1 Vector Case

Figure 2: Example of 2-state fits in the vector case on the 32I lattice at the unitary point. Points with error bars are from lattice results and the curves are from the 2-state fits. For , it is too flat to use 2-state fit, so we use a constant fit instead. The black lines are the final fit values and the gray bands indicate the fit errors.

As we mentioned above in section (III), we will compute to check whether the same determined from the pion also applies to the proton case. We use a 2-state fit to handle the ratios


where is the energy difference between the first excited-state and the ground state, is the source-sink separation and is the time slice with current insertion. Constant is the desired matrix element, coefficients and are related to the transition between the ground state and the first excited-state, while accounts for the excited-state to excited-state contribution. So in different channels their significance can be different; we need to pick out the significant terms in order to get a stable fit. As an example, the fitted results on 32I at the unitary point can be found in Fig. 2. For , it is too flat to apply the 2-state fit, so we use a constant fit instead. For , there is no plateau on the plot, which is presumably due to the fact that we used rather than in our 3-point function contraction code. However, a 2-state fit can handle this case very well, the lattice data points almost all lay on the fit curves. (We also tested in the clover case using sink-sequential methods and the final results are not affected; the figures will be shown in the next section.) For , the final error band is much larger than the error of the lattice data, which is because the fitted is small and the excited-state contributions cannot be accurately fixed by the data. We also tried to use a constant fit for and the results are consistent with the 2-state fit, but the is around 2 which is unacceptable. For , since it is noisier than the other 3 cases, a 2-state fit gives a result consistent with the constant fit; we choose to use the 2-state fit as shown in Fig. 2.

Figure 3: Results of the vector cases as a function of pion mass squared for 32I. and are rescaled by a factor of 20 for clarity.

The pion mass dependence is shown in Fig. 3. keeps constant when changes and at each is consistent with within errors; decreases a little as the pion mass decreases but the values are all consistent with . Given this situation, the equation is already satisfied within errors so we cannot determine a unique factor for the improvement term. In other words, there is no obvious need for this kind of improvement in the vector channel, since no evident deviation between the spatial and temporal parts is observed.

The results from 24I are similar. Using the bare value of , we find that within errors with determined from pion using Ward identity Liu et al. (2014), which means the pion also applies to the proton case. So we will use the normalization constant provided in the above reference, which is, for 24I and for 32I.

iv.2 Axial Vector Case

Figure 4: 2-state fits of 32I at the unitary point for the axial case. (Points with error bars are from lattice data and the curves are from the 2-state fits). For different cases we keep different terms in the 2-state fit. The gray bands indicate the final fit results.
Figure 5: The same as in Fig. 4 but for lattice 24I at the unitary point.

For the axial vector, we also try to use the standard 2-state fit to handle the excited-state contamination. However, for , the lattice data points show no obvious curvature at either the source side or the sink side on both lattices 24I and 32I within errors (the fit examples at the unitary points are shown in Fig. 4 and 5, respectively), which means the transition terms are heavily suppressed. If we still force a full 2-state fit on the data, the final uncertainty will be uncontrollable since the data have no constraint on the and terms. But the data at different source-sink separations have some discrepancy, so we need to retain the term in addition to the constant . For , we utilize and both the and terms to fit the data. In the practical fit procedure, we combine and into a joint fit with shared parameter , where a wide prior is used in some channels to ensure stable fit results. For , we merely use a constant fit on the combination of three time separations. For , although there is no evident plateau (for the same reason as for ), a 2-state fit fits the data well. We can see that is larger than by a factor of .

Figure 6: Results of and as a function of squared pion mass, for both 24I and 32I.

The valence pion mass dependence of the axial vector charges can be found in Fig. 6, the values are all normalized. has hardly any pion mass dependence, which is consistent with other calculations Alexandrou et al. (2013); Bhattacharya et al. (2014), and the results from 24I and 32I are consistent with each other. At the unitary point, we have and , which are lower than the experimental value, and the deviation is around 7 percent. This is also consistent with other calculations Alexandrou et al. (2013); Owen et al. (2013); Bhattacharya et al. (2016) for pion mass MeV. On the other hand, is smaller than by , and it decreases with decreasing pion mass. At the unitary point, the gap between and is around . This deviation gets smaller with increasing pion mass. In fact, in the strange quark region, this difference vanishes. This behavior is evidence that the deviation has a physical origin and is not a mere statistical fluctuation.

Figure 7: Results of the improvement currents and as a function of pion mass squared, for 24I and 32I.

It is very interesting that the currents with derivatives (see Fig. 7) manifest exactly the opposite behavior to those of the point currents. This time has strong dependence and the value increases as pion mass decreases, while for , the central values are tiny and almost constant within errors. These opposite behaviors between the dimension-3 local currents and the dimension-4 derivative currents is exactly the pattern that is needed to implement the improvement and makes larger and thus closer to the experimental value. Using Eq. (21), we can indeed obtain the coefficient at each pion mass and, furthermore, the value of is constant ( for 32I and for 24I) at different pion masses within errors, which is what we expect (since mass dependence entails a higher dimensional correction) and is shown in Fig. 8. The errors of come from bootstrap resampling: we do 2-state fit (or constant fit) and solve Eq. (21) in each bootstrap sample.

Figure 8: The coefficient solved on 32I. The errors of come from bootstrap.

Since the coefficient is positive and the charges related to the currents with derivative operators are also positive, the final improved axial charges will be larger than the original values. The improved results are presented in Fig. 9. At the unitary point, the improved axial charges are and , which represents a 3.4% and 1.7% increment towards the experimental value, respectively. For other pion masses, there are improvements as well.

Figure 9: The final improved axial charge values of 24I and 32I compared with the unimproved .

Although the improvement is still not enough to fill the gap of between lattice and experiments, it at least is in the correct direction. The results here are all from lattices with MeV; for the results around the physical pion mass, this several percent of improvement may be more significant. We believe that the long-standing deviation of from the experimental value is probably not due to one single source. It may be a combined effect of finite lattice volume, finite lattice spacing, and excited-state contamination. Although the improvement currents contributes only two or three percent, they should be taken into account at finite cut off before approaching the continuum, in addition to the finite normalization factors and .

V a test of clover case

As a benchmark test of our results and the implementation of our low-mode substitution sandwich method, we carry out a test calculation using clover fermions as valence on the same DWF sea on lattice 24I. We also want to know whether the discrepancy observed above is only an artifact of overlap fermion or a more common phenomenon. We use the standard sink-sequential method for constructing the 3-point functions without any low-mode substitution. Therefore, the forward-backward derivative operator can be implemented easily; we can therefore check the difference between using and . Moreover, we also try to figure out whether our improvement is still valid for the and flavor separately (connected insertion part only). In this test, we choose the clover coefficient and mass parameter to make the pion mass to be around MeV. We only consider one source-sink separation