Isovector charges of the nucleon from 2+1flavor QCD with clover fermions
Abstract
We present highstatistics estimates of the isovector charges of the nucleon from four 2+1flavor ensembles generated using Wilsonclover fermions with stout smearing and treelevel tadpole improved Symanzik gauge action at lattice spacings and fm and with and 170 MeV. The truncated solver method with bias correction and the coherent source sequential propagator construction are used to costeffectively achieve measurements on each ensemble. Using these data, the analysis of twopoint correlation functions is extended to include four states in the fits, and of threepoint functions to three states. Control over excitedstate contamination in the calculation of the nucleon mass, the mass gaps between excited states, and in the matrix elements is demonstrated by the consistency of estimates using this multistate analysis of the spectral decomposition of the correlation functions and from simulations of the threepoint functions at multiple values of the sourcesink separation. The results for all three charges, , and , are in good agreement with calculations done using the cloveronHISQ lattice formulation with similar values of the lattice parameters.
pacs:
11.15.Ha, 12.38.GcNucleon Matrix Elements (NME) Collaboration
I Introduction
This work presents highstatistics estimates of isovector charges of the nucleon, , and , on four ensembles of (2+1)flavor lattice QCD using cloverWilson fermions and a stout smeared treelevel tadpoleimproved Symanzik gauge action Edwards et al. (2016). With increased precision, we demonstrate control over excitedstate contamination using a multistate analysis of the spectral decomposition of the correlation functions.
Nucleon charges play an important role in the analysis of standard model (SM) and beyond the standard model (BSM) physics. The nucleon axial charge is an important parameter that encapsulates the strength of weak interactions of nucleons. The ratio, , is best determined from the experimental measurement of neutron beta decay using polarized ultracold neutrons by the UCNA Collaboration, Mendenhall et al. (2013), and by PERKEO II, Mund et al. (2013). In the SM, up to second order corrections in isospin breaking Ademollo and Gatto (1964); Donoghue and Wyler (1990) as a result of the conservation of the vector current. Since is so well measured, it serves to benchmark lattice QCD calculations and our goal is to provide estimates with total uncertainty.
The isovector charges and , combined with the helicityflip parameters and extracted from the measurements of the neutron decay distribution, probe novel scalar and tensor interactions at the TeV scale Bhattacharya et al. (2012). Assuming that and are measured at the precision level Alarcon et al. (2007); Wilburn et al. (2009); Pocanic et al. (2009), one requires that and be calculated with a precision of 10%–15% Bhattacharya et al. (2012). This precision has recently been reached using the cloveronHISQ lattice formulation Bhattacharya et al. (2016) and the current cloveronclover analysis is a necessary independent check using a unitary lattice QCD formulation. The tensor charge is also given by the zeroth moment of the transversity distributions that are measured in many experiments including DrellYan and semiinclusive deep inelastic scattering (SIDIS). Accurate calculations of the contributions of the up and down quarks to the tensor charges will continue to help elucidate the structure of the nucleon in terms of quarks and gluons and provide a benchmark against which phenomenological estimates utilizing a new generation of experiments at Jefferson Lab (JLab) can be compared Dudek et al. (2012). We also use the conserved vector current relation to determine the neutronproton mass difference in QCD by combining the estimates of with the difference of light quarks masses MeV obtained from independent lattice QCD calculations GonzálezAlonso and Martin Camalich (2014); Bhattacharya et al. (2016).
Most extensions of the Standard Model designed to explain nature at the TeV scale have new sources of violation, and the neutron electric dipole moment (EDM) is a very sensitive probe of these. Planned experiments aim to reduce the current bound on the neutron EDM of cm Baker et al. (2006) to around cm. To put stringent constraints on many BSM theories, one requires that matrix elements of novel violating interactions, of which the quark EDM is one, are calculated with the required precision. The contributions of the quark EDM to the neutron EDM Bhattacharya et al. (2015); Pospelov and Ritz (2005) are given by the flavor diagonal tensor charges. Precise results for the connected contributions to these charges from 2+1+1flavor cloveronHISQ lattice formulation have been reported in Bhattacharya et al. (2016). Here we present results from a similar high statistics study using the cloveronclover formulation. The needed contributions of disconnected diagrams are being done in a separate study Gambhir et al. (2016).
The methodology for the lattice QCD calculations of the matrix elements of quark bilinear operators within the nucleon state is welldeveloped Lin (2012); Syritsyn (2014); Green (2016); Constantinou (2014); Bhattacharya et al. (2015, 2016). Our goal is to first calculate the charges as functions of the lattice spacing , the quark mass characterized by the pion mass and the lattice size expressed in dimensionless units of . After renormalization of these lattice estimates, physical results will be obtained by taking the continuum limit (), the physical pion mass limit (set by MeV) and the infinite volume limit () using a combined fit in these three variables Bhattacharya et al. (2016, 2015). Here, we present results for four ensembles at lattice spacings and fm with and MeV. These ensembles are labeled , , , and and described in Table 1.
In this work, we demonstrate that precise estimates for matrix
elements within nucleon states can be obtained by combining the
allmodeaveraging (AMA) errorreduction
technique Bali et al. (2010); Blum et al. (2013) (Sec. II.4) and the
coherent source sequential propagator
method Bratt et al. (2010); Yoon et al. (2016). A detailed analysis of
excitedstate contamination, comparing the variational method
and the 2state fit to data at multiple sourcesink separations
, was presented in Yoon et al. (2016) using the
ensemble.
The high statistics data allow us to perform a first analysis including up to four states in fits to the twopoint correlators and three states in fits to the threepoint functions. To obtain results for the charges in the limit that the sourcesink separation , we generate data at 4–5 values of on each ensemble. Using these data we perform a detailed comparison of results obtained using 2state versus 3state fits. Our final estimates of the charges are from 3state fits.
The renormalization constants of the various quark bilinear operators are calculated on three ensembles , , in the RIsMOM scheme and then converted to the scheme at GeV using 2loop matching and 3loop running. Our final estimates of the renormalized isovector charges of the nucleon in the scheme at GeV are given in Table 15. Results for the connected part of the flavor diagonal charges are given in Table 16.
Estimates of all three isovector charges, , and and of the flavor diagonal charges are in very good agreement with similar high precision calculations done using a 2+1+1flavor cloveronHISQ formulation Bhattacharya et al. (2016). Our estimates of obtained with heavy and quark masses corresponding to and MeV are within of the experimental result, , from neutron beta decay Mendenhall et al. (2013); Mund et al. (2013).
This paper is organized as follows. In Sec. II, we describe the parameters of the gauge ensembles analyzed and the various methods used to obtain high precision results. Twostate fits to two and threepoint functions to extract the unrenormalized charges are presented in Sec. III along with a discussion of our understanding of, and control over, excitedstate contamination. In Sec. III.2, we extend the analysis to include up to four states in fits to twopoint functions and three states in threepoint correlation functions. The calculation of the renormalization constants in the RIsMOM scheme is discussed in Sec. IV. Our final renormalized estimates are given in Sec. V and compared to previous results obtained using a 2+1+1flavor cloveronHISQ lattice formulation but with similar statistics and lattice parameters Bhattacharya et al. (2016, 2015) in Sec. VI. We end with conclusions in Sec. VII.
Ii Lattice Methodology
A detailed description of the lattice methodology and our approach has been presented in Refs Yoon et al. (2016); Bhattacharya et al. (2015, 2016). Here we reproduce the discussion necessary to establish the notation and give details relevant to the analysis and the results.
The four ensembles of 2+1flavor lattice QCD analyzed in this work were generated by the JLab/W&M collaboration Edwards et al. (2016) using clover Wilson fermions and a treelevel tadpoleimproved Symanzik gauge action. The update is carried out using the rational hybrid Monte Carlo (RHMC) algorithm Duane et al. (1987). One iteration of stout smearing with the weight for the staples is used in the fermion action. A consequence of the stout smearing is that the tadpole corrected treelevel clover coefficient used is very close to the nonperturbative value determined, a posteriori, using the Schrödinger functional method Edwards et al. (2016).
The lattice parameters of the four ensembles are summarized in
Table 1. Estimates for the lattice spacing were
obtained using the Wilsonflow scale following the prescription
given in Ref. Borsanyi et al. (2012). We caution the reader that an
alternate estimate of for the ensemble we label with
and MeV, has been quoted in
Ref. Leskovec et al. (2016) to be (and
MeV since is unchanged) using the
mass difference. Thus, different
estimates of from this coarse ensemble may vary by
depending on the observable used to set them.
The strange quark mass is first tuned in the 3flavor theory by requiring the quantity to equal its physical value 0.1678. We choose this quantity since it is independent of the light quark masses to lowest order in chiral perturbation theory, i.e., the ratio depends only on the value of the strange quark mass Lin et al. (2009) and can, therefore, be tuned in the SU(3) symmetric limit. The resulting value of is then kept fixed as the lightquark masses are decreased in the (2+1)flavor theory to their physical values. Further details involving the generation of these gauge configurations will be presented in a separate publication Edwards et al. (2016).
The parameters used in the calculations of the two and threepoint functions carried out on the four ensembles are given in Table 2. Analyzed configurations are separated by 6, alternating 4 and 6, 4, and 4 trajectories on the , , and ensembles, respectively. Note that the ensemble has been analyzed in 5 different ways labeled as runs R1–R5 in Table 2 to understand and control excitedstate contamination in nucleon matrix elements. As discussed in Ref. Yoon et al. (2016), and analyzed further here, the five calculations give consistent results. Relevant details of the lattice methods used and of the analyses carried out are summarized next.
Ensemble ID  (fm)  (MeV)  

