Computing Nucleon Charges with Highly Improved Staggered Quarks

This work continues our program of lattice-QCD baryon physics using staggered fermions for both the sea and valence quarks. We present a proof-of-concept study that demonstrates, for the first time, how to calculate baryon matrix elements using staggered quarks for the valence sector. We show how to relate the representations of the continuum staggered flavor-taste group $\text{SU}(8)_{FT}$ to those of the discrete lattice symmetry group. The resulting calculations yield the normalization factors relating staggered baryon matrix elements to their physical counterparts. We verify this methodology by calculating the isovector vector and axial-vector charges $g_V$ and $g_A$. We use a single ensemble from the MILC Collaboration with 2+1+1 flavors of sea quark, lattice spacing $a\approx 0.12$ fm, and a pion mass $M_\pi\approx305$ MeV. On this ensemble, we find results consistent with expectations from current conservation and neutron beta decay. Thus, this work demonstrates how highly-improved staggered quarks can be used for precision calculations of baryon properties, and, in particular, the isovector nucleon charges.


I. INTRODUCTION
Accurate first-principles calculations of nuclear cross sections are an important objective in the particle physics community. In particular, heavy nuclei, such as 12 C and 40 Ar, are used as targets in neutrino-scattering and darkmatter detection experiments. In calculations of cross sections, a necessary component is the modeling of nuclei as a collection of nucleons, opening up an opportunity for lattice QCD [1]. At the quasielastic peak, for example, the electromagnetic and axial-vector form factors of the nucleon, which characterize the electric charge and spin distribution within the nucleon, are key ingredients. Such form factors can be obtained from the first-principles lattice-QCD framework. However, these hadronic inputs remain one of the largest sources of systematic error as the experimental precision on these cross-sections continues to improve [2][3][4].
The electromagnetic form factors have been extracted precisely from high statistics electron-nucleon scattering experiments [5,6]. At zero momentum transfer, the proton's electric form factor becomes the total electric charge g V = 1, and the slope at the origin is related to the charge radius. Recently, experiments that make use of the Lamb shift of muonic hydrogen report significantly smaller proton radii than those measured via scattering [7]. (For recent reviews of the proton radius puzzle, see Refs. [8,9].) In addition, a recent reanalysis has demonstrated that the vector form factors at intermediate Q 2 also exhibit tensions outside of their quoted uncertainties [10]. These disagreements could benefit from better knowledge of the Standard Model predictions, which necessitates using lattice QCD to calculate the form factor.
In comparison, the nucleon axial-vector form factor is much less constrained from experimental data. A recent re-analysis [11] of the deuterium bubble-chamber data found greater uncertainties than previously assumed. Again, lattice QCD can be illuminating here, computing the axial-vector form factor from first principles as an independent check on the form factor extracted from experimental data. At zero momentum transfer, the axialvector form factor gives the so-called nucleon axial charge g A = 1.2756 (13), which has been measured precisely in neutron beta decay [12]. Thus, the axial-charge can be used to validate lattice-QCD calculations before studying the momentum dependence of the form factor. In addition, a percent-level first principles calculation of g A could shed light on the neutron lifetime puzzle [13].
Lattice-QCD calculations of baryonic observables are hindered by the well-known exponential growth of the noise relative to the signal, which sets in at large times [14,15]. At early times, where the signal-to-noise ratio is favorable, the lattice-QCD correlator data contain significant contributions from several states in an infinite tower. When using a fit to disentangle the higherlying states from those of interest, some residual, unwanted contamination remains in the parameters of interest. It is imperative, therefore, to demonstrate control over both the noise and the excited-state contamination.
In this work, we use an ensemble generated by the MILC Collaboration [16], which incorporates a sea with equal-mass up and down quarks, the strange quark, and the charm quark. MILC uses the highly improved staggered-quark (HISQ) action [17] for the sea quarks; here we use the HISQ action for the valence quarks too. Because staggered fermions have only one component per site and retain a remnant chiral symmetry, they are computationally efficient. Nevertheless, staggered fermions are complicated by the fermion doubling problem, leading to four species, known as tastes, for each fermion field. The four tastes become identical in the continuum limit, leading to an SU(4n f ) flavor-taste symmetry for n f flavors. Consequently, the spectrum of staggered lattice baryons is rich and intricate. For nucleons, the spectrum has been classified [18][19][20], finding many states that have the same properties as the physical nucleon.
In a recent paper, we used staggered baryons to calculate the nucleon mass [20]. Computing nucleon charges is the next step and a necessary one en route to the full momentum dependence of the form factors. As discussed in Ref. [20], it can be advantageous to use unphysical nucleon-like states to carry out the calculation. These states obtain the same properties as the physical nucleon in the continuum limit, where the full SU(8) F T flavor(isospin)-taste symmetry emerges. For matrix elements such as charges and form factors, however, one must find the correct group-theoretic normalization factors relating nucleon-like matrix elements to their physical counterparts. This exercise is a straightforward if complicated application of the generalization of the Wigner-Eckart theorem to SU (8).
To demonstrate this approach, we compute the nucleon vector and axial-vector charges on a single MILC HISQ ensemble with lattice spacing a ≈ 0.12 fm and pion mass M π ≈ 305 MeV. We employ local vector and axial-vector currents. We also outline the steps needed to apply this method to matrix elements of other baryons, with an eye to future studies including staggered baryons, such as N → ∆ transition form factors. This paper is organized as follows. In Sec. II, we discuss staggered-baryon correlators, starting with a brief review of the two-point correlator methodology [20]. We then present an overview of our three-point correlators.
Here, we also present one of the key results of this paper: the correct normalization of the nucleon-like matrix elements. In Sec. III, we describe strategies for removing excited-state contamination. Section IV provides the details of our simulation, while Sec. V describes Bayesian fits to the correlator data. Our computational results are presented in Sec. VI, including the robustness of our results under variations of our fitting procedure, the renormalization of the bare lattice charge to the physical charges, and the final values for g V and g A on the single ensemble being used. Finally, we compare our results to mixed-action results on the same ensemble and provide our conclusions in Sec. VII. Appendices A and B present the group theory relating the nucleon-like matrix elements to their physical counterparts, including a numerical demonstration that these derivations are correct.

