Isovector Axial Charge and Form Factors of Nucleons from Lattice QCD

Abstract

A survey of the calculations of the isovector axial vector form factor of the nucleon using lattice QCD is presented. Attention is paid to statistical and systematic uncertainties, in particular those due to excited state contributions. Based on a comparison of results from various collaborations, a case is made that lattice results are consistent within 10%. A similar level of uncertainty is in the axial charge $g_A^{u-d}$, the mean squared axial charge radius $⟨r_A^2⟩$, the induced pseudoscalar charge $g_P^{\ast }$, and the pion–nucleon coupling $g_{\pi NN}$. Even with the current methodology, a significant reduction in errors is expected over the next few years with higher statistics data on more ensembles closer to the physical point. Lattice QCD results for the form factor $G_A\left(Q^2\right)$ are compatible with those obtained from the recent MINER$\nu$A experiment but lie 2–3$\sigma$ higher than the phenomenological extraction from the old $\nu$–deuterium bubble chamber scattering data for $Q^2>0.3$ GeV². Current data show that the dipole ansatz does not have enough parameters to fit the form factor over the range $0\le Q^2\le 1$ GeV², whereas even a $z^2$ truncation of the z expansion or a low order Padé are sufficient. Looking ahead, lattice QCD calculations will provide increasingly precise results over the range $0\le Q^2\le 1$ GeV², and MINER$\nu$A-like experiments will extend the range to $Q^2\sim 2$ GeV² or higher. Nevertheless, improvements in lattice methods to (i) further control excited state contributions and (ii) extend the range of $Q^2$ are needed.

1. Introduction

The axial charge, $g_A^{u-d}$, gives the strength of the coupling of the weak current to the nucleons. It has been determined very accurately from the asymmetry parameter A (relative to the plane defined by the directions of the neutron spin and the emitted electron) in the decay distribution of the neutron, $n\to p+e^{-}+{\overline{\nu }}_e$. The best determination of the ratio of the axial to the vector charge, $g_A/g_V$, comes from using (i) polarized ultracold neutrons (UCN) using the the UCNA collaboration, $1.2772\left(20\right)$ [1,2], and (ii) cold neutron beam using the PERKEO III, $1.27641\left(45\right)\left(33\right)$ [3,4]. Note that, in the SM, $g_V=1$ up to second-order corrections in isospin breaking [5,6] as a result of the conservation of the vector current.

The axial charge enters in many analyses of nucleon structure and of the Standard Model (SM) and probes of beyond-the-SM (BSM) physics [7,8]. For example, it enters in the relation between the Cabibbo–Kobayashi–Maskawa (CKM) matrix element $V_{ud}$ and the neutron lifetime, ${\tau }_n$. High precision extraction of $V_{ud}$, knowing ${\tau }_n$ and $g_A$, is important for the test of the unitarity of the first row of the CKM matrix [9,10,11]. It is needed in the analysis of neutrinoless double-beta decay [12] and in the rate of proton–proton fusion [13], which is the first step in the thermonuclear reaction chains that power low-mass hydrogen-burning stars like the sun.

The axial vector form factor (AVFF) gives the dependence of this coupling on the momentum squared transferred by the weak current to the nucleon. It is an input in the theoretical calculation of the neutrino–nuclei scattering cross-section needed for the analysis of neutrino oscillation experiments [14,15,16]. The cleanest experimental measurement would be from scattering neutrinos off liquid hydrogen targets; however, these are not being carried out due to safety concerns. Extractions from ongoing neutrino scattering experiments (T2K, NOvA, MINERvA, MicroBooNE, SBN) have uncertainty due to the lack of precise knowledge of the incoming neutrino energy and the reconstruction of the final state of the struck nucleus and thus of the cross-section.

The MINER$\nu$A experiment [17] has recently shown that the axial vector form factor of the nucleon can be extracted from the charged current elastic scattering process ${\overline{\nu }}_{\mu }H\to {\mu }^{+}n$ in which the free proton in hydrogen (H) (part of the hydrocarbon in the target) is 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 that have unresolved systematics. On the theoretical front, lattice QCD provides the best method for first principal nonperturbative predictions with control over all sources of uncertainty [14,15].

A recent comparison [18] of results for the AVFF from lattice QCD [19], the MINER$\nu$A experiment [17], and the phenomenological extraction from neutrino–deuterium data from 1980s and 1990s [20] showed that, in the near term, the best prospects for determining the AVFF will be a combination of lattice QCD calculations and MINER$\nu$A-like experiments. Lattice QCD will provide the best estimates for $Q^2\lesssim 0.5$ GeV² and be competitive with MINER$\nu$A for $0.5\lesssim Q^2\lesssim 1$ GeV². For $Q^2\gtrsim 1$ GeV²; new ideas are needed for robust predictions using lattice QCD.

The goal of theory efforts in support of neutrino oscillation experiments is robust calculations of the cross-section for targets, such as ${}^{12}C$, ${}^{16}O$, and ${}^{40}\mathit{Ar}$, being used in experiments. This involves a four step process: a precise determination of the AVFF, nuclear models of the ground state of the nuclei from which the neutrino scatters, the intranucleus evolution of the struck nucleon using many-body theory to include complex nuclear effects up to $\approx 5$ GeV for the DUNE experiment, and the evolution of the final state particles to the detectors. The overall program requires complete implementation of these within Monte Carlo neutrino event generators [14,15,16]. The output of the generators provides the input essential to experimentalists for determining neutrino oscillation parameters from current and future experiments.

