{'fraction_non_alphanumeric': 0.03590277655199195, 'fraction_numerical': 0.017068822480715937, 'mean_word_length': 4.365470214982479, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nucleon Charges and Electromagnetic Form Factors from 2+1+1-Flavor Lattice QCD

We present lattice-QCD results on the nucleon isovector axial, scalar and tensor charges, the isovector electromagnetic Dirac and Pauli form factors, and the connected parts of the isoscalar charges. The calculations have been done using two ensembles of HISQ lattices generated by the MILC Collaboration with 2+1+1 dynamical flavors at a lattice spacing of 0.12 fm and with light-quark masses corresponding to pions with masses 310 and 220 MeV. We perform a systematic study including excited-state degrees of freedom and examine the dependence of the extracted nucleon matrix elements on source-sink separation. This study demonstrates with high-statistics data that including excited-state contributions and generating data at multiple separations is necessary to remove contamination that would otherwise lead to systematic error. We also determine the renormalization constants of the associated quark bilinear operators in the RI-sMOM scheme and make comparisons of our renormalized results with previous dynamical-lattice calculations.


I. INTRODUCTION
Precision measurements of the properties of neutrons provide an opportunity to determine a fundamental parameter of the Standard Model (SM), to compute many parameters of effective theories in nuclear physics and to explore novel physics at the TeV scale. Decays of neutrons provide one of the best determinations of the CKM quark-mixing parameter V ud and the nucleon axial charge g A . Calculations of nucleon matrix elements of bilinear quark operators yield a variety of interesting physical quantities, including the nucleon charges such as g A,S,T , the nucleon σ-term and strangeness, and the electromagnetic form factors. In addition to these quantities, isoscalar bilinear matrix elements can be related to measurements of the neutron electric dipole moment (nEDM), which will shed light on CP violation in and beyond the Standard Model (BSM). In this paper we report a lattice-QCD (LQCD) calculation of the isovector charges g A , g S and g T , the connected-diagram part of the isoscalar charges and the isovector electric and magnetic radii, extracted from the electromagnetic form factors. The axial charge g A is a key parameter in nuclear physics, while estimates of g S and g T are needed to constrain possible scalar and tensor interactions at the TeV scale [1].
The nucleon isovector axial charge g A is a key parameter in the description of nucleon structure, since it encapsulates the interaction of the charged axial current with the nucleon. For example, it affects the rate of protonproton fusion, which is the first step in the thermonuclear reaction chains that power low-mass hydrogen-burning stars like the Sun; and g A is central to the extraction and phenomenology of the CKM matrix element V ud . Presently, g A is best known from the experimental measurement of neutron beta decay using polarized ultracold neutrons by the UCNA collaboration [2], which dom-inates the PDG average with uncertainty at the 0.2% level [3]. (See the figure in Sec. VI for the collected experimental g A measurements used in PDG2012, the recently updated UCNA number 1.2756 (30) [4] and a recent result from Perkeo II 1.2761 +14 −17 [5].) The nucleon scalar and tensor charges have not, until very recently, been well studied since the contributions of effective scalar and tensor interactions in the SM are small, at the 10 −3 level, and still below the current experimental limits. In many extensions to the SM (e.g. supersymmetry), novel scalar and/or tensor interactions arise via exchanges in either the s or t channels or through loop effects, and these can also contribute at the 10 −3 level to neutron decay. Since the SM contributions are known to high precision, 10 −5 , one has the opportunity to measure these scalar and tensor couplings in neutron beta decay experiments with sufficient precision to isolate the BSM from the SM contributions. The current status of the theory and experimental measurements is summarized in Ref. [1], in which it was shown that, assuming that planned experiments achieve 10 −3 sensitivity, to constrain new physics at the TeV scale estimates of g S and g T with 10-15% accuracy are needed. A number of such experiments are being developed at Los Alamos (UCNB [6] and UCNb [7]) and Oak Ridge National Lab (Nab [8][9][10]), and they aim to make measurements in the coming years.
Lattice QCD is the best method to obtain g S and g T with the desired precision. The calculation steps, conceptually and procedurally, are the same as for the vector and axial charges. The main hurdles are the statistical error in g S , which are a factor of about five larger than those in g A and g T .
Lattice calculations of the properties of nucleons are more difficult than those for mesons for a number of reasons. First, the lowest excited state, the Roper N (1440), lies relatively close to the nucleon mass. Therefore, a creation operator that has substantial overlap with excited states will require large values of Euclidean time t to reduce contributions of excited states to correlation functions. Alternatively, the nucleon operators must be carefully crafted to maximize overlap with the ground state. Second, there is a signal-to-noise problem. In Euclidean space, the signal-to-noise ratio in nucleon correlation functions scales like exp((−M N + 3M π /2)t). Since M N > 1.5M π , high statistics are needed to obtain a signal at large t, where excited-state contributions have become negligible, in order to extract the properties of the ground state. With finite computer resources we must, therefore, trade off between statistical and systematic error due to contamination from excited states. As a result, much larger computational resources are needed for precision studies of nucleons on the lattice than would be required for mesons. Third, heavy-baryon chiral perturbation theory (HBXPT) is more difficult because of the nearby ∆ resonance [11]. As a result, there are several different HBXPT expansions, and it is not a priori clear how to tell which will converge well for a given observable. Consequently, although chiral perturbation theory has been a very important tool for understanding the dependence of meson observables on the light-quark masses in a lattice calculation, it is not as useful for baryons. Thus, extrapolations of lattice calculations from unphysically large light-quark (u and d) masses to the physical point have significant uncertainty. This problem is being addressed gradually as simulations are being performed closer to physical values of the light-quark masses. Lastly, finite-volume effects are observed to be larger in baryon correlation functions as compared to mesons. Historically, studies of finite-volume effects have been carried out using mesonic observables, for which it has been established empirically that finite-volume effects are negligible for M π L > 4. Since the generation of gauge ensembles is expensive and driven by physics goals in the meson sector, only a single ensemble of lattices with M π L > 4 is usually generated at a given lattice spacing. As a result, finite-volume effects are less well understood in the baryon sector. A few studies suggest even larger values of M π L are needed to have this systematic under control [12][13][14].
In this work we present a detailed analysis of excitedstate contamination in matrix elements and results on the three isovector charges g A , g S and g T , corresponding to the matrix elements of isovector quark bilinear operators within the nucleon, and the Dirac and Pauli (or the related electric and magnetic) charge radii extracted from the corresponding electromagnetic form factors. We also present the connected part of the isoscalar charges. All calculations have been done at one lattice spacing, a ≈ 0.12 fm, so no extrapolation to the continuum limit is possible. Ensembles of gauge configurations generated at two values of the light quark mass, corresponding to M π ≈ 310 and 220 MeV, have been analyzed, so we will make some comments on quark-mass dependence. Our main focus will be on understanding three issues: statis-tical errors, excited-state contributions, and nonperturbative calculations of the renormalization constants in the RI-sMOM scheme.
The paper is organized as follows. In Sec. II we describe the gauge ensembles and lattice parameters used in this study. Details of the calculations of the two-and three-point functions are given in Sec. III. Analysis of the statistical signal is presented in Sec. III C and of the excited-state contamination in Sec. IV. We discuss the calculation of the renormalization constants in the RI-sMOM scheme in Sec. V. The results for the charges and comparison with other published estimates are given in Secs. VI and VII. The electromagnetic form factors and charge radii are discussed in Sec. VIII. Our conclusions are given in Sec. IX.

II. LATTICE PARAMETERS AND SETUP
We analyze two ensembles of gauge configurations generated by the MILC collaboration [15] with N f = 2+1+1 flavors of highly improved staggered quarks (HISQ) [16][17][18][19][20][21] as described in Table I. The HISQ action, proposed by HPQCD/UKQCD Collaboration [16,17], has, among existing variations of staggered fermions, at nonzero a the smallest splittings between the staggered "tastes" that become four degenerate flavors in the continuum limit [19,22]; this leads to a significant reduction in the discretization errors associated with the staggered action. The work presented here is the first step in our analysis of HISQ sea-quark ensembles of about 1000 configurations generated at three lattice spacings a ∈ {0.12, 0.09, 0.06} fm, two light-quark masses corresponding to M π ≈ 310 and 220 MeV and the strange and charm-quark masses set to their physical values. (The actual values of the lattice spacings we use in this study, based on the data presented in Ref. [15], are a = 0.120(1), 0.088(1) and 0.058(1) fm.) Our goal is to perform a simultaneous continuum and chiral extrapolation of physical quantities using these six ensembles. In this paper, we focus on demonstrating control over systematic errors associated with our lattice approximations using the two ensembles at a ≈ 0.12 fm.
Staggered-type fermions are notorious for their complications in calculations involving baryons, especially of matrix elements. Therefore, we use clover (O(a)improved Wilson) fermion action in the valence sector for our calculation of nucleon matrix elements. Strictly speaking, such a mixed-action approach with HISQ fermions for sea quarks and clover for valence quarks, results in a non-unitary formulation. One consequence of this mixed-action approach is the possibility of exceptional configurations. These are configurations in which the spectrum of the clover Dirac matrix has nearzero modes. Such configurations would have been suppressed if the lattices had been generated with the same clover action, so their presence is an artifact of using the mixed-action approach. Signatures of such configura-tions, which manifest at sufficiently small quark masses, include: (i) correlation functions calculated on them have anomalously large values, thus biasing the ensemble average, and/or (ii) the calculation of the inverse of the clover Dirac matrix fails to converge due to poor condition number. Our approach to this problem is empirical. Based on the two signatures listed above we have determined that exceptional configurations are absent in the data on M π ∈ {220, 310} MeV MILC ensembles [15] at 0.12 fm, but the same is not the case for the ensemble at M π ≈ 135 MeV. On these we find signatures of exceptional configurations, and therefore do not analyze them. We expect that exceptional configurations will be suppressed at finer lattice spacing, and our ongoing analysis of an ensemble with a = 0.09 fm and M π = 130 MeV shows no exceptional configurations on the full set (883 configurations). Thus, the mixed-action approach can be used at the physical quark mass for lattice spacings of 0.09 fm and smaller.
The mixed-action approach is used for two reasons. First, to calculate matrix elements within baryon states, we need high-statistics analyses on lattices with large volume (M π L > 4). Second, to take the continuum limit and to elucidate the dependence on pion mass requires large ensembles at multiple pion masses and lattice spacings. Because the generation of ensembles requires very large computing resources sustained over many years, few collaborations can meet these requirements. The MILC 2+1+1-flavor ensembles are the only ones that satisfy these requirements that are available to us. Since we are able to avoid exceptional configurations, the mixedaction allows us to test new-physics ideas and computational methods. In our mixed action approach, lattice discretization errors are dominated by our Wilson-clover action, which has not been fully O(a) improved; the the clover term has not been nonperturbatively tuned, but merely set to the tadpole-improved perturbative coefficient. Also, the operators used to calculate the matrix elements and renormalization constants have not been improved. Thus, our discretization errors start at O(a), and we must include a linear term in continuum extrapolation.
We use hypercubic (HYP) smearing [23] of the gauge links before inverting the clover Dirac matrix needed to construct correlation functions [24,25]. Using gauge fields averaged over a hypercube reduces short-distance noise (lattice artifacts) without changing long-distance physics. We observe this improvement in the calculations of two and three-point correlation functions. HYP smearing also modifies the discretization artifacts appearing at high momentum in the calculation of the renormalization constants [26]. We describe our strategy for estimating the associated systematic uncertainty in Sec. V. We find that even using conservative error estimates, since HYP smearing drives the renormalization constants close to the tree-level value (unity), the uncertainty due to renormalization constants in our preliminary study [1] is reduced.
Further details regarding the tuning of the valence clover action to match the HISQ sea-quark action and issues regarding the mixed action are discussed in Refs. [27][28][29][30][31][32]. In Table I, we show the level of agreement between the pion and nucleon masses calculated with the two actions. Similar parameter choices for the same valence and sea-quark actions in the light-quark sector are also used in a study of charmed-hadron physics in Refs. [27,28].

III. LATTICE METHODOLOGY
In this section we describe the lattice calculation of two-and three-point correlation functions. After establishing the notation and methodology in Secs. III A and III B, we discuss the statistical errors in Sec. III C and our understanding and mitigation of excited-state contamination in the extraction of the ground-state matrix elements in Sec. IV.

A. Two-Point Correlators
The correlation functions with the quantum numbers of the spin-1/2 nucleon are constructed using the baryonic interpolating operator where C is the charge-conjugation matrix iγ 4 γ 2 , {a, b, c} are color indices, is the antisymmetric tensor and q 1 and q 2 are one of the two quarks {u, d}. For example, in the case of the proton, we want q 1 = u and q 2 = d. Two-point correlators are derived from these interpolating fields as where p is the baryon momentum and P = 1 4 (1 + γ 4 )(1 + iγ 5 γ 3 ) is the spin projection. A and B label the smearing parameters used for the source and sink operators as discussed below. Eq. 2 can be decomposed in terms of energy eigenstates: where n runs over all energy eigenstates that couple to the operator defined in Eq. 1 with amplitude A n,A for smearing parameter A. The normalization of these states is defined as 0|(χ N ) † |p, s = Au N ( p, s) with the spinors in Euclidean space satisfying To construct correlation functions, we generate valence clover quark propagators with a gauge-invariant β L 3 × T N cfgs Nprops (aMπ)sea (aMπ) val (MπL) val (MπT ) val (aMN )sea (aMN ) val (aMR) val 6.00 24 3  TABLE I. Details of the two ensembles analyzed and lattice parameters used in this study. The subscript "sea" labels the masses of the Goldstone pseudoscalar meson and nucleon calculated using the HISQ on HISQ action [15], while the subscript "val" labels the masses calculated with the valence clover fermions on the HISQ lattices. The "sea" masses have a single statistical uncertainty, while the valence masses include statistical and systematic uncertainty due to fitting-window selection added in quadrature. We also list the spatial (L) and temporal lattice extents (T ) in lattice units, the value of MπL, MπT , the number of configurations analyzed, and the total number of measurements (Nprops) performed on each ensemble. Note that although we call our fitted excited-state nucleon mass MR the Roper mass, it requires study of higher excited states and the volume dependence of the correlator to properly distinguish a true Roper resonance from a scattering state.
Gaussian-smeared source centered at x, the point at which the nucleon operator is defined in Eq. 1. Smearing is done using a fixed number, n KG , of applications of the Klein-Gordon operator with coefficient σ. The ideal smearing for calculating the ground-state zeromomentum nucleon mass and matrix elements leading to g A,S,T is one that minimizes overlap with excited states. An effective-mass plot with the results of one-and two-state fits for the zero-momentum nucleon is given on the top of Fig. 1 with diagonal Gaussian smearing A = B. For Gaussian smearing parameters of {σ, n KG } = {5.5, 70} (using the conventions found in Chroma [33]), we find that the excited-state signal dies out around t = 6 (≈ 0.7 fm), giving enough data points to extract groundstate and first-excited masses and amplitudes as shown.
The bottom of Fig. 1 shows that the extracted groundstate mass agrees between the two-state and one-state fits when one-state fits are constrained to t ≥ 6, and the twostate fit gives a consistent ground-state mass for all fitting windows. Our estimate of the masses of the ground and first-excited states for each ensemble are given in Table I. One challenge in the selection of smearing parameters is the need to simultaneously improve the signal in states at nonzero momentum, as high as say | p| 2 = 5 in units of 2π La , since these are needed to study form factors. We find the signal in our smeared correlator deteriorates significantly compared to smearings with smaller radii. This is not surprising, since high-momentum states are expected to have a smaller overlap with a broadly smeared source. Thus, to improve the quality of the signal at higher momenta, we need to produce multiple correlators with different smearings and in each case explicitly subtract any excited-state contributions to obtain results for the ground state. The choice of smearing in this study was driven by improving results for g A,S,T . We, therefore, used diagonal Gaussian-smeared sources A = B with a single smearing {5.5, 70}, optimized to improve the signal in the zero-momentum two-point nucleon correlator.
To simplify notation, we drop the subscripts defining the smearing, A and B, in all further discussions. The masses of the ground and first-excited state will be labeled as M 0 and M 1 and the corresponding amplitudes with which they couple to the operator defined in Eq. 1 as A 0 and A 1 . These masses and amplitudes are needed as inputs to extract the charges and form factors from three-point correlators. Our final results for these quantities are obtained by applying a fit to the smeared-smeared zero-momentum correlators, keeping only two states in Eq. 3:

B. Three-Point Correlators
To calculate the nucleon matrix elements (such as isovector charges or electromagnetic form factors), we first calculate the matrix element of general form for the isovector vector current, A µ = uγ µ γ 5 d for isovector axial current, etc., and p {i,f } are the initial and final nucleon momenta. Such a matrix element is extracted from an appropriate three-point correlation function after Fourier transforming out the spatial dependence and projecting on the baryonic spin, leaving a time-dependent threepoint correlator of the form where f n,n contains kinematic factors involving the energies E n and amplitudes A n between the creation and annihilation operators and the corresponding states. The latter are obtained from the analysis of the two-point correlators with n and n labeling the different energy states. Z Γ is the operator renormalization constant which is determined nonperturbatively in this work; see Sec. V. The projection T used is T mix = 1 4 (1 + γ 4 )(1 + iγ 5 γ 3 ). In this work we are interested in only the ground-state matrix element with n = n = 0. The parameter of interest in quantifying excited-state contamination, discussed in Sec. IV, is the source-sink separation (t f − t i ). In this study it is varied between 8 and 12 time slices in lattice units, which in physical units corresponds to source-sink separations between about 0.96 and 1.44 fm. By fitting the time dependence of the three-point correlators to the  Table I. The data are for momenta p = 0 (squares) and | p| 2 = 5 2π  Table I. The two fits agree for tmin ≥ 6, and the two-state fit yields a consistent ground-state mass for all tmin > 0.
form of Eq. 6 with n and n restricted to 0 and 1, we isolate the matrix elements in the ground state from those in the first excited state as shown in Sec. IV.
To study electromagnetic form factors, we use one final momentum ( p f = 2π La {0, 0, 0}) and vary the initial momenta over all p i = 2π La {n x , n y , n z } with integer n x,y,z such that n 2 x + n 2 y + n 2 z ≤ 5. For all other charges we project onto p i = 0 by inserting the operator at zero momentum.

C. Statistics
The MILC Collaboration has produced ensembles of roughly 5500 trajectories of 2+1+1-flavor HISQ lattices at the two quark masses. The ensembles at a ≈ 0.12 fm are described in Table I. Five hundred trajectories are discarded for thermalization, which is somewhat more conservative than the 300 discarded by MILC. Correlators are then calculated on configurations separated by five trajectories. On each configuration, we use four smeared sources, displaced both in time and space directions to reduce correlations. Furthermore, two sets of these four source points, again maximally separated in space and time directions, are used on each alternate configuration to reduce correlations. To evaluate the statistical significance of the data, in addition to the full set, we also analyze the roughly 500 configurations with each of these two sets of sources. We verify that the two sets give compatible results and the errors are roughly √ 2 larger compared to the full set.

IV. EXCITED-STATE CONTAMINATION
All observables reported in this paper (charges, charge radii, form factors) need to be calculated between ground-state nucleons. The operators used to create and annihilate the states, defined in Eq. 1, however, couple to the nucleon and all its radially excited states. There are two possible ways to reduce contributions from excited states: by reducing the overlap of the interpolating operator with the excited states and by increasing the time separation t sep = t f − t i between the source and sink to exponentially suppress excited-state contamination. As discussed above, we use sources and sinks with one fixed smearing size that improves overlap with the ground state; nevertheless, the two-point correlator shows significant excited-state contribution extending to t = 5. Thus, using single-state analysis on the threepoint correlator is likely to be problematic for small t sep . Statistics limit the upper value of t sep that can be explored, and we find that the signal degrades very significantly by t sep = 12. We, therefore, investigate up to five time separations between t sep = 8 and 12 to quantify the excited-state contamination as discussed below.
We consider the leading excited-state contamination mass M 1 and its coupling to our operator with amplitude A 1 . We can write the three-point function with source shifted to t i = 0, operator insertion at t = t and sink at where n |O Γ |n is an abbreviation for N n ( p f , s ) |O Γ | N n ( p i , s) . To extract 0|O Γ |0 from the two-and three-point functions we make the following different kinds of fits. In each case, we apply a nonlinear least-square fitter that automatically selects a fit range appropriate to the form used. For each form on each correlator, the fit range is expanded as long as the quality of the fit (in terms of uncorrelated χ 2 /dof) does not sharply decline.
• one-one method assumes a single state dominates the two-point and three-point functions. A 0 and M 0 are extracted from a fit to the two-point function given in Eq. 5 and 0|O Γ |0 is estimated from the three-point functions keeping only the first term in Eq. 7.
• Ratio method also assumes a single state dominates the three-point function. 0|O Γ |0 are estimated from the ratio of three-point to two-point functions, which for large t sep is expected to be a constant, the desired matrix element. Some statistical noise may cancel in the ratio as long as the source and sink operators are identical between the two-and three-point functions, but this relies on there being a good signal in both at separation t sep .
• two-two method: A 0 , A 1 , M 0 and M 1 are extracted from a fit to the two-point function. These amplitudes and masses are used in a two-parameter fit to the three-point function to estimate 0|O Γ |0 and 1|O Γ |0 . In the case of charges where both initial and final nucleon operators are at rest, we can assume 0|O Γ |1 and 1|O Γ |0 are equal, and we analyze only the real part of the three-point function. However, in the case of the form factors, the initial and final states are not the same and both matrix elements must be retained.
• The two-sim method is a simultaneous fit to all t sep , and the same as the two-two method for extracting A 0 , A 1 , M 0 and M 1 . The fit to the three-point function is made to the expression in Eq. 7 using data from all investigated values of t sep simultaneously.
• The two-simRR method -this simultaneous fit to all t sep is the same as the two-sim method but includes the 1|O Γ |1 term. The 1|O Γ |1 term cannot be distinguished from the 0|O Γ |0 term if simulations are done at a single t sep since it depends only on t sep , not t. This contamination is, however, exponentially suppressed as it is proportional to e (M1−M0)tsep .
Note that when fitting the form factors, each mass M n should be replaced by the appropriate energy E n for the momentum p i used.
In this section, we will briefly demonstrate our analysis method on the isovector charges of g A,S,T . We will leave the source-sink dependence and various analysis methods used for the isoscalar charges and form factors to the following sections.
We study two versions of the simultaneous fit including excited-state degree of freedom, with and without the higher-order 1|O Γ |1 term in Eq. 7. Figure 2 shows the fit with the worst quality among all our data: unrenormalized g S from the 220-MeV ensemble with (upper) and without (lower) the 1|O Γ |1 contribution. Both fits capture the data, and the final fit keeping all the terms is marginally better. The two-simRR fit is about factor of 2 noisier than two-sim one. We calculate the difference between the two fits within the jackknife process, and find that for these ensembles the difference is consistent with zero for all charges at zero momentum with our current statistics. To completely illuminate the systematic error, we will use two-simRR fit for the better determined isovector and isoscalar charges. The form factors are noisier than charges, and there we use only the two-sim fit analysis.
The results of fits for the unrenormalized isovector charges are shown in Fig. 3. Estimates from the two-simRR method are shown by the horizontal bands. Based on these two ensembles of roughly 1000 configurations at a ≈ 0.12 fm and the tuned Gaussian-smeared sources used in calculating the quark propagators, we note the following features for isovector g A,S,T : • The statistical errors increase by about 40% with each unit increase in t sep .
• Only data for g A on the 310-MeV ensemble show a small increase (by about 1σ) with t sep between 8 and 12; g S shows a decrease of similar magnitude.
• Based on the trends first seen in the 310-MeV ensemble, we considered it sufficient to investigate the 220-MeV ensemble using only t sep ∈ {8, 10, 12}.
• The two-simRR fit estimates of the central values and errors in the isovector charges are consistent with data from other fits for all values of t sep within statistical errors.
• The errors increase by about 20% on lowering the light (u and d) quark masses by a factor of two, going from 310-to 220-MeV ensemble.
• The signal in g S is the noisiest. Nevertheless, on the 220-MeV ensembles, the error estimate is about 15%, reasonably close to our desired accuracy.
Our conclusion, based on the analysis of the a ≈ 0.12 fm lattices, is that with a well tuned smeared operator the central values and error estimates from the two-sim fit agree with those from the other fits for separation t sep = 10 which corresponds to t sep ≈ 1.2 fm in physical units for pion mass as light as 220 MeV in our case. If restricted to simulations at a single t sep , we consider the data at t sep = 10 the best compromise between reducing excited-state contamination and having a good statistical signal with O(1000) lattices.
Noting the analysis of the bare axial charges presented by the ETMC [34,35], CSSM [36] and LHPC [37] collaborations, we conclude that excited-state contamination becomes comparable to (or smaller than) statistical errors for t sep ≥ 1.2 fm. These collaborations in their g A analyses have explored using summation and variational methods at different values of the lattice spacing and quark masses and with different number of flavors (see Table III). The summation method implemented by the CLS-Mainz [38][39][40] collaboration sums over the full range t f -t i for multiple t sep , including time-slices close to the source and the sink where excited-state contributions are the largest. Also, their fit ansatz does not take into account contributions from the transition terms such as 0|O Γ |1 . They conclude that the summation method gives estimates 1-2σ larger than the "plateau" (ratio) method with t sep ≈ 1 fm. However, examining their data in detail at each t sep , the various estimates are consistent within errorbars.

V. NONPERTURBATIVE RENORMALIZATION IN RI-SMOM SCHEME
This section describes the lattice calculation of the renormalization constants Z A,S,T in the RI-sMOM scheme (regularization-independent symmetric momentum-subtraction) [41,42] and their conversion to MS scheme at 2 GeV. For this calculation, the HISQ configurations are fixed to Landau gauge after hypercubic (HYP) smearing and the clover propagators S(0, x) are calculated using point sources. From these propagators, we extract the wavefunction renormalization constant Z ψ and calculate the truncated three-point correlators as functions of renormalization scale µ to estimate the three Z A,S,T .
The clover action we use is improved to O(a) only in tadpole-improved perturbation theory. Our renormalized operators, defined as O R Γ = Z Γ O Γ , also do not include any O(a) improvements. The errors in the results, therefore, start at O(a). The three-point function we calculate is defined by Λ Γ (x, 0, y) = ψ(x)O Γ (0)ψ(y) = S(x, 0)ΓS(0, y), (8) with Γ representing the Dirac matrices I (scalar), γ µ γ 5 (axial-vector) and σ µν (tensor). In momentum space, this three-point function is where S(p) is the Fourier transform of S(x, 0), and we have applied γ 5 -Hermiticity to the right quark leg. From this we construct the amputated three-point correlator and the projected amputated three-point function where the projector P Γ for the RI-sMOM scheme is I (scalar), (q µ /q 2 )γ 5 / q (axial-vector) and (i/12)γ [µ γ ν] (tensor). In the RI-sMOM scheme, the allowed momenta satisfy the relations The renormalized projected amputated three-point function is defined as where Z Γ is the operator renormalization constant. In the RI-sMOM scheme, this is set equal to one, its treelevel value, for all tensor structures. This condition fixes the value Z −1 ψ Z Γ at the subtraction point. Similarly, the wavefunction renormalization constant Z ψ is defined by Having extracted the renormalization constants in the RI-sMOM scheme at scale µ = q 2 , we convert them to the MS scheme at µ = 2 GeV using the one-loop conversion factors in Landau gauge [42,43]: where α s (µ) in this horizontal matching can, to O(α 2 s ), be taken to be the coupling in any scheme, RI-sMOM or MS scheme or from the plaquette using the BLM procedure [44]. We use the coupling α MS . Note that the above conversion factors are computed in the chiral limit, and possible O(ma) corrections are ignored.
These Z, now in MS scheme defined at scale µ = q 2 , are then run to µ = 2 GeV using where the evolution function E Γ (α s (µ)/π) at two-loop is [45] β 0 and β 1 are the universal coefficients of the β-function, with C A = 3, C F = 4/3 and T F = 1/2; γ 0 Γ and γ 1 Γ are the first two coefficients in the anomalous dimension of the operators in the MS scheme, and α s (µ) at two-loop has the following expression where We have analyzed 101 (60) lattices on the M π = 310 MeV (220 MeV) ensemble and plot the results for the ratios Z Γ /Z V in Fig. 4. In extracting the final estimates, we incorporated the following observations: • The data for the ratios Z Γ /Z V are more smooth in all cases, and the fits are more stable. Some sources of systematic uncertainty cancel in the ratios.
• We find no discernible difference as a function of the pion mass, and the fits to the two data sets are indistinguishable as shown in Fig. 4 where both sets are plotted together. We, therefore, neglect possible mass dependence and make a single fit to the combined data.
• The estimates of Z Γ are sensitive to the different possible combinations of momentum components p µ for a given value of p 2 = q 2 . These differences are indicative of lattice discretization effects due to reduction of the continuum Lorentz group to the hypercubic group on the lattice. Since these effects are smallest in momenta with symmetric components, for example the combination {1, 1, 1, 1} versus {2, 0, 0, 0}, we choose the most symmetric combinations for our final analysis, shown in Fig. 4. That is, for a given p 2 , we choose momenta that minimize p [4] 1,2 /(p 2 ) 2 with p [4] = µ p 4 µ .  (2) TABLE II. The results for Z ψ , ZV , and the three ratios of renormalization constants ZA,S,T /ZV in the MS scheme at 2 GeV. The lattice calculation is done in the RI-sMOM scheme. The error estimates cover the spread in values from different methods used as discussed in the text. We do not find significant variation with the pion mass and quote a single result that is used for both ensembles.
At weak enough coupling the final results for the renormalization constants in the MS scheme at 2 GeV should be independent of the q 2 value (within a window Λ QCD q 2 π/a) selected to define them in the RI-sMOM scheme. This requires that the errors due to lattice discretization and the use of truncated perturbation theory to convert from RI-sMOM to MS scheme are negligible. The data for Z S /Z V , however, show very significant q 2 dependence. This calls into question whether on 0.12-fm lattices there exists a window in which the nonperturbative effects at low q 2 and lattice-discretization effects at high q 2 are both small.
Recent analyses have shown that smeared lattices alter the window in which lattice data are consistent with perturbation theory [26]. HYP smearing smooths the gauge fields in a 2 4 hypercube, so gluons with momentum of order 1/a are suppressed, modifying the high-q 2 behavior of Z Γ . Since our data for Z S /Z V shows large q 2 dependence, we investigate two prescriptions, each defined such that the renormalized charges, Z Γ /Z V × g bare Γ /g bare V , have a well-defined continuum limit. In the first case we pick the value of Z Γ /Z V at a fixed physical value, q 2 = 5 GeV 2 . The error on this estimate is taken to be half of the total variation in the range 4 < q 2 < 6 GeV 2 . In the second case we assume an ansatz for the lattice artifacts and fit the q 2 dependence of the data [26]. We try several variations of c/q 2 + Z + d 1 q + d 2 q 2 and of the fit range. We do not include dependence on the pion mass, since the data do not show any significant difference between the two ensembles. We find that the ansatz c/q 2 + Z + d 1 q is a good fit to all the data for q 2 > 1 GeV 2 , as shown in Fig. 4. We take the estimates of Z from these fits as our central values. In addition to fits to the ratios Z Γ /Z V , we also constructed the ratios, within a jackknife process, from fits to individual Z Γ .
It turns out that the three estimates are consistent. This is easy to see, from data shown in Fig. 4, for Z A /Z V and Z T /Z V since the q 2 dependence above q 2 = 4 GeV 2 is small.
To assign an overall error to Z Γ /Z V we took into account the different estimates. Our final error estimate covers (i) the variation with q 2 between 4 ≤ q 2 ≤ 6 GeV 2 in method one, (ii) the variation with the fit ansatz in method two and (iii) the difference in the three estimates. Note that this conservative error estimate, given in Table II, is much larger than the error from the fits as shown in Fig. 4. It also captures the spread in data cor- responding to points with different breaking of rotational symmetry, i.e. data with larger p [4] 1,2 /(p 2 ) 2 values. Our final estimate for the renormalized charges is obtained by multiplying the ratios g bare Γ /g bare V with the corresponding ratios Z Γ /Z V since the vector Ward identity implies Z V g bare V = 1.

VI. NUCLEON ISOVECTOR CHARGES
We present results for the three unrenormalized charges from two-simRR fit in Table IV and the final  renormalized values in Table V. To facilitate comparison with previous work with improved actions, we give a compilation of lattice parameters used by other collaborations in Table III, and selected results are shown in Fig. 5.
We employ two strategies to extract renormalized charges g A,S,T , and their difference is used as an estimate of systematic errors. In the first method, we extract, under two separate jackknife analyses due to the different numbers of configurations analyzed, the unrenormalized charges g bare A,S,T and the renormalization constants Z A,S,T in the MS scheme at 2 GeV. These are multiplied together with relative errors added in quadrature. In the second method, we extract the ratios g bare A,S,T /g bare V and Z A,S,T /Z V and use the fact the Z V g bare V = 1. Our data, however, give Z V g V = 0.95(3) and 0.96(4) for the 310and 220-MeV ensembles, respectively; this leads to a difference of up to 0.04 between the two estimates.

A. Axial Charge gA
The best lattice-QCD calculations, involving multiple lattice spacings (including continuum extrapolation of estimates) and high statistics, yield estimates with about 5% statistical error and are a few standard deviations lower than the experimental values. As shown on the right-hand side of Fig. 5, most lattice estimates lie between 1.05 and 1.20, which is 5-15% below the experimental values.
Most collaborations, such as ETMC [34,35], CSSM [36], CLS-Mainz [38][39][40], and LHPC [37], find that estimates of central value on g A increase with nucleon source-sink separation t sep in three-point function and that the statistical errors also grow. Our two-simRR fit result agrees with estimates from the larger values of t sep ≥ 10 and after multiplication by Z A gives 1.193(68) as shown in Table V. This estimate, without continuum extrapolation, is about 1σ below the experimental value.
It has previously been shown that estimates of g A can be underestimated due to insufficiently large spatial volumes (see Refs. [48,58] for example), especially those with M π L < 4. Finite-volume corrections based on HBXPT [48,57] or a simple parametrization formula [58] are often used to correct this systematics. However, a recent global survey on lattice g A [13] suggested that there might exist ambiguities in the HBXPT correction and that larger volume M π L ≈ 6 might be needed to reproduce the experimental value of g A . Thus, a more detailed finite-volume study with pion masses below 250 MeV is needed to understand this systematic better.
In their most recent work [48], QCDSF collaboration find g A = 1.24(4) based on a new data point at the physical pion mass. Data at and above their next lightest M π = 253 MeV ensemble lie in the range 1.05 ≤ g A ≤ 1.10. Their chiral fit, consequently, suggests that g A increases significantly between 140 and 250 MeV. Recent ETMC results [57] do not show a significant increase between 210 and 250 MeV. To clarify the chiral behavior, therefore, requires data below 210 MeV.
To summarize, current data suggest that to obtain a precise value of g A will require simulations close to the physical light-quark masses, large lattices, high statistical precision and a careful study of excited-state contamination.

B. Scalar and Tensor Charges gS and gT
Our final estimates, given in Table V and shown in Fig. 6, are g S = 0.72 (32) and g T = 1.047 (61). LHPC      [24,25,34,39,46,49,50,52,56]. Note the change in scale between the experimental and theory plots. Most of the errorbars here are statistical only. In data from the few calculations that also quote systematic errors, we add these to the statistical ones as outer errorbar bands, marked with dashed lines. (right) Calculations of gA using at least 2+1 flavors O(a)-improved dynamical fermions, plotted as a function of M 2 π , with mπL > 4 to avoid systematics due to small spatial extent.  ) and using ZV g bare V = 1 as discussed in the text. The error quoted is obtained by combining the statistical and systematic errors in the ratios of Z's and g's in quadrature under the assumption that they are independent. The last row gives estimates extrapolated to the physical pion mass Mπ = 140 MeV using a fit linear in M 2 π .
has recently published lattice calculations giving g S = 1.08(28)(16) and g T = 1.038(11)(12) [37]. Their main result exploits a single ensemble created with tree-level clover-improved Wilson fermions on 48 4 lattices at a very similar value of the lattice spacing a ≈ 0.116 fm but at M π ≈ 150 MeV. Their result [56] for g A on the same lattices is near 1 (1.00 (8)), which, combined with current experimental values for the neutron lifetime, implies V ud > 1. It seems evident that not all systematic uncertainties are under control for this ensemble. One should also note that both calculations lack continuum extrapolation. Our future update will include 0.09-fm and 0.06fm ensembles and thus reduce the systematic uncertainty due to the continuum extrapolation. The phenomenological implications of the estimates of g S and g T are given in section IX.

