Isovector charges of the nucleon from 2+1-flavor QCD with clover fermions

Isovector charges of the nucleon from 2+1-flavor QCD with clover fermions

Abstract

We present high-statistics estimates of the isovector charges of the nucleon from four 2+1-flavor ensembles generated using Wilson-clover fermions with stout smearing and tree-level 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 cost-effectively achieve measurements on each ensemble. Using these data, the analysis of two-point correlation functions is extended to include four states in the fits, and of three-point functions to three states. Control over excited-state 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 three-point functions at multiple values of the source-sink separation. The results for all three charges, , and , are in good agreement with calculations done using the clover-on-HISQ lattice formulation with similar values of the lattice parameters.

Nucleon matrix elements, lattice QCD, isovector charges
pacs:
11.15.Ha, 12.38.Gc
12

Nucleon Matrix Elements (NME) Collaboration

I Introduction

This work presents high-statistics estimates of isovector charges of the nucleon, , and , on four ensembles of (2+1)-flavor lattice QCD using clover-Wilson fermions and a stout smeared tree-level tadpole-improved Symanzik gauge action Edwards et al. (2016). With increased precision, we demonstrate control over excited-state 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 helicity-flip 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 clover-on-HISQ lattice formulation Bhattacharya et al. (2016) and the current clover-on-clover 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 Drell-Yan and semi-inclusive 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 neutron-proton mass difference in QCD by combining the estimates of with the difference of light quarks masses  MeV obtained from independent lattice QCD calculations González-Alonso 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+1-flavor clover-on-HISQ lattice formulation have been reported in Bhattacharya et al. (2016). Here we present results from a similar high statistics study using the clover-on-clover 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 well-developed 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 all-mode-averaging (AMA) error-reduction 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 excited-state contamination, comparing the variational method and the 2-state fit to data at multiple source-sink separations , was presented in Yoon et al. (2016) using the ensemble.3 In this work, we extend the 2-state fit results presented there by doing the calculation at an additional value of the lattice spacing () and at a lighter pion mass on two ensembles with different volumes ( and ).

The high statistics data allow us to perform a first analysis including up to four states in fits to the two-point correlators and three states in fits to the three-point functions. To obtain results for the charges in the limit that the source-sink separation , we generate data at 4–5 values of on each ensemble. Using these data we perform a detailed comparison of results obtained using 2-state versus 3-state fits. Our final estimates of the charges are from 3-state fits.

The renormalization constants of the various quark bilinear operators are calculated on three ensembles , , in the RI-sMOM scheme and then converted to the scheme at  GeV using 2-loop matching and 3-loop 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+1-flavor clover-on-HISQ 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. Two-state fits to two- and three-point functions to extract the unrenormalized charges are presented in Sec. III along with a discussion of our understanding of, and control over, excited-state contamination. In Sec. III.2, we extend the analysis to include up to four states in fits to two-point functions and three states in three-point correlation functions. The calculation of the renormalization constants in the RI-sMOM 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+1-flavor clover-on-HISQ 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+1-flavor lattice QCD analyzed in this work were generated by the JLab/W&M collaboration Edwards et al. (2016) using clover Wilson fermions and a tree-level tadpole-improved 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 tree-level clover coefficient used is very close to the non-perturbative 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 Wilson-flow 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.4 Similar but smaller differences in obtained using different observables are expected for the other three ensembles. Also, note that the ensemble labeled here was labeled in Ref. Yoon et al. (2016). In this paper, we use estimates of and primarily to label the ensembles and for comparing against previous results with similar lattice parameters in Sec. VI. For this reason, we postpone a more detailed study of scale setting on these ensembles to future works.