Here, I review the status of lattice QCD calculations of the axial charge, $g_A^{u-d}$, and the AVFF. In addition, note that the flavor diagonal axial charges $g_A^{u,d,s,c,b}$ provide the contribution of each quark flavor to the spin of the nucleon, whose calculation is computationally more expensive due to the additional disconnected contributions. The current status of the results for these nucleon charges has been reviewed by the Flavor Lattice Averaging Group (FLAG) in 2019 and 2021 [21,22]). Including results post FLAG 2021 [23,24,25,26,27], the values from the various calculations with 2+1- and 2+1+1-flavors of sea quarks lie in the ranges $1.22\lesssim g_A^{u-d}\lesssim 1.32$, $0.74\lesssim g_A^u\lesssim 0.89$, $-0.48\lesssim g_A^d\lesssim -0.38$, and $-0.06\lesssim g_A^s\lesssim -0.025$. There have been no substantial new results for flavor diagonal charges since the FLAG reports, so I will not discuss them further in this work.

Based on the results in Refs. [19,23,26,27,28,29], I present the case that lattice results for AVFF, over the range $0<Q^2\le 1$ GeV², are also available with $\lesssim 10\%$ uncertainty and agree with MINER$\nu$A results to within a combined sigma as discussed in Ref. [18] but disagree with the neutrino–deuterium results for $Q^2>0.3$ GeV². At the same time, I also highlight the need for much higher statistics and better control over excited state contributions to nucleon correlators in lattice calculations for the uncertainty to be reduced to the few percent level.

The outline of this review is as follows. I will summarize the methodology and steps in the calculation of the axial and pseudoscalar form factors in Section 2. This includes a discussion of the nucleon three-point correlation functions calculated in Section 2.1, how form factors are obtained from them in Section 2.2, and possible excited state contributions (ESC) that must be removed in Section 2.3. I then review the operator constraint imposed on the three form factors, the axial, $G_A\left(Q^2\right)$, the induced pseudoscalar, ${\tilde{G}}_P\left(Q^2\right)$, and the pesudoscalar $G_P\left(Q^2\right)$ by the axial Ward–Takahashi (also referred to in the literature as the partially conserved axial current (PCAC)) identity in Section 2.4, and how it provides a data-driven method for validating the enhanced contributions of multihadron, $N\pi$, excited states. These enhanced excited state contributions are due to the coupling of the axial and pseudoscalar currents to a pion, i.e., the pion pole dominance hypothesis. Extrapolation of the lattice results in the physical point defined by the continuum ($a=0$) and infinite volume ($M_{\pi }L\to \infty$) limits at physical light quark masses in the isospin symmetric limit, i.e., $m_u=m_d$ set using the neutral pion mass ($M_{{\pi }^0}=135$ MeV) is discussed in Section 2.5. A consistency check on the extraction of the axial charge is discussed in Section 2.6. I will then review the results for the AVFF obtained by the various lattice collaborations after extrapolation to the physical point in Section 3, the comparison of the lattice QCD result, the recent MINER$\nu$A data, and the phenomenological extraction from the old neutrino–deuterium scattering data, along with my perspective on future improvements in Section 4. I end with a few concluding remarks in Section 5.

2. Calculation of the Axial Vector form Factors Using Lattice QCD

The quark line diagrams for the two-point $C^{\mathrm{2pt}}$ and the three-point $C_J^{\mathrm{3pt}}(\mathit{q};t,\tau )$ (with the insertion of the axial, $A_{\mu }$, and pseudoscalar, P, currents) correlators are shown in Figure 1. The methodology for calculating these is the same in all ongoing calculations. For $C_J^{\mathrm{3pt}}(\mathit{q};t,\tau )$, two kinds of quark propagators are calculated—the original moving forward from the source point (say from the left blob), and a sequential propagator moving backwards from a nucleon source with definite momentum ${\mathbf{p}}_f$ at the sink (the right blob). This nucleon source is constructed by tying together two original propagators shown by the two bottom quark lines. The insertion of the current with three-momentum $\mathbf{q}$ between the source and sink nucleons then reduces to that between the original propagator and the sequential propagator, as shown by the top line. By momentum conservation, the source nucleon is projected to momentum ${\mathbf{p}}_i=\mathbf{q}+{\mathbf{p}}_f$. The Euclidean four-momentum transfer squared is given by $Q^2={\mathbf{q}}^2-{(E_N-M_N)}^2$.

In current calculations (standard method) the nucleon interpolating operator, $\mathcal{N}$, used is

$$\begin{array}{cc}\hfill \mathcal{N}\left(x\right)=& {ϵ}^{abc}\left[q_1^{aT}\left(x\right)C{\gamma }_5\frac{1\pm {\gamma }_4}{2}{q}_2^b\left(x\right)\right]q_1^c\left(x\right),\hfill \end{array}$$

(1)

where $C={\gamma }_4{\gamma }_2$, and the optional factor $1\pm {\gamma }_4$ projects on to positive parity nucleon states propagating forward/backward in time for zero momentum correlators. Developing a variational basis of interpolating operators that includes all significant $N\pi$ states, the holy grail of taming the ESC, is still a work in progress [30,31].

