Nucleon Isovector Axial Form Factors

We present results for the isovector axial vector form factors obtained using thirteen 2+1+1-flavor highly improved staggered quark (HISQ) ensembles generated by the MILC collaboration. The calculation of nucleon two- and three-point correlation functions has been done using Wilson-clover fermions. In the analysis of these data, we quantify the sensitivity of the results to strategies used for removing excited state contamination and invoke the partially conserved axial current relation between the form factors to choose between them. Our data driven analysis includes removing contributions from multihadron $N \pi$ states that make significant contributions. Our final results are: $g_A = 1.292 (53)_\text{stat}\,(24)_\text{sys}$ for the axial charge; $g_S = 1.085 (50)_\text{stat}\, (103)_\text{sys}$ and $g_T = 0.991 (21)_\text{stat}\, (10)_\text{sys}$ for the scalar and tensor charges; $\langle r_A^2 \rangle = 0.439 (56)_\text{stat} (34)_\text{sys}$ fm${}^2$ for the mean squared axial charge radius, $g_P^\ast = 9.03(47)_\text{stat}(42)_\text{sys} $ for the induced pseudoscalar charge; and $g_{\pi NN} = 14.14(81)_\text{stat}(85)_\text{sys}$ for the pion-nucleon coupling. We also provide a parameterization of the axial form factor $G_A(Q^2)$ over the range $0 \le Q^2 \le 1$ GeV${}^2$ for use in phenomenology and a comparison with other lattice determinations. We find that the various lattice data agree within 10\% but are significantly different from the extraction of $G_A(Q^2)$ from the $\nu$-deuterium scattering data.


I. INTRODUCTION
In ongoing neutrino scattering experiments (T2K, NOvA, MINERvA, MicroBooNE, SBN), the lack of precise reconstruction of the final state of the struck nucleus gives rise to uncertainty in the cross-section.Theoretical calculations of the cross-section for targets, such as 12 C, 16 O, and 40 Ar, being used in experiments take as input axial-vector form factor of the nucleon and build in nuclear effects using nuclear many body theory [1][2][3].Both of these steps, calculating nucleon axial form factors using lattice QCD and including nuclear effects using many-body theory, have uncertainties.Incorporating nuclear effects involves modeling of the complex physical phenomena (quasi-elastic, resonance, deep inelastic scattering) that contribute when considering incoming neutrino energies up to 5 GeV relevant for ongoing and future (DUNE) experiments.These complex phenomena make it hard to reconstruct the incident neutrino energy or the cross-section from the imprecise knowledge of the final state of the struck nucleus.On the other hand, the cleanest experimental measurements of the nucleon axial-vector form factor from scattering neutrinos off liquid hydrogen targets are not being carried out due to safety concerns.
The MINERνA experiment [4] has recently shown that the axial-vector form factor of the nucleon can be extracted from the charged current elastic scattering process ν µ H → µ + n in which the free proton in hydrogen (H) gets converted into a neutron.This opens the door to direct measurements of the nucleon axial-vector form factor without the need for extraction from scattering off nuclei, whose analysis involves nuclear corrections which have unresolved systematics, and making detailed comparisons with predictions from lattice QCD.For example, our result for the axial charge radius, ⟨r 2  A ⟩ = 0.439(56) stat (34) sys fm 2 , given in Eq. (36), is consistent, within one combined sigma, with the MINERνA result, ⟨r 2  A ⟩ = 0.73 ± 0.17 fm.Similarly, recent advances in simulations of lattice QCD have enabled robust results for the nucleon charges that have been reviewed by the Flavor Lattice Averaging Group (FLAG) in their 2019 and 2021 reports [5,6]).Results for axial-vector form factors [7] are now available with ≲ 10% uncertainty as we show in this work.At the same time, there continues to be progress in nuclear many-body theory for the calculation of the neutrinonucleus cross-section [1].
In this work, we present lattice QCD results for the isovector axial, G A (Q 2 ), induced pseudoscalar, G P (Q 2 ), and pesudoscalar G P (Q 2 ) form factors, the axial, scalar and tensor isovector charges g u−d A , g u−d S and g u−d T , the axial charge radius squared ⟨r 2  A ⟩, the induced pseudoscalar coupling g * P , and the pion-nucleon coupling g πN N .The calculation has been done using thirteen ensembles generated with 2+1+1-flavors of highly improved stag-gered quarks (HISQ) by the MILC collaboration [8].The construction of nucleon two-and three-point correlation functions has been done using Wilson-clover fermions as described in [9].The analysis of the data generated using this clover-on-HISQ formulation includes a study of excited state contributions (ESC) in the extraction of ground state matrix elements (GSME) and a simultaneous chiral-continuum-finite-volume (CCFV) fit to obtain results at the physical point, which throughout the paper will be defined as taking the continuum (a = 0) and infinite volume (M π L → ∞) limits at physical light quark masses in the isospin symmetric limit, m u = m d , which are set using the neutral pion mass (M π 0 = 135 MeV).The masses of the strange and charm quarks in the lattice generation have been tuned to be close to their physical values in each of the thirteen ensembles [8].
The three form factors G A (Q 2 ), G P (Q 2 ) and G P (Q 2 ) must, up to discretization errors, satisfy the constraint in Eq. ( 16) imposed by the partially conserved axial current (PCAC) relation ∂ µ A µ = 2mP between the axial and pseudoscalar currents.The decomposition of the matrix elements (ME) into form factors, given in Eqs. ( 1) and (2), assumes that they are GSME.Post-facto, deviations from the PCAC relation larger than those expected due to lattice discretization artifacts are indicative of residual ESC in the extraction of ME from the spectral decomposition of the three-point correlation functions.They point to the need for reevaluation of the key inputs in this analysis-the number and energies of the excited states that contribute significantly to the threepoint functions.The strategies used to remove ESC are described in Secs.II B and V, and the use of the PCAC relation to evaluate how well ESC have been controlled is discussed in Sec.II C.
In Ref. [10], we showed that the standard method of taking the excited-state spectrum from fits to the nucleon two-point correlation function to analyze the three-point functions lead to form factors that fail the PCAC test by almost a factor of two on the physical pion mass ensemble a09m130W , and identified the cause as enhanced contributions to ME from multihadron, N π, excited states that have mass gaps smaller than of radial excitations [11,12].These contributions had been missed in all prior lattice calculations.Including N π excited states in the analysis reduces the disagreement to within 10%, an amount that can be attributed to discretization effects.In this paper, we include N π states in the analysis of all thirteen ensembles described in Table I.Data from various analyses discussed in Secs.III A, III B, and IV B are then extrapolated to the physical point using simultaneous CCFV fits and results compared to understand systematics..
In order to extract g A and ⟨r 2 A ⟩, we parameterize the Q 2 behavior of G A (Q 2 ) using the dipole and the model independent z-expansion.We find that the dipole ansatz does not provide a good fit and our final results are obtained using the model independent z-expansion.We show that the pion-pole dominance (PPD) hypothesis, Eq. (20), tracks the improvement observed in satisfying the PCAC relation when N π states are included in the analysis.We, therefore, use it to parameterize G P (Q 2 ) and extract g * P and g πN N in Sec.IV.Similarly, the analysis of the ESC in isovector charges extracted from the forward matrix elements is carried out using information from both the 2-and 3-point correlation functions and the noninteracting energy of the lowest N π state.
Our final result for the axial form factor, parameterized using the z 2 truncation, is given in Eq. (34); the axial charge obtained from extrapolating it to Q 2 = 0 in Eq. (30), and the charge radius in Eq. (31).The results for the induced pseudoscalar charge g * P and g πN N are given in Eqs. ( 44) and (45).Lastly, the results for the three isovector charges g u−d A,S,T from the forward matrix elements are given in Eq. ( 50).
This paper is organized as follows.In Sec.II, we briefly review the notation and the methodology for the extraction of the three form factors: the axial, G A , the induced pseudoscalar G P , and the pseudoscalar, G P , from matrix elements of the axial and pseudoscalar currents within ground state nucleons.In Sec.II B, we explain the three strategies used to remove the ESC to the three-point functions.The analysis of the form factors with respect to how well they satisfy the relations imposed between them by PCAC relation and the PPD hypothesis is presented in Sec.II C. Based on this analysis, we present our understanding of the excited states that contribute in Sec.II D. The parameterization of the axial form factors as a function of Q 2 and the extraction of the axial charge g A and the charge radius squared ⟨r 2  A ⟩ is carried out in Sec.III.Parameterization of the induced pseudoscalar form factor, G P , and the extraction of the induced pseudoscalar coupling g * P and the pion-nucleon coupling g πN N is carried out in Sec.IV.The calculation of the isovector charges g u−d A,S,T from forward matrix elements is described in Sec.V. A summary of our results and a comparison with previous lattice calculations is presented in the concluding Sec.VI.Six appendices give further details of the analysis and the data.

II. METHODOLOGY FOR EXTRACTING THE FORM FACTORS
The matrix elements of the axial A µ = ūγ µ γ 5 d and pseudoscalar P = ūγ 5 d currents between the ground state of the nucleon can be decomposed, in the isospin symmetric limit, into the axial G A , induced pseudoscalar G P , and pseudoscalar G P form factors as where u(⃗ p i ) is the nucleon spinor with momentum ⃗ p i , q = p f − p i is the momentum transferred by the current, 0.1510 (20) 306.9 (5) I.The parameters of the 2+1+1-flavor HISQ ensembles generated by the MILC collaboration and analyzed in this study are quoted from Ref. [8].On two ensembles, a06m310 and a06m220, a second set of calculations labeled a06m310W and a06m220W , have been done with a larger smearing size σ as described in Ref. [13].In this clover-on-HISQ study, all fits are made versus M val π , which is tuned to be close to the Goldstone pion mass M sea π .The finite-size effects are analyzed in terms of M val π L. Columns 7-10 give the values of the source-sink separation τ used in the calculation of the three-point functions, the number of configurations analyzed, and the number of measurements made using the high precision (HP) and the low precision (LP) truncation of the inversion of the Wilson-clover operator [14].The last column gives the value of Q 2 |max for two cuts, n 2 ≤ 6 (n 2 ≤ 5 for the four a06m310 and a06m220 ensembles) and n 2 ≤ 10 used in the analysis.The full set of Q 2 values simulated on each ensemble are given in Tables X-XXII. 2 is the spacelike four momentum squared transferred.The spinor normalization used is The process of obtaining the GSME needed in Eqs. ( 1) and ( 2) from fits to 2-and 3-point correlation functions is described next.
A. Two-and three-point correlation functions The lattice calculation starts with the measurement and analysis of the two-and three-point correlation functions C 2pt (p; τ ) and C J (q, t; τ ) constructed using the nucleon interpolating operator N , where the ± sign give positive parity states propagating forward/backward in time.The spectral decompositions of the two time-ordered correlation functions are and where J Γ = A µ or P is the quark bilinear current inserted at time t with momentum q, and |Ω⟩ is the vacuum state.
In our set up, the nucleon state |j⟩ is, by construction, projected to zero momentum, i.e., p j = (M, 0), whereas ⟨i ′ | is projected onto definite momentum p i = (E, p) with p = −q by momentum conservation.Consequently, the states on the two sides of the inserted operator J are different for all q ̸ = 0.The prime in ⟨i ′ | indicates that this state can have nonzero momentum.For large time separations, τ and τ −t, only the ground state contributes and the GSME, ⟨0 ′ |J|0⟩, whose Lorentz covariant decomposition is given in Eqs. ( 1) and (2), can be extracted reliably.Assuming this is the case, and choosing the nucleon spin projection to be in the "3" direction, the decompositions become where i ∈ 1, 2, 3 and the kinematic factor K −1 ≡ 2E(E + M ).These correlation functions are complex valued, and the signal, for the CP symmetric theory, is in Im C Ai , Re C A4 , and Re C P .For a given q 2 ̸ = 0, the correlators with momentum combinations q = (2π/L)n ≡ (2π/L)(n 1 , n 2 , n 3 ) related by cubic symmetry can be averaged to increase the statistics before making fits.We construct the following averages A µ and P : where sgn(x) = x/|x| is a sign function with sgn(0) = 0, /n 2 = N/3, q = (2π/L)n, β ≡ q3=0 1, and N ≡ q 1 is the number of equivalent (under the cubic group) momenta averaged.
The pseudoscalar form factor, G P , is given uniquely by Eq. (14).For a subset of momenta, G A and G P are determined uniquely from Eqs. (10) and (12).In general, we solve the over-determined system of equations, Eqs. ( 10)- (13).Of these, correlators A 3,L and A 4 are nonvanishing for all q, and are thus sufficient to solve for G A and G P .In practice, the A 4 correlator has a poor signal and is dominated by excited states contributions, which we exploit to determine the relevant low-lying excited states.These turn out to be towers of multihadron N π and N ππ states.We find that including these states in fits to the spectral decompositions given in Eqs. ( 5) and ( 6) is essential for extracting the GSME.With the GSME in hand, the form factors G A and G P are determined using Eqs.( 10)- (12).

