Nucleon Structure from Lattice QCD Using a Nearly Physical Pion Mass

We report the first Lattice QCD calculation using the almost physical pion mass mpi=149 MeV that agrees with experiment for four fundamental isovector observables characterizing the gross structure of the nucleon: the Dirac and Pauli radii, the magnetic moment, and the quark momentum fraction. The key to this success is the combination of using a nearly physical pion mass and excluding the contributions of excited states. An analogous calculation of the nucleon axial charge governing beta decay has inconsistencies indicating a source of bias at low pion masses not present for the other observables and yields a result that disagrees with experiment.


Introduction
Lattice QCD is the only known rigorous framework for ab-initio calculation of the structure of protons and neutrons with controllable uncertainties. It can provide quantitative answers to both fundamental questions such as the quark and gluon composition of the nucleon spin and phenomenological questions such as the sensitivity of modern detectors to physics beyond the Standard Model (BSM), to fundamental symmetry violations, and to hypothetical dark matter particles [1,2,3]. It also has the potential to help resolve discrepancies between inconsistent experimental results. However, with current computer resources, its predictive power is limited by uncertainties arising from heavier than physical quark masses, finite lattice spacing and volume, incomplete removal of excited states, and omission of disconnected contractions. Progress in lattice calculations of nucleon structure hinges on identifying and removing the most significant among these uncertainties.
Significant effort has been focused on lattice calculations of the isovector Dirac and Pauli radii (r 2 1,2 ) v , anomalous magnetic moment κ v , axial charge g A , and quark momentum fraction x u−d : p|qγ µ γ 5 q|p = g Aūp γ µ γ 5 u p , where Q 2 = −q 2 = −(p − p) 2 and u p , u p are nucleon spinors, and we note that computationally expensive disconnected contractions cancel in isovector observables. Until recently, although some success has been achieved [4,5,6,7,8,9], lattice results relied heavily on large extrapolations using Chiral Perturbation Theory (ChPT) yielding potentially uncontrollable corrections. This is particularly problematic for some observables, e.g., (r 2 1,2 ) v and x u−d , for which ChPT predicts rapid change towards the chiral regime. For example, in typical lattice calculations with pion masses 250 MeV, prior to extrapolation to m phys π ≈ 135 MeV, (r 2 1 ) v is underestimated by ≈ 50% [10,6,7,8], x u−d overestimated by 30 − 60% [5,11,12], and g A underestimated by ≈ 10% [4,13,14], compared with experiment. These glaring discrepancies and the dependence on large extrapolations clearly indicate the need for calculations near the physical pion mass. Moreover, it has recently been found that excited-state effects become worse with decreasing pion mass [15], and their careful analysis is required before even attempting extrapolations in the pion mass towards the physical point using ChPT.
In this paper, we report the first Lattice QCD calculation of five nucleon structure observables using pion masses as light as m π = 149 MeV. This is so close to the physical value that chiral extrapolation of these observables yields changes between the lowest pion mass point, m π = 149 MeV, and the physical point, m phys π ≈ 135 MeV, that are less than or comparable to the statistical uncertainty of the m π = 149 MeV data. For each ensemble, we remove excited-state contaminations by varying the source-sink separation in the range ≈ 0.9 . . . 1.4 fm and apply the summation method [16] to extract the ground state matrix elements. We observe remarkable agreement with experiment, well within the statistical uncertainties, for the isovector Dirac and Pauli radii, the anomalous magnetic moment, and the quark momentum fraction, all computed with the same methodology. However, as discussed below, there are inconsis-tencies in the results for the axial charge, g A , calculated at the lowest two pion masses that do not arise for the other observables, that make the results suspect, and that indeed lead to disagreement with experiment. There appears to be a source of bias significantly affecting the evaluation of the axial charge at low pion masses, leading to a qualitative difference in behavior between g A and the other observables.