The strange quark mass is first tuned in the 3-flavor 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 light-quark 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 three-point 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 excited-state 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
Table 1: Parameters of the 2+1 flavor lattices generated by the JLab/W&M collaboration Edwards et al. (2016) using clover-Wilson fermions and a tree-level tadpole-improved Symanzik gauge action. The lattice spacing is obtained using the Wilson-flow scale and is the dominant source of error in estimates of . The bare quark masses are defined as . Note that the ensemble labeled here was labeled in Ref. Yoon et al. (2016).
ID Method Analysis Smearing Parameters
C1: a127m285 AMA 2-state 1000 4000 128,000
C2: a094m280 (R1) AMA 2-state 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 2-state 1005 3015 96,480
C6: a094m280 (R5) AMA 2-state 1005 3015 96,480
C7: a091m170 AMA 2-state 629 2516 80,512
C8: a091m170L AMA 2-state 467 2335 74,720
Table 2: Description of the ensembles and the lattice parameters used in the analyses. Results from the four runs, R1–R4, on the ensemble were presented in Ref. Yoon et al. (2016). We have extended the statistics in runs R1 and R4 and added R5 to further understand the dependence of the estimates on the smearing size, the efficacy of the variational method and the 2-state fit to data at multiple source-sink separations . The smearing parameters are defined in the text. AMA indicates that the bias in the low-precision measurements (labeled LP) was corrected using high-precision measurements as described in Eq. (10). VAR indicates that the full matrix of correlation functions with smearing sizes listed was calculated and a variational analysis performed to extract the ground state eigenvector as described in Sec. II.3. Analysis of data with multiple to obtain the estimate is carried out using Eqs. (7) and (8).

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 non-relativistic 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 two-point and three-point 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 zero-momentum 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 three-point 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 three-point contributions, is used. Note that, at zero-momentum, the defined in Eq. (6) becomes zero unless , , and .

The two- and three-point correlation functions defined in Eq. (2) are constructed using quark propagators obtained by inverting the clover Dirac matrix with gauge-invariant Gaussian smeared sources. These smeared sources are generated by applying to a unit point source. Here is the three-dimensional 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 three-point 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 two-point functions and from the first two excited states for the three-point functions. The analysis, using Eqs. (7) and (8), is called a “2-state fit” or “3-state fit” or “4-state fit” depending on the number of intermediate states included. The 2-state analysis (keeping one excited state) requires extracting seven parameters (, , , , , and ) from fits to the two- and three-point functions. The 3-state 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 3-state 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 excited-state 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 excited-state 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 three-point function with operator insertion at zero momentum. In the 2-state fits discussed in Sec. III.1, we first estimate the four parameters, , , and from the two-point 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 three-point 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 source-sink 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 three-point data are done within the same single elimination jackknife process to estimate the errors. The same procedure is then followed in the 3-state analysis described in Sec. III.2.