VII. NUCLEON CONNECTED ISOSCALAR CHARGES
In this section we summarize results for the connected diagrams contributing to the isoscalar charges g s A,S,T . One motivation for their study is that g s T probes novel contributions to the quark electric dipole moment inside the nucleon [59], as discussed below.
The neutron electric dipole moment (nEDM) d n is a measure of the distribution of positive and negative charge inside the neutron. To generate a nEDM, the theory must include processes that violate CP-symmetry. There are two sources of CP-violation in the Standard Model: the phase in the CKM matrix and the Θ-term in the Lagrangian. The CKM phase gives rise to a nEDM d n ∼ 10 −32 e · cm [60] that is too small to account for the observed baryon asymmetry of the universe. The current upper limit Θ < 10 −10 , an unnaturally small number, is obtained from the current experimental limit d n < 2.9 × 10 −26 e · cm [61]. Possible new interactions at the TeV scale (supersymmetry, left-right models, extra dimensions) are a rich source of additional CP violation that could give rise to a large nEDM in the range 10 −28 -10 −26 e · cm, enough to explain baryogenesis. This is an exciting scenario, since the next generation of EDM experiments are targeting 10 −27 e · cm. Independent of the details of the candidate theories at the TeV scale, in the effective field theory language, there are two CP-violating operators at dimension five that give rise to CP-violating interactions of the electric field with the neutron. These are the quark EDM (qEDM) and quark chromoelectric dipole moment (CEDM) operators [59,62] Here, F µν is the electromagnetic field, G µν is the gluon field and e and g 3 are their respective couplings. The couplings, {d γ,G u , d γ,G d }, encapsulate the interaction of the quarks with the photon and gluon, and v H = 246 GeV is the vacuum expectation value of the Higgs field.
The matrix elements of these operators are very poorly known [63] and are needed in order to use future measurements of nEDM to tighten constraints on the allowed parameter space of BSM theories.
The matrix element of the qEDM operator is an extension of the lattice-QCD calculation of g T ; one needs to calculate terms of the form which can be expressed in terms of the isoscalar and isovector tensor charges of the neutron. We have already discussed the calculation of the isovector tensor charge g T , and present first results for the connected part of isoscalar tensor charge. The remaining disconnected part is beyond the scope of this study.
In Fig. 7, we show estimates of the unrenormalized connected parts of the isoscalar axial, scalar and tensor charges, g s A,S,T , extracted using the same analysis methods described in Sec. IV. Again, we find consistency between the methods for our tuned smearing-parameter choices. The one exception is the isoscalar scalar charge, for which estimates based on one-one or ratio methods show an increase with the source-sink separation for t < 10. Overall, the two-simRR method gives our best estimates, and these agree with those from the other methods for t sep ≥ 10.
The renormalization constants for the isoscalar charges also receive contributions from disconnected diagrams. In the approximation that disconnected diagrams are neglected, the renormalization constants for the isoscalar and isovector charges are the same. We, therefore, use results in Table II to renormalize the connected isoscalar charges at 2 GeV in MS scheme. These estimates of the renormalized charges are given in Table VI. We also include estimates for the value at the physical pion mass using a linear extrapolation in M 2 π .