A short description of the six steps in the calculation of the AVFF that are common to all fermion discretization schemes and independent of the selection of input simulation parameters is given next in Section 2.1, Section 2.2, Section 2.3, Section 2.4, Section 2.5 and Section 2.6.

2.1. Correlation functions $C^{\mathrm{2pt}}$ and $C_J^{\mathrm{3pt}}(\mathit{q};t,\tau )$

Two kinds of smeared sources have been used to generate the original and sequential quark propagators in most lattice calculations: (i) the Wuppertal source [32] and (ii) the exponential source [25]. These quark propagators are stitched together to construct the gauge-invariant time-ordered correlation functions $C^{\mathrm{2pt}}(\mathit{p};\tau )$ and $C_J^{\mathrm{3pt}}(\mathit{q};t,\tau )$ shown in Figure 1, whose spectral decompositions are

$$\begin{array}{cc}\hfill C^{\mathrm{2pt}}(\mathit{p};\tau )\equiv & \langle \Omega \left|\mathcal{T}\left(\mathcal{N}\left(\tau \right)\overline{\mathcal{N}}\left(0\right)\right)\right|\Omega \rangle =\sum _{i=0}{\left|A_i^{\prime }\right|}^2{e}^{-E_i\tau },\hfill \end{array}$$

(2)

and

$$\begin{array}{cc}\hfill C_J^{\mathrm{3pt}}(\mathit{q};t,\tau )\equiv & \langle \Omega \left|\mathcal{T}\left(\mathcal{N}\left(\tau \right)J_{\Gamma }\left(t\right)\overline{\mathcal{N}}\left(0\right)\right)\right|\Omega \rangle ,\hfill \\ \hfill =& \sum _{i,j=0}{A_i^{\prime }}^{\ast }{A}_j\langle i^{\prime }\left|J_{\Gamma }\right|j\rangle e^{-E_{i}t-M_j(\tau -t)},\hfill \end{array}$$

(3)

where $J_{\Gamma }=A_{\mu }=\overline{\psi }{\gamma }_{\mu }{\gamma }_5\psi$, $J_{\Gamma }=P=\overline{\psi }{\gamma }_5\psi$ is the quark bilinear current inserted at time t with momentum $\mathit{q}$, and $|\Omega \rangle$ is the vacuum state. In the $C_J^{\mathrm{3pt}}(\mathit{q};t,\tau )$, the nucleon state $|j\rangle$ is, by construction, projected to zero momentum, i.e., $p_j=(M,\mathbf{0})$, whereas $\langle i^{\prime }|$ is projected onto definite momentum $p_i=(E,\mathit{p})$, with $\mathit{p}=-\mathit{q}$ by momentum conservation. The prime in $\langle i^{\prime }|$ indicates that this state can have nonzero momentum. Consequently, the states on the two sides of the inserted operator J are different for all $\mathbf{q}\ne 0$. The goal is to extract the ground-state matrix elements (GSME), $\langle 0^{\prime }\left|J\right|0\rangle$, from fits to Equation (3).

A major challenge in the analysis of all nucleon correlators is that the signal-to-noise ratio decays exponentially as $e^{-(M_N-1.5M_{\pi })\tau }$ with the source–sink separation $\tau$ [33,34]. In current calculations ($O\left({10}^5\right)$ measurements), a good signal in $C^{\mathrm{2pt}}(\mathit{p};\tau )$ and $C_J^{\mathrm{3pt}}(\mathit{q};t,\tau )$ extends to ≈2 and ≈1.5 fm, respectively.

At these $\tau$, the residual contribution of many theoretically allowed radial and multihadron excited states are observed to be significant. These states arise because the standard nucleon interpolating operator $\mathcal{N}$, defined in Equation (1), used to construct the correlation functions in Equations (2) and (3), couples to a nucleon and all its excitations with positive parity including multihadron states, the lowest of which are $N\left(\mathbf{p}\right)\pi (-\mathbf{p})$ and $N\left(\mathbf{0}\right)\pi \left(\mathbf{0}\right)\pi \left(\mathbf{0}\right)$. Since it is not known a priori which excited states contribute significantly to a given $C_J^{\mathrm{3pt}}(\mathit{q};t,\tau )$, the first goal is to develop methods to identify these and remove their contributions. Operationally, this boils down to determining the energies $E_i$ to put in fits to data using the theoretically rigorous spectral decomposition given in Equation (3).

2.2. Extracting the form Factors

Once the GSME, $\langle 0^{\prime }\left|J\right|0\rangle$, have been extracted, their Lorentz covariant decomposition into the axial $G_A$, induced pseudoscalar ${\tilde{G}}_P$, and pseudoscalar $G_P$ form factors is

$$\begin{array}{cc}\hfill & \langle N\left({\mathit{p}}_f\right)\left|A_{\mu }\left(\mathit{q}\right)\right|N\left({\mathit{p}}_i\right)\rangle \hfill \\ & =\overline{u}\left({\mathit{p}}_f\right)\left[G_A\left(Q^2\right){\gamma }_{\mu }{\gamma }_5+q_{\mu }{\gamma }_5\frac{{\tilde{G}}_P\left(Q^2\right)}{2M}\right]u\left({\mathit{p}}_i\right),\hfill \end{array}$$