Lattice Results
We perform calculations using ten ensembles of Lattice QCD gauge fields generated with O(a 2 )-improved Symanzik gauge action and tree-level cloverimproved Wilson fermion action using 2-HEX stout gauge links [17]. Two light u and d (with m u = m d ) and one heavier strange (m s m u,d , tuned to be close to physical) quark flavors are simulated fully dynamically, while effects of heavier quarks are neglected. In addition to varying the pion mass in the range 149 . . . 357 MeV, we include different spatial volumes, time extents of the lattice, and one ensemble with a smaller lattice spacing in order to estimate the size of the corresponding systematic effects; see Tab. 1. Nucleon matrix elements are extracted from nucleon 3-point functions that are computed with the standard sequential source method. Nucleon field operators are optimized to overlap as much as possible with the single-nucleon ground state at rest by tuning the spatial width of Gaussian smeared quark sources.
In order to discriminate between the ground and excited-state matrix elements on a Euclidean lattice, we vary the timelike distance ∆t between nucleon sources and sinks in the 3-point functions. With increasing ∆t, excited states in 3-point functions are suppressed as ∼ e −∆E∆t/2 and disappear in the ∆t → ∞ limit, where ∆E is the (potentially m π -dependent) energy gap to the closest contributing state. However, using a large source-sink separation is impractical, since statistical noise grows rapidly with ∆t. Instead, we combine calculations with three values of ∆t ≈ 0.9, 1.2, 1.4 fm using the summation method [16], which benefits from improved asymptotic behavior [18,19] and which we find the most reliable and robust method with the presently existing statistics. In a separate publication [20] we report detailed studies and comparisons of different methods, including also the generalized pencil-of-function (GPoF) method [21].
We present and discuss our results for all the observables in Figs. 1-5 below. To emphasize the importance of controlling excited states, in the same figures we also show the standard plateau method results at given ∆t (open symbols) for the three lightest m π . As this separation is increased, the data points consistently approach our final results, while their uncertainties increase as expected. As illustrated by the data, the effect of the excited states may be very dramatic, especially for (r 2 1 ) v , κ v (r 2 2 ) v and x u−d . Our lightest pion mass value m π = 149(1) MeV is only ≈ 10% higher than the physical pion mass in the isospin limit m π ≈ 134.8 MeV [22]. To make the small extrapolation of our data in m π to the physical point, we add the 202(1) MeV and the large-volume 254(1) MeV ensembles to minimize Table 1: Lattice QCD ensembles. The strange quark mass ms is tuned to be close to physical. dependence on the functional form and volume and perform ChPT fits using standard parameters from ChPT phenomenology. Table 2 gives details of the fits, results and comparisons with experiment. Note that due to the short extrapolation distance, the change in the values of the chiral fits between m π = 149 MeV and m π = 135 MeV is smaller than the statistical uncertainties of the m π = 149 MeV data in all cases except for (r 2 1 ) v , for which the two are comparable. Although they are not shown, again because of the short extrapolation distance linear fits to the data yield extrapolated values statistically consistent with the chiral fits. Hence, the validity of our final extrapolated results is not dependent on the convergence of ChPT and in contrast to all previous calculations, the dominant uncertainties in the final results are the statistical uncertainties in the lowest mass data points rather than the chiral extrapolations. To explore the consistency of our data with the low-energy theory, we repeat the chiral fits using the full available range of m π = 149 . . . 357 MeV and show the error bands for both ranges in Figs. 1-5. For (r 2 1 ) v , (r 2 2 ) v , κ v , and x u−d results from m π 350 MeV and m π 250 MeV fits agree within statistical uncertainty whereas for g A , both the error bands and the extrapolated result are inconsistent, providing one indication of an unphysical bias in the calculation of g A to be discussed below.
To explore nucleon electromagnetic structure, we compute matrix elements of the quark-nonsinglet vector current p | ūγ µ u −dγ µ d |p between polarized proton (uud) states with different momenta p, p . We extract isovector Dirac and Pauli form factors F v 1,2 (Q 2 ) from these matrix elements, and then the "radii" The Dirac radius (r 2 1 ) v is shown in Fig. 1. We compare it to the experimental value (r 2 , with the error bar dominated by the uncertainty in (r 2 E ) p , the proton electric charge radius. We show two experimental values for (r 2 1 ) v in Fig. 1, which correspond to two inconsistent values for (r 2 E ) p : the CODATA recommended value [36] (also cited Table 2: Extrapolations to the isospin-limit m phys π = 134.8 MeV [22] using data points 149 ≤ mπ ≤ 254 MeV. The footnotes describe details of the ChPT formulas and fixed parameters used in the fits. The right column shows the deviation from experiment, relative to the lattice uncertainty.  [28,4,29,30], c A = 1.5 [10] [31]. c Includes higher-order "core" term [32]. [28,4,29,30], [24]. Using Ref. [24] for both (r 2 E ) p,n results in (r 2 1 ) v = 0.640(9) (higher exp. point in Fig. 1).
by PDG [24]) based on hydrogen spectroscopy and e − p scattering and the recent and controversial result from measurement of the µ − p Lamb shift [23] 1 . Relative to the lattice uncertainty, the extrapolated value deviates from the µp Lamb shift value by −0.07σ lat and from the PDG value by −2.17σ lat , where σ lat is the uncertainty of the present calculation.
In a similar fashion, we extract the isovector anomalous magnetic moment κ v = κ p − κ n and the Pauli radius (r 2 2 ) v = (κ p (r 2 2 ) p − κ n (r 2 2 ) n )/(κ p − κ n ) from the Pauli form factor F 2 (Q 2 ). The results are less precise than (r 2 1 ) v because the forward values F 2 (0) and dF2(0) dQ 2 must be extrapolated, and we use , the Q 2 fit is less precise in smaller spatial volumes. This explains the significant increase of error bars in Figs. 2, 3 going from 32 3 to 24 3 lattices (e.g., at m π ≈ 250 MeV), compared to the corresponding error bars of (r 2 1 ) v in Fig. 1. Fortunately, our calculation at the lightest pion mass uses the largest V s ≈ (5.6 fm) 3 yielding the smallest Q 2 min ≈ 0.05 GeV 2 and is the least affected by this problem.
We compute the isovector quark momentum fraction 2 x u−d from the for- 1 It is notable that the new experimental value from the µ − p Lamb shift is in better agreement with phenomenological fits based on dispersion relations [37]. ward matrix element of the operator in Eq. (2) and renormalize the lattice value to the standard MS(2 GeV) scheme using the RI /M OM method [38]. In Fig. 4 we show our results together with the CTEQ6 value [25]. Achieving agreement between Lattice QCD and phenomenology for x u−d is one of the most important accomplishments of this paper: previous lattice calculations [35,39,12] consistently overestimated 3 the phenomenological value. The combination of our use of a nearly physical pion mass and removal of excited state contamination has eliminated the discrepancy. The presence of substantial excited-state contaminations in x u−d was first hinted at in Ref. [40] for m π = 493 MeV, later substantiated by an early report from the present study [15] and at m π = 373 MeV [41], and recently further confirmed at m π ≥ 195 MeV [42].
In contrast to the observables discussed above, the computed value of the nucleon axial charge g A disagrees significantly with g exp A = 1.2701(25) [24]. The most striking feature of the results in Fig. 5 is the dramatic decrease in g A as m π decreases from 254 MeV to 202 MeV and on to 149 MeV. Whereas the m π 250 MeV and m π 350 MeV chiral fits for the other observables in Figs. 1-4 are consistent, this decrease makes the two fits qualitatively different for g A . In looking for the origin of the strong decrease for g A at pion masses antiquarks, i.e. x q = 1 0 dx x fq(x) + fq(x) , where f q(q) is a parton distribution function (PDF). 3 There is a hint in a recent lattice calculation in Ref. [12] that x u−d decreases and becomes closer to phenomenology as the pion mass approaches the physical value.  202 and 149 MeV we note that when removing excited state contaminants by increasing ∆t, g A decreases strongly for both, whereas at 254 MeV it increases by a comparable amount. None of the other observables show this kind of reversal in behavior between 254 MeV and lower pion masses. Excluding our two lightest m π ensembles, our results are consistent with a recent calculation [9,42] that uses O(a) improved Wilson fermions and the summation method for 195 MeV m π 649 MeV. In that calculation, chiral extrapolation to the physical mass is statistically consistent with experiment and g A (∆t) typically increases 10% when ∆t increases by 0.5 fm, which is the behavior for our 254 MeV ensemble. All this evidence suggests that there is an unphysical source of bias that affects g A at the small pion masses 202 and 149 MeV more significantly than the other observables we calculate.
One possible source of this bias could be the presence of thermal states. The lightest state in QCD is a pion at rest, and its contribution to Green's functions is suppressed as e −mπLt , where L t is the time extent of the lattice. The ensembles available to us in Tab. 1 have values of m π L t much smaller than the values L t = 2L s normally used, which may make our calculations susceptible to thermal effects, especially with lighter pion masses. Note that, in contrast to finite-L s effects (FVE) discussed in the next paragraph, finite-L t effects can add contributions to the spectrum that lie below the groundstate nucleon, significantly changing the time dependence of our two-point and three-point functions and making it difficult to isolate the desired single-particle hadron state. Although we find that the behavior of g A (∆t) at m π ∼ 250 MeV does not depend significantly on m π L t , the qualitative difference between the behavior of g A (∆t) at m π = 254 MeV and at 149 and 202 MeV, which both have very small m π L t , is suggestive of thermal effects for these lower pion masses. To assess the effect of the bias in these lowest two masses, we show the result of a chiral extrapolation of all the other ensembles in Fig. 6. It is noteworthy that this extrapolation is consistent both with experiment and Ref. [9], suggesting that the bias is indeed concentrated in the two ensembles with small m π and m π L t that are expected to be particularly susceptible to thermal state effects.
In other works, discrepancies in g A have been attributed to finite volume effects (FVE). According to Ref. [13], FVE may lead to a 9% bias in g A at m π L s ≈ 4.5 in calculations with domain wall fermions, and as much as ≈ 25% in calculations with N f = 2 flavors of Wilson fermions. However, our data at m π ≈ 250 MeV suggest that FVE are negligible: if we assume that the FVE correction scales as δ FVE g A ∼ e −mπLs as in Ref. [13], and use the g A values from the two lattices with m π ≈ 250 MeV that differ only in spatial volume (m π L s = 3.6 and 4.8), we obtain the estimate δ FVE g A = −0.02 (6) for m π = 149 MeV, L s = 5.6 fm, which is much smaller than the observed discrepancy δg A ∼ −0.30.

Discussion
The crucial advances in this work are the careful and uniform control of excited-state contamination and extending the pion mass range down to the nearly physical value m π = 149 MeV. Because of this, we have achieved agreement with experiment for the first time for the Dirac and Pauli radii, the magnetic moment, and the quark momentum fraction and observed indications of internal inconsistency in the calculation for the one observable that disagrees with experiment, the axial charge. This is strong evidence that this combination of methods is sufficient to correctly estimate and remove systematic uncertainties in most nucleon structure calculations and diagnose cases in which some additional bias renders the calculation unreliable. Although the lightest pion mass used in our calculation is still above the physical one, only a very short-range extrapolation is required. As Figs. 1-4 display, the extrapolated results are strongly influenced by the lowest pion mass data, which are themselves already within 1σ of experiment for three observables and nearly so for (r 2 1 ) v . Moreover, the change in the values of our fits from the lowest pion mass, m π = 149 MeV, to the physical point, m π = 134.8 MeV, is smaller than, or in the case of (r 2 1 ) v comparable to, the statistical uncertainty of the m π = 149 MeV lattice data. Thus, in contrast to all previous calculations, the dominant uncertainties in the final results are the statistical uncertainties in the lowest mass data points rather than the extrapolation in m π .
Although additional statistics are necessary to reduce the statistical uncertainties of the results to match experimental ones, the current precision is a major advance. We have achieved 5% precision for the Dirac radius (r 2 1 ) v , 10% precision for the anomalous magnetic moment κ v and the Pauli radius κ v (r 2 2 ) v , and 15% precision for the isovector quark momentum fraction x u−d (see Tab. 2). Remarkably, the precision of (r 2 1 ) v is already comparable to the discrepancy between the two contradictory experimental values, and future improved calculations will contribute to a resolution of this puzzle.
It is useful to comment on the role of ChPT in this work. As emphasized above, because of the small extrapolation distance to the physical pion mass, our primary results for observables are independent of ChPT and its range of convergence. However, our data provide interesting input to the study of it. We use standard parameters from ChPT phenomenology given in Tab. 2 and do not attempt a global fit to determine them. When fitting data in the two ranges m π 250 MeV and m π 350 MeV, we find consistency for all observables except g A . In the case of (r 2 1 ) v , we find that ChPT fits in the two ranges become inconsistent if contributions of the ∆ are removed (not shown in    1). Both observations suggest that ChPT as we are using it is working well in this range of pion masses. Hence, it is one useful diagnostic for the presence of bias in the calculation of g A .
Lattice results are also subject to other systematic errors such as dependence on finite lattice spacing a and finite volume. However, we expect that these errors are negligible at the present level of statistical precision. The finite lattice spacing effects, including chiral symmetry breaking, should vanish in the limit a → 0. Discretization effects of this QCD action in the hadron spectrum have been carefully studied in Ref. [17] and they were found to be quite mild. In addition, our results with two values of a near m π ∼ 300 MeV shown in Figs. 1-4 agree within statistics.
Finite volume effects are widely believed to be negligible if the spatial box size L s 4m −1 π . Most of our ensembles satisfy this criterion, including the lightest data point m π = 149 MeV. Two ensembles at m π ≈ 250 MeV serve as a direct check that finite volume effects are insignificant: we have compared data from lattices with m π L s = 3.6 (24 3 × 48) and m π L s = 4.8 (32 3 × 48), otherwise completely identical, and find no statistically significant discrepancies in any of the reported quantities.
In the case of g A , we have noted that internal inconsistencies such as the lack of agreement between the two ranges of chiral fits and the qualitative difference between the behavior of g A (∆t) at m π = 254 MeV and at lower pion masses indicate some source of bias that is not present for the other observables. One possible source of bias could be the presence of thermal states, as discussed in connection with Fig. 6. Another possibility is the breaking of chiral symmetry by Wilson fermions. Although we have verified that finite a effects are smaller than the statistical uncertainty at m π = 300 MeV, it is possible that chiral symmetry breaking leads to a progressively stronger systematic bias at a given a as the pion  mass is lowered, and that this disproportionately affects the axial charge. Yet another possibility is very long-range autocorrelations as recently reported for a calculation [43] of g A using domain wall fermions with m π down to 170 MeV, where values for g A averaged over the first and second halves of the ensemble differed by almost four standard deviations while no other observables showed similar autocorrelations. While the physical cause has not yet been determined, slow equilibration across topological sectors could conceivably affect the axial charge. Although we do not observe such anomalies in our calculation, they may still lead to underestimated stochastic error of g A and apparent deviation from experiment. We also note that a recent preliminary report on work with twisted mass fermions using a 48 3 × 96 ensemble at the physical pion mass [44] showed agreement with the experimental value for g A . Although it had no study or removal of excitedstate effects, a small volume corresponding to m π L = 3.0, N f = 2, and included no data on electromagnetic form factors or radii, we expect that further results from that work will be useful in helping to resolve the problem with g A .
More lattice spacings and further volumes and temporal extents at pion masses below 200 MeV are required to obtain strong bounds on the systematic uncertainties of our calculations, and in particular to understand the unique behavior of the axial charge. With appropriately increased statistics, our control over excited states could be cross-checked by further varying the nucleon source-sink separation and by comparing against other excited state removal techniques.
It is a significant advance that Lattice QCD, for four key properties of nucleon structure, now is in good agreement with experimental results.