Axial Vector Form Factors from Lattice QCD that Satisfy the PCAC Relation

Previous lattice QCD calculations of axial vector and pseudoscalar form factors show significant deviation from the partially conserved axial current (PCAC) relation between them. Since the original correlation functions satisfy PCAC, the observed deviations from the operator identity cast doubt on whether all the systematics in the extraction of form factors from the correlation functions are under control. We identify the problematic systematic as a missed excited state, whose energy as a function of the momentum transfer squared, $Q^2$, is determined from the analysis of the 3-point functions themselves. Its mass is much smaller than those of the excited states previously considered and including it impacts the extraction of all the ground state matrix elements. The form factors extracted using these mass/energy gaps satisfy PCAC and other consistency conditions, and validate the pion-pole dominance hypothesis. We also show that the extraction of the axial charge $g_A$ is very sensitive to the value of the mass gaps of the excited states used and current lattice data do not provide an unambiguous determination of these, unlike the $Q^2 \neq 0$ case. To highlight the differences and improvement between the conventional versus the new analysis strategy, we present a comparison of results obtained on a physical pion mass ensemble at $a\approx 0.0871\,\mathrm{fm}$. With the new strategy, we find $g_A = 1.30(6)$. A very significant improvement over previous lattice results is found for the axial charge radius $r_A = 0.74(6)\,\mathrm{fm}$, extracted using the $z$-expansion to parameterize the $Q^2$ behavior of $G_A(Q^2)$, and $g_P^\ast = 8.06(44)$ obtained using the pion pole-dominance ansatz to fit the $Q^2$ behavior of the induced pseudoscalar form factor $\widetilde{G}_P(Q^2)$.