VIII. ISOVECTOR ELECTROMAGNETIC FORM FACTORS OF THE NUCLEON
The Dirac and Pauli form factors (F 1,2 ) are extracted from the matrix elements of the isovector vector current in the nucleon state N through the relation where the momentum transfer q = p f −p i . Another common set of definitions of these form factors, widely used in experiments, are the Sachs (electric and magnetic) form factors; which can be related to the Dirac and Pauli form factors through A compilation of the lattice parameters used by various collaborations performing simulations of the electromagnetic form factors is given in Table VII.
The vector-current matrix elements, N |V µ |N in the ground state (with n = n = 0) at different momenta are obtained by using the same projection matrices T as in Eq. 6. This overdetermined system of linear equations allows solution for the Dirac and Pauli form factors F 1,2 with various Lorentz indices µ and momenta p i for a particular Q 2 . To minimize the excited-state contribution to the ground-state matrix element, we again employ a number of fits, using the notation established in Sec. IV. The upper two plots of Fig. 8 show 310-MeV Dirac and Pauli form factor at each Q 2 as a function of t sep . Once again, we observe that fitted values from two-two method are consistent with those from two-sim fit, while one-one method is less consistent. The lower two plots of Fig. 8 show examples from 220-MeV ensemble for all values of t sep investigated. We adopt as our preferred value the two-sim fit, which takes into account excited-state systematics. Relative to our adopted values, we find the central values for the one-one and two-two method shift by no more than 2σ with our statistics with the smallest separation. We find that the differences with and without the 1|O Γ |1 term in Eq. 6 are consistent with zero within errors, and including the RR term increases the errors. Therefore, the final results presented use the two-sim fit neglecting this term.
We normalize the Dirac and Pauli form factors F 1,2 by the value F v 1 (Q 2 = 0) determined directly in the calculation; thus, the renormalization factor Z V cancels, and the O(a)-systematics are reduced. In the case of F 1 , we explore two functional forms to characterize the Q 2 behavior: a conventional dipole (1 + Q 2 /M 2 ) −2 with one free parameter and a more general quadratic in Q 2 , (1 + b 1 Q 2 + b 2 Q 4 ) −1 with two free parameters. We find that dipole form does not work for all of our form factor data at any source-sink separation; we have to cut the data as low as Q 2 ≤ 0.4 GeV 2 to make it work. Since there is no fundamental physics reason for using this form, we take the central value from the general quadratic form, which gives a much better fit. Unfortunately, including more free parameters in the fit results in the final extrapolated value of the charge radii having larger uncertainty. In the case of F 2 , we also investigated multiple fit ansätze: with the anomalous magnetic moment κ v ≡ F v 2 (0). We find that all three ansätze capture the data reasonably, since the errorbars are larger in Pauli form factor. However, we choose to use the general ansatz for the final fit.  The renormalization constants used to convert these to the MS scheme at 2 GeV are the same as given in Table II. The error quoted on the renormalized charges are obtained by combining the statistical and systematic errors in the ratios of Z's and g's in quadrature under the assumption that they are independent. The extrapolation of the renormalized charges to 140 MeV is carried out using a fit linear in M 2 π .
Collaboration    Figure 9 shows both 310 and 220-MeV Dirac and Pauli form-factor results with the general ansatz, and the extrapolation to the physical pion mass point. We found small pion-mass dependence on these ensembles. For Dirac form factors, our form factors are larger than the experimentally reconstructed values; consequently, the smaller slope around the Q 2 = 0 point gives smaller charge radius as defined below. This feature has been observed in the past with pion mass larger than 300 MeV. The disagreement in the Pauli form factors is less severe, but there is also misalignment in the small-Q 2 region.
The size of the nucleon characterized by the effective Dirac and Pauli radii can be determined from the electromagnetic form factors. These are determined from the slope of the corresponding form factor in the zero-Q 2 limit.
Since the value of the smallest momenta allowed in typical lattice simulations is large, to extract the radii it is important to develop ansätze that capture the low-Q 2 behavior well. We use the general ansatz defined above to extract the radii, since they give the best fit. Attempts are being made to obtain data at smaller momenta to improve the determination [69]. Figure 10 shows our results for the Dirac and Pauli radii for three fit methods and the two ensembles. In the case of the Dirac radii, the one-one method on the 310-MeV ensemble data becomes reliable only at larger separation; the central values increase (within errorbar) with separation. We find that using the two-two method, which includes the leading Roper-nucleon contribution, gives consistent radii for all source-sink separations. The 220-MeV data do not show any significant trend with t sep . In the case of Pauli radii, we do not observe significant dependence on the separation nor on the analysis method. Similar conclusions also apply to the anomalous magnetic moment. In all cases, the estimates from the two-two method agree well with those from the two-sim method, which we adopt as our preferred results, collected in Table VIII. Recent analysis by the LHPC [56] shows a significant increase in r 2 1 with t sep , especially for their M π = 150 MeV ensemble. We do not observe a statistically significant effect and need data on lower-M π ensembles to check the trend. The CLS-Mainz Collaboration [38] find a dependence on t sep , however their four values of t sep are smaller than 1 fm, within the range of separations where we find excited-state contamination.
A summary of all the N f = 2 + 1 and N f = 2 + 1 + 1 lattice calculations of the isovector Dirac and Pauli meansquared radii are summarized in Fig. 11 along with the lowest-order heavy-baryon chiral perturbation theory (HBXPT) using experimental inputs [70,71]. Note that most groups only report statistical errors, which are shown in this figure; ours also includes the systematics due to the choice of fit-form ansatz.