II. STAGGERED BARYON CORRELATORS
For simplicity, we focus here on two flavors, up and down, with isospin symmetry. With staggered fermions, instead of the usual SU(2) F isospin symmetry, an enlarged SU(8) F T flavor-taste symmetry group emerges in the continuum limit. It is important to note that the irreducible flavor-taste representations contain components with non-trivial taste and unexpected isospin. For example, Bailey has shown [19] that nucleon-like states exist with unphysical isospin yet masses equal in the continuum limit to the physical tasteless nucleon. In fact, all physics of such nucleon-like states can be related to that of the physical nucleon. In particular, here we show how to relate nucleon-like matrix elements to their physical counterparts. As such, we are allowed to choose any nucleon-like representation, for example, one that reduces the computational complexity.
We use the isospin-3 2 operators that transform in the 16 irrep of the geometric timeslice group (GTS) [18,21], as presented in Ref. [20]. They are less complicated to analyze because only a single nucleon-like taste appears in the spectrum. On the other hand, this irrep contains contributions from three ∆-like tastes.
A. Two-point correlators Using the same notation as in Ref. [20], the two-point correlators read s D (0) defined in Ref. [20]. To increase the statistical precision, we average over the eigenvalues s = ± of the staggered rotation in the x-y plane, and also the eight corners of the cube D; together, s and D label the components of the 16 irrep. Here, r 1 , r 2 = 2, 3, 4, 6 represent four different operator constructions, or "classes" [18][19][20], as well as other possible properties, such as smearing.

B. Staggered-baryon matrix elements
In this work, we are specifically interested in the isovector nucleon vector and axial-vector charges, namely g V and g A , respectively. These are defined through the nucleon matrix elements where Γ A = γ z γ 5 or Γ V = γ 4 , u and d are continuum-QCD up-and down-quark fields, and u N is the nucleon spinor at zero momentum. We calculate these nucleon matrix elements using (highly improved) staggered quarks. To achieve this, we must extend the mass relations of Bailey [19] to matrix elements. The baryon-like matrix elements and the physical matrix elements are related through symmetry transformations in the continuum. In the appendices, we find the appropriate Clebsch-Gordan coefficients that relate the single-taste baryon matrix elements to the physical tasteless QCD matrix elements by applying the generalized Wigner-Eckart theorem of SU(8) F T .
The correctly normalized three-point correlators for our baryon-like operators are then where t is the source-sink separation time and τ is the current insertion time. The factor −1 in front of C V and the factor −3 in front of the second term of C A come from the group theory just described.
where χ f is the field in the HISQ action of flavor f . The local vector and axial-vector currents, V and A, have spin-taste γ 4 ⊗ ξ 4 and γ z γ 5 ⊗ ξ z ξ 5 [21], respectively. The opposite-parity partners then arise from spin-taste γ 5 ⊗ ξ 5 and γ z γ 4 ⊗ ξ z ξ 4 , respectively. Because these local currents are not derived from Noether's theorem, they require a finite renormalization, that is, Z V V and Z A A have the same matrix elements as the continuum isovector currents in Eq. (2.2).
In the limit τ → ∞ and t − τ → ∞, the ratio of the three-point to the two-point correlators approaches the desired nucleon charge where g J is now the bare lattice charge, that is g J = Z J g J . In practice, of course, we compute the correlators for several values of t and τ and fit the t and τ dependence to extract the charges. The finite renormalization factors Z J are determined first by noting that the remnant chiral symmetry requires Z A = Z V + O(m q a) 2 . At zero momentum transfer, the vector current simply counts the number of up quarks minus the number of down quarks, for the proton 2 − 1 = 1. One could, thus, define Z V by demanding Z V g V = g V = 1. Here, however, we prefer to define Z V via a similar relation obtained from a pseudoscalar-meson matrix element [22], then use the result to renormalize our nucleon matrix elements. With this choice our result for g V is a genuine test of our methodology.

III. EXCITED-STATE CONTAMINATION
Excited-state contamination is one of the most difficult challenges when accurately estimating nucleon ma-trix elements from lattice QCD. The problem is even more complicated with staggered nucleons because of the presence of negative-parity and low-lying ∆-like states in the spectrum, which non-staggered formulations do not contain. These are both significant sources of excitedstate contamination in the present calculation. We have, however, demonstrated control of excited-state contamination when extracting nucleon physics from two-point staggered-baryon correlators [20]. Here we describe extensions of those techniques to the three-point correlators of the present work. In particular, we show how to suppress contributions from the lowest-lying negative-parity states and the lowest three ∆-like tastes.

A. Negative parity states
Let C 2pt (t) and C 3pt (t, τ ) be any staggered-baryon correlators. The source-sink separation is denoted t and the current insertion time is denoted τ . Any staggered operator that is local in time will create negative parity states, which in turn causes the characteristic oscillations in time. This is obvious from the correlators spectral decomposition where M ± are the lowest-lying ± parity masses,z ± and z ± are, respectively, the source and sink overlap factors for states of parity ±, and A ±± and A ±∓ are the transition matrix elements. For simplicity, we have ignored backward propagating terms proportional to e −M±(T −t) , which are assumed to contribute negligibly in the following. Equation (3.2) shows that the terms involving negative parity states change sign when either t/a or τ /a change by one unit. With this in mind, a time-averaging procedure can be applied to suppress the negative parity contributions to the correlator. A similar scheme was deployed in Ref. [23]. The first ingredient is where we call aM snk the time-averaging parameter. Substituting this expression into the spectral decomposition in Eq. (3.2), one sees that the functional forms of primed correlators are unchanged except that the sink overlap factors becomes If one chooses aM snk = aM − , then terms with the M − state at the sink will vanish, while the overlap factors for the positive parity states become slightly larger. In practice, the time-averaging parameter does not need to be exact to suppress the negative-parity states. A similar time-averaging parameter, aM src , can be introduced to reduce the negative parity contributions at the source via C 3pt (t, τ ) = e −aMsrc C 3pt (t, τ ) + C 3pt (t + a, τ + a). If several negative parity states contribute significantly to the data, successive applications of this procedure, with suitable parameters [aM (1) src , aM (2) src , . . .] and [aM (1) snk , aM (2) snk , . . .], can appreciably suppress them. On the other hand, because the relative error in the correlators becomes larger with time, too much time-averaging renders the data statistically less precise. Moreover, time-averaging reduces the available τ range in the modified correlators, thereby producing fewer data for the fit. For each data set, some study is necessary to strike an optimal balance.