B. Strategies to extract ground state matrix elements
Calculations of nucleon correlation functions face two key challenges.First, the statistical signal-to-noise ratio decays exponentially with the source-sink separation τ as e −(M N −1.5Mπ)τ [15,16].This limits current measurements of two-point (three-point) functions to ≲ 2 (≲ 1.5) fm.Second, at these τ , the residual contribution of many theoretically allowed radial and multihadron excited states can be significant.These states arise because the standard nucleon interpolating operator N , defined in Eq. ( 4), used to construct the correlation functions in Eqs. ( 5) and ( 6), couples to nucleons and all its excitations with positive parity including multihadron states, the lowest of which are N (p)π(−p) and N (0)π(0)π(0).The goal is to remove the contributions of all these excited states to three-point functions to obtain the GSME, ⟨0 ′ |J|0⟩, which we do by fitting the averaged correlators A µ and P using Eq. ( 6).
An important note applicable to all fits used to remove excited states.For all our ensembles, the energy of the two lowest in these tower of positive parity states, N (p = 1)π(p = −1) and N (0)π(0)π(0), are approximately the same.Since fits to Eqs. ( 5) and ( 6) depend only on E i and not on the nature of the states, the contribution of both states is taken into account using the E 1 calculated for either state.Thus, the reader should understand that the contribution of both states are being included when, for brevity, we say N π state.
To extract the GSME, we need to address two questions: (i) which excited states make large contributions and (ii) how large is this contribution to various observables.The most direct (and statistically the best motivated assuming that a common set of states dominate the ESC in all correlators) way to get the ground-state matrix element that addresses these questions is to simultaneously fit, with the full covariance matrix, all five nucleon three-point functions including two or more excited states, and then solve the linear set of Eqs. ( 10)- (14).
In general, simultaneous fits to the current data (3point or a combination of 2-and 3-point) do not resolve the excited states, i.e., there are large regions of parameter space where fits give similar χ 2 /dof .The one exception is fits to the correlation functions ⟨N A 4 N ⟩ that, as discussed below, play a central role in our analysis as they expose the large contribution of N π states.Unfortunately, even 2-state simultaneous fits to all five nucleon three-point functions are not stable for all momentum channels and ensembles.We, therefore, resort to taking the energies and amplitudes, especially those of the ground state, from separate, but within a single overall jackknife process, fits to the 2-point function.
The analysis of the nucleon two-point functions and the extraction of the spectrum is presented in Appendix A and the extrapolation of the data for the nucleon mass, M N , to the continuum limit in Appendix B. Results for the excited state parameters, i.e., the energies E i (q), the masses M i , and the amplitudes A i , have large uncertainty.For example, in a four state fit, there is a large region of parameter space where fits have similar χ 2 /dof .
In short, statistical precision of current data does not allow simultaneous fits to all five nucleon three-point functions using a 3-state (or higher) fit with the full covariance matrix.Even a robust determination of the energies, E i , and amplitudes, A i , of excited states that make significant contributions from fits to 2-point functions is lacking.Our best approach is a hybrid of using information from fits to 2-and 3-point functions.To this end, we construct three strategies with different estimates of excited-state parameters to fit the three-point data using Eq. ( 6).These are described next.
In the standard approach, labeled S 2pt , we take E i , M j , A ′ 0 and A 0 from 4-state fits to C 2pt , and input them into an m-state truncation (m ≤ n) of Eq. ( 6) to extract the matrix element ⟨0 ′ |J|0⟩.In this paper, we truncate the spectral decompositions given in Eqs. ( 6) and ( 5) at m = 3 and n = 4, respectively.
The second strategy, labeled S A4 , was proposed in Ref. [10].Again E 0 , M 0 , A ′ 0 and A 0 are taken from 4-state fits to C 2pt , however, E 1 and M 1 are determined from two-state fits to the three-point correlator A 4 .The output E 1 and M 1 are then fed into the fits to the other four correlation functions defined in Eqs. ( 10)-( 12), ( 14).This strategy assumes that the same [first] excited state parameters apply to all five correlation functions, and these are given by fits to A 4 .
The third strategy S sim is similar to S A4 except that E 1 and M 1 are outputs of simultaneous two-state fits to all five three-point correlators defined in Eqs. ( 10)- (14).It is, from a statistical point of view, better motivated than S A4 because the underlying assumption in both cases is that the same excited states contribute to all five correlators.It avoids the two-step procedure used in S A4 , i.e., first obtain E 1 and M 1 from fits to A 4 and then use them in fits to the other four correlators.In S sim , we used the averaged correlator A xy = (A 1 + A 2 )/2, since these two correlators are equivalent under cubic rotational symmetry, thus reducing the number of correlators fit simultaneously to four.
We used the full covariance matrix for all fits to the 2point and 3-point functions with the S 2pt and S A4 strategies.In the S sim strategy, the covariance matrix was restricted to be block diagonal in each correlation function.In the O(1000) fits made to remove ESC (ensembles ⊗ Q 2 values ⊗ correlation functions ⊗ strategies) the selection of parameters was done individually due to the large differences in ESC behavior versus ensembles, Q 2 values, and correlation functions.
We emphasize from the very outset that in all fits with each of the three strategies, the excited state amplitudes, A (′) i and A j , are not needed since these arise only in the combinations |A ′ i ||A j |⟨i ′ |J|j⟩, which are fit parameters but are not used thereafter in the analysis.Second, the ground state parameters, M 0 , E 0 , A ′ 0 and A 0 are common for all three strategies and are taken from four-state fits to the two-point correlators.
The unrenormalized values of the three form factors at various values of Q 2 simulated, for each of the three strategies and for the 13 ensembles are given in Tables X-XXII in Appendix C. The size of the effect of N π state can already be inferred from the difference between the S 2pt and S sim data even though S sim includes only the lowest N (−1)π(1) state in the fit.Overall, this comparison shows that the contribution of the N π state to G P and G P is enhanced, reaching ∼ 45% at the physical pion mass.The roughly 5% effect observed in G A is important phenomenologically and needs to be made more precise.Later, in Sec.II C, we choose the S sim strategy to present the final results based on the three form factors satisfying the PCAC relation.
A comparison of G A (Q 2 ) and the combination N ), which should be proportional to G A according to the PPD hypothesis, obtained using the three strategies S 2pt , S A4 , and S sim , is shown in Fig. 1. Results for both form factors are consistent between S A4 and S sim for each of thirteen ensembles with errors from S sim being slightly larger.On the other hand G P (and G P ) from strategy S 2pt show noticeable differences that increase as Q 2 → 0 and M π → 135 MeV (see also the data in Tables X-XXII in Appendix C).This effect is correlated with the increase in the difference between ∆E 2pt 1 (used in S 2pt fits) compared to ∆E A4 1 and ∆M A4 1 (output of S sim fits) in the same two limits as shown later in Fig. 6.Also, from Eq. ( 6) it is obvious that a smaller ∆E 1 implies a larger ESC.
The same data for N ) from the 13 ensembles with the S sim strategy are plotted in Fig. 2. Remarkably, they show no significant variation with respect to the lattice spacing a or M π except for a 1σ lower values on the a06m135 ensemble, which we identify to be statistics limited.
In Appendix C we summarize why, with the methodology for momentum insertion through the operator used in this work, improving the lattice calculations (M π → 135 MeV, increasing M π L > 4, and reducing a) will increasingly give data at Q 2 < 0.5 GeV 2 .Even in this work, most of the data for Q 2 > 0.7 GeV 2 comes from the M π ≈ 310 MeV ensembles.If the optimistic scenario presented by the current data, mild dependence on {a, M π } as shown in Fig. 2 and in Ref. [7], holds then one will have confidence in the final result for G A also for 0.5 ≲ Q 2 ≲ 1.5 GeV 2 even though its extraction in that region will be from data on M π > 200 MeV ensembles.To push precision to Q 2 ∼ 5 GeV 2 to meet the needs of the DUNE experiment will need further algorithmic developments and much higher statistics data.
A detailed analysis of the extrapolation of these G A to Q 2 → 0 to get the axial charge g A is presented in Sec.III.The marked improvement in satisfying the PCAC relation and the PPD hypothesis shown by the G A , G P , and G P obtained with the S A4 and S sim strategies (which include the N π state) is discussed next in Sec.II C.