a127m285  0.127(2)  285(3)  6.1  1.24930971  0.2850  0.2450  5.85  
a094m280  0.094(1)  278(3)  6.3  1.20536588  0.2390  0.2050  4.11  
a091m170  0.091(1)  166(2)  6.3  1.20536588  0.2416  0.2050  3.7  
a091m170L  0.091(1)  172(6)  6.3  1.20536588  0.2416  0.2050  5.08 
ID  Method  Analysis  Smearing Parameters  
C1: a127m285  AMA  2state  1000  4000  128,000  
C2: a094m280 (R1)  AMA  2state  1005  3015  96,480  
C3: a094m280 (R2)  LP  VAR  , ,  443  0  42,528  
C4: a094m280 (R3)  AMA  VAR  , ,  443  1329  42,528  
C5: a094m280 (R4)  AMA  2state  1005  3015  96,480  
C6: a094m280 (R5)  AMA  2state  1005  3015  96,480  
C7: a091m170  AMA  2state  629  2516  80,512  
C8: a091m170L  AMA  2state  467  2335  74,720 
ii.1 Correlation Functions
The interpolating operator used to createannihilate the nucleon state is taken to be
(1) 
with color indices denoted by , charge conjugation matrix , and and the two different flavors of light quarks. The nonrelativistic projection is inserted to improve the signal, with the plus and minus sign applied to the forward and backward propagation in Euclidean time, respectively.
The twopoint and threepoint nucleon correlation functions at zero momentum are defined as
(2) 
where and are the spinor indices. In writing Eq. (2), the source time slice has been translated to time ; the sink time slice, at which a zeromomentum nucleon insertion is made using the sequential source method Bratt et al. (2010); Yoon et al. (2016), is at for forward propagation; and is the time slice at which the bilinear operator is inserted. The Dirac matrix is , , and for scalar (S), vector (V), axial (A) and tensor (T) operators, respectively, with . In this work, subscripts and on gamma matrices run over , with .
The nucleon charges are defined as
(3) 
where the normalization of the spinors in Euclidean space is
(4) 
To analyze the data, we construct the projected two and threepoint correlation functions
(5)  
(6) 
The projection operator is used to project on to the positive parity contribution for the nucleon propagating in the forward direction. For the connected threepoint contributions, is used. Note that, at zeromomentum, the defined in Eq. (6) becomes zero unless , , and .
The two and threepoint correlation functions defined in Eq. (2) are constructed using quark propagators obtained by inverting the clover Dirac matrix with gaugeinvariant Gaussian smeared sources. These smeared sources are generated by applying to a unit point source. Here is the threedimensional Laplacian operator and and are smearing parameters that are given in Table 2 for each calculation. Throughout this paper, the notation will be used to denote a calculation with source smearing and sink smearing . Variations of the parameter over a large range does not impact any of the results Yoon et al. (2016), and it is dropped from further discussions since our choice lies within this range. Before constructing the smeared sources, the spatial gauge links on the source time slice are smoothed by 20 hits of the stout smearing procedure with weight . A more detailed discussion of the efficacy of source smearing used in this study is given in Ref. Yoon et al. (2016).
ii.2 Behavior of the Correlation Functions
Our goal is to extract the matrix elements of various bilinear quark operators between ground state nucleons. The lattice operator , given in Eq. (1), couples to the nucleon, all its excitations and multiparticle states with the same quantum numbers. The correlation functions, therefore, get contributions from all these intermediate states. Using spectral decomposition, the behavior of two and threepoint functions is given by the expansion:
(7)  
(8) 
where we have shown all contributions from the ground state and the first three excited states , and with masses , and to the twopoint functions and from the first two excited states for the threepoint functions. The analysis, using Eqs. (7) and (8), is called a “2state fit” or “3state fit” or “4state fit” depending on the number of intermediate states included. The 2state analysis (keeping one excited state) requires extracting seven parameters (, , , , , and ) from fits to the two and threepoint functions. The 3state analysis introduces five additional parameters: , , , and . These simultaneous fits to data at multiple values of provide estimates of the charges in the limit . Throughout this paper, values of and are in lattice units unless explicitly stated.
Nine of the twelve parameters in the 3state analysis—the three masses , and and the six matrix elements —are independent of the details of the interpolating operator. Our goal is to obtain their values by removing the discretization errors and the higher excitedstate contaminations. The amplitudes depend on the choice of the interpolating nucleon operator and/or the smearing parameters used to generate the smeared sources. It is evident from Eqs. (7) and (8) that the ratio of the amplitudes, , is the quantity to minimize in order to reduce excitedstate contamination as it determines the relative size of the overlap of the nucleon operator with the excited states. A detailed analysis of how it can be reduced by tuning the smearing size and a comparison of the efficacy with a variational analysis (run R2 and R3), described in Sec. II.3, was presented in Ref. Yoon et al. (2016) using the ensemble. We present an update on the comparison using renormalized charges obtained from fits with the full covariance matrix in Sec. VI.
We extract the charges and ( and ) from the real (imaginary) part of the threepoint function with operator insertion at zero momentum. In the 2state fits discussed in Sec. III.1, we first estimate the four parameters, , , and from the twopoint function data. The results for all four ensembles and for three selected fit ranges investigated are collected in Table 4. These are then used as inputs in the extraction of matrix elements from fits to the threepoint data. For the insertion of each operator , extraction of the three matrix elements , and is done by making one overall fit to the data versus the operator insertion time and the various sourcesink separations using Eq. (8). In these fits, we neglect the data on time slices on either end adjacent to the source and the sink for each to reduce the contributions of the neglected higher excited states. Fits to both the two and threepoint data are done within the same single elimination jackknife process to estimate the errors. The same procedure is then followed in the 3state analysis described in Sec. III.2.
In this study, we demonstrate that stable estimates for the masses, massgaps and the charges can be obtained with measurements. The errors in the other matrix elements are large, nevertheless certain qualitative features can be established.
ii.3 The variational Method
One can also reduce excitedstate contamination by constructing the two and threepoint correlation functions incorporating a variational analysis (see Dragos et al. (2015, 2016) and references therein for previous use of the variational method for calculating nucleon matrix elements). To implement this method on the ensemble, we constructed correlation functions using quark propagators with three different smearing sizes that are summarized under runs R2 and R3 in Table 2 but with a single fm. The twopoint correlation function for the nucleon at any given time separation is then a matrix made up of correlation functions with source smearing and sink smearing . The best overlap with the ground state is given by the eigenvector corresponding to the largest eigenvalue obtained from the generalized eigenvalue relation Fox et al. (1982):
(9) 
where are the eigenvectors with eigenvalues . The matrix should be symmetric up to statistical fluctuations, so we symmetrize it by averaging the offdiagonal matrix elements. Our final analysis for the calculation of the ground state eigenvector was done with and as discussed in Yoon et al. (2016).
Similarly, in our variational analysis, the threepoint function data , from which various charges are extracted, are matrices . The ground state estimate is obtained by projecting these matrices using the ground state vector estimated from the twopoint variational analysis, i.e., . This projected correlation function is expected to have smaller excited state contamination compared to the correlation function with single smearings. Since the variational correlation function has been calculated at a single , we analyze it using only 2state fits. Note that the contribution of the matrix element cannot be isolated from from fits to data with a single value of . Consequently, our variational estimates of the charges, collected in Table 13, include the contamination from the transition unlike results from multistate fits to data at a number of values of .
ii.4 The AMA Method for High Statistics
The high statistics calculation on the four ensembles was carried out using the allmodeaveraging (AMA) technique Bali et al. (2010); Blum et al. (2013) and the coherent sequential source method Bratt et al. (2010); Yoon et al. (2016). To implement these methods, we choose at random four time slices separated by on each configuration of the and ensembles and on five time slices separated by twentyfive time slices on the ensemble. On the lattices, we choose three time slices separated by time slices on each configuration and staggered these by time slices between successive configurations to reduce correlations.
On each of these time slices, we choose randomly selected source locations from which lowprecision (LP) evaluation of the quark propagator is carried out. The resulting LP estimates for two and threepoint functions from these sources are, a priori, biased due to the lowprecision inversion of the Dirac matrix. To remove this bias, we selected an additional source point on each of the time slices from which a highprecision (HP) and LP measurement of the correlation functions was carried out. The total number of measurements made on each ensemble are given in Table 2.
On each configuration, the bias corrected two and threepoint function data are constructed first using the HP and the LP correlators as
(10) 
where and are the two and threepoint correlation functions calculated in LP and HP, respectively. Correlators from the two kinds of source positions and , are assumed to be translated to a common point when defining Eq. (10). The bias in the LP estimate (first term) is corrected by the second term provided the LP approximation is covariant under lattice translations, which is true for the two and threepoint functions. The contribution to the overall error by the second term is small provided the HP and LP calculations from the same source point are correlated. To estimate errors, the measurements on each configuration are first averaged and then the single elimination Jackknife procedure is carried out over these configuration averages.
We used the adaptive multigrid algorithm for inverting the Dirac matrix Babich et al. (2010) and set the lowaccuracy stopping criterion and the HP criterion to the analogous . We have compared the AMA and LP estimates for both the two and threepoint correlation functions themselves and for the fit parameters , , and the matrix elements . In all cases, we find the difference between the AMA and LP estimates is a tiny fraction (few percent) of the error in either measurement Yoon et al. (2016). In short, based on all the calculations we have carried out, possible bias in the LP measurements with as the stopping criteria in the adaptive multigrid inverter is much smaller than the statistical errors estimated from measurements.
ii.5 Statistics
The total number of LP and HP measurements and the values of sourcesink separations analyzed are given in Table 2. Our statistical tests show that correlations between measurements are reduced by choosing the source points randomly within and between configurations Yoon et al. (2016). Also, using the coherent source method for constructing the sequential propagators from the sink points to reduce computational cost does not significantly increase the errors Bratt et al. (2010); Yoon et al. (2016).
On all the ensembles, we first estimate the masses and the amplitudes using the 2, 3 or 4state fits to the twopoint function data and then use these as inputs in the extraction of matrix elements from fits to the threepoint data. Both of these fits, to two and threepoint data, are done within the same Jackknife process to take into account correlations in the estimation of errors. We performed both correlated and uncorrelated fits to the nucleon two and threepoint function data. In all cases, correlated and uncorrelated fits gave overlapping estimates. For the final results, we use fits minimizing the correlated .
We find that the central values from the 3state fits are consistent with those from the 2state fits, and the error estimates are comparable. Our final quoted estimates are from 4state fits to the twopoint data and a 3state fit to the threepoint data with the matrix element set to zero as discussed in Sec. III.2. Our overall conclusion is that to obtain the isovector charges and with 1% uncertainty (or 2% uncertainty at the physical pion mass and after extrapolation to the continuum limit) will require measurements on each ensemble. Approximately five times larger statistics are needed to extract with the same precision.




