B. ∆-like states
Another source of excited-state contamination arises from the presence of the three ∆-like states in the 16-irrep correlators. With four different classes of interpolators at both the source and the sink, we adopt the strategy from Ref. [24] and solve the generalized eigenvalue problem (GEVP) [25]. In Ref. [20], we applied the GEVP to our two-point correlators and successfully disentangled the nucleon-like state from the ∆-like states. We extend that strategy to the three-point functions here.
Given a matrix two-point correlator, C 2pt (t), the left and right nucleon eigenvectors, u(t 1 , t 0 ) and v(t 1 , t 0 ), are the solutions of Here, we focus on the eigenvectors for the nucleon-like state, the ones with the lowest eigenvalues, and put the others aside. These eigenvectors optimize the projection onto the nucleon-like state in both the two-and threepoint correlators via (3.14) One has to decide which t 1 and t 0 to use in Eqs. (3.11) and (3.12). The stability of our results under such variations will be discussed in Sec. IV. Below we call the correlators in Eqs. (3.13) and (3.14) the nucleon-optimized two-and three-point correlators.
To summarize our strategy, we start with the correlators in Eqs. (2.1), (2.4), and (2.3), and apply two iterations of time-averaging at both the source and sink, and then project the time-averaged correlation matrix as in Eq. (3.14). The time-averaging suppresses the negative parity states contributions, and the projection suppresses the ∆-like baryons contributions.