(4)

and

$$\begin{array}{cc}\hfill \langle N\left({\mathit{p}}_f\right)\left|P\left(\mathit{q}\right)\right|N\left({\mathit{p}}_i\right)\rangle & =\overline{u}\left({\mathit{p}}_f\right)\left[G_P\left(Q^2\right){\gamma }_5\right]u\left({\mathit{p}}_i\right),\hfill \end{array}$$

(5)

where $u\left({\mathit{p}}_i\right)$ is the nucleon spinor with momentum ${\mathit{p}}_i$, $q=p_f-p_i$ is the four-momentum transferred by the current, and $Q^2=-q^2={\mathit{q}}^2-{(E\left({\mathit{p}}_f\right)-E\left({\mathit{p}}_i\right))}^2$ is the space-like four-momentum squared transferred. On the lattice, the discrete momenta are $\mathit{p}=2\pi \mathit{n}/La=2\pi (n_x,n_y,n_z)/La$ with $n_i\in \pm \{0\dots L\}$. The spinor normalization used is

$$\begin{array}{cc}\hfill \sum _{s}u(\mathit{p},s)\overline{u}(\mathit{p},s)& =\frac{E\left(\mathit{p}\right){\gamma }_4-i\gamma ·\mathit{p}+M}{2E\left(\mathit{p}\right)}.\hfill \end{array}$$

(6)

It is important to note that the excited states have to be removed from the C, which have the spectral decomposition given in Equation (3), and not from the form factors, i.e., after the decompositions. Equations (4) and (5) are only valid for the GSME. If there are residual ESC, then additional “transition” form factors have to be included in the rhs of Equations (4) and (5).

Assuming the GSME and choosing the nucleon spin projection to be in the “3” direction, the explicit forms of the decompositions in Equations (4) and (5) become

$$\begin{array}{cc}\hfill C_{A_{\{1,2\}}}\left(\mathit{q}\right)\to & K^{-1}\left[-q_{\{1,2\}}{q}_3\frac{{\tilde{G}}_P}{2M}\right],\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill C_{A_3}\left(\mathit{q}\right)\to & K^{-1}\left[-q_3^2\frac{{\tilde{G}}_P}{2M}+(M+E)G_A\right],\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill C_{A_4}\left(\mathit{q}\right)\to & K^{-1}{q}_3\left[(M-E)\frac{{\tilde{G}}_P}{2M}+G_A\right],\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill C_P\left(\mathit{q}\right)\to & K^{-1}{q}_3{G}_P,\hfill \end{array}$$

(10)

where the kinematic factor $K\equiv \sqrt{2E(E+M)}$. In each case, data with all equivalent momenta that have the same ${\mathit{q}}^2$ are averaged to improve the statistical signal. These correlation functions are complex valued, and the signal for the CP symmetric theory is in $\mathrm{Im}{C}_{A_i}$, $\mathrm{Re}{C}_{A_4}$, and $\mathrm{Re}{C}_P$.

It is clear that $G_P$ is determined uniquely from $C_P$ (Equation (10)) and for certain momenta ${\tilde{G}}_P$ from $C_{A_{\{1,2\}}}$ using Equation (7). The $C_{A_3}\left(\mathit{q}\right)$ and $C_{A_4}\left(\mathit{q}\right)$ give linear combinations of $G_A$ and ${\tilde{G}}_P$, and Equation (8) gives only $G_A$ when $q_3=0$.

2.3. Extracting the Ground State Matrix Elements: Exposing and Incorporating $N\pi$ States

The most direct way to extract $\langle 0^{\prime }\left|J\right|0\rangle$ is to make fits to Equation (3) keeping as many intermediate states as allowed by the data’s precision that demonstrate convergence. The problem is that even unconstrained two-state fits are numerically ill behaved. The next option is to take the $A_i$ and $E_i$ from $C^{\mathrm{2pt}}(\mathit{p};\tau )$, as $\mathcal{N}$ creates the same set of states in $C^{\mathrm{2pt}}$ and $C_J^{\mathrm{3pt}}$ and input these in fits to $C_J^{\mathrm{3pt}}(\mathit{q};t,\tau )$, either by doing simultaneous fits or via priors within say a bootstrap process to correctly propagate the errors. Of these, the ground states $A_0^{\prime }$, $A_0$, and $E_0$ are well determined from fits to the two-point function. Similarly, one would expect that $E_1$ can also be taken from $C^{\mathrm{2pt}}$. This was the strategy used until 2017 when it was shown in Ref. [35] that the resulting form factors do not satisfy the constraint imposed on them by the PCAC. Deviations from the PCAC due to discretization effects of about $\approx 5\%$ were expected; however, almost a factor of two was found on the physical pion mass ensembles.

The reason was provided by Bär [36,37] using $\chi$PT: enhanced contributions to ME from multihadron, $N\pi ,\dots$, which are excited states that have much smaller mass gaps than those of radial excitations, with the lowest being $\approx 1230$ versus N(1440). These states were not evident in fits to $C^{\mathrm{2pt}}(\mathit{p};\tau )$, as they have small amplitudes. A different approach to analysis to include the $N\pi$ states was needed.