In this study, we demonstrate that stable estimates for the masses, mass-gaps 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 excited-state contamination by constructing the two- and three-point 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 two-point 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 off-diagonal 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 three-point 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 two-point 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 2-state 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 all-mode-averaging (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 twenty-five 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 low-precision (LP) evaluation of the quark propagator is carried out. The resulting LP estimates for two- and three-point functions from these sources are, a priori, biased due to the low-precision inversion of the Dirac matrix. To remove this bias, we selected an additional source point on each of the time slices from which a high-precision (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 three-point function data are constructed first using the HP and the LP correlators as

(10)

where and are the two- and three-point 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 three-point 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 low-accuracy stopping criterion and the HP criterion to the analogous . We have compared the AMA and LP estimates for both the two- and three-point 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 source-sink 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 4-state fits to the two-point function data and then use these as inputs in the extraction of matrix elements from fits to the three-point data. Both of these fits, to two- and three-point 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 three-point 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 3-state fits are consistent with those from the 2-state fits, and the error estimates are comparable. Our final quoted estimates are from 4-state fits to the two-point data and a 3-state fit to the three-point 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.

Figure 1: Data for the the effective mass of the pion (left) and the nucleon (nucleon) obtained using 2-state fits to the zero-momentum two-point correlation functions. The ensemble ID and smearing size are specified in the labels. The 2-state fits to the nucleon are made with our “best” value of described in the text.
Figure 2: Comparison of the effective mass of the pion (left) and the nucleon (right) obtained using 2-state fits to the zero-momentum two-point correlation functions on the ensemble for three different smearings, , and, . The 2-state fits to the nucleon are made with our “best” value of described in the text.
Figure 3: Two- (left) and three-state (right) fits to from the 6 simulations on the 4 ensembles as described in the text. (Left) Data not included in the fits based on are shown in grey. (Right) Lines showing the fits are limited to points fit.
Figure 4: Two- (left) and three-state (right) fits to from the 6 simulations on the 4 ensembles as described in the text. (Left) Data not included in the fits based on are shown in grey. (Right) Lines showing the fits are limited to points fit.
Figure 5: Two- (left) and three-state (right) fits to from the 6 simulations on the 4 ensembles as described in the text. (Left) Data not included in the fits based on are shown in grey. (Right) Lines showing the fits are limited to points fit.
Figure 6: Two- (left) and three-state (right) fits to from the 6 simulations on the 4 ensembles as described in the text. (Left) Data not included in the fits based on are shown in grey. (Right) Lines showing the fits are limited to points fit.
Figure 7: Comparison of the 2-state fits to the data for for the three calculations on the ensemble with different smearing parameter . Each row shows results for the (left) (middle) and (right) calculations.

Iii Fits and Excited-state contamination

To understand and control excited-state contamination we present analyses using 2-, 3- and 4-state fits to the two-point functions and 2- and 3-state fits to the three-point 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 two-point 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 two-point (three-point) functions. In each case, the methodology employed in the analysis is the same except that when using three and higher state fits to two-point functions we introduce non-uniform priors.

In each fit, to understand and quantify the excited-state contamination, there are three parameters that we optimize: (i) the starting time slice used in fits to the two-point data; (ii) the number of time slices , adjacent to the source and sink, skipped in fits to the three-point 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 excited-state 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 4-state fits to the two-point functions. Even then, in the case of 4-state fits, only about eight points contribute to determining the six excited-state parameters since the plateau in the effective-mass 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 auto-correlations 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 2-state fits

The selection of the best combination of , and for the quoted results using 2-state fits was carried out as follows:

  • Step 1: Fits to the two-point 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, 2-state fits to the three-point data were performed for the three sets of , labeled A, B and C, in Table 3. Since the pattern of excited-state 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 2-state results, we use as it removes the most points adjacent to the source and sink that have the largest excited-state 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
,
,
Table 3: The three fits, defined by the values of used in the 2-state fit, investigated to quantify the stability of the estimate.
ID Type Fit Range d.o.f.
3–20 0.6206(20) 1.099(45) 3.51(7)e-08 2.19(11)e-08 0.624(28) 1.26
4–20 0.6193(27) 1.048(71) 3.46(10)e-08 1.98(23)e-08 0.572(55) 1.31
5–20 0.6181(31) 0.980(85) 3.40(13)e-08 1.66(29)e-08 0.487(77) 1.36
4–20 0.4721(25) 0.851(28) 2.86(9)e-08 3.41(13)e-08 1.195(41) 1.33
5–20 0.4676(37) 0.776(39) 2.67(15)e-08 2.92(16)e-08 1.096(57) 0.96
6–20 0.4643(51) 0.724(50) 2.51(22)e-08 2.62(19)e-08 1.042(90) 0.91
4–20 0.4711(30) 0.864(65) 5.55(21)e-10 3.91(40)e-10 0.705(59) 1.08
5–20 0.4670(48) 0.739(77) 5.19(40)e-10 3.02(29)e-10 0.583(66) 0.9
6–20 0.4654(65) 0.70(10) 5.04(58)e-10 2.86(37)e-10 0.568(92) 0.97
3–20 0.4672(25) 0.925(47) 4.52(13)e-12 3.83(24)e-12 0.847(45) 0.94
4–20 0.4635(37) 0.809(65) 4.29(21)e-12 3.02(28)e-12 0.705(58) 0.73
5–20 0.4642(40) 0.83(10) 4.33(24)e-12 3.21(75)e-12 0.74(15) 0.78
3–22 0.4209(24) 0.859(29) 4.67(12)e-10 4.90(14)e-10 1.050(27) 1.42
4–22 0.4180(30) 0.802(37) 4.50(16)e-10 4.41(22)e-10 0.981(43) 1.25
5–22 0.4183(32) 0.808(55) 4.51(18)e-10 4.49(55)e-10 1.00(10) 1.34
4–22 0.4252(19) 0.895(44) 4.87(10)e-10 4.80(44)e-10 0.985(76) 1.75
5–22 0.4210(36) 0.755(72) 4.59(24)e-10 3.35(40)e-10 0.730(62) 1.46
6–22 0.410(11) 0.584(79) 3.73(85)e-10 2.82(50)e-10 0.76(30) 1.1
Table 4: Estimates of the nucleon masses and and the amplitudes and extracted from fits to the two-point correlation functions data using the 2-state ansatz given in (7). For the ensemble, we give the estimates from the three runs with different smearing parameters described in Table 2. The notation labels a nucleon correlation function with source and sink constructed using smearing parameter . We give estimates for three different fit ranges and , expressed in lattice units, as described in the text.
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)
Table 5: Estimates of the three matrix elements , , for the three isovector operators obtained using the 2-state fit to the three-point correlators with our “best” choices of , and .
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)
Table 6: Estimates of the ratios of the unrenormalized isovector charges with our “best” choices of , and in the 2-state fits.