C. The PCAC relation and pion-pole dominance
In this section, we evaluate how well the form factors from the three strategies, tabulated in Appendix C, 1. Data from the 13 ensembles with a ≈ 0.15, 0.12, 0.09, and 0.06 fm for the unrenormalized axial form factor GA(Q 2 ) (first and third columns) and N ) (second and fourth columns) plotted versus Q 2 .Each panel compares the data obtained using the three strategies S2pt, SA4, and Ssim for controlling ESC.
N ) (right) from the 13 HISQ ensembles analyzed.No significant variation with respect to the lattice spacing a or Mπ is observed except for the a06m135 ensemble, which we consider statistics limited.satisfy the PCAC relation, which in terms of the bare axial, A µ (x), and pseudoscalar, P (x), currents is: where the quark mass parameter m ≡ Z m m ud Z P Z −1 A includes all the renormalization factors, and m ud = (m u + m d )/2 = m l is the light quark mass in the isospin symmetric limit.Using the decomposition in Eqs. ( 1) and (2) of GSME, the PCAC relation requires that the three form factors G A , G P , and G P satisfy, up to discretization errors, the relation which we rewrite as with The PPD hypothesis relates G P to G A as Tests of whether the form factors satisfy the PCAC (R 1 + R 2 = 1) and PPD (R 3 = 1) relations are presented in Fig. 3 and Fig. 4, respectively.Data with the S 2pt strategy show about 10% deviation for both the PPD and PCAC relations for Q 2 > 0.3 GeV 2 .Below it, the deviation grows to about 40% at the lowest Q 2 point on the two physical pion mass ensembles.See also the discussion in Appendix C on the differences in data for the form factors obtained using the three ESC strategies.
There is a very significant reduction in the deviations for both the S A4 and S sim strategies for Q 2 < 0.3 GeV 2 .In fact, except for three M π ≈ 220 MeV ensembles, data below Q 2 = 1 GeV 2 is essentially independent of Q 2 and the deviations from unity and the variations between ensembles is in most cases within about 5%, which can be due to possible discretization errors.The differences between data with S sim and S A4 are much smaller.Also, the improvement in the PPD relation, Eq. ( 20), tracks that in PCAC, Eq. (17).We point out a caveat in our clover-on-HISQ calculation of the quark mass m used in Eq. ( 16).For four ensembles, a12m310, a09m130W , a06m220, and a06m135 we have calculated m 2pt using the following ratio of pion two-point correlators, For the other ensembles, the data for these twopoint functions were not collected, so we use the HISQ sea quark mass am sea ud for m since for staggered fermions, in fact all lattice fermions with chiral symme-try, Z m Z P Z −1 A = 1.These quark masses are given in Table II and we find that m 2pt is 5 -20% larger than am sea ud , which is not unexpected for our clover-on-HISQ calculation.Noting that R 2 ≈ 0.5R 1 (see Fig. [15] in Ref. [9]), such a 20% systematic error would increase R 1 + R 2 by about 7%.This would bring the data from the physical mass ensembles, a09m130W and a06m135, in better agreement but would not alter our conclusion that form factors obtained with S A4 and S sim strategies show better agreement with the PCAC relation compared to S 2pt .Also, m does not enter in the PPD relation, Eq. ( 20), and the deviation from unity of the PPD relation with S A4 and S sim data is observed to be smaller than seen in the PCAC relation as shown in Fig. 4. Equally important, this caveat does not impact the extraction of individual form factors or their subsequent analysis since m only enters in the test of how well the three form form factors satisfy the PCAC relation, Eq (16).
We further examine whether the deviation from unity in Fig. 3 at small Q 2 is a discretization error.The O(a) improvement affects only the axial current, A µ → A µ + c A a∂ µ P , and adds to the left hand side in Eq. ( 16) the term −Q 2 ac A G P , i.e., under improvement, Eq. ( 16) can be written as where the improvement coefficient, c A , is typically O(10 −2 ) and negative.Thus, this effect is expected to be small for Q 2 < 1 GeV 2 , and will not change our conclusions.On the other hand, effects due to possible mistuning of the clover coefficient, c SW , and c A , and O(a 2 ) corrections are likely to increase with Q 2 .Similarly, artifacts are expected to increase with the quark mass since the improvement coefficient b m is not included The PPD relation (Eq.( 20)) can be derived from PCAC (Eq.( 16)) provided In this case, R 1 + R 2 = 1 would also imply R 3 = 1.In Fig. 5, we compare R 4 from the three strategies for all ensembles except a06m220W , a06m310, and a06m310W where G P is not available.We note a roughly linear increase in R 4 with Q 2 , which is consistent with the behavior observed in Ref. [7] and with the analysis of the Goldberger-Trieman discrepancy using χPT in Ref. [17].Lastly, we note that the data for R 4 from all three strategies, S 2pt , S A4 and S sim , overlap implying that the changes in G P and G P , both of which have a pion pole, between different treatments of ESC (S 2pt versus S A4 or S sim ) cancel in the ratio R 4 within our statistics.This observation supports our hypothesis that the same excited states contribute to all five correlation functions.15) and ( 16), requires R1 +R2 = 1 up to discretization errors.The dashed lines give the ±5% deviation band.

D. Excited States Spectrum
In Fig. 6 we show data for the energy gaps, ∆E 1 and ∆M 1 on the two sides of the operator insertion for the various ensembles, including the two physical pion mass ones, a09m130W and a06m135.The results for ∆E A 1 and ∆M A 1 , outputs of the simultaneous fits to all five correlators (insertions of A µ and P ) at a given momentum transfer p = 2πn/L overlap with the results ∆E A4 1 and ∆M A4 1 obtained from fits to just C A4 .This indicates that the energy gaps in the S sim fits are essentially controlled by C A4 .The momentum dependence of the data is consistent with the expectation that the relevant excited states on the two sides are N (n) + π(−n) and N (0) + π(−n).This is based on the rough agreement between the data and the corresponding noninteracting energies of these states, ∆M 1 and ∆E 1 , shown by the dashed red and blue lines, respectively, and consistent with the PPD hypothesis that the currents inject a pion with momentum ⃗ q.
The data with open circles in Figure 6 are the energy gaps ∆E 2pt 1 obtained from the nucleon two-point correlators.These are roughly independent of momentum and larger than those from S A4 or S sim fits, especially Results for the ratio R4 = (GP / GP ) × (4 mMN /M 2 π ).For the pion-pole dominance hypothesis to be exact (derivable from the PCAC relation), R4 should be unity independent of Q 2 .The data show an approximate linear increase with Q 2 , which is consistent with the Goldberger-Trieman discrepancey as discussed in Ref. [7,17].
for the smaller Q 2 points.The difference increases as Q 2 → 0 and M π → 0. This behavior is consistent with ∆E 2pt 1 corresponding to a mixture of radial and higher multiparticle excitation whereas the energy of the intermediate excited states identified by the S A4 and S sim fits, N (n) + π(−n) and N (0) + π(−n), decreases with decreasing n and M π .A word of caution when making these identifications: it is very important to qualify that the ∆E 1 and ∆M 1 from the two-state fits in S A4 and S sim strategies are effectively trying to account for all the intermediate states that make significant contributions and not just the lowest or the most intuitive ones.Given the size of the effect, identifying and improving control over all the excited states that make significant contribution to these correlation functions will be key to obtaining, in the future, higher precision results for the form factors.

E. Renormalization Constant ZA
The renormalization constant Z A for the axial current needed for the form factors G A and G P and the charges g A , g * P and g πN N was determined nonperturbatively using the RI-sMOM intermediate scheme in Ref. [13].We use the results given in Table V  In this section, we present the analysis of G A (Q 2 ) without including the values of the axial charge g A obtained from the forward matrix element.Its extraction is discussed separately in Sec.V.This is done to keep the two extractions of g A -from the forward matrix element ((50)) and by extrapolating 36))-separate.The final result, given in Table VI, is taken to be the average of the two.
The axial-vector form factor G A (Q 2 ) can be parameterized, near Q 2 = 0, by the axial charge g A and the axial charge radius squared ⟨r 2 A ⟩: where g A ≡ G A (0) and To extract these from lattice data obtained at ).Among the various parameterizations, we study the dipole ansatz and the model-independent zexpansion.The dipole ansatz has two free parameters, the axial charge g A and the axial mass M A .The z-expansion is the series in terms of the variable that maps the kinematically allowed analytical region Q 2 ≥ 0 to that within a unit circle, |z| < 1 [18].The parameter t 0 is discussed later.For sufficiently small z, fits with the first few terms should suffice.In practice, to stabilize the fits we impose the condition |a k | ≤ 5 for all z k≥1 truncations [18].With increasingly precise data over a sufficiently large range of Q 2 , our goal is to demonstrate that a data-driven choice can be made between the various parameterizations.
In the data presented here, the statistical signal is good for momentum transfer with n 2 ≤ 6 but often poor in the four points with 8 ≤ n 2 ≤ 10.To test the stability of the dipole and z k fits with such few points, we compare the output of the fits to the lowest six versus all ten Q 2 points on nine ensembles where data on all ten Q 2 values exist.Observing consistency, the final results are taken from fits to six (five in 4 cases) points, i.e., to data up to the max given in Table I.This implies that results for form factors presented for Q 2 ≳ 0.5 GeV 2 come mainly from the M π ∼ 310 MeV ensembles, and, a priori, could have an large systematic uncertainty.In practice, however, the observed weak dependence of form factors on M π (see Fig. 2 and similar result from a clover-on-clover calculation presented in Ref. [7]) suggests that the uncertainty is likely within the quoted errors.Readers should, nevertheless, keep in mind that results at Q 2 ≳ 0.5 GeV 2 are based mainly on data from the M π ≈ 310 MeV ensembles.
Estimates of G A from the dipole fit to data from the three ESC strategies are consistent, however, on six ensembles these dipole fits to S sim and S A4 data have poor p-values.Our evaluation of the failure is that the dipole ansatz does not have enough parameters to capture the change in the curvature over the range of Q 2 studied.A consequence is that estimates for g A and ⟨r 2 A ⟩ are smaller than those from z 2 fits to the same S sim and S A4 data.Since agreement with PCAC is essential and S sim data do the best job, while the dipole fits are poor, we rule out the dipole ansatz.Henceforth, our final results are obtained using the z-expansion, and the dipole results are given only for comparison.
In the z-expansion fits, the free parameter t 0 in Eq. ( 28) is used to adjust the maximum value of z within the unit circle |z| ≤ 1.We take t 0 = 0.4, 0.2, 0.12 for M π ≈ 310, 220, 130 MeV ensembles, which gives |z| ≲ 0.2.We have checked that using t 0 = 0 does not change the fits or the values significantly.
To ensure that the form factors satisfy the expected 1/Q 4 perturbative behavior in the limit Q 2 → ∞, sum rules can be imposed as was done in Ref. [9].However, to obtain the behavior near Q 2 = 0 from six or ten data points with Q 2 max ≈ 1 GeV 2 , we choose to make fits without the sum rules [19], i.e., to not increase the weight of the larger error high Q 2 points by imposing the sum rules.The z 1 and z 2 fits to G A (Q 2 ) from the S sim strategy are shown in Fig. 13.The resulting bare axial charge g A ≡ G A (0) and the charge radius squared ⟨r 2 A ⟩ from the z 2 fits are shown in Fig. 7, and the data are summarized in Table XXIII in Appendix D. From these data and the z-expansion fits, we conclude the following: • There is agreement in results between z 2 and z 3 fits in all cases.To account for the small curvature observed in the data shown in Fig. 13 and yet avoid ∆X/Mπ ∆X/Mπ ∆X/Mπ ∆X/Mπ ∆X/Mπ ∆X/Mπ ∆X/Mπ ∆X/Mπ FIG. 6. Results for the energy (mass) gaps ∆E1 (∆M1) for the first excited state extracted from (i) simultaneous fits to axial three-point correlators C[Aµ] and the pseudoscalar correlator CP (Ssim strategy and labeled ∆E A 1 and ∆M A 1 ), and (ii) from fits to the C[A4] correlator (SA4 strategy and labeled ∆E A4 1 and ∆M A4 1 ).These mass gaps are compared with the first excited state energy ∆E 2pt 1 from four-state fits to the nucleon two-point correlator.Note that the difference between them (black circles versus blue triangles), and consequently the difference between the form factors extracted, increases as Mπ → 135 MeV and n 2 → 0 (equivalently Q 2 → 0).overparameterization, evaluated using the Akaika Information Criteria (AIC) [20], we will present final results with the z 2 truncation.
• The errors in the data from the two physical mass ensembles a09m130W and especially a06m135 are large and underscore the need for higher statistics.
Results for bare gA from the strategy Ssim and the differences ∆(SA4) = gA|S A4 − gA|S sim and ∆(S2pt) = gA|S 2pt − gA|S sim .To facilitate visualization of the spread, the errors plotted for ∆(X) are those in gA(X).Results are shown for the dipole fit labeled "D" and z 1,2,3 -truncations.(Right) Analogous results for ⟨r 2  A ⟩.
• Results for both g A and ⟨r 2 A ⟩ from both S A4 and S sim analyses overlap and increase in value as M π → 135 MeV.This increase is correlated with the increasing ESC of the N π state.
We take the final results from the S sim strategy in which a simultaneous fit is made to all five correlators and the form factors come closest to satisfying the PCAC relation as shown in Fig. 3.This is a 2-state fit and we find that stable 3-state fits require higher statistics.Thus, with the current data, we do not have a reliable way of estimating the systematic uncertainty associated with possible residual ESC.
The analysis of g A obtained from the forward matrix elements is postponed to Sec.V.