Iii Fits and Excitedstate contamination
To understand and control excitedstate contamination we present analyses using 2, 3 and 4state fits to the twopoint functions and 2 and 3state fits to the threepoint functions. We find that to get reliable estimates of the masses and amplitudes for the first states we need to include states in the fit to the twopoint function. For this reason, we analyze the correlation functions with the following combinations , , , and where the first (second) value is the number of states included in fits to the twopoint (threepoint) functions. In each case, the methodology employed in the analysis is the same except that when using three and higher state fits to twopoint functions we introduce nonuniform priors.
In each fit, to understand and quantify the excitedstate contamination, there are three parameters that we optimize: (i) the starting time slice used in fits to the twopoint data; (ii) the number of time slices , adjacent to the source and sink, skipped in fits to the threepoint functions; and (iii) the values of at which data are collected used in the fits. The final values of these parameters, chosen on the basis of the and the stability of the fits, represent a compromise between statistical precision and reducing excitedstate contamination. In general, we reduce the value of and and enlarge the number of values included when increasing the number of states in the fit ansatz. For example, we set in 4state fits to the twopoint functions. Even then, in the case of 4state fits, only about eight points contribute to determining the six excitedstate parameters since the plateau in the effectivemass plot starts at as shown in Figs. 1 and 2.
Our focus is on obtaining estimates for the charges in the limit for the six calculations labeled as in Table 2. Two overall caveats that will be made explicit at appropriate places are: the statistics in the case of the ensemble () are insufficient as the autocorrelations between configurations are large. Similarly, the errors in the data for are much larger than for , consequently the fits used to extract are much less stable. In all cases, the fits and the error analysis presented here are based on using the full covariance matrix.
iii.1 Analysis using 2state fits
The selection of the best combination of , and for the quoted results using 2state fits was carried out as follows:

Step 1: Fits to the twopoint correlators using the full covariance matrix for different values of were made. The best value minimizing the correlated was determined to be for the six calculations.

Step 2: Using these values of , we determined the four parameters , , and . Then, 2state fits to the threepoint data were performed for the three sets of , labeled A, B and C, in Table 3. Since the pattern of excitedstate contamination is different in the various charges, the best set was determined separately for each charge. Our final results are based on Fit B for the scalar and vector charges, and Fit C for the axial and tensor charges as discussed below.

Step 3: Repeat Step 2 for the three cases and for the ensembles, respectively. All three cases gave overlapping results. For our quoted 2state results, we use as it removes the most points adjacent to the source and sink that have the largest excitedstate contamination.

Step 4: To quantify the dependence on , we repeat Steps 1–3 for .
Results for , , and for and are given in Table 4. Results for the three matrix elements for our best choice of , and are given in Table 5. Estimates of the ratios of the unrenormalized charges, , are given in Table 6.
Fit A  Fit B  Fit C  

,  
,  
ID  Type  Fit Range  d.o.f.  

3–20  0.6206(20)  1.099(45)  3.51(7)e08  2.19(11)e08  0.624(28)  1.26  
4–20  0.6193(27)  1.048(71)  3.46(10)e08  1.98(23)e08  0.572(55)  1.31  
5–20  0.6181(31)  0.980(85)  3.40(13)e08  1.66(29)e08  0.487(77)  1.36  
4–20  0.4721(25)  0.851(28)  2.86(9)e08  3.41(13)e08  1.195(41)  1.33  
5–20  0.4676(37)  0.776(39)  2.67(15)e08  2.92(16)e08  1.096(57)  0.96  
6–20  0.4643(51)  0.724(50)  2.51(22)e08  2.62(19)e08  1.042(90)  0.91  
4–20  0.4711(30)  0.864(65)  5.55(21)e10  3.91(40)e10  0.705(59)  1.08  
5–20  0.4670(48)  0.739(77)  5.19(40)e10  3.02(29)e10  0.583(66)  0.9  
6–20  0.4654(65)  0.70(10)  5.04(58)e10  2.86(37)e10  0.568(92)  0.97  
3–20  0.4672(25)  0.925(47)  4.52(13)e12  3.83(24)e12  0.847(45)  0.94  
4–20  0.4635(37)  0.809(65)  4.29(21)e12  3.02(28)e12  0.705(58)  0.73  
5–20  0.4642(40)  0.83(10)  4.33(24)e12  3.21(75)e12  0.74(15)  0.78  
3–22  0.4209(24)  0.859(29)  4.67(12)e10  4.90(14)e10  1.050(27)  1.42  
4–22  0.4180(30)  0.802(37)  4.50(16)e10  4.41(22)e10  0.981(43)  1.25  
5–22  0.4183(32)  0.808(55)  4.51(18)e10  4.49(55)e10  1.00(10)  1.34  
4–22  0.4252(19)  0.895(44)  4.87(10)e10  4.80(44)e10  0.985(76)  1.75  
5–22  0.4210(36)  0.755(72)  4.59(24)e10  3.35(40)e10  0.730(62)  1.46  
6–22  0.410(11)  0.584(79)  3.73(85)e10  2.82(50)e10  0.76(30)  1.1 
Axial  Scalar  Tensor  