IV. SIMULATION DETAILS
To demonstrate the feasibility of nucleon matrix elements with staggered quarks, we use a single gauge ensemble, which was generated by the MILC collaboration [16]. MILC implemented the one-loop, tadpoleimproved Lüscher-Weisz gauge action [26], as well as the HISQ action [17] for the sea, which contains equal-mass up and down quarks, the strange quark, and the charm quark. In this work, we also employ the HISQ action for the valence quarks, with the same mass as the up-down sea quarks.
The ensemble has dimension L 3 × T = 24 3 × 64, a lattice spacing a = 0.1222(3) fm (determined from the F p4s mass-independent scheme [27]), a pion mass M π ≈ 305 MeV, and a light-to-strange-quark mass ratio of 1/5. Other parameters of this ensemble are listed in Ref. [27]. Note that the CalLat [28] and the PNDME [29] collaborations have both used this same ensemble to calculate g A , albeit with either the Möbius domain wall or Wilson-clover valence fermion actions, respectively.
We generate the two-and three-point correlators according to Eqs. (2.1), (2.3), and (2.4). We measure each correlator on 872 configurations, and further increase the statistics by randomly placing the corner-wall sources on eight maximally separated timeslices to give a total of 6976 measurements per correlator. (The τ /a = 7 correlators have only four time sources per configuration. ) We block all measurements in a single gauge configuration and every four consecutive gauge trajectories to avoid autocorrelations. The covariance matrix between different correlator components are estimated with the non-linear shrinkage method [30] to avoid ill-conditioning from finite sample sizes.
As described in Ref. [20], we use corner-wall sources to optimize the signal-to-noise ratio and point sinks. In the present work, we remove the Coulomb-gauge fixed links, as we have empirically observed that leaving out the links has little effect on correlators but with the added advantage of a simpler code. Here we also incorporate the Wuppertal smearing [31,32] at the sink by applying in order to reduce excited state contamination. In Eq. (4.1), n is the n th iteration of N total iterations; all shifts are stride 2 to preserve the staggered symmetries. We include the appropriate gauge transporters to make the smearing gauge covariant [20], but for succinctness they are omitted from Eq. (4.2). We generate data with two different root-mean-squared (rms) smearing radii, σ, which are about 0.2 and 0.6 fm. We label the two smearings as Gr2.0N30 and Gr6.0N70. For the three-point correlators, we invert the propagators from the current insertion to obtain all operator classes at the sink. Five current insertion times, τ /a = [3,4,5,6,7], are generated for both the vector and axial-vector current insertions.
To eliminate the unwanted negative parity states, we then pass all the correlators through two iterations of time-averaging using Eqs. (3.4) and (3.7) with [aM (1) src , aM (2) src ] = [aM (1) snk , aM These two numbers are based on an observation in Ref. [20] that the lowest-lying negative parity state seems to have energy around the S-wave N π state, which in this ensemble is about 0.9 in lattice units. We then set the second averaging parameters about aM π ∼ 0.2 higher than the first ones, which again is consistent with our findings in Ref. [20]. As the goal is to suppress the negative-parity states, the accuracy of these parameters is not crucial. Note that each iteration of the source timeaveraging in Eq. (3.8) reduces the current insertion timeslices by one, so the time-averaged three-point correlators have only τ /a = [3,4,5]. The smearing in Eq. (3.4), on the other hand, does not reduce the range for τ /a. After time-averaging, we solve for the left and right eigenvectors using Eq. (3.12) in order to optimize our correlators as in Eq. (3.14). To ensure the robustness of our fitting methodology, we test the stability of our results under variations of the choice of t 0 . To do so, we compute the effective mass of the optimized two-point correlators, which we define as The t 0 stability plots are shown in Fig. 1. All choices produce similar results, and so we choose t 0 = 6a for the subsequent analyses. Similarly, we vary (t 1 − t 0 )/a from 2 to 6 and find, again, that the differences are negligible, so we fix t 1 − t 0 = 2a.
Since we have normalized the nucleon-like three-point correlators correctly, and each of the correlator transformations that we perform preserve the functional form of the spectral decomposition, the optimized three-totwo-point correlator ratios converge to the desired nucleon charges in the large-time limits. In Figs. 2 and 3, we plot the ratio of the nucleon-optimized three-to-twopoint correlators, with and without time-averaging at the source and sink. The left column shows the optimized correlators without time-averaging, and the right column shows them time-averaged with the parameters given in Eq. (4.3). The two smearing are shown in the top (Gr2.0N30) and bottom (Gr6.0N70) rows. Significant oscillations are clearly present in the unaveraged correlators, particularly for the vector current in Fig. 2. This is expected, because the parity partner of the vector current is the pseudoscalar current P , and the N π|P |N matrix element gives a large contribution to the vector-current data and, thus, causes large oscillations. For the g A data, on the other hand, the parity partner of the axial current is the tensor current T z4 , and when the nucleon is at rest N |T z4 |N = 0. Consequently, the first non-zero contribution in the axial-vector parity partner channel will likely be from N π|T z4 |N , leaving small oscillations.

V. CORRELATOR FITTING
We apply the Bayesian fitting methodology implemented in corrfitter [33] to extract the nucleon mass and matrix elements. We observe in Fig. 2 that the vector correlators have noticeable oscillatory contributions, whereas the axial-vector correlators shown in Fig. 3 do not. Further, the vector correlators seem relatively insensitive to our choice of Wuppertal smearing. We perform separate fits to the vector and axial-vector correlators, but include their correlations through bootstrapping.
It should be stressed that, after applying the excitedstate suppression techniques from Sec. III, the interpretation of the higher exponentials in the correlators is ambiguous. For the positive-parity channel, the first "excited state" could be a mixture of any leftover ∆-like states, the P-wave N π states, or other finite volume energy levels higher up in the spectrum that are related to resonances. For the negative-parity channel, we found in Ref. [20] that the ground state is likely to contain S-wave N π states. The time-averaging procedure to cancel out

Quantity
Prior value ± width M+0 = nucleon mass 1100 the negative-parity states makes identification of these states even more ambiguous. Regardless of the origin, we can treat the excited states as nuisance parameters and fit them away with an exponential fit function. In this case, each excited exponential mass parameter describes a conglomeration of several eigenstates of the Hamiltonian. Still, we will refer to each exponential in the fit function as a state without necessarily identifying it with any single eigenstate. As discussed in Ref. [20], the stability of the extracted fit parameters as a function t min indicates lack of excited-state contamination, as long as they are modeled accurately. This t min stability plot is shown in Fig. 5, and is discussed in Sec. V B.

A. Functional forms for fitting
For the g A analysis, we perform simultaneous fits to the optimized two-and three-point correlators, and include both the Gr2.0N30 and Gr6.0N70 sink smearings. Observation of the strong suppression of excited states in Fig. 2 leads us to use a fit ansatz that contains two positive-parity states and one negative-parity state: Here sink overlap factors. The terms A ±i,±j are the unrenormalized axial-vector matrix elements, with A +0,+0 = g A the desired bare axial charge. Note that the two-point correlator terms involving finite temporal T extent are not included here since we average our data symmetrically around the T /2 point as described in Ref. [20].
In the Bayesian fit, we choose Gaussian priors for the ground-state masses, overlap factors, and matrix elements. We choose log-normal priors for the mass differences between adjacent states to enforce the ordering of states. It has been observed empirically that the nucleon mass has an approximate linear dependence on the pion mass (see, for example, Ref. [34]), so we choose a prior of 1100 ± 200 MeV for the nucleon mass on our ensemble with M π = 305 MeV. We put a wide prior of 300 ± 200 MeV centered at the ∆-like mass for the mass splitting M +1 − M +0 to accommodate for potential mixing of many physical states. For the same reason, we also impose a wide mass prior of 1600 ± 300 MeV for the negative-parity mass M −0 centered at the S-wave N π state. All prior choices are summarized in Table I. A priori, we have no knowledge of the sign or magnitude of the overlap factors. Consequently, all overlap factors are effectively unconstrained. Very wide priors of 0 ± 5 are chosen for all matrix elements, apart from g A = A +0,+0 , for which we choose a wide prior of 1.2±0.3 centered near the PDG [12] value of g A . As discussed below in Sec. VI, we know from other work with pseudoscalar mesons that Z A is close enough to unity not to influence the choice of prior.
For the g V analysis, we use the same two-point functional form as Eq. (5.1). However, for the three-point correlators we use where the notation is identical to that of Eq. (5.2). The V i,j are the unrenormalized vector matrix elements. The V +i,+j with i = j are omitted on the first line of Eq. (5.3), because they are forbidden by vector charge conservation, up to small discretization effects. The priors are also identical to the g A fits except for the bare vector charge, g V . Given that the renormalization constant is close to unity, we choose the g V prior to be 1.0 ± 0.3.

B. Fit Stability
The most important part of the nucleon matrix element fitting procedure is separating the nucleon observables of interest from the excited-state contributions. To demonstrate the lack of excited-state contamination, we examine the stability of the observables as choices in the TABLE II. Summary of the nominal fit range parameters. t is the source-sink separation time, and τ is the current insertion time. tmin and tmax are the minimum and maximum source-sink separation included in the nominal fits; ∆τmin is the minimum time after the current insertion time that we include in the three-point fits.

Correlator
Fit Parameter Nominal value Two-point tmin/a 5 tmax/a 13 Three-point ∆τmin/a 3 tmax/a 13 fit are varied. Specifically, we vary t min , ∆τ min , and t max where t min is the minimum source-sink separation time that we include in our two-point correlator fits, ∆τ min is the minimum source-sink separation time after the current insertion time, τ , that we include in our three-point correlator fits, and t max is the maximum source-sink separation time. The nominal parameters for the nominal fits are given in Table. II. We plot the stability of the extracted M N ( g V and g A ) as a function of t min and ∆τ min in Fig. 4 (Fig. 5). The x-axes are different choices of t min , and the y-axes are the corresponding observables. The four different choices of ∆τ min are also shown slightly displaced for each t min . The solid squares are the nominal fits with parameters given in Table II. As can be seen in Fig. 4, the extracted nucleon mass is stable as a function of t min /a and ∆τ min /a, which illustrates the lack of excited-state contamination in these posteriors. Similar behavior is seen for g V in Fig. 5. The only noticeable structure in the stability plots is for g A , where the observable is stable for t min /a ≥ 4. Note that as we increase t min /a or ∆τ min /a, fewer data are available to fit, and, consequently, the results become less precise. Thus, as in Ref. [20], we have demonstrated control over excited-state contamination when extracting matrix elements from staggered-baryon correlators.

VI. RESULTS
In this section, we present our Bayesian fitting results and our final renormalized values for the nucleon charges g V and g A . All fitting errors are estimated from 1000 bootstrap samples. We take correlations into account by using the same bootstrap samples for both g V and g A .

A. Nucleon Mass
In Fig. 6, we plot the extracted posterior fitted value for the nucleon mass from simultaneous fits of both smearings of the optimized two-point correlator and three-point correlator of a given current. We also plot the nucleon-optimized effective masses. This effectivemass data is identical to the t 0 /a = 6 data shown in Fig. 1. The green-shaded bands are the posterior estimates with the g A three-point correlators, while the yellow-shaded bands are with the g V three-point correlators. We obtain aM N = 0.707(6) from the g A fit, and aM N = 0.704(9) from the g V fit.
There are some notable features in our fits. First, the g V fit has larger posterior uncertainties than the g A fit.
Both fits include the same information from the twopoint correlators, so the difference must arise from the three-point correlators. As one can see in Fig. 2, the g V three-point correlators are less sensitive to the Wuppertal smearing than the g A correlators. On the other hand, the g V three-to-two-point correlator ratios show remarkably little curvature, even at the early times. This behavior implies that the vector three-point correlators become quickly saturated by the ground state, and therefore provide limited additional information about the overlap factors and masses than what is contained in the two-point correlators. The g A data does not share these features, and thus contains additional information about the twopoint posteriors. This explains why the g V fit has a less precise nucleon mass than the g A fit.
For these reasons, we quote the posterior nucleon mass from the g A fits as the nominal result, which has value aM N = 0.707 (6), M N = 1141(10) MeV, (6.1) where the error shown is statistical only. It is crucial to bear in mind that this result is for a lattice spacing of a = 0.1222(3) fm and pion mass of M π = 305 MeV [27]. For comparison, a fit including only the two-point correlators yields aM N = 0.704 (9), which is identical to the posterior of the fit with g V . In Ref. [20], we computed the nucleon mass at the same lattice spacing but with a physical pion mass, obtaining M N = 960(9) MeV. The difference between these two masses is ∆M N = 181(13) MeV, assuming uncorrelated statistical errors. Given that the pion mass difference between these two ensembles is about 170 MeV, ∆M N agrees within 1σ with the empirical observation that M N = 800 MeV + M π within a few per cent [34].

B. Nucleon gV and gA charges
In Figs. 7 and 8, we plot the optimized g V and g A threeto-two-point correlator ratios as a function of source-sink separation t. The raw data are identical to the righthand plots of Figs. 2 and 3. The posterior fit results are superimposed as gray bands. In the limits τ, t − τ → ∞, the data points are seen to converge to these posteriors. It should be emphasized, however, that the ratio data points are shown only for illustration: we perform direct fits to the optimized correlators, as discussed in Sec. V, in order to obtain results, namely g V = 1.03(2), (6.2) g A = 1.24 (5). (6.3) It should be mentioned that the g V and g A fits have some different features. First, the residual oscillations from the parity partner matrix element are noticeable in the g V fit. Second, the g V data turn out to be relatively insensitive to the Wuppertal smearing radius. Both of these features can be observed in Fig. 7. This highlights that there is less uncorrelated data available with which to extract g V as compared to g A . In contrast, we observe that the vector correlators in Fig. 7 contain less positive parity excited state contamination at early times than the axial-vector correlators in Fig. 8. As such, since the oscillations turn out to be easier to constrain and there is less contribution from the same parity excited states, we obtain a more precise estimate for g V than for g A .
As discussed in Sec. II B, the remnant chiral symmetry enforces Z A = Z V + O(am q ) 2 . Therefore, the ratio of bare charges is renormalized, and we obtain a value of (5). where the correlation between g A and g V is taken into account via bootstrapping. We can also obtain Z V by imposing current conversation on a pseudoscalar meson vector-current matrix element [22]. Then the renormalized charges are

VII. DISCUSSION AND CONCLUSIONS
We have presented two key results in this work. First, we have shown how to analytically relate the staggered nucleon-like matrix elements with non-trivial tastes to the physical nucleon matrix elements. This step is crucial for our on-going program of extracting high-precision nucleon results from staggered fermions. The general procedure, which can be applied to any staggered baryon matrix element, is outlined in Appendices A and B. Specifically, for the nucleon charges g V and g A , we summarize our key results for the zero-momentum isovector (axial) vector three-point correlators in Eqs. (2.3) and (2.4). These equations explicitly show the non-trivial normalizations needed to relate the nucleon-like matrix elements to their physical counterparts. Our successful computation of g V and g A shows that continued use of the 16 irrep of the staggered symmetry group GTS is feasible, which is convenient because the 16 contains a single nucleon-like taste in the spectrum [20].
This finding is encouraging, because the additional complexity of staggered baryons, compared with staggered mesons is probably the reason staggered-baryon matrix elements have not been explored until now. There are as many meson tastes (16) as bosonic irreps of GTS. As such, each staggered meson interpolating operator excites only a single taste of meson. In contrast, there are 64 = 4 3 different tastes of a staggered baryon, yet only three unique irreps of GTS, denoted 8, 8 , and 16 after their dimension. Consequently, there are not enough unique components of these irreps to accommodate all 64 tastes of baryons, and more than one taste of the same baryon can appear in each irrep's tower of states. Choosing an irrep with only one nucleon taste simplifies the correlator analysis and, as we have shown in this paper, allows for accurate and precise results for nucleon matrix elements.
The second key result of this work is demonstrating the practicality of staggered baryons by computing the isovector nucleon vector and axial-vector charges. For this purpose, we choose a single ensemble with a ≈ 0.12 fm, 2+1+1 flavors in the sea, and, when using identical sea and valence HISQ quarks, M π = 305 MeV. With approximately 7000 measurements and techniques designed to handle staggered correlators, we find fewpercent statistical uncertainty. Our final values for g V , g A , and g A /g V on this ensemble are 3) The conservation of the vector charge, g V = 1, is a nontrivial verification of our methodology. As discussed in Sec. V, we include two positive-parity states and one negative-parity state in our fit function. The number of matrix elements included in the fit grows quadratically as a function of the states included. With more precise data, we could constrain more matrix elements. Alternatively, we could also impose tighter priors on the transition matrix elements and overlap factors, for example with the empirical Bayes method [29]. This proof-of-concept study does not attempt a full calculation with all errors included, so we leave exploration of those options for future work.
The same ensemble has been used by both the CalLat [28] and PNDME [29] collaborations in their calculations of g V and g A . CalLat uses Möbius domainwall fermions for the valence quarks, while PNDME uses Wilson fermions with the clover action. CalLat defines Z V by demanding g V = Z V g V = 1 and uses the remnant chiral symmetry to set Z A = Z V . They then quote g V = 1.021(2) and g A = 1.21 (1). PNDME determines Z V and Z A independently via the regularizationindependent symmetric momentum-subtraction scheme, commonly known as RI-sMOM, and quote g V = 0.97(2), g A = 1.21(4), and g A /g V = 1.25 (2). Our result is consistent with both, despite the different choices of valencequark formulation. Other calculations in the literature, which have 2 + 1 flavors in the sea, have not been per-formed at values of a, M π , and physical volume close enough to ours to allow a straightforward comparison.
With an eye towards sub-percent determinations of the axial charge, it is instructive to compare how the precision on g A is influenced by each collaboration's data and methodology. Presumably influenced by the common ensemble, the three analyses share a few common aspects. First is the use of eight sources (with high-precision solutions of the Dirac equation) per gauge-field configuration, so the raw statistics are about the same. Second, the time range of the central fits for the two-point correlators turns out to be the same: t max + 1 − t min = 8. Third, all three collaborations simultaneously fit a correlator containing the matrix element with the two-point correlators. Last, PNDME and we use time ranges in the central fits of the three-point correlators, such that there are 21 data points in the fit.
In addition, each collaboration employs techniques to improve the signal. We have two smeared sinks and start with 4 × 4 matrix correlators, which is natural and necessary with our choice of staggered irrep. We apply the GEVP to the 4 × 4 matrix for each smearing radius to find the optimal source and sink operators for the nucleon. PNDME increases statistics via the truncatedsolver method with bias correction [35,36]. CalLat reduces noise with an a-independent number of steps of a gradient flow [37]. In the future, we could easily take advantage of the truncated-solver method, while the gradient flow would prevent us from using numerous technical results from the Fermilab Lattice and MILC collaborations, such as lattice-spacing and renormalization-factor determinations.
A more striking difference is CalLat's introduction of the currents into a propagator in a way inspired by the Feynman-Hellmann theorem [38]. A key feature of the technique is that instead of a three-point function, the matrix element lies within another two-point function. Thus, the CalLat method requires a fit to a single time variable instead of two; indeed the matrix element pops out of a fit to the ratio of the two two-point correlators.
In the end, the relative precision on g A is quoted as 1%, 3%, and 4% for CalLat [28], PNDME [29], and this work, respectively. One should bear in mind, however, the effective number of components per site, which are four for Wilson fermions, eight for staggered fermions (corresponding to the corners of the unit cube), and 4L 5 for domain-wall fermions (where L 5 is the extent of the fifth dimension; L 5 = 8 in Ref. [28]). Taking the number of components into account but ignoring algorithmic speedups from the code or specific features of each action, the cost for given precision is roughly the same. It would, therefore, be interesting to explore the truncated-solver and Feynman-Hellmann-inspired methods with staggered fermions.
This work sets the foundation needed to continue a program of precise nucleon form-factor calculations. Calculations of the vector and axial-vector form factors at nonzero momentum transfer are indeed underway on the same ensemble as used here. Further we, have started computing g V and g A on the same ensembles used in Ref. [20]. These ensembles have physical pion masses, and a range of lattice spacings to enable a continuum extrapolation.