A. Extrapolation of gA and ⟨r 2
A ⟩ to the Physical Point Extrapolation of the renormalized axial charge g A and the axial charge radius squared ⟨r 2  A ⟩ to the physical point (a → 0, M π → 135 MeV, L → ∞) is performed using a simultaneous CCFV fit keeping only the lowest order corrections in the ansatz where Y = ⟨r 2 A ⟩ or g A and {b Y i } denote the corresponding set of fit parameters.The discretization artifacts are taken to start at O(a) since the clover action used is only tadpole improved and the axial and pseudoscalar currents are unimproved [21].Similarly, only the lowest order term in the chiral expansion is kept to avoid overparameterization as data at only three values of the pion mass have been simulated.
We have performed four CCFV fits to (i) the full set of thirteen ensembles (13-pt); and three "12-pt" fits that exclude (ii) the coarsest lattice point a15m310, (iii) the smallest volume point a12m220S that also has large errors, and (iv) the point a06m135 that has large statistical errors and shows the largest difference from the other 12 points.The three 12-pt fits are used to estimate systematics due to discretization and finite volume effects, and the impact of the a06m135 point.Results of the 13-point CCFV extrapolation for g A and ⟨r 2  A ⟩ are summarized in Table III for six cases: the three strategies used for removing ESC, S sim , S A4 , and S 2pt , and the two Q 2 parameterizations, z 2 and dipole.The parameters of the 13-point CCFV extrapolation of the z 2 fit to the S sim data, used to get the final central values, are given in Table XXVI in Appendix E for both g A and ⟨r 2 A ⟩. Results for all the other cases can be constructed using the data for the form factors given in Tables X-XXII in Appendix C.

gA
The result for g A , taken from the 13-point CCFV extrapolation of z 2 fits to the S sim data, shown in Fig. 8, is The first error is the total analysis uncertainty given by the overall bootstrap process, and the next two are additional systematic uncertainties: (i) the difference between using z 2 and z 3 fits and (ii) the difference of this central value from the average of the three 12-point CCFV fits.The two systematics are added in quadrature to get the total systematic error given in the second  sim FIG. 8.The axial charge gA given by the 13-pt CCFV fit to Ssim data using the z 2 fit to GA(Q 2 ̸ = 0).The pink band in each panel gives the result of the CCFV fit (Eq.( 29)) versus the x-axis variable with the other two variables set to their physical values.The data points in each panel have been shifted in the other two variables using the same CCFV fit, however, the size of errors are not changed.The final result at the physical point is shown by the red cross.A ⟩ given by the 13-pt CCFV fit to data obtained using the z 2 fit to GA(Q 2 ̸ = 0) with the Ssim strategy.(Bottom panels) The 12-pt fit without the a06m135 point (open blue square in the three panels on the top) that has large errors.The pink band in each panel gives the result of the CCFV fit (Eq.( 29)) versus the x-axis variable with the other two variables set to their physical values.The data points in each panel have been shifted using the same CCFV fit, however, the size of errors are not changed.The final result at the physical point is shown by the red cross.line in Eq. (30).In Sec.V, this result is compared with an independent analysis of g A obtained from the forward matrix element, i.e., from the zero momentum correlator, C A3 (p = 0), as defined in Eq. (7).