ID  Type  
a127m285  1.423(14)  0.179(21)  0.9(2.4)  1.07(4)  0.35(4)  0.6(1.1)  1.166(13)  0.182(16)  0.2(1.2)  
a094m280  1.349(19)  0.130(20)  0.6(0.7)  1.18(6)  0.42(5)  1.0(0.8)  1.071(17)  0.157(15)  0.7(4)  
a094m280  1.384(28)  0.111(36)  0.3(1.3)  1.23(12)  0.52(12)  1.4(1.7)  1.085(30)  0.221(36)  0.0(0.8)  
a094m280  1.372(25)  0.026(39)  0.4(2.5)  1.28(9)  0.42(9)  0.6(3.3)  1.067(25)  0.276(28)  0.6(1.6)  
a091m170  1.388(23)  0.133(33)  2.1(2.6)  1.17(11)  0.48(7)  0.1(3.9)  1.091(20)  0.154(22)  0.2(1.7)  
a091m170L  1.401(20)  0.118(26)  1.0(2.4)  1.15(8)  0.44(8)  1.4(2.2)  1.067(25)  0.235(23)  0.5(8) 
ID  Type  

a127m285  1.125(11)  0.848(27)  0.922(12)  
a094m280  1.130(17)  0.987(50)  0.897(15)  
a094m280  1.154(24)  1.030(95)  0.906(27)  
a094m280  1.143(22)  1.068(76)  0.889(22)  
a091m170  1.146(21)  0.963(84)  0.901(16)  
a091m170L  1.166(18)  0.960(69)  0.888(21) 
The quality of the data for the twopoint functions on the four ensembles is illustrated by plotting the effective mass,
(11) 
for the pion and the nucleon in Fig. 1. As expected, the signal in the pion does not degrade with , whereas that for the proton becomes noisy by , with 1–2 fluctuations in apparent already at . The onset of a plateau indicates that the groundstate pion mass can be extracted using 1state fits to data at . In practice, the ground state mass is largely determined from the region , while the excited state masses and amplitudes are determined from the region . The value of is, therefore, adjusted depending on the number of states included in the fit.
To assess the statistical quality of the data, the autocorrelation function was calculated using two quantities that have reasonable estimates on each configuration: (i) the pion twopoint correlator at and (ii) the threepoint correlation function at the midpoint in for . Autocorrelations increase as or is decreased or is increased. In particular, the data from the ensemble showed significant autocorrelations. In this case, the 467 configurations consist of four streams with roughly 170, 100, 100 and 100 configurations. These are too few to even determine the autocorrelation time reliably. For the other ensembles, the autocorrelation function falls to by 1–2 configurations, and binning the data by a factor of two did not change the Jackknife error estimates. Our overall conclusion is that much larger statistics are needed to get reliable error estimates on the ensemble and it is very likely that the quoted errors for this ensemble, evaluated without taking into account autocorrelations, are underestimates.
To exhibit the dependence of the twopoint correlation functions on the smearing size given by , we show a comparison of the effective mass for the pion and the nucleon for the , and correlators on the ensemble in Fig. 2. We find that the errors in the effectivemass data (and in the raw twopoint functions) increase with the smearing size for both the pion and the proton. The onset of the plateau in both states, however, occurs at earlier times with larger . Thus, the relative reduction in the excitedstate contamination in the correlation functions with larger smearing has to be balanced against the increase in statistical noise. Based on these trends on the ensemble, our compromise choice for the three ensembles at fm is and for the ensemble it is . In physical units, this choice corresponds to setting the size of the smearing parameter fm.
Our final estimates for the four parameters , , , and the ratio for three values of are given in Table 4. In addition to minimizing , we required stability in the value of under the variation as criteria for choosing our best . With the selected , we find that is also consistent within except on the ensemble which, as stated above, requires much higher statistics.
To illustrate the threepoint function data and the size of excitedstate contamination, we plot an “effective” charge,
(12) 
i.e., the ratio of the threepoint function to the nstate fit that describes the twopoint function. This ratio converges to as the time separations and become large provided the fit to the twopoint function, , gives the ground state. Our methodology for taking into account excitedstate contamination and obtaining estimates of the charges from data with in the limited range 1–1.5 fm is described next.
The data and the 2state fits to the ratio of the three to twopoint functions using our best choice of , and are shown in Figs. 3, 4, 5 and 6 for the four isovector charges. In the right panels of these figures, we show the 3state fits, discussed in Sec. III.2, to facilitate comparison.
The magnitude of the excitedstate contamination as a function of and the smearing parameter is different for the four charges. The dependence on is exhibited in Fig. 7 for the , and calculations on the ensemble. In , the magnitude of the excitedstate contamination, measured as the difference between the data at the central values of for (about 1 fm) and the estimate, is about 10%, 5% and 3% for the , and calculations, respectively. The pattern in is similar, however, the reduction in the contamination with is smaller. For , the overall variation with and between the three estimates is . The vector charge shows insignificant excitedstate contamination and no detectable dependence on . On the other hand, the errors in individual data points increase with the smearing for all four charges.
We use the data with the three values of the smearing on the ensemble to test whether the 2state fit gives equally reliable estimates in spite of differences in the excitedstate contamination. We find that the three estimates are consistent within for all four charges as shown in Fig. 7. However, because the magnitude of the excitedstate effect is different in the four charges , we do not uniformly use the same set of values of in our final 2state fit, but tune them for each case.
Based on the data shown in these figures and on the results of the fits with our best choices of the input fit parameters given in Tables 5 and 6, we evaluate below the excitedstate effect in each of the charges and the efficacy of the 2state fit in providing estimates.