The quality of the data for the two-point 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 ground-state pion mass can be extracted using 1-state 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 auto-correlation function was calculated using two quantities that have reasonable estimates on each configuration: (i) the pion two-point correlator at and (ii) the three-point 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 auto-correlations. 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 auto-correlation time reliably. For the other ensembles, the auto-correlation 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 auto-correlations, are underestimates.

To exhibit the dependence of the two-point 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 effective-mass data (and in the raw two-point 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 excited-state 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 three-point function data and the size of excited-state contamination, we plot an “effective” charge,

(12)

i.e., the ratio of the three-point function to the n-state fit that describes the two-point function. This ratio converges to as the time separations and become large provided the fit to the two-point function, , gives the ground state. Our methodology for taking into account excited-state 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 2-state fits to the ratio of the three- to two-point 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 3-state fits, discussed in Sec. III.2, to facilitate comparison.

The magnitude of the excited-state 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 excited-state 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 excited-state 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 2-state fit gives equally reliable estimates in spite of differences in the excited-state 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 excited-state effect is different in the four charges , we do not uniformly use the same set of values of in our final 2-state 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 excited-state effect in each of the charges and the efficacy of the 2-state 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 excited-state 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 excited-state 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 2-state 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 excited-state 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 excited-state 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 4-states in fits to the two-point function data and 3-states in fits to the three-point functions to evaluate the stability of estimates obtained using 2-state fits.

iii.2 Analysis using 3-state fits

In this section, we investigate the stability of the estimates from 2-state fits by increasing the number of states kept in the fits to the two- and three-point function data. The additional features introduced into the analysis, over and above those discussed in Sec III.1 for the 2-state fits, are:

  • The two-point function data were analyzed using 3- and 4-state fits. In fits with more than two states, the excited state masses and amplitudes are, in many cases, ill-determined. 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).

  • With more fit parameters, the values of and were reduced to increase the number of data points included in the fits to both two- and three-point functions. The results with the final choices of these parameters for various fits are listed in Tables 7 and 8.

  • Data with and 14 were used for all four charges in fits to the three-point functions on the ensemble and with and 16 for the three  fm ensembles.

The priors for the 3- and 4-state fits to the two-point 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 non-trivial priors were needed or used for determining and .

  • Results for and obtained from 2-state fits without priors were used to guide the selection of their priors in 3-state 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