⟨r 2
A ⟩ The CCFV fit to ⟨r 2 A ⟩, obtained from the S sim data with z 2 fit, is shown in Fig. 9 (top panels).It gives with the errors derived in the same way as for g A .For both g A and ⟨r 2 A ⟩, the largest dependence in the CCFV fit is on M 2 π for the S sim and S A4 strategies.This is a consequence of the increasing influence of the N π state as its mass gap decreases well below the N(1440) as M π → 135 MeV.In contrast, the S 2pt data, which do not include the N π state in the analysis, show mild dependence on all three variables {a, M π , M π L}.
Estimates from the 12-pt CCFV fit excluding the a06m135 point, shown in the bottom panels of Fig. 9, are consistent with the 13-point results.This is expected A ⟩ from the 13-point CCFV fit.Results are given for the z 2 and dipole fits to GA(Q 2 ̸ = 0), and for the three strategies used to control ESC.In each case, in addition to the central value and the total analysis error, the two systematic errors are the difference between the z 2 and z 3 estimates, and the difference from the 12-pt CCFV fits explained in the Sec.III A.
since the errors in the a06m135 point are large.Clearly, to further improve the estimates of both g A and ⟨r 2 A ⟩ requires much higher statistics data at small Q 2 on the physical pion-mass ensembles.
The Q 2 dependence of the axial form factor up to 1 GeV 2 , obtained at the physical point, is shown in Fig. 10 for the three strategies S 2pt , S A4 , and S sim .The pink band in these figures was obtained using the following three step process starting with the renormalized lattice data for G A (Q 2 )/g A , which on each of the thirteen ensembles are at different discrete values of Q 2 .First the data on each ensemble were fit using the z 2 -ansatz (see Eq. ( 27)) and the result is taken to specify G A (Q 2 )/g A for 0 < Q 2 ≤ 1 GeV 2 .Second, we chose a set of eleven Q 2 values evenly distributed over this range, and for the thirteen data points at each of these Q 2 values carry out a 13-point CCFV extrapolation using Eq. ( 29).The result was taken to be the physical point value of G A /g A at that Q 2 .In each of these CCFV fits, the thirteen points from the thirteen ensembles are uncorrelated as these are independent calculations.Lastly, these eleven extrapolated points were fit by the z 2 ansatz to obtain the final parameterization valid in the interval 0 ≤ Q 2 ≤ 1.0 GeV 2 and shown by the pink band in Fig. 10.The errors in the original lattice data were fully propagated through this three step process carried out within a single bootstrap setup.This gives the central value and error.Possible uncertainty due to incomplete removal of ESC or due to using only the leading order CCFV fit ansatz is to be estimated separately.
Figure 10 also shows the experimental bubble chamber data and the dipole ansatz with M A = 1.026 (21) GeV extracted from it (green band) [17].A recent z-expansion analysis of the ν-deuterium data [22] finds a ≈10X larger uncertainty.In our analysis, only the S 2pt data are roughly consistent with a dipole ansatz with Exp. pheno: 1.026 (21) 10. Results for GA/gA at the physical point for the three strategies S2pt, SA4 and Ssim (labeled "2pt", "A4", and "sim", respectively) used to control the excited-state contamination.
The three step process used to get these results shown by the pink band is described in the text.In each case, the error band represents the full analysis error for that strategy but with the value at Q 2 = 0 fixed to unity.The label 1⊕a⊕M 2 π ⊕ F V specifies that all 4 terms in the CCFV ansatz, Eq. ( 29), were kept.The experimental ν-deuterium data (gray crosses labeled Exp.) were provided by Ulf Meissner and the dipole result MA = 1.026 (21) GeV is taken from Ref. [17].This and the two other dipole fit with MA = 1.20 and 1.35 GeV are shown only for comparison.
1.30 GeV, however, the three form factors extracted using S 2pt fail to satisfy the PCAC relation.We, therefore, reemphasize that the dipole curves with M A = 1.026, 1.2 and 1.35 GeV are shown only for comparison.
The data from the S A4 and S sim strategies are consis- tent, and show a more rapid fall until Q 2 ≈ 0.3 GeV 2 , and give results roughly consistent with the dipole values M A = 1.026( 21) GeV (and ⟨r 2 A ⟩ = 0.444(28) fm 2 )) in Ref. [17].They then level out falling more slowly, however, note that for Q 2 > 0.5 GeV 2 our results are mainly from the heavier M π ≈ 310 MeV ensembles.On the heavier M π ensembles, the mass gaps in the N π analyses (blue triangles and red squares in Fig. 6) increase rapidly towards those from the S 2pt fit (black circles).Consequently, as shown in Fig. 10, the G A /g A from the S sim analysis moves towards the S 2pt result for Q 2 ≳ 0.5 GeV 2 .
To obtain data for Q 2 > 0.5 GeV 2 on physical pion mass ensembles with M π L > 4 requires simulations at much larger values of q where statistical and discretization errors are large with the methodology used in this work.A more promising method for generating data at large Q 2 is momentum smearing [23].Also, when including points with larger Q 2 , the z-expansion fits with and without sum-rules should be compared since it is not known, a priori, when the expected 1/Q 4 asymptotic behavior becomes significant.Alternately, as shown in Ref. [7], one can analyze the data using a Padé parameterization.Our fits to the final G A using the Padé ansatz ), which has the 1/Q 4 behavior built in, gave estimates consistent with Eqs. 30 and 31.
The coefficients b i in the CCFV fit using Eq. ( 29) are shown in Fig. 11 for the three strategies S 2pt , S A4 , and S sim .The coefficients b 1 (Q 2 ) and b 2 (Q 2 ) are similar within errors for S sim and S A4 , significantly different from zero, and qualitatively different from the case S 2pt .The b i for S sim and S A4 show a change in behavior at Q 2 ≈ 0.3 GeV 2 , coincident with the region where the curvature in G A changes as shown in Fig. 10.This could be due to the fact that most of the raw data controlling the parameterization at Q 2 ≳ 0.3 GeV 2 comes from the M π ≈ 220 and 310 MeV ensembles.On the other hand, support for the parameterization comes from the observation that the data, plotted in Fig. 2, do not show a significant variation versus {a, M π }.Cutting out the data with Q 2 > 0.3 GeV 2 to see if fits change, unfortunately eliminates most of the M π ≈ 220 and 310 MeV ensemblesthey do not have enough points to perform even a z 2 fit.We are, therefore, not able to resolve the reason for the change in parameterization around Q 2 ≈ 0.3 GeV 2 .
To provide our best parameterization of G A (Q 2 ) for phenomenology, we repeated the above 3-step procedure using the S sim strategy data.Again, the data after extrapolation to the continuum limit (at Q 2 = 0, 0.1, 0.2, ..., 1.0 GeV 2 ) were fit with a z 2 ansatz, t c = 9M result, shown in Fig. 12, has the parameterization with the correlation matrix: This fit gives which are consistent with the estimates in Eqs. ( 30) and ( 31), albeit with a ⟨r 2 A ⟩ larger by roughly 1σ.We also carried out this final z 2 fit setting t 0 = 0 in the definition of z.The results are g A = 1.282(54) While consistent, the coefficient c 2 in this fit is essentially undetermined.We, therefore, choose the results given in Eq. (34).
For our final results from the analysis of G A , we take the average weighted by the "stat" errors of values given in Eqs.(30), (31) and (34) For errors, we take the larger of the "stat" error and keep the "sys" errors given in Eqs. ( 30) and (31).
We now return to two related issues that arise because the majority of the data used to get the parameterization in Eq. ( 32) are at Q 2 ≲ 0.5 GeV 2 as shown in Fig. 1.The first is whether this parameterization is compatible with G A (Q 2 = ∞) = 0 since the sum rule constraints needed to build in the perturbative 1/Q 4 behavior have not been imposed?And the second is-how reliable is this G A for Q 2 > 0.5 GeV 2 since most of the data in this region are from the M π ≈ 310 MeV ensembles?Regarding the first issue, the parameterization in Eq. ( 32) at z = 1 These are consistent with zero within one sigma.For the second issue, data in Fig. 2 show that G A (Q 2 ) extracted on 12 ensembles show little dependence on {a, M π }.The one exception, a06m135, where the data lie about one sigma lower, has already been identified as statistics limited.As stated before, if future data continue to show little dependence on a and M π , then even data from the M π ≈ 310 MeV ensembles would provide a good approximation to the continuum G A and increase confidence in the result for Q 2 > 0.5 GeV 2 .Lastly, we made a z 2 fit to the same continuum extrapolated data but imposed a prior, G A (Q 2 = ∞) = 0 with width 0.3 based on the value in Eq. (37).The result is with G A = −0.07(10) and dG A /dz = −0.17(17) at z = 1.The main difference from the result in Eq. ( 32) is the tightening of the estimate of the z 2 term.Overall, in all these fits, the coefficients a 0 and a 1 , defined in Eq. ( 32) are stable, whereas higher precision data are needed to improve a 2 .
IV. COUPLINGS g * P AND gπNN FROM THE INDUCED PSEUDOSCALAR FORM FACTOR The induced pseudoscalar coupling g * P is defined as where m µ is the muon mass and Q * 2 = 0.88m 2 µ is the energy scale of muon capture.Similarly, the pion-nucleon coupling g πN N is obtained from the residue of G P (Q 2 ) at the pion pole, i.e., through the relation where F π = 92.9MeV is the pion decay constant.The function F P defined as is G P without the pion pole and should equal G A if PPD were exact.This requires R 3 ≡ F P /G A , plotted in Fig. 4, to be unity.Deviations for S 2pt are significant, while those for S sim and S A4 are one within the size of discretization errors and/or violations of PPD expected.
To extract g * P and g πN N from the lattice data, a parameterization of the Q 2 behavior of G P and F P was carried out.A comparison of the z 1 and z 2 fits to G A and F P from the S sim strategy is shown in Fig. 13 for the thirteen ensembles.Results from z 2 and z 3 fits are consistent, indicating convergence, while z 1 fits miss the small curvature seen.To avoid over parameterization, we again take the z 2 results as the central values.
A. Parameterization of GP (Q 2 ) and FP Based on the analysis of PCAC (see Figs. 3), we focus on the G P (Q 2 ) data from the S sim and S A4 strategies and again give the S 2pt results only for comparison.
We consider two ways to parameterize G P (Q 2 ), both of which build in the pion-pole dominance hypothesis.The first is the expansion where the c i (i = 0, 1, 2) are fit parameters, and we have kept as many terms in the polynomial as can be resolved by the data.Results for g * P , g πN N F π and g πN N F π /M N using this fit (labeled PD) to the S sim data are given in Table XXIV along with the χ 2 /DOF and p-value for the fits.
In the second way, we treat F P (Q 2 ) as an analytic function that can be fit using either the dipole ansatz (with free parameters F P (0) and M P ) or the z-expansion, Eq. ( 27), with z again defined by Eq. ( 28).Results for g * P , g πN N F π and g πN N F π /M N , from z 2 fits to F P (Q 2 ) obtained with the S sim strategy are given in Table XXV and agree with those in Table XXIV.
Results for the bare g * P from the S sim strategy for five Q 2 parameterizations of G P are shown in Figure 14 (left) along with their differences from results obtained using the S A4 and S 2pt strategies.The analogous results for unrenormalized g πN N F π are shown in Fig. 14 (right).
B. Extrapolation of g * P and gπNN to the Physical Point 1. g * P Renormalized g * P is extrapolated to the physical point in two ways.In the first method 2m µ M N F P (Q * 2 ) is extrapolated using the CCFV fit function given in Eq. ( 29) and multiplied by the the pion-pole factor at the physical point: In the second method, extrapolation of g * P is carried out by adding a pion-pole term, b π ), to the CCFV fit function in Eq. (29).The two methods give consistent estimates and their unweighted average is used to get the final results summarized in Table IV for each of the three strategies, S sim , S A4 , and S 2pt .
The error obtained from the overall analysis is quoted as the first "stat" uncertainty.The systematical errors associated with truncation of the z-expansion, and the largest difference of the central value from the three 12pt CCFV fits are quoted as the second and third errors.The difference between the two extrapolation methods described above is quoted as the fourth error.For the final result, we take the S sim with z 2 fits value: =9.03(47) stat (42) sys .
In the second line, the three systematic errors are combined in quadrature.The 13-pt CCFV fit to the S sim data on each ensemble fit with z 2 , is shown in Fig. 15.