Appendix A: Relating Staggered-QCD Matrix Elements to QCD Matrix Elements
Lattice gauge theory with staggered fermions can be thought of as an extension of QCD with four degenerate flavors, called tastes, for each quark. The associated taste symmetry allows for many more composite states which can have non-trivial taste structures. We call states that have non-trivial taste "baryon-like" states, to distinguish from the physical single-taste baryons. In this work, we focus on the nucleon and restrict ourselves to that case going forward. The nucleon-like states can be mapped onto the physical nucleon states through appropriate flavor-taste symmetry transformations. This allows the freedom to choose which nucleon-like state to study in order to extract observables. As highlighted in Ref. [20], the two-point correlator data constructed from nucleon-like states are easier to analyze than their physical counterparts due to the smaller multiplicity of tastes in the spectrum. However, one needs the mapping from the specific nucleon-like state to the physical state.
We use isospin-3 2 , GTS-16 nucleon-like interpolating operators to extract nucleon observables, since the spectrum contains only a single nucleon-like state. The relationship between the nucleon-like matrix elements and the physical nucleon matrix elements is, unfortunately, not at all transparent. In this and the following appen-dices, we will establish the relationship between the 16irrep nucleon-like matrix elements and the single-taste physical nucleon matrix elements.
Bailey [19] inferred the spectrum of staggered baryons by subducing nucleon-like representations of the full SU(8) F T flavor-taste symmetry of the continuum limit into GTS. We expand that work to matrix elements. Specifically, we will demonstrate how one can apply the generalized Wigner-Eckart theorem to SU(4) and relate the lattice nucleon-like matrix elements to the physical tasteless nucleon matrix elements through appropriate normalization factors, which are generalized Clebsch-Gordan coefficients. The procedure outlined here can be applied to any staggered baryon matrix elements in any SU(n f ) × GTS flavor-taste irrep.
Following the notation from Ref. [20], we first determine the continuum quantum numbers of the nucleonlike states that subduce into the 16 irrep of GTS. This step is needed for the generalized Wigner-Eckart theorem. We focus on the continuum symmetry group SU(2) S × SU(8) F T , where SU(2) S is the spin symmetry and SU(8) F T is the flavor (F ) and taste (T ) symmetry for two equal-mass flavors. This group breaks on a discrete lattice to the unbroken flavor symmetry subgroup SU(2) F and the "geometric timeslice group" (GTS) [18,19]. GTS can be decomposed into [20,39] where Q 8 is generated by the discrete taste transformations {Ξ 12 , Ξ 23 }, SW 3 by the cubic rotations {R 12 , R 23 }, and D 4 by the discrete taste and spatial inversion transformations {Ξ 123 , I S }. (These symbols are all defined in the appendix of Ref. [20].) The subgroup chain we work with is 1 where P = I S Ξ 4 becomes the usual parity operation in the continuum limit [18]. The factor SU(2) D4 on the second line arises from decomposing the SU(4) T taste symmetry onto a discrete lattice, which leads to the factor D 4 in Eq. (A1), combined with the U(1) D4 phase factor. Note that in Ref. [20] we omitted the U(1) D4 factor, but here we make it explicit. The other the groups are defined and explained in Ref. [20].