The data for the axial charge shown in Fig. 3 converges to the value from below and at the central values of show up to 10% variation with due to excitedstate contamination. We, therefore, use Fit C based on data with the larger values of . In all but the case, the data at lies above the result of the fit. The errors in these data, however, are large.
^{5} Thus, to confirm the estimate, requires additional high precision data for .The data with correlators on the ensemble show the least excitedstate effect: the estimates at the central values of show only a tiny increase with and their error bands overlap. As a result, the matrix elements and , given in Table 5, are poorly determined.

The data for the scalar charge have larger uncertainty so we choose Fit B to include all the data except that with and on the ensemble. As shown in Fig. 4, the 2state fits again converge from below, the three estimates from the ensemble are consistent within as shown in Fig. 7, and the estimates from the and ensembles overlap.

The data for the tensor charge show small excitedstate contamination and converge to the estimate from above. Fits to the the data for are the most stable and all the fits give consistent estimates. We choose Fit C for the final estimates. These estimates agree with those from Fit B, and the of both fits are also consistent.

The data for the vector charge show little variation with or . The excitedstate contamination is highly suppressed because is associated with a conserved charge at . As a result, statistical fluctuations in the data stand out. Note that the error estimates in from the fits are comparable to the 1–2% variance in individual data points at the largest . Only the data with on the ensemble deviate by about from the result of the fit.
In the next section, we extend the analysis to include up to 4states in fits to the twopoint function data and 3states in fits to the threepoint functions to evaluate the stability of estimates obtained using 2state fits.
iii.2 Analysis using 3state fits
In this section, we investigate the stability of the estimates from 2state fits by increasing the number of states kept in the fits to the two and threepoint function data. The additional features introduced into the analysis, over and above those discussed in Sec III.1 for the 2state fits, are:

The twopoint function data were analyzed using 3 and 4state fits. In fits with more than two states, the excited state masses and amplitudes are, in many cases, illdetermined. The fits were stabilized by carrying out an empirical Bayesian analysis with Gaussian priors for both the mass gaps and the amplitudes of the excited states Lepage et al. (2002); Chen et al. (2004).

Data with and 14 were used for all four charges in fits to the threepoint functions on the ensemble and with and 16 for the three fm ensembles.
The priors for the 3 and 4state fits to the twopoint function data were determined empirically as follows:

The ground state mass and amplitude are very well constrained by the plateau in the effective mass for . Thus, no nontrivial priors were needed or used for determining and .

Results for and obtained from 2state fits without priors were used to guide the selection of their priors in 3state fits. The width was chosen to be large but consistent with the requirement that the bands for and the are positive. These priors did not need any subsequent changes.

Results for