gπNN
The CCFV extrapolation to obtain g πN N is carried out using Eq. ( 29) for (i) the product g πN N F π = M N F P (−M 2 π ), and the result, in the continuum, divided by F π = 92.9MeV; and (ii) F P (−M 2 π ) and the result multiplied by M N (= 939 MeV)/F π (= 92.9 MeV).It turns out that these two extrapolations have different systematics: the slopes with respect to M 2 π of g πN N F π and F P (−M 2 π ) are different as shown in Fig. 15.The two estimates are, however, consistent and we take their average to get the g πN N for the nine cases summarized in Table V: the three strategies and the three types of Q 2 fits.
is taken from the S sim data with z 2 fits and the errors are estimated as for g * P .
C. FP (Q 2 ) at the Physical Point The F P (Q 2 ) at the physical point was obtained following the same three step procedure used for extrapolating G A (Q 2 ) that is described in Sec.III B. This F P (Q 2 ) is compared with the G A (Q 2 ) already shown in Fig. 12  Similar to G A (Q 2 ), the two physical mass ensembles impact the coefficients b i (Q 2 ) shown in Fig. 17 only for Q 2 ≲ 0.4 GeV 2 .The plots show some pion mass dependence for Q 2 < 0.2 GeV 2 , i.e., b 2 (Q 2 ) ̸ = 0.The coefficients for the lattice spacing dependence, b 1 (Q 2 ), and for finite volume, b 3 (Q 2 ), have large uncertainty.Also, neglecting the finite volume term does not change b 1 (Q 2 ) and b 2 (Q 2 ) significantly.Overall, the shape of these coefficients versus Q 2 is somewhat different from those for G A shown in Fig. 11.
The z 2 fit to the physical point F P , shown in Fig. 16, with t c = 9M 2 π and t 0 = 0.25 GeV 2 has the parameterization The agreement, within errors, with the parameterization of G A (Q 2 ) given in Eqs. ( 32) and ( 33) is very good.This is not unexpected based on the overlap between the two shown in Fig. 16, nevertheless, one should keep in mind the Q 2 dependence shown in Fig. 5, and the Goldberger-Trieman discrepancey [7,17].pole from strategy Ssim for a ≈ 0.15, 0.12 0.09, 0.06 fm lattices.The z 2 (magenta lines with error band) is compared to the z 1 fit (blue dash dot).The four largest z points are excluded from the fits to the Mπ ≈ 310 MeV and a12m220S data.The vertical black dotted line corresponds to Q 2 = 0.The value of t0 and ensemble name are given in the labels.
FIG. 14. (Left) Results for the bare g * P from the strategy Ssim and the differences ∆(X) = g * P (X) − g * P (Ssim).To facilitate visualization of the spread, the errors plotted for ∆(X) are those in g * P (X).The fits used to parameterize the Q 2 behavior are labeled "PD" defined in Eq. (42); "D" for the dipole fit, and z k for various truncations of the z-expansion.(Right) Results for bare values of gπNN Fπ obtained with the strategy Ssim and the differences ∆(X) = gπNN Fπ(X) − gπNN Fπ(Ssim).) × g * P (top row), g * P (middle row), and gπNN Fπ (bottom row) data using the 13-point fit.In each case the data were obtained using the z 2 parameterization of FP (Q 2 ̸ = 0) with strategy Ssim.The black solid line with the pink error band represents a hyperplane obtained by taking the physical limit of all CCFV fit variables except the one shown on the x-axis.The plotted data points have been shifted by using the fit coefficients, while the errors are unchanged.FIG. 17.The behavior of the four CCFV fit coefficients, b0, b1, b2 and b3 defined in Eq. ( 29) versus Q 2 obtained in the extrapolation of FP (Q 2 ) to the physical point.The rest is the same as in Fig. 16.
The spectral decomposition of the forward, q = 0, three-point function truncated at three states, |i⟩ with i = 0, 1, 2, can be written as where The prefactors in terms involving the excited states are products such as r 2 2 b 22 .These products are simply parameters in the fits and are not used subsequently.Thus the ratios r i for the excited states are, by themselves, not needed.
To remove excited states contributions, we made three kinds of fits to the 3-point functions using Eq. ( 48): • 3 * : This is a 3-state fit with b 22 set to zero.The four parameters A 0 , M 0 , M 1 , M 2 are taken from four-state fits to the two-point function, leaving only ⟨0|O Γ |0⟩, and products such as r 2 1 b 11 as free parameters.This strategy (along with its two-state version) was used to get the results presented in Ref. [13] that are reproduced in Eq. (49). to zero, otherwise the fits become unstable.The three parameters A 0 , M 0 , M 1 are again taken from four-state fits to the two-point function.The value of the second mass gap, ∆M 2 , is left as a free parameter in the fit.The sign of ∆M 2 for a given charge determines whether |1⟩ lies above or below |2⟩ as shown pictorially in Fig. 18.
• 3-RD-N π: In this fit, M 1 is fixed to the noninteracting energy of the (N (n)π(−n)) state with n = (1, 0, 0)).For the value of M 2 , we use a Bayesian prior with a narrow width centered about the first excited state mass determined from the two-point correlator as given in Table XXX in Appendix F.
We also tried two-state fits with ∆M 1 left as a free parameter.For the axial charge, we found large fluctuations in ∆M 1 between the jackknife samples leading to unreliable values.So we do not present these estimates.Results for unrenormalized isovector nucleon charges, g A , g T , and g S , using the 3 * , 3-RD, and the 3-RD-N π fits are given in Table XXXI, and the other parameters of the 3-RD fits are given in Table XXXII in Appendix F.
The final renormalized charges are presented in the MS scheme at 2 GeV.We carry out the renormalization using the RI-sMOM intermediate scheme as described in Ref. [13].To understand systematics, we use three methods: , where X = A, T, S; and (ii) Empirically, some of the systematics cancel in each ratio in the second method, giving slightly smaller overall errors.In method three, we take the average of the first two within the jackknife process, and use it to get the final estimates.The renormalization factors Z X and Z X /Z V used in this study are given in Table V in Ref. [13].
We use the same leading order CCFV ansatz, given in Eq. ( 29), for extrapolating results to the physical point for all three strategies: 3 * , 3-RD, and 3-RD-N π.
We now focus on the 3-RD analysis and make three CCFV extrapolation with the following cuts on the data • "13-point" CCFV fit uses all thirteen points.
• "11-point" CCFV fit: This fit excludes the a06m310W and a06m220W points obtained with larger smearing radius for sources used to calculate the quark propagators [13].Larger smearing radius reduces the ESC at smaller values of τ but gives larger statistical errors at the values of τ used in our excited-state fits to get the τ → ∞ values as discussed in Ref. [14].Also, we expect significant correlations between the two pairs, (a06m310, a06m310W ) and (a06m220, a06m220W ), since they use the same gauge configurations.Comparing the two pairs, the results for the three charges agree except for g S between a06m220 and a06m220W , which can be accounted for by the larger statistical errors, especially in the a06m220W data.Consequently, the "11-point" CCFV fit is used to get the central value for g S , which is shown in Fig. 19.
• "10-point" CCFV fit excludes a15m310, a06m310W and a06m220W points.Since the variation with the lattice spacing is the dominant systematic, removing the a15m310 point (coarsest lattice with a ∼ 0.15 fm) aims to provide a handle on higher order, O(a 2 ), corrections neglected in Eq. ( 29).
Results from these three CCFV extrapolations, different truncations of the CCFV ansatz, and the three renormalization methods are given in Tables XXXIII, XXXIV, and XXXV in Appendix F and used to assess the various systematics.
The central values are taken from the "13-point fit" for g A and g T and the "11-pt fit" for g S with the 3-RD data renormalized using the third (average of the first two) method.Note that we find a systematic shift of ≈ 0.03, 0.02 and 0.03 between the first two renormalization methods for the three charges, g A , g T and g S , respectively.
These CCFV fits are shown in Fig. 19.Each panel in a given row shows the fit result versus one of the three variables with the other two set to their physical point values.In the left two panels, we show two fits: (i) using the full ansatz given in Eq. ( 29) (pink band), and (ii) assuming there is dependence only on the x-axis variable (gray band).For example, in the left panels the gray band corresponds to a fit with b g X 2 = b g X 3 = 0.The data show that the discretization errors are the dominant systematic, i.e., there is an almost complete overlap of the two fits (pink and gray bands) for g A and a significant overlap for g S and g T .The variation with a over the range 0 < a ≤ 0.15 fm is about 10%, 5% and 30% for g A , g T , and g S , respectively.The large variation with a in g S is similar to that found in the clover-on-clover calculation [7].
The final results of the 3-RD analysis are: The first error quoted (labeled stat) is the total uncertainty from the central analysis.The second error is an estimate of the uncertainty in the CCFV extrapolation.For g A and g T , this is taken to be the average of the differences |11-pt -13-pt| and |10-pt -13-pt|.For g S , it is the difference |10-pt -11-pt|.The third error is half the difference in estimates between the first two renormalization methods.The g A from the 3-RD fit is in good agreement with the result obtained from the extrapolation of the axial form factor G A (Q 2 ) to Q 2 = 0 given in Eq. (30).It is also consistent with the experimental value g A = 1.2766 (20) but has much larger errors.The difference between the 3 * (PNDME [13]) value reproduced in Eq. ( 49), and the 3-RD is due to different excited state energies used in the fits to the spectral decomposition in Eq. 48.The data in Table XXXII show that the fit parameter M 2 when left free satisfies M π ≲ M 2 − M 0 ≲ 3M π for all but the a ≈ 0.12 fm lattices.In [10], we showed evidence that the N (p 1 ) + π(−p 1 ) with p 1 = (1, 0, 0)2 π/La state makes a significant contribution on the zero momentum side of the operator insertion in the calculation of the form factors, and the noninteracting energy of this state is M N π − M 0 ≈ 2M π .In short, the M 2 output by the 3-RD fit has a mass lower than M 1 obtained from the two-point correlator and broadly consistent with the hypothesis that the N π states contribute.We again caution the reader that these excited state masses should only be regarded as effective fit parameters that encapsulate the effect of the full tower of N (p) + π(−p) states with momenta p = (2π/L)n as well as other multihadron and radial excitations that can contribute.
For g S and g T , the 3-RD fit reduces to a 2-state fit if ∆M 2 (= M 2 −M 1 ) = 0, i.e., M 2 ≃ M 1 .This is the case for many of the ensembles as shown in Table XXXII.Results given in Eq. ( 50) are consistent with those in Eq. ( 49) indicating that sensitivity to excited state energies input into the analysis is small.
Based on the 3-RD fits, which indicate that the data for g A prefer a low-mass excited state with ∆M ≈ 2M π , the 3-RD-N π fit defined above, with the mass gaps summarized in Table XXX, were performed.Charges from this fit are compared with the 3-RD and the 3 * -state fits (or the 2-state fit for the g S ) in Table XXXI for the thirteen ensembles.The p-value for many of the 3-RD-N π fits are low.To stabilize the 3-RD-N π fits, we increased the width of the priors for M 2 , however, this still did not lead to stable fits for several ensembles.The 3-RD-N π results are, therefore, not included in the final analysis.
Results for g A,S,T in Eq. ( 50) are compared with those from other collaborations in Table VI.This comparison complements the latest FLAG review that considered results published prior to 2021 [5,6].Overall, results for g A and g T from all calculations are consistent within five percent and for g S at ten percent.Their precision will continue to improve steadily as higher statistics data are generated at additional {a, M π } points with M π L ≳ 4.
Our conclusion is that, with current statistics, fits for the axial charge are more stable with input of M 0 and M 1 from the 4-state fit to the 2-point function and letting the 3-point function determine M 2 (corresponding roughly to the N π state), i.e., the 3-RD fit.In future works with higher statistics, we expect results from 3-RD and 3-RD-N π fits to come together as the same three states provide the dominant contributions.Only the details of their inclusion are different.