The nucleon axial form factor G A (Q 2 ) is an important input needed to calculate the cross-section of neutrinos off nuclear targets. It is not well-determined experimentally [1], and the most direct measurements using liquid hydrogen targets are unlikely to be performed due to safety concerns. At present, these form factors are typically extracted from measurements of scattering off nuclear targets and involves modeling of nuclear effects [2,3], which introduces uncertainties [4]. Lattice QCD provides the best approach to calculate these from first principles, however, one has to demonstrate that all systematics are under control.
The axial, G A , and the induced pseudoscalar, G P , form factors are extracted from the matrix elements of the four components of the isovector axial current A µ ≡ uγ 5 γ µ d between the ground state of the nucleon: (1) and the pseudoscalar form factor G P from where P = uγ 5 d is the pseudoscalar density, N (p) is the nucleon state with mass M and lattice momentum p ≡ 2πn/La with n ≡ (n 1 , n 2 , n 3 ). We neglect the induced tensor form factorG T in Eq. (1) since we assume isospin symmetry, m u = m d , throughout [5]. All the form factors will be presented as functions of the spacelike four-momentum transfer Q 2 ≡ p 2 −(E −M ) 2 = −q 2 .
In our previous work [6], we showed that form factors with good statistical precision can be obtained from lattice simulations, however, these data do not satisfy the partially conserved axial current (PCAC) relation: where m is the PCAC quark mass. Such a failure has also been observed in all other lattice calculations [7][8][9][10][11][12].
Since PCAC is an operator relation, it is important to find the systematic responsible for the deviation, and remove it prior to comparing lattice data with phenomenology.
In this work we show that the problematic systematic is a missed lower energy excited state. Using data from a physical pion mass ensemble, a09m130W [13], we show how the mass and energy gap of this state can be determined from the analysis of nucleon 3-point correlation functions. We then demonstrate that form factors extracted including these parameters satisfy PCAC and other consistency conditions. With these improvements, we claim that the combined uncertainty in the lattice data is reduced to below ten percent level.
All lattice data presented here are from our calculations using the clover-on-HISQ formulation [6,13]. The gauge configurations are from the physical-mass 2 + 1 + 1-flavor HISQ ensemble a09m130W generated by the MILC collaboration [14] with lattice spacing a ≈ 0.0871 fm and M π = 130 MeV. The pion mass on these configurations with the clover valence quark action is M π ≈ 138 MeV. Further details of the lattice param-arXiv:1905.06470v2 [hep-lat] 19 Sep 2019 eters and methodology, statistics, the interpolating operator used to construct the 2-and 3-point correlation functions can be found in Refs. [6,13].
The nucleon operator used to create and annihilate the nucleon state couples to the ground and all the excited and multiparticle states with appropriate quantum numbers. To isolate the ground state matrix elements, we fit the data for the 2-and 3-point functions, C 2pt and C 3pt Γ , using their spectral decompositions. For the 2-point functions, the four states truncation is where A i are the amplitudes and E i are the energies with momentum p. The data and fits using Eq. (4) are shown in Fig. 1 (left). There is a reasonable plateau at large τ in M eff (τ ) ≡ log C 2pt (τ ) C 2pt (τ +1) for all momenta up to n 2 = 6. The right panel shows M eff , M 0 , and the first two mass gaps, determined using a variant of the Prony's method [15], that are consistent with those obtained from 4-state fits [13].
The two-state truncation of the 3-point functions C where the source point is translated to t = 0, the operator is inserted at time t, and the nucleon is annihilated at the sink time slice τ . In this relation, |0 and |n are the ground and n th excited state. The superscript denotes that the state could have nonzero momentum p . The momentum transfer q = p − p = p since p at the sink is fixed to zero. The M i , E i and A i |A i are the masses, energies and the amplitudes for the creation|annihilation of these states by the nucleon interpolating operator. To display and discuss the data, it is much more convenient and common to consider the five ratios, R 5µ and R 5 , of the 3-point correlation functions of A µ and P to the 2-point correlator as defined in Ref. [6]: where q i is the momentum transferred by the current in the "i" spatial direction. The direction "3" is singled out since it is chosen to be the direction of the spin projection of the Dirac spinors in the construction of the 3-point functions. In the limit of large source-sink separation, τ → ∞, these ratios give the combination of the form factors shown on the right hand side. We have explicitly displayed the kinematic factors to show which momentum combinations have a signal in each case. Data with equivalent momenta are averaged in the analysis. Equations (6)-(9) form an over-complete set. R 51 and R 52 can be averaged as they are related by the lattice cubic symmetry and give G P . For q 3 = 0, R 53 gives G A . For the other momentum combinations, one gets a linear combination of G A and G P . Thus, the three A i correlators give results for G A andG P for all values of momentum transfer. Consequently, data from A 4 correlators have traditionally [6][7][8][9][10][11][12] been neglected because they exhibits very large excited-state contamination (ESC) as shown in Fig. 2. Lastly, G P is obtained uniquely from Eq. (10).
In our previous work [6], the energies of the excited state used to isolate the ground state matrix elements in fits to the 3-point functions were taken from fourstate fits to the 2-point correlation function defined in Eq. (4). The resulting form factors violated PCAC. Furthermore, the violation increased as Q 2 → 0, a → 0 and M π → M physical π . Correcting for O(a) lattice artifacts in the axial current showed no significant reduction in the violation [6].
In this paper we show that by using these values of M i and E i to remove the ESC we missed a lower excited state. Furthermore, the energy of this state can be determined from the analysis of the A 4 correlator, ie, the channel with the largest ESC is the most sensitive to it.
In Fig. 2, we compare the conventional 3 * -state fit to the A 4 correlator with masses and energies, M i and E i , taken from the 4-state fit to the 2-point function [6] versus the new two-state fit with M 1 and E 1 left as free parameters. The χ 2 /DOF and p-value of the fits for all ten momentum cases are given in Table I. Note that for n = (0, 0, 1), χ 2 /DOF reduces from 21.8 to 0.8. The values of the mass|energy gaps of the "first" excited state shown in Fig. 3 are much smaller, and close to those expected for non-interacting N (p)π(−p) and N (0)π(p) states [16]. By n 2 6, the mass gaps become similar and, correspondingly, the violation of PCAC at larger momentum-transfer are observed to be small (see Fig. 6 and Ref. [6]). We hypothesize that this lower energy excited state provides the dominant contamination in all five 3-point correlation functions. On the basis of consistency checks including PCAC, we make the case that it is essential to include this lower energy excited  To highlight the differences and improvements, we define two analysis strategies, "conventional", S 2pt , and "new", S A4 : • S 2pt : All the ground and excited state M i and E i , are taken from 4-state fits to the nucleon 2-point function and used in the 3 * -state analysis of all the 3-point functions as detailed in Ref. [13].
• S A4 : The ground state parameters M 0 and E 0 are taken from the 4-state fits to the nucleon 2-point function. These are considered reliable based on the observed plateau in the effective-mass data at large τ as shown in Fig. 1. The parameters for the first excited state, M 1 and E 1 , are taken from fits to the A 4 3-point correlator as discussed above. These are then used in a two-state analysis of all other 3-point functions.
In both cases, it is important to note that residual ESC may still be present in the M i and E i . Future higher precision calculations will improve the precision of the calculations by steadily including more states in the fits. In Fig. 3, we show three sets of data for the energy gaps of the first excited state: ∆E 2pt 1 ≡ E 2pt 1 − E 0 obtained from fits to the 2-point correlator. These are compared with the two values on either side of the A 4 operator insertion, which are expected to be different since the correlator is projected to zero-momentum at the sink: and ∆M 2pt 1 , whereas S A4 corresponds to using ∆E A4 1 and ∆M A4 1 . In Fig. 3, we also show, using dotted lines, the expected values for ∆E and ∆M if we assume that the leading contribution of the current A 4 (q) is to insert or remove a pion with momentum q. Thus the plotted ∆E corresponds to the values for a non-interacting state N (p = 0)π(p), while ∆M 1 to N (p)π(−p). In calculating these values, we have used the lattice values for M 0 and M π and the relativistic dispersion relation, which is consistent with the data from the 2-point function. The values and variation of ∆E A4 1 and ∆M A4 1 with n 2 are roughly consistent with this picture.
Using the excited-state parameters extracted from the analysis of the A 4 correlator, and following the strategy S A4 gives very different values for the ground state matrix elements extracted from the three spatial, A i , and the P correlators. A comparison of the 2-state fits using S A4 and the 3 * -state fit using S 2pt is shown in Fig. 4 for the lowest non-zero momentum channels. Based on  the χ 2 /DOF of the fits, we cannot distinguish between the two strategies except for the P channel in spite of having high statistics data (165K measurements on 1290 configurations) [13]. The key point in each of the four channels is the convergence-it is from below and including the "new" lower excited state (S A4 ) gives significantly larger values of the matrix elements and thus the form factors. This pattern persists for n 2 5, above which the difference in the mass gap does not have a significant effect.
The results for the three form factors G A , G P and G P are compared in Fig. 5. The effect of using S A4 is clear and largest for n = (1, 0, 0). Also, the change in G A (Q 2 ) is only apparent for n = (1, 0, 0), consequently data at smaller Q 2 are needed to quantify its Q 2 → 0 limit.
The pattern, that the effect increases as Q 2 → 0, a → 0 and M π → M physical π , is confirmed by the analysis of the 11 ensembles described in Ref. [13], and these detailed results will be presented in a separate longer paper [19].
With G A , G P and G P in hand, we present the test of the PCAC relation, Eq. (3), in Fig. 6. The figure also shows data for the pion-pole dominance (PPD) hypothesis that relates It is clear that both relations are satisfied to within 5% at all Q 2 with S A4 , whereas the deviation grows to about 40% with S 2pt at n = (1, 0, 0) as first pointed out in Ref. [6]. What is also remarkable is that the PPD relation with the expected proportionality factor 4M 2 provided by the Goldberger-Treiman relation [20] tracks the improvement in PCAC. In fact, the data for the two tests overlap at all Q 2 .
The last test we perform is the relation ∂ 4 A 4 = (M − E)A 4 that should be satisfied by the ground state matrix element. The data and fits for ∂ 4 A 4 are shown in Fig. 7. The values of (M − E)A 4 are essentially zero in both cases; for S 2pt because M − E is small. Again, it is clear that the relation is only satisfied for S A4 .
The bottom line is that the two relations, PCAC and ∂ 4 A 4 = (M − E)A 4 , and the pion-pole dominance hypothesis are all satisfied using S A4 but not with S 2pt . The data shown in Fig. 3 is consistent with the picture that the "new" lower energy state is mainly due to the current A µ (q) injecting a pion with momentum q. There are two open questions: (i) how do we extract g A , ie, what is the analogous lowest excited state at zero momentum since we cannot determine its parameters from the A 4 correlator, and (ii) why it was not clear from the data shown in Figs. 2 and 4 that the mass gaps used in S 2pt were too large. These points are addressed below.
Results for g A have been obtained from the A 3 correlator at zero-momentum in all previous calculations because it has the best signal. The states with the lowest energy that are candidates for the 1 2  ππ state at zero momentum, whereas in the other case it would need to insert a pion with (p = (1, 0, 0)) and at the same time cause the transition N (0) → N (p = (−1, 0, 0)) to ensure zero total momentum. In any case, since the only quantity that enters in our analysis is the mass gap and not the specifics of the excited state, we take the common value, ∆M = 1230 MeV, in the reanalysis of A 3 In what follows, all results for the renormalized axial current are presented using Z A = 0.95(4) taken from Ref. [18]. Fits to the zero-momentum A 3 correlator with prior ∆M a = 0.12 (4) give g A in the range 1.29-1.31 depending on the values of τ used in the fit compared to g A = 1.25(2) using S 2pt given in Ref. [18] (column for PCAC (open symbols) and   (7) 4.67 (24) 2 in Table II). However, fits with priors in the range 0.1 ∆M 1 0.4 are not distinguished on the basis of χ 2 /DOF. The output ∆M 1 tracks the input prior, and the value of g A increases as the prior value is decreased. Thus, we regard this method as giving g A with uncontrolled systematics-the relevant ∆M 1 has to be determined first. Parenthetically, similar fits to extract the scalar and tensor charges g S and g T are much more stable, the value of the output ∆M 1 is far less senstitive to the prior and the results vary by 2σ as will be shown in Ref. [19].
Our current best estimate for ∆M 1 on a given ensemble is to take the lower of the N ππ or N (n 1 = 1)π(n 2 = 1) states. Assuming they are roughly degenerate, one can use the value of ∆M A4 1 a ≈ 0.1 shown in Fig. 3 at n 2 = 1 that, as we have argued above, corresponds to the latter state. Using this ∆M A4 1 a, our analysis of the A 3 data gives g A = 1.30 (6).
The second way we extract g A is to parameterize the Q 2 dependence of G A (Q 2 ) using the z-expansion and the dipole ansatz. The z-expansion fits using the process defined in Ref. [6,13] give g A = 1.30(7) for S A4 compared to g A = 1.19(5) using S 2pt . These results are independent of k for k > 2 in the z k -expansion. The dipole fit gives g A = 1.20(6) with a large χ 2 /DOF = 1.97 and the results are essentially the same for S A4 and S 2pt as can be seen in Fig 5. One can fix the dipole fit to not miss the crucial low Q 2 points by putting a cut on Q 2 , however, for this study we choose to neglect it.
The root-mean-squared charge radius extracted using the z-expansion fits gives r A = 0.74(6) fm with S A4 and r A = 0.45(7) fm with S 2pt . Once the lattice data have been extrapolated to the continuum limit, they can be compared with (i) a weighted world average of (quasi)elastic neutrino and antineutrino scattering data [1], (ii) charged pion electroproduction experiments [1], and (iii) a reanalysis of the deuterium target data [21]: The induced pseudoscalar charge g * P , defined as is obtained by fittingG P (Q 2 ) using the small Q 2 expansion of the PPD ansatz: Our result using S A4 is g * P = 8.06(44), while the MuCap experiment gave g * P | MuCap = 8.06(55) [22,23]. We caution the reader that all the results summarized in Table II are at fixed a ≈ 0.0871 fm. Comparison to the phenomenological values should be made only after extrapolation to the continuum limit. The goal of this paper is to highlight the changes on using S A4 .
Second, we comment on why the lower-energy state is missed when following S 2pt . It is well known that extracting E i from n-state fits to C 2pt gives E i with ESC since the number of pre-plateau data point that are sensitive to excited states are typically 8-12 as shown in Fig. 1. While we find an ≈ 15% change between a 2and 4-state fit, we did not anticipate E A4 1 ∼ E 2pt 1 /4 at small Q 2 as shown in Fig. 8. The known methodology to getting a more realistic excited state spectrum in a finite box with nucleon quantum numbers is to construct a large basis of interpolating operators, including operators overlapping primarily with multiparticle states, and solve the generalized eigenvalue problem (GEVP) [24] in a variational approach [25][26][27][28][29]. One should then compare the energies with lattice data, for example in the axial case with E A4 1 , to determine which states contribute to a given 3-point function. This option will be explored in future calculations.
To get a rough picture of the impact of the choice of E 1 on the ESC in 3-point functions, assume that the prefactors in the two 0 ↔ 1 transition terms are equal and unity, and the 1 ↔ 1 term can be neglected in Eq. 5. Then the ESC should fall off exponentially as e −(M1−M0)(τ /2) + e −(E1−E0)(τ /2) . In Fig. 8, we plot this function for three typical values of (M 1 − M 0 ) and (E 1 − E 0 ) with n = (0, 0, 1). These values are from the 2-and 4-state fits to the 2-point function and those extracted from a 2-state fit to the A 4 3-point function. Over the interval 10 < τ /a < 16, corresponding to 0.9-1.4 fm in which lattice data are typically collected, the exponential fall-off is approximately linear. Furthermore, the three curves in this range can be roughly aligned by a constant shift in their magnitude, ie, by just a change in the prefactors we have set to unity in Eq. 5, and which are free parameters in the actual fits. Thus, over a limited range of τ , the expected exponential convergence can be masked to look linear. On the other hand, the size of the ESC is very sensitive to E 1 and large even at τ /a = 25 for E A4 1 . The lesson is that while the excited state energy gaps impact the magnitude of the ESC at any given τ , checks on the E i using 2-state fits and the convergence to τ → ∞ of ground state matrix elements is hard to judge from a limited range of τ even for very different energy gaps. As discussed above, the extraction of g A is plagued by this problem since we are not able to extract M A4 1 from a 3-point function. In short, it is very important to determine E 1 reliably. Once this is done, even 2-state fits give reasonable estimates of the τ → ∞ value based on the consistency checks discussed above.
An attempt to resolve the PCAC conundrum has been presented in [30]. We contend that it missed resolving the lower energy state and did not solve the problem. The projected currents A ⊥ µ and P ⊥ introduced in their work consist of a rotation in the basis of the five currents A µ and P . For the lattice ensembles and parameters explored in their [30] and our [6,13] calculations, the three A ⊥ i essentially remain within the space of the A i . Thus, the S 2pt strategy with A ⊥ i gives G A (Q 2 ) and G P that are essentially unchanged, and one continues to get a low value for g * P [30]. The operator A ⊥ 4 is mostly rotated into the A i . Thus A ⊥ 4 no longer shows the large ESC, and the "sinh" behavior illustrated in Fig. 2 becomes "cosh" like. Their "fix" to PCAC comes from P ⊥ , which now gets its dominant contribution from A 4 and ∂ 4 A 4 . Analysis of our a09m130W ensemble shows that the contribution of the ∂ 4 A 4 part is roughly three times that of P due to the small value of the PCAC mass m in the definition of P ⊥ . Also, note that, by construction, the total contribution of (P ⊥ − P ) is supposed to be zero in the ground state. On the other hand, we contend that the solution to the PCAC problem lies in the identification of the lower energy excited state[s] that, as we have presented, should be used to remove the ESC in all 3-point axial/pseudoscalar correlators. Using S A4 changes the results for all three form factors, especially at low Q 2 .