It turned out that two-state fits to $C_{A_4}^{\mathrm{3pt}}$ provided a data-driven method [38] that exposed these states and confirmed that the lowest of the tower of $N\left(\mathbf{p}\right)\pi (-\mathbf{p})$ states makes a very significant contribution. By itself, $C_{A_4}^{\mathrm{3pt}}$ is dominated by excited states and fits to it using the $E_i$ from $C^{\mathrm{2pt}}$, which gave very poor ${\chi }^2/dof$. Making fits leaving $E_1$ a free parameter dramatically improved the ${\chi }^2/dof$ (compare left two panels in top row of Figure 2), and the resulting output values of $E_1$ on either side of operator insertion were roughly consistent with $N\left(\mathbf{p}\right)\pi \left(\mathbf{p}\right)$ as shown in Figure 3 (reproduced from Ref. [38]). An illustration of the current understanding of the transitions contributing to the GSME and of the lowest-contributing excited states is shown in Figure 4 (right).

The fact that there is an enhancement of the ME in the axial channel has been understood for over 60 years as the “pion pole dominance” (PPD) hypothesis. On the lattice, the creation of a $N\pi$ state by $\mathcal{N}$ is suppressed by the 3D volume compared to just the nucleon, as each state has a normalization factor of $1/V$ for a point (local) source $\mathcal{N}$. The axial current can, however, couple to this pion, and because the pion is light, this coupling can occur anywhere in the time slice at which the current is inserted with momentum $\mathit{q}$ (see Figure 4 (left)). This gives a factor of V enhancement, thus approximately canceling the normalization factor $1/V$ [36,37]. Thus, the ME obtains an enhanced contribution that is an artifact to be removed when the the pion comes on shell. Note that since energy is not conserved on the lattice, both the neutron and the pion can come on shell, however, since momentum is conserved, and possible excited states must have the same total momentum as the neutron state. PPD tells us that the axial current with momentum $\mathit{q}$ can be viewed as the insertion of a pion with $\mathit{q}$, and this has a large coupling to the nucleon. These processes are illustrated in Figure 4.

Having identified large contributions from the $N(\mathbf{p}=1)\pi (\mathbf{p}=-1)$ state, certainly in the extraction of ${\tilde{G}}_P$ and $G_P$, the question is—do we need to include other multihadron and radial excited state contributions if we want results with percent level precision? What about in $G_A$? Note that, in addition to the enhanced contribution shown in Figure 4, $\chi$PT also indicates that the one-loop contributions due to the diagram shown in Figure 5 (again a $N\pi$ contribution) could be $O(5\%)$ in all the five $C_J^{\mathrm{3pt}}$. Thus, the $N\pi$ state could be significant for extracting $G_A$ (the GSME $\langle 0^{\prime }\left|A_3\right|0\rangle$) from $C_{A_3}$ at the percent precision desired. Based on these arguments, it is clear that one needs at least three-state fits to the five $C_J^{\mathrm{3pt}}$ in Equation (3)—the ground state, the $N\pi$ state, and one other that effectively accounts for all other excited state contributions.

The caption of Figure 2 points out some of the features of the ESC observed in current data and the efficacy of fits to the spectral decomposition of $C_J^{\mathrm{3pt}}$ with and without including the lowest $N\pi$ state to remove the ESC.

In my evaluation, the details of the fits made to remove ESC are the most significant differences between the calculations performed by the different collaborations. With the current methodology, higher statistics data are needed to improve these fits and reduce the dependence on exactly how the analyses are done.

A very important point to keep in mind is that fit parameters in a truncated ansatz (the $A_i$ and $E_i$ in say a three-state fit in our case) try to incorporate the effects of all contributions. Thus, the connection between parameters coming out of fits to a truncated Equation (3) and physical states made in Figure 3 and Figure 4 are very approximate at best.

2.4. Satisfying PCAC

The nonsinglet PCAC relation between the bare axial, $A_{\mu }\left(x\right)$, and pseudoscalar, $P\left(x\right)$, currents is the following:

$${\partial }_{\mu }{A}_{\mu }=2\widehat{m}P,$$

(11)

where the quark mass parameter $\widehat{m}\equiv Z_m{m}_{ud}{Z}_P{Z}_A^{-1}$ 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 Equations (4) and (5) of the GSME, the PCAC relation requires that the three form factors $G_A$, ${\tilde{G}}_P$, and $G_P$ must satisfy on each ensemble up to discretization errors; thus, the relation is

$$2M_N{G}_A\left(Q^2\right)-\frac{Q^2}{2M_N}{\tilde{G}}_P\left(Q^2\right)=2\widehat{m}{G}_P\left(Q^2\right).$$

(12)

All prior Ref. [35] calculations did not check this relation and missed observing that the data showed large deviations. Calculations subsequent to Ref. [38] that include the lowest mass gap state $N(\mathbf{p}=1)\pi (\mathbf{p}=-1)$ in the analysis obtain form factors that already satisfy the PCAC to within ≈10% at a lattice spacing of $a\approx 0.09$ fm. (The ETMC result is an exception, as explained in Ref. [23] ). An illustration of the size of the deviation from unity of $R_1+R_2\equiv \frac{2\widehat{m}{G}_P\left(Q^2\right)}{M_N{G}_A\left(Q^2\right)}+\frac{Q^2{\tilde{G}}_P\left(Q^2\right)}{4M_N^2{G)}_A\left(Q^2\right)}$, without and with the lowest $N\pi$ state included, is shown in Figure 6 taken from Ref. [19].

To summarize, satisfying the PCAC relation in Equation (12) provides a strong and necessary constraint on the extraction of the three axial form factors. $\chi$PT analysis by Bär [36,37] and data-driven validation in Refs. [29,30,38] show that the lowest $N(\mathbf{p}=1)$$\pi (\mathbf{p}=-1)$ state makes a large contribution and needs to be included in the analysis. For percent level precision, the next question is—what other states need to be included? Current analyses include up to three states, where the third state, if the parameters are left free, is effectively trying to account for all of the residual ESC. Such fits have been implemented in different ways. For example, in Ref. [29], the $N\pi$ state is hardwired, and the third state is taken to be the lowest excited state in fits to $C^{\mathrm{2pt}}$. In Refs. [19,27,38], a simultaneous fit to all five $J=A_{\mu }$ and P correlators is made wherein the $A_4$ correlator fixes $E_1$ to being close to the $N\pi$ state. Over time, with much higher statistics data, results from various methods and collaborations should converge.

2.5. Extrapolating Lattice AVFF to The Physical Point for Use in Phenomenology

The next step, once ESC have been removed and form factors have been extracted from GSME on a given ensemble, is to extrapolate these data to the physical point to provide a parameterized form for $G_A$ and ${\tilde{G}}_P$ that can be used in phenomenology. The challenge is that the discrete set of $Q_i^2$ values at which the data are obtained are different on each ensemble.

One simple way to implement this consists of the following three steps:

• Parameterize the $G_A\left(Q_i^2\right)$ data on each ensemble. Depending on the number of $Q_i^2$ values, it could be a suitably truncated z expansion or a Padé. The expected $1/Q^4$ asymptotic behavior can be built in by using sum rules in the z expansion [39] or through an $\{n,n+2\}$ Padé in $Q^2$. The Mainz collaboration [26] combines the removal of the ESC at various values of $Q_i^2$ and the $Q^2$ parameterization on a given ensemble to include correlations.
• Pick n values, $Q_k^2$, over a range of say $0\le Q^2\le 1$ GeV². Extrapolate the data at each of these values of $Q_k^2$ using a simultaneous fit in $\{m_q,a,M_{\pi }L$} to the physical point. (Note that I use $M_{\pi }$ or $m_q$ to define the light quark mass interchangeably.) A typical ansatz used for such chiral continuum finite volume (CCFV) extrapolations is

  $$\begin{array}{cc}\hfill g(M_{\pi },a,M_{\pi }L)& =g_0+c_{1}a+c_2{M}_{\pi }^2\hfill \\ \hfill & +\frac{c_3{M}_{\pi }^{2}exp(-M_{\pi }L)}{\sqrt{M_{\pi }L}},\hfill \end{array}$$

  (13)

  where I have kept only the lowest order corrections in each of the $\{m_q,a,M_{\pi }L\}$ variables and assumed that discretization errors start at $O\left(a\right)$.
• Having obtained the form factor in the continuum limit at $Q_k^2$ points with $k=1\dots n$, carry out the final parameterization again using a truncated z expansion or a Padé.

This three-step process can be done within a single bootstrap procedure to propagate errors as has been done in Refs. [19,27] to produce the NME and PNDME results shown in Figure 7. Or, these steps can be combined, especially if there are correlations between them. For example, as can be done in step (ii) to account for correlations between the coefficients of the CCFV fits for different values of $Q_k^2$.

The plots in Figure 7 provide two comparisons. In the panels on the left, the physical point results from the RQCD [24,29], ETMC [23], NME [27], and Mainz [26] collaborations are compared against those from the PNDME [19]. On the right, they are overlaid and compared to the phenomenological extraction from the old neutrino–deuterium bubble chamber data [20]. The PNDME, RQCD, and NME data mostly overlap, whereas the ETMC and Mainz data overlap and fall off slower for $Q^2\gtrsim 0.3$ GeV². On the other hand, the neutrino–deuterium ($\nu$D) data [20] fall off much faster for $Q^2\gtrsim 0.2$ GeV². Overall, as shown in the right plot, the five lattice QCD estimates are consistent within $1\sigma$ and lie about $2\sigma$ above the $\nu$D band for $Q^2\gtrsim 0.3$ GeV².

There also are results from the CalLAT [16], PACS [25,40], and LHP+RBC+UKQCD [41] collaborations, which have not been included in the comparison because they have not been extrapolated to the physical point. The Fermilab collaboration [31] has embarked on the much harder problem of calculating transition matrix elements as well, e.g., $N\to \Delta$ or $N\to N\pi$.

From the analysis of the NME and PNDME data, my understanding is that the differences in exactly how the ESC values are handled by the various collaborations and the uncertainty in the final results should be considered as a work in progress. The uncertainty from the differences in the overall procedure for parameterization and CCFV extrapolation is, I believe, smaller, especially since the data do not show a large dependence on any of the three parameters $\{m_q,a,M_{\pi }L\}$, especially for $a\lesssim 0.1$ fm and $M_{\pi }L\gtrsim 4$, as illustrated in Figure 8 [19,27]. Hopefully, the next generation calculations will shed light on and possibly resolve the various differences.

Other findings in Refs. [19,27] are that (i) the dipole ansatz $G_A\left(Q^2\right)=\frac{g_A}{{(1+c{Q}^2/M_N^2)}^2}$ gives poor fits (very low p values) to data on many ensembles. My conclusion, therefore, is that the lattice data already show that the dipole ansatz does not have enough parameters to capture the $Q^2$ behavior over the range $0\le Q^2\le 1$ GeV². The second finding is that (ii) The PPD relation between $G_A$ and ${\tilde{G}}_P$ works very well.

2.6. Consistency Check in the Extraction of the Axial Charge $g_A^{u-d}$

There are two ways in which one can extract the axial charge $g_A^{u-d}$. The first is from the forward matrix element using $C_{A_3}$ in Equation (8) with $\mathit{q}=0$, and the second is by extrapolating the form factor $G_A(Q^2\ne 0)$ to $Q^2=0$. I am considering them as separate because the extraction from the forward matrix element is computationally clean: $C_{A_3}(\mathit{q}=0)$ has the smallest errors, and the verification of the symmetry of the data about $\tau /2$ is a good test. The errors grow with $\mathit{q}$, as shown in Figure 2. On the other hand, $G_A\left(Q^2\right)$ is constrained by being part of the PCAC relation, Equation (12), that has to be satisfied. The two results must agree after CCFV extrapolation. Based on the data in Ref. [19], I conclude the following:

• The difference between $g_A^{u-d}$ extracted without and with including $N\pi$ states is $O(\approx 5\%)$, i.e., $1.218\left(39\right)\to 1.294\left(48\right)$ on including one (the lowest) $N\pi$ state in the analysis. Note that the errors in each result are $O(\approx 3\%)$.
• The difference between $g_A^{u-d}=G_A(Q^2\to 0)$ extracted by extrapolating $G_A\left(Q^2\right)$ data obtained without and with including the lowest $N\pi$ state is also $O(\approx 5\%)$, i.e., $1.213\left(39\right)\to 1.289\left(56\right)$. Again, the errors in each are $O(\approx 4\%)$.

• Thus, for each of the two cases, without and with including $N\pi$ states, we obtain consistent estimates for the charge from the two methods; however, the results including the $N\pi$ state are about $6\%$ larger. This difference is consistent with the expected $\sim 5\%$ one-loop correction to the charge in $\chi$PT; however, it is roughly a one combined $\sigma$ effect and, therefore, needs validation. My pick for the final result is the analysis including the $N\pi$ state, since it gives form factors that satisfy the PCAC relation. Higher precision data are needed to further clarify the other significant ESCs and how to include them.

3. Comparison of Charges Obtained by Various Lattice Collaborations

The results for the axial charge, $g_A^{u-d}$, the charge radius squared, $\langle r_A^2\rangle$, the induced pseudoscalar coupling $g_P^{\ast }$, and the pion–nucleon coupling $g_{\pi NN}$ extracted from $G_A$ and ${\tilde{G}}_P$ by various collaborations using the relations result in the following:

$$\begin{array}{cc}\hfill G_A\left(Q^2\right)& =g_A(1-\frac{\langle r_A^2\rangle }{6}{Q}^2+\cdots ),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill g_P^{\ast }& \equiv \frac{m_{\mu }}{2M_N}{\tilde{G}}_P\left(Q^{\ast 2}\right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill g_{\pi NN}& \equiv \underset{Q^2\to -M_{\pi }^2}{lim}\frac{M_{\pi }^2+Q^2}{4M_N{F}_{\pi }}{\tilde{G}}_P\left(Q^2\right).\hfill \end{array}$$

(16)

which are summarized in

Table 1. Comparison of $g_A$, $\langle r_A^2\rangle$, $g_P^{\ast }$, and $g_{\pi NN}$ from recent calculations: PNDME 23 [19], RQCD [24,29], ETMC [28], PACS [25,40], Mainz [26], and NME [27]. Results for $g_A^{u-d}$ are in the $\overline{\mathrm{MS}}$ scheme at scale 2 GeV. Earlier results for $g_A^{u-d}$ have been reviewed in the FLAG reports [21,22].

Collaboration	${\mathit{g}}_{\mathit{A}}^{\mathit{u}-\mathit{d}}$	$\langle {\mathit{r}}_{\mathit{A}}^2\rangle$ fm2	${\mathit{g}}_{\mathit{P}}^{\ast }$	${\mathit{g}}_{\mathit{\pi }\mathbf{NN}}$
PNDME 23	1.292(53)(24)	0.439(56)(34)	9.03(47)(42)	14.14(81)(85)
RQCD 19/23	$1.{284}_{27}^{28}$	0.449(88)	8.68(45)	12.93(80)
ETMC 23	1.283(22)	0.339(46)(6)	8.99(39)(49)	13.25(67)(69)
PACS 23	1.264(14)(1)	0.316(67)
Mainz 22	1.225(39)(25)	0.370(63)(16)
NME 21	1.32(6)(5)	0.428(53)(30)	7.9(7)(9)	12.4.(1.2)

. Here, $m_{\mu }$ is the muon mass, $Q^{\ast 2}=0.88m_{\mu }^2$ is the energy scale of muon capture, and $F_{\pi }=92.9$ MeV is the pion decay constant.

The results show about $10\%$ variation in $g_A^{u-d}$, $g_P^{\ast }$, and $g_{\pi NN}$ and about $25\%$ in $\langle r_A^2\rangle$. Part of this is likely due to different methodologies used in the analysis, in particular how and if the lowest $N\pi$ state is included in the analysis. These results will improve over time.

4. Comparison of the Differential Cross-Section Using Lattice AVFF with MINER$\mathbf{\nu }$A Data

A comparison of the antineutrino–nucleon charged current elastic cross sections calculated using predictions of AVFF from lattice (PNDME 23 [19]) and neutrino–deuterium analysis [20] with MINERvA measurement [17,42] is presented in Figure 9, which is reproduced from Ref. [18]. A ${\chi }^2$ test was performed to determine the significance of the differences between the three. No significant difference was found between MINERvA–lattice QCD (PNDME) and between MINERvA–deuterium results. A $\approx 2.5\sigma$ tension was, however, found between the PNDME and the deuterium results. Based on data shown in Figure 7, the deviation in the deuterium–Mainz and deuterium–ETMC will be even larger. To assess the scope for future progress, three regions of $Q^2$ with different prospects for the extraction of AVFF from lattice QCD and MINERvA-like experiments were identified.

For $Q^2\lesssim 0.2{\mathrm{GeV}}^2$, LQCD predictions and fits to the deuterium bubble chamber data are in good agreement. In this region, the experimental errors in the measurement on hydrogen by MINERvA are large, whereas the errors in the parameterization of the deuterium bubble chamber data are smaller. The $\nu$D result has often been used as a benchmark; however, note that there is unresolved uncertainty in the deuterium data as discussed in Ref. [20]. Also, no new deuterium data are expected in the near-term, so I do not comment on its future prospects. Lattice QCD data are competitive and will improve steadily. This region will be well characterized by the axial charge, the axial charge radius, and well parameterized by a low-order z expansion or a Padé.

For $0.2{\mathrm{GeV}}^2\lesssim Q^2\lesssim 1{\mathrm{GeV}}^2$, the AVFF from PNDME has the smallest errors, and the predicted differential cross-section lies above the hydrogen and $\nu$D values, i.e., the same ordering as for the AVFF shown in Figure 9 (left). Future improvements in both the hydrogen data and lattice calculations will provide robust crosschecks in this region.

The region $Q^2\lesssim 0.5{\mathrm{GeV}}^2$ is where lattice QCD data, even with current methodology, will improve rapidly as more simulations are done closer to $M_{\pi }=135$ MeV, $a\to 0$ and on larger volumes, because on a given ensemble and for given statistics, the value of $Q^2{|}_{\max }$ with a good signal (usually specified by the lattice momentum ${\mathbf{n}}^2$) decreases in these limits.

For $Q^2\gtrsim 1{\mathrm{GeV}}^2$, current LQCD data have large statistical errors and systematic uncertainties—these are discretization and residual excited state contributions. With the current methodology, the lattice AVFF comes mostly from simulations with $M_{\pi }\gtrsim 300$ MeV [19]. New methods are needed to obtain data from simulations with physical pion masses, $M_{\pi }\approx 135$ MeV. Similarly, improvements further MINERvA and follow-up experiments are needed to cover the full range of $Q^2$ relevant for DUNE.

5. Concluding Remarks

Extensive calculations of the AVFF are being carried out by at least the following nine lattice QCD collaborations: PNDME [19], RQCD [24,29], ETMC [23], NME [27], Mainz [26], CalLAT [16], PACS [25,40], LHP+RBC+UKQCD [41], and Fermilab [31]. As shown in Section 3 and Section 2.5, we now have results to within 10% precision. The major of the uncertainty comes from resolving and removing excited state contributions.

The good news is that the methodology for the calculation of the correlation functions, $C^{2\mathrm{pt}}$ and $C_J^{3\mathrm{pt}}$, is robust. The bad news is that the exponentially falling signal-to-noise ratio in them means that ESC are large at source–sink separations possible in today’s calculations. Second, it is also clear that multihadron, $N\pi$, excited states give large contributions and must be included in the analysis to remove them. Unfortunately, it is not yet known how many of these states will need to be included in the analysis for percent level precision. The operator constraint that the form factors satisfy the PCAC relation in Equation (12) provides a valuable check that must be carried out in all calculations. The third challenge is obtaining data at large $Q^2$, because the discretization and statistical errors grow with $Q^2$ on a given ensemble, and the $Q^2{|}_{\max }$, i.e., the lattice momenta $2\pi \mathit{n}/La$ with the largest ${\mathit{n}}^{\mathbf{2}}$ that has a good signal-to-noise ratio, decreases as simulations are done closer to the physical point. Thus, to obtain data for $Q^2\gtrsim 1$ GeV² on physical pion mass ensembles will need/benefit from new methodology.

I estimate that a factor of ten increase in statistics will reduce the statistical errors to a level that will provide much more clarity in removing the ESC. Similary new developments, including variational methods [30] with multihadron states and momentum smearing [43], will improve the calculations and extend the range of $Q^2$. I anticipate improvements in both statistics and methods will provide LQCD predictions of AVFF for nucleons in the range $Q^2\lesssim 1{\mathrm{GeV}}^2$ (hopefully higher) with percent level precision by about 2030 in concert with DUNE producing data.