VI. CONCLUSIONS AND COMPARISON WITH PREVIOUS CALCULATIONS
We have presented results for the axial, G A (Q 2 ), the induced pseudoscalar, G P (Q 2 ), and the pseudoscalar, G P (Q 2 ), form factors of nucleons using thirteen ensembles of 2+1+1-flavors of HISQ ensembles generated by the MILC collaboration [8].A large part of the focus of this work is on understanding the nature of the excited states that contribute significantly to the relevant correlation functions and removing their contributions.The analysis presented here strengthens the case for including multihadron excited states, such as N π, made in Ref. [10].Our data driven analysis strategy, labeled S sim , identifies the contributions from N π state in the extraction of the GSME.The three form factors obtained including N π state satisfy the PCAC relation to within 10% as opposed to a ∼ 50% deviation without them.For the final results, we therefore choose the data obtained with the S sim strategy for removing ESC, parameterize the Q 2 behavior using the model-independent z 2 fit; and extrapolate the data to the physical point using a simultaneous chiral-continuum-finite-volume (CCFV) ansatz including the leading order corrections given in Eq. (29).For errors, we quote two estimates: the first labeled "stat" is the total error obtained from the analysis used to produce the central value, and the second, labeled "sys", includes the various systematic uncertainties combined in quadrature as discussed in the appropriate sections.
Our final results are: • The axial charge is g A = 1.292(53) stat (24) sys .This is the unweighted average of the value from the extrapolation of G A (Q 2 ) to Q 2 = 0 (Eqs.(36)) and from the forward matrix element (Eq.( 50)).The "stat" and "sys" errors quoted are the larger of those from the two determinations.This result is consistent with the experimental value but has much larger errors.
• The extraction of the axial charge radius squared is discussed in Sec.III B, and the result from Eq. ( 36) is ⟨r 2 A ⟩ = 0.439(56) stat (34) sys fm 2 .
• The extraction of the induced pseudoscalar charge is discussed in Sec.IV B 1 with the result g * P = 9.03(47) stat (42) sys given in Eq. (44).
• Our procedure for obtaining the axial form factor, G A (Q 2 ), in the continuum limit is discussed in Sec.III B. The final parameterization is given in Eq. ( 32), the covariance matrix of the fit in Eq. ( 33), and the corresponding values of g A = 1.281(53) and ⟨r 2 A ⟩ = 0.498(56) fm 2 in Eq. ( 34).The overall final values from the analysis of G A are given in Eq. (36).
A comparison of lattice results from various collaborations for all the above quantities was presented recently in Ref. [7].The charges g A,S,T have also been reviewed by FLAG [5,6].Since then, new results have been presented in Refs.[27][28][29]32].The full list of relevant publications that have included N π states in the analysis of ESC and checked whether form-factors satisfy the PCAC relation are [7,10,24,26,28,29].We first summarize the results and the important features in each of these calculations, and then show a comparison of G A (Q 2 ) obtained by the various collaborations in Fig. 20.Results for the charges are compared in Table VI.
The observation that the form factors extracted using the spectrum from the nucleon 2-point function fail to satisfy the PCAC relation Eq. ( 17) was made in Ref. [9].The possible cause, enhanced contributions of multihadron (N π) excited states in the axial channel was proposed by Bär [12] using a χP T analysis.This was confirmed using the data for the three-point function with the insertion of the A 4 current in Ref. [10].This datadriven analysis, including only the lowest N π excited state, found that the ESC to the G P (Q 2 ) and G P (Q 2 ) form factors were about 35%, while that in G A (Q 2 ) could be O(5%) as the latter is affected only at one-loop in χP T [11,12].The smallest Q 2 data in Tables XVII and XXII in Appendix C for the two physical mass ensembles a09m130W and a06m135 show ∼ 5%, ∼ 45% and ∼ 45% difference between the S 2pt and S sim values.A ∼ 5% level of effect in G A (Q 2 ) is also consistent with what is observed in the axial charge g A extracted from the forward matrix element as shown in Table XXXI in Appendix F.
The RQCD collaboration [24] has extracted G A (Q 2 ) from a two-state fit to thirty-six 2+1-flavor Wilson-clover ensembles generated by the coordinated lattice simulations (CLS) collaboration.The G P (Q 2 ) and G P (Q 2 ) are, on the other hand, extracted using a 3-state fit in which the first excited state energies are fixed to be the noninteracting energies of the lowest N π state and the second excited state energies are taken to be the first excited state (values higher than N (1440)) given by fits to the 2-point nucleon correlators.While their form factors satisfy the PCAC relation, they are equally well fit by the dipole ansatz and z 4+3 (i.e., z 3 and ZX /ZV × g where the gV is the vector charge.In each panel, the pink bend with black solid line represents the full CCFV fit.In the left (middle) panels, the gray band shows the fit to the date keeping only the a (M 2 π ) dependent term in Eq. (29).The value at the physical point is marked by the red star.The data in each panel have not been shifted to the physical point in the other two fit variables.
with sum rule constraints).The axial charge g A = 1.302(86) from z-expansion (fit labeled !z 4+3 ) is larger than g A = 1.229 (30) from dipole (labeled !2P ) with the latter agreeing with that from the forward matrix element.The corresponding difference in ⟨r 2  A ⟩ is 0.449(88) versus 0.272(33) fm 2 .Results for g * P (8.68(45) versus 8.30 (24)) and for g πN N (12.93(80)versus 14.78(1.81))are consistent.They have recently [25] updated their results for g A,S,T based on the analysis of 47 ensembles.We quote these new values in Table VI.
The preferred estimates from the ETM collaboration [26] are from a single 2+1+1-flavor physical mass 64 3 × 128 ensemble at a ≈ 0.8 fm.For the analysis of G A (Q 2 ), they take excited-state energies from the 2-point function and find ⟨r 2 A ⟩ = 0.343(42)(16) fm 2 .Their result for g A = 1.283 (22) is obtained from the forward matrix element extracted without including possible contamination from N π states.When results from the direct calculations of G P (Q 2 ) and G P (Q 2 ) are used, the three form factors show large deviations from the PCAC relation which they attribute partially to large discretization errors in their twisted mass formulation [33].Consequently, they quote final estimates of G P (Q 2 ) derived from G A (Q 2 ) using the pion-pole dominance hypothesis, i.e., the quoted G P is not independently determined.
The NME collaboration [7] analyzed seven ensembles generated with 2+1-flavors of Wilson-clover fermions.They make a simultaneous fit to all five correlation functions with insertion of the axial, A µ , and pseudoscalar, P , currents, i.e., same as the S sim strategy used in this work.The A 4 correlator provides the dominant contribution to fixing the excited-state energies which turn out to be close to the lowest N π states as also discussed in this paper and in Ref. [10].The resulting form factors satisfy the PCAC relation to within ten percent.Observing only a small dependence of G A (Q 2 ) on a and M π , they provide a continuum parameterization of G A (Q 2 ) neglecting these effects, and thus underestimate the uncertainty.This G A (Q 2 ) is reproduced in Fig. 20.The value of the axial charge without including N π state is g A = 1.242(46)(42) and including it gives 1.32(6) (5).Their other results are ⟨r 2  A ⟩ = 0.428(53)(30) fm 2 , g * P = 7.9(7)(9) and g πN N = 12.4.(1.2).
The Mainz Collaboration [28] analyze fourteen 2+1flavor Wilson-clover ensembles also generated by the coordinated lattice simulations (CLS) collaboration.They TABLE VI.Comparison of gA,S,T , ⟨r 2 A ⟩, g * P and gπNN from recent calculations labeled as: PNDME 23 is this work, RQCD [24] (here we list values obtained with the !z 4+3 fit, and take gA,S,T from their recent work [25]), ETMC [26,27], NME [7], Mainz [28], and PACS [29].All results for gA,S,T are in the MS scheme at scale 2 GeV.For completeness, we also give results for gS,T from the Mainz collaboration [30] and from the ETMC collaboration [27,31].These and earlier results are reviewed by the FLAG [5,6].) at the physical point obtained using a z-expansion fit by the RQCD [24] (light faun band), ETMC [26] (light tan band), NME [7] (tan band), Mainz [28] (brown band) collaborations and this work (turquoise band).The GA extracted using the ν-deuterium bubble chamber scattering experiments data [22] is shown by the gray band and labeled νD in the lower panel.
obtain a parameterization of G A (Q 2 ) in the continuum from a single combined fit-summation method for dealing with ESC and the z 2 fit for the Q 2 behavior.This result is shown in Fig. 20 and from it they get g A = 1.225(39)(25) and ⟨r 2 A ⟩ = 0.370(63)(16) fm 2 .The PACS collaboration [29] has analyzed one physi-cal pion mass ensemble with a large volume (128 4 , i.e., (10.9fm) 4 ) at a = 0.085 fm and get g A = 1.288 (14) (9).Remarkably, they find that sources with exponential smearing for the generation of quark propagators, in contrast to Gaussian smearing used by all other calculations, leads to essentially no excited-state effects.The limitation of this calculation is only 20 configurations, each separated by 10 molecular dynamics trajectories, were analyzed.Most likely, this total of 200 trajectories represents less than one unit of autocorrelation time.Consequently, the errors are likely underestimated.Since no parameterization of G A (Q 2 ) was presented, we do not include their results in the comparison.
Phenomenologically, the most important quantity needed for the analysis of neutrino oscillation experiments is G A (Q 2 ), and we show a comparison of results from various lattice collaborations in Fig. 20 along with the extraction from the ν-deuterium bubble chamber scattering experiments [22].In all cases, except ETM, the data are extrapolated to the physical point and then fit using a truncated z-expansion.The bands in Fig. 20 overlap indicating that the lattice results are consistent within one sigma and the envelope of the bands suggests a roughly 10% uncertainty throughout the range 0 < Q 2 < 1.0 GeV 2 .The other significant observation is that the lattice results fall slower than the phenomenological extraction (the νN band) for Q 2> ∼ 0.3 GeV 2 .When comparing these lattice data, it is important to note that the various collaborations handle various systematics differently.These systematic effects will become clearer and the analysis more robust as the precision of the data increases.Similarly, recent calculations including N π states in a variational basis of interpolating operators [34] is a step forward in providing further insight into the excited states contributions and developing better methods for removing them.
What has become clear is that N π states need to be included in the analysis for the three form factors, and satisfying the PCAC relation ( 17) is an essential and necessary condition.The questions that remain for higher precision are how many multihadron states need to be kept in the spectral decomposition for a given precision and the size of their contributions.The roughly 10% spread in the lattice results compared in Fig. 20 will be reduced with much higher statistics data that will be available over the next few years, and better analyses of excited states contributions.

VII. ACKNOWLEDGEMENT
We thank the MILC collaboration for providing the 2+1+1-flavor HISQ lattices used in this study.The calculations used the Chroma software suite [35] Here ns is the number of states kept in the fits with ns = 2 or 3 implying a frequentist fit and ns = 4 implying an empirical Bayesian fit.The speed of light c 2 , the χ 2 /dof and p-value are for the fits to the dispersion relation, which are shown for the a09m130W and a06m135 ensembles in Fig. 22.
and making an association with physical states.
eff reaches the first excited state energy E1 given by the four-state fit.The time interval used in the four-(three-) state fits to C 2pt (t) is given in the third (fourth) column.These fit ranges were chosen individually for each case balancing between keeping the maximum number of points with signal in the effective mass plot and the χ 2 /dof .
1 , from the frequentist 3-state (or 2-state) fit.The mass differences M   Here we revisit the extrapolation of the nucleon mass M N (a, M 2 π , M π L) given in Table VII to the physical point and extend the discussion in the Appendix B in Ref. [19].We use the following CCFV ansatz: Results and the fit parameters c i for various truncations of this ansatz are given in Table IX.The AIC score is defined as where k is the number of parameters and L is the likelihood function.When quoting AIC scores, we drop the irrelevant constant.The CCFV fits F1 and B1 are shown in Fig. 24.Our analysis indicates • The CCFV fits, F1-F4, to the M • The F3 and B3 fits, which include the higher order M 3 π term give a c 4 that is roughly consistent with the χPT prediction c 4 = 3g 2 A /(32πF 2 π ) = −5.716.On including a 2 and/or finite volume correction terms in addition to the M 3 π term, c 4 remains consistent with the χPT prediction for F1, F2 and F4 fits but becomes smaller for B1, B2 and B4.
• The finite volume coefficient, c 5 , is not well determined in any of the fits.Without it, fits F1 and B1 have small p-value but give results consistent with the experimental value.Including it, the p-value of F2 and B2 fits improves to an acceptable level, but the coefficients of the lattice spacing dependence, c 1 and c 2 , become less well determined.Neglecting the c 2 term (F4 and B4 fits), the c 1 becomes well determined, while the other c i are essentially unchanged.In these cases, the ∼ 25 MeV shifts in the M N from the F1 or B1 fits persist.
Overall, with the current data, we are not able to determine whether M (3) 0 or M (4) 0 give better estimates of the ground state nucleon mass.Also, we can at best make four-state fits to the two-point function and three-state fits to the three-point functions.after shifting them in a to a = 0 using the CC fits.The CC fit is also shown versus M 2 π with a set to a = 0.06 fm (blue dashed line), 0.09 fm (orange dotted line), 0.12 fm (green dotted line), and 0.15 fm (purple dash-dot line).In a perfect fit, these curves should pass through points with the same color, i.e., with the same lattice spacing a. in the text, and B1-B4 are to the 4-state empirical Bayesian fit data and labeled M (4) 0 .To make the interpretation of coefficients ci defined in Eq. (B1) easier, we give both the functional dependance within square parentheses and the units.The results for M (4) 0 are the same as in Ref. [19], except for a small change in a06m135 value due to increased statistics.Fits corresponding to the B2 and B4 were given in Table XV in Ref. [19] (labeled B1 and B2 there) and led to MN = 0.976 (20) GeV and 0.972(6) GeV, respectively.The table also gives the AIC score and the p value of the fits.41).The χ 2 /dof and p-value of the fits are also given, and Fπ and MN are in units of GeV.