CONCLUSIONS
All previous lattice calculations of the three form factors G A , G P and G P [6][7][8][9][10][11][12], showed significant violations of the PCAC relation, Eq. (3). This failure had cast doubts on the lattice methodology for extracting these form factors. In this work, we show that the systematic responsible for the violation is a lower energy excited state missed in previous analyses. Furthermore, its energy can be extracted from fits to the A 4 3-point function. Detailed analysis of the A 4 correlator had, so far, been neglected as it is dominated by ESC and is not needed to extract the form factors. Using the mass/energy gaps of this lower excited-state, we show that lattice data satisfy PCAC to within 5%, the level expected with reasonable estimates of the current level of statistical and systematic errors. An additional consistency check is that the ground state matrix elements now satisfy the relation ∂A 4 = (M − E)A 4 . We also show that pion-pole dominance works to the same level as PCAC with the proportionality constant 4M 2 suggested by the Goldberger-Treiman relation.
We show that the direct extraction of g A from the A 3 correlator at zero-momentum requires knowing the energies of the excited states that give the dominant contamination, ie, the result for g A is particularly sensitive to the input value of the mass gap ∆M 1 . We show that the ∆M 1 obtained from the 2-point function is much larger than what is expected, so alternate methods for determining it are needed because fits to the A 3 correlator data, while precise, are not able to distinguish between ∆M 1 in a wide range. Our new analysis using two plausible estimates of ∆M 1 gives g A = 1.30 (6).
We provide heuristic reasons for why previous fits to remove ESC with a large ∆M 1 did not exhibit large χ 2 /DOF, and why the smaller values of the mass gaps that impact the extraction of the form factors were missed. For the form factors at Q 2 = 0, the good news is that implementing this improvement, in the axial channels (and an analogous procedure for the vector current for extracting electromagnetic form factors), does not require the generation of new lattice data but only a reanalysis.
We demonstrate the improvement in G A (Q 2 ), G P (Q 2 ) and G P (Q 2 ) by analysing a physical mass ensemble with a ≈ 0.0871 fm, M π ≈ 138 MeV [6,13]. We perform both the dipole and z-expansion fits to G A (Q 2 ) to parameterize the Q 2 behavior and extract the axial charge radius squared, r 2 A . The dipole ansatz does not fit the data well and is dropped. The z-expansion fit gives r A = 0.74(6) fm. We fit G P (Q 2 ) using the pion-pole dominance ansatz and find g * P = 8.06(44). To obtain results for these quantities in the continuum limit, a full analysis of the 11 ensembles described in Ref. [13] is in progress.