IX. CONCLUSIONS
In this paper we demonstrate that g S and g T can be calculated to a precision of 20% or better on ensembles with O(1000) configurations; we plan to increase the statistics in the future by doubling the number of source points simulated on each lattice. This is significant, since this level of precision is needed to leverage experimental measurements of b and b ν at the 10 −3 level to constrain novel scalar and tensor interactions at the TeV scale.
We show that contamination from excited states can be understood and taken into account by doing the calculations at multiple values of t sep and performing a simultaneous fit to all the data using Eq. 7 while keeping one excited state in the analysis. We also find that with O(1000) lattices, a consistent estimate is obtained with t sep ≈ 1.2 fm. In cases where sufficient computer resources are not available to carry out studies with multiple t sep , this separation should be sufficiently large for well tuned nucleon creation/annihilation operators at pion masses above 220 MeV.
We find that renormalization constants Z A,S,T,V can be calculated with better than 5% accuracy using the RI-sMOM scheme. We do not find a window in which the data for Z S and Z T versus q 2 matches perturbative behavior on the a = 0.12 fm ensembles, however, different estimates of individual Z's and the ratios Z A,S,T /Z V lie within 10% of unity. We find that the data for the ratios Z A,S,T /Z V have less scatter than individual Z's. We assign a conservative estimate for the systematic errors to cover the spread. Within these error estimates, the final values of the renormalized charges obtained using different analysis strategies give consistent results. The substantial artifacts in the calculation of Z's due to the coarse discretization of the lattice underscore the need for a controlled continuum extrapolation of the renormalized charges, which we plan to investigate in the future using data at three lattice spacings (a = 0.12, 0.09 and 0.06 fm).
The estimates of charges at the two values of quark masses corresponding to M π = 310 and 220 MeV, are in most cases within their respective 1σ errors. A simple linear chiral extrapolation to the physical pion mass introduces additional uncertainty. To reduce systematic errors due to chiral (and continuum) extrapolation will require higher statistics simulations for at least three values of the quark mass or simulations at the physical mass. Since the number of configurations in the MILC ensembles are fixed, our current strategy to increase statistics is to double the number of source points simulated on each lattice.   11. Summary of the isovector Dirac and Pauli meansquared radii from all currently existing N f = 2 + 1 and 2 + 1 + 1 nucleon electromagnetic form-factor calculations [24, 56-58, 65, 72]. The dashed curve indicates the leading-order HBXPT prediction.