Appendix E: Summary of CCFV Fits
This appendix presents results after chiral-continuum-finite-volume extrapolation for g A and ⟨r 2 A ⟩ in Table XXVI; g * P in Tables XXVII and XXVIII; and g πN N in Table XXIX.All the data for G A (Q 2 ) and G P (Q 2 ) used were obtained with the S sim strategy to remove ESC, and the Q 2 behavior was fit using either the z 2 truncation (G A (Q 2 )) or the PD fit ( G P (Q 2 )) given in Eq. ( 42).The four parameters, b 0,1,2,3 , define the CCFV ansatz given in Eq. ( 29).The tables also give the χ 2 /DOF, the p-value, and the Akaika Information Criteria (AIC and AICc) [20] scores for the CCFV fit.The definition of AIC is given in Appendix B, and including correction for small sample sizes, AICc is defined as rmAICc = AIC + (2k 2 + 2k)/(n − k − 1) where n is the number of data points and k is the number of parameters.

AQ 2 [GeV 2 ]FIG. 12 .
FIG.12.The final estimate of GA(Q 2 ) at the physical point.The eleven fiducial points used to make the fit are shown in blue and their errors are obtained from the one overall bootstrap analysis covering the three step process described in the text.The parameterization versus z is given in Eq.(32).

in Fig. 16 .
If PPD were exact, then F P should equal G A .The overlap of the two bands turns out to be surprisingly good over the whole Q 2 interval.Results for the four fit parameters, b i (Q 2 ), versus Q 2 obtained in the CCFV extrapolation process are shown up to 1 GeV 2 in Fig 17 for data obtained with the S sim strategy.

× f − 1 poleFIG. 13 .
FIG.13.Comparison of GA and FP ≡ GP × f −1 pole from strategy Ssim for a ≈ 0.15, 0.12 0.09, 0.06 fm lattices.The z 2 (magenta lines with error band) is compared to the z 1 fit (blue dash dot).The four largest z points are excluded from the fits to the Mπ ≈ 310 MeV and a12m220S data.The vertical black dotted line corresponds to Q 2 = 0.The value of t0 and ensemble name are given in the labels.

FIG. 15 .
FIG.15.The chiral-continuum-finite-volume extrapolation of the (Q * 2 + M 2 π ) × g * P (top row), g * P (middle row), and gπNN Fπ (bottom row) data using the 13-point fit.In each case the data were obtained using the z 2 parameterization of FP (Q 2 ̸ = 0) with strategy Ssim.The black solid line with the pink error band represents a hyperplane obtained by taking the physical limit of all CCFV fit variables except the one shown on the x-axis.The plotted data points have been shifted by using the fit coefficients, while the errors are unchanged.

FIG. 16 .
FIG.16.The close overlap in the physical point results for FP (Q 2 ), defined in Eq. (41), with GA(Q 2 ) (black lines) reproduced from Fig.12.Both were obtained using the three step process described in Sec.III B and the Ssim strategy data.The full CCFV fits to the FP data are shown with the solid red line and pink error band and the fit without the FV term with dashed blue line and hatched error band.These error bands show only the central analysis uncertainty. b

• 3 -
RD: This is a 3-state fit with b 01 , b 11 and b 22 set

FIG. 19 .
FIG.19.The simultaneous chiral-continuum-finite-volume (CCFV) fit to the axial gA (top, 13-point), tensor gT (middle, 13-point), and scalar gS (bottom, 11-point) charges.The data are extracted using the 3-RD fit described in the text and are the average over the two renormalization methods ZX × g (bare) X

0 FIG. 20 .
FIG. 20.A comparison of the isovector axial form factorGA(Q 2 ) at the physical point obtained using a z-expansion fit by the RQCD[24] (light faun band), ETMC[26] (light tan band), NME[7] (tan band), Mainz[28] (brown band) collaborations and this work (turquoise band).The GA extracted using the ν-deuterium bubble chamber scattering experiments data[22] is shown by the gray band and labeled νD in the lower panel.

0 ,
the second and third panels.The difference in the ground state mass, M is given in the bottom panel.

FIG. 23 .
FIG. 23.Data for the effective masses m (n) eff , defined in Eq. (A4), from the a06m135 two-point correlators.Top panel shows results for p = 0 with the En−1 in Eq. (A4) taken from the four-state fit (left) and three-state fit (right).These input energy levels En are shown by the dashed lines with yellow error bands.The red plus, green square, and blue circle symbols correspond to m (n) eff with n = 0, 1, 2, respectively.The bottom panel shows m (n)eff for p = 2πn/L, n = (1, 0, 0).In the bottom right panel, E0 is taken from the four-state fit and E1 = E0 + 2Mπ (solid black line) is assumed.
smaller continuum M N than fits to M shown in Fig.21(bottom panel) for each of the thirteen ensembles.•Only F1 (M N = 0.939(12) GeV) and B1 (M N = 0.945(16) GeV) fits give estimates consistent with the physical value of M N = 939 GeV.The other fits give ≈ 25 MeV higher values.

FIG. 24 . 4 ) 0 () 0 (
FIG. 24.The result of the chiral-continuum (CC) fit (no finite volume term) to the nucleon mass M(4) 0 (B1 fit in Table IX)(top panel) and M (3) 0 (F1 fit in Table IX) (bottom panel) is shown by the red line with the error band.The data for B1 and F1 fits given in Table VIII are plotted versus M 2 π there.

TABLE IV .
g * P from the z 2 -expansion, dipole, and pion-pole dominance (PD) fits.The first column gives the strategy used for extracting the matrix elements.In each value, the first error is the total analysis error and the rest are systematic errors explained in the text.

TABLE V .
Results for gπNN from the z 2 , dipole, and pion-pole dominance (PD) fits.The first column gives the ESC strategy used to extract the matrix elements.The first error is statistical and the rest are systematic as explained in the text.
⟨0|O Γ |0⟩ is the bare charge g Γ , the transition matrix elements are b ij ≡ ⟨i|O Γ |j⟩/⟨0|O Γ |0⟩, the ratios of amplitudes are r i = |A i |/|A 0 |, and the successive mass gaps are ∆M . This research used resources at (i) the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract 222.The plot of E 2 obtained from the Bayesian four-state fit versus p 2 for the a09m130W (squares) and a06m135 (circles) ensembles.The slope gives the speed of light, c 2 listed in TableVII.It shows significant deviation from unity for the a06m135 ensemble.Note that the blue line lying above most square points is a consequence of including the full covariance matrix. FIG

TABLE XI .
The bare form factors GA, GP , and GP versus Q 2 for the 3 strategies Ssim, SA4, and S2pt on ensemble a12m310.

TABLE XII .
The bare form factors GA, GP , and GP versus Q 2 for the 3 strategies Ssim, SA4, and S2pt on ensemble a12m220S.

TABLE XIII .
The bare form factors GA, GP , and GP versus Q 2 for the 3 strategies Ssim, SA4, and S2pt on ensemble a12m220.

TABLE XVII .
The bare form factors GA, GP , and GP versus Q 2 for the 3 strategies Ssim, SA4, and S2pt on ensemble a09m130W .

TABLE XVIII .
The bare form factors GA, GP , and GP versus Q 2 for the 3 strategies Ssim, SA4, and S2pt on ensemble a06m310W .Data for GP were, by accident, not saved.

TABLE XIX .
The bare form factors GA, GP , and GP versus Q 2 for the 3 strategies Ssim, SA4, and S2pt on ensemble a06m310.Data for GP were, by accident, not saved.

TABLE XX .
The bare form factors GA, GP , and GP versus Q 2 for the 3 strategies Ssim, SA4, and S2pt on ensemble a06m220W .Data for GP were, by accident, not saved.

TABLE XXIII .
Results for gA and ⟨r 2 A ⟩ given by z 2 fits to the axial form factor, GA(Q 2 ), obtained with the Ssim strategy.The χ 2 /dof and p-value of the fits are also given.

TABLE XXIV .
Results for g * P , gπNN Fπ and gπNN Fπ/MN given by the "PD" fits (defined in Eq. (42)) to GP obtained using Ssim strategy.The χ 2 /dof and p-value of the fits are also given, and Fπ and MN are in units of GeV.

TABLE XXV .
The values of g * P , gπNN Fπ and gπNN Fπ/MN given by z 2 fits to Ssim strategy data for FP defined in Eq. (

TABLE XXVI .
Summary of the parameters in the 13-point CCFV fit (Eq.(29)) to gA and ⟨r 2 A ⟩.The data used are given in TableXXIII.These were obtained by fitting the Q 2 behavior of GA, obtained with the Ssim strategy, using the z 2 truncation.Details are given in Sec.III A.
A ⟩ (obtained from GA with Ssim and z 2 fit) extrapolated using the 13-point CCFV

TABLE XXVII .
Summary of parameters values in the 13-point CCFV fit (see Eq. (29)) for obtaining g * P .The data used are given in Tables XXIV and XXV.In the top half, the quantity(Q * 2 + M 2 π )g * P = 2mµMN FP (Q * 2), with FP is defined in Eq. (41) and fit using z 2 , is extrapolated, while in the bottom half(Q * 2 + M 2 π )g * P = (Q * 2 + M 2 π )(mµ/2MN ) GP (Q * 2) is used.The extrapolated results are then converted to g * P by dividing by the physical value of (Q * 2 + M 2 π ).Details are given in Sec.IV.(obtained from GP with Ssim and "PD" fit, Eq. (42)) extrapolated using the 13-point CCFV plus pole fit * P (obtained from FP (Q * 2 ) using Ssim data and fit using z 2 ) extrapolated using the 13-point CCFV plus pole

TABLE XXIX .
Summary of the 13-point CCFV fit parameters for the extraction of gπNN as described in Sec.IV B 2. The data used are given in Tables XXIV and XXV.In the top table, the product gπNN Fπ = MN FP (−M 2 π ) is extrapolated, and the result, in the continuum, is divided by Fπ = 92.9MeV.In the bottom table, FP (−M 2 π ) is extrapolated and the result in the continuum multiplied by MN /Fπ.

TABLE XXX .
Mass gaps of excited state for the 3-RD-N π fit in units of the lattice pion mass for each ensemble.The M1 is fixed to the noninteracting energy of the N (n) + π(−n) state with n = (1, 0, 0).The M2 is constrained to be near the first excited state mass given by the two-point correlator by using the narrow prior shown in the last column.These mass gaps can be compared with the 3-RD fit results given in the TableXXXII.

TABLE XXXI .
Summary of bare charges gA, gS, gT , and gV obtained from forward matrix elements along with the p-value of the three fits used to remove ESC: 3-RD (first row), 3 * (or 2-state for gS) (second row), and 3-RD-N π (third row) described in the text.The mass gap M2 − M0 output by the 3-RD-N π fits is summarized in TableXXX.

TABLE XXXIV .
Summary of CCFV fits to gT using Eq.(29).The rest is same as in TableXXXIII.