Using shift symmetries to relate staggered correlators
The goal is to assign continuum quantum numbers of SU(2) S ×SU(2) F ×SU(4) T to each nucleon-like state cre- ated by every component of the 16 irrep. We begin by investigating the continuum quantum numbers of the simplest nucleon-like states created by the 16 irrep. Afterwards, we can use the lattice symmetry transformations to obtain the remaining components.
We can form nonvanishing two-point correlation functions is by contracting any one of the 16 irrep components with the same component on a later timeslice. One can then apply lattice rotations and shifts 2 to show that these 16 two-point correlators are identical in the ensemble average.
The 16 irrep components split into two sets of 8 different components that reside on the eight corners of a cube (see the appendix of Ref. [20] for explicit constructions). The construction of nonvanishing three-point correlator data also depends on the current insertion. For the local vector and axial-vector currents we use in this work, the zero-momentum three-point correlators do not vanish if and only if the source and sink interpolators are identical. Correlators constructed from the same set of 8 components can be related to each other with the lattice shift symmetries.
To summarize, this means that the nonvanishing twopoint correlators satisfy where the superscript denotes the 16 irrep operators, M and N are equal to any one of the eight corners of the cube, and s = ±1 are the eigenvalues of the lattice rotation R 12 for M = N = 0. The notation is defined in detail in Ref. [20]. We are using local currents J = V, A in this work, so the nonvanishing three-point correlators satisfy Going forward, it is sufficient to study the correlator with component N = 0 located at the origin of the staggered unit cube, x B 16 ± 0 ( x, t). Then, owing to Eq. (A4), the other seven components follow immediately.
The quantum numbers of the nucleon-like states created by x B 16 ± 0 ( x, t) will be denoted as [ 3 2 , 3 2 ] F [16, ± 0] GTS . The first bracket gives the unbroken SU(2) F flavor quantum numbers, which here has total and z-component isospins 3 2 , and the second bracket denotes the 16 irrep with the eigenvalues of R 12 .

Quantum numbers of nucleon-like states
Next, we must find a convenient basis for the continuum nucleon-like states and then subduce them down to the [ 3 2 , 3 2 ] F [16, ± 0] GTS lattice states. From Eq. (A2), we want to track the quantum numbers of SU(2) S ×SU(2) F × SU(4) T and may ignore the passive phase U(1) D4 and parity P = +1 factors. From the group subduction presented in Ref. [19,20], the 16 irrep is subduced from the continuum spin-flavor-taste irrep via Here we have adopted a convention that labels non-SU(2) group irreps by their dimensions and subscript M (mixed), S (symmetric), or A (antisymmetric). The irreps of SU(2) are denoted with standard spin notation. The task of classifying a general irrep of SU(4) amounts to finding the maximal set of commuting operators and uniquely labeling the states by their eigenvalues; for a general SU(4) irrep, there are 6 eigenvalues to classify [41]. Because there are no degenerate irreps when decomposing any of the irreps in this work from SU(4) into SU(2) × SU(2), we can use the eigenvalues of the pair of SU(2) factors to identify SU(4) states. Therefore, only 4 of those 6 eigenvalues are necessary to completely characterize the states. As such, the 4 eigenvalues of each state can be uniquely identified with two pairs of the |L 2 , L z quantum numbers.

Matching the continuum and lattice nucleon-like states
Now that we have established an appropriate basis for the nucleon-like states, both on the lattice and in the continuum, we are ready match the two sets. In particular, we are interested in which linear combination of states from Eq. (A7) combine to subduce into the lattice states For the 16 irrep nucleonlike states, we have shown in Ref. [20] that j Q8 = 3 2 and j D4 = 1 2 . Consequently, we only need to determine m S , m Q8 , and m D4 .
We start with determining m D4 of SU(2) D4 from Eq. (A7). To do so, it is illuminating to study the decomposition where {I S } is the group generated by the lattice spatial inversion. As Eq. (A8) shows, I S receives contributions from three different factors: the taste factor SU(2) D4 , a phase factor e −iπ/2 = −i from U(1) D4 to match the eigenvalues of I S , and the continuum-limit parity P = I S Ξ 4 . For the spin-1 2 irreps of SU(2) D4 , which include the 16 irrep nucleons [20], the matrix representation of I S is the tensor product of those three factors where σ 3 is the third Pauli matrix. The representation in Eq. (A9) can be mapped onto the groups in Eq. (A8). The first factor arises from the 180 degrees rotation in the "x-y plane" of the spin-1 2 representation of SU(2) D4 , the second e −iπ/2 phase is from U(1) D4 , and the +1 from parity. As can be seen from Eq. (A9), for the spin-1 2 irrep of SU(2) D4 , the I S matrix admits ±1 eigenvalues which arise from the m D4 = ± 1 2 components of σ 3 . Since the nucleon is a positive-parity state with I S = 1, we assign m D4 = 1 2 to the [ 3 2 , 3 2 ] F [16, ± 0] GTS lattice states. We now consider the quantum numbers of m S and m Q8 . The 16 irrep components can be labeled by the irreps of W 3 = SW 3 × {1, I S }, where SW 3 is the cubic rotation group, as [18] 16 where E is the two-dimensional irrep of SW 3 , T 1 and T 2 are the different three-dimensional irreps of SW 3 , and the superscripts show the eigenvalues of I S . By applying lattice rotations to [ 3 2 , 3 2 ] F [16, ± 0] GTS , we can show they belong to the two-dimensional E + irrep of W 3 .
Subducing SU(2) SW3 ⊂ SU(2) Q8 × SU(2) S to the lattice angular momentum of SW 3 is a problem common to all fermion formulations [42]. We can write the irrep components of SU(2) SW3 that subduce into E as [42] 2, 0 of SU(2) S × SU(2) Q8 , and matching Q 8 factors, V can only transform as (0, 1) irrep of SU(2) Q8 × SU(2) D4 from Eq. (A17). We have just found that V is a triplet of SU(2) D4 , and so we need to determine its z-component quantum number. With positive parity, the three m D4 components of the (0, 0, 1) irrep from SU(2) S × SU(2) Q8 × SU(2) D4 subduce into the lattice currents γ 4 ⊗ γ 4 , γ 4 ⊗ ξ 4 ξ 5 , and γ 4 ⊗ ξ 5 . Each transforms trivially in RF. The first lattice current is local and the other two are non-local with multi-link connections between the quarks and antiquarks. The eigenvalues of I S are +1 for the local current, and −1 for the other two. As discussed in Eq. (A8), the matrix representation of I S in the continuum can be constructed from the tensor product of representations of SU(2) D4 , U (1) D4 , and P to give where the SU(2) D4 factor is in a spin triplet as discussed, U (1) D4 is a trivial factor to give the correct I S eigenvalues, and the parity is also trivial by construction. Consequently, to get the correct I S = 1 eigenvalue on the lattice, the local γ 4 ⊗γ 4 current must have zero z-component in the triplet irrep of SU(2) D4 in the continuum limit. This completes the subduction of V into V .
In the continuum, A is a spin-1 operator of SU(2) S . The A 1 irrep subduces from the spin-0 irrep of SU(2) SW3 and the E irrep subduces from the spin-2 irrep of SU(2) SW3 . With the rules for the addition of angular momentum, this requires A to be in the irrep (1,1) of SU(2) S × SU(2) Q8 with zero z-component spins in both SU(2) factors. Now, according to Eq. (A17), A can either be a spin-0 or 1 operator of SU(2) D4 . Recall that on the lattice, D 4 is generated by the transformations I S and Ξ 123 [20]. A is an eigenvector of both these symmetries with respective eigenvalues 1 and −1. Because SU(2) D4 subduces into the D 4 factor of the GTS group, these non-trivial eigenvalues mean that A cannot transform trivially under SU(2) D4 . As such, A can only belong to spin-1 irrep of SU(2) D4 . Further, it has zero z-component following the same argument in Eq. (A19).
In summary, we have determined the continuum quantum numbers of A and V, which subduce into the desired lattice current operators, A and V respectively. Using the same notation as in Eq. (A7), the continuum currents transform as The spin and flavor quantum numbers of the tensor operators are denoted by the superscripts, whereas the taste quantum numbers are given in the subscripts. n t = 4 is the number of tastes and √ n t = 2 is required to properly normalize tensor operators. The minus sign in front of the axial current is a convention that we follow according to Table I of Ref. [41].
As an aside, there is an easy way to obtain the continuum taste quantum numbers of an arbitrary quark bilinear without explicit group subduction. Table I of Ref. [41] outlines the SU(4) generators and their corresponding tensor operators. Once we adopt the Euclidean Dirac representation for the taste gamma matrices ξ 4 = σ 3 ⊗ I, ξ j = σ 2 ⊗ σ j (where σ j are the usual Pauli matrices), those generators give the components of the continuum taste matrices. For example, the local axial and vector currents we use have taste gamma matrices of  Table I of Ref. [41], and similarly T as Q 8 , we can recognize the tensor product S ⊗ T = σ 3 ⊗ σ 3 and σ 3 ⊗ I, indicating a spin-1 representation whenever a σ 3 appears in the tensor product. This yields the continuum taste quantum numbers of these states as (1, 0 x B 16 + 0 ( x, t) ("16+" with eigenvalue of +1) interpolating operators, as a function of source-sink separation time t and current insertion time τ . Both interpolators are used in Eq. (2.4) to compute gA. Each plot represents a different Wuppertal smearing at the sink, with parameters 0.2 (Gr2.0N30) and 0.6fm (Gr6.0N70) rms radii. The group theory requires that in the continuum and t, τ → ∞ limits that the ratio is equal to −3, which is shown as a dashed lines.
x B 16 − 0 ( x, t) and x B 16 + 0 ( x, t) interpolators in Fig. 9. In the limits τ, t − τ → ∞ and a → 0, the ratio should converge to the dashed lines at −3 as predicted by the above group theory. The small deviation is caused by a combination of excited state contamination, discretization effects, and taste-breaking effects. The same ratio for the vector current is consistent with one to high precision, as enforced by the lattice relation in Eq. (B6). Fig. 9 is therefore a non-trivial verification of our group theory understanding of staggered baryon matrix elements.
To relate the two remaining staggered matrix elements to their counterparts in QCD without tastes we observe that