Nucleon axial, scalar, and tensor charges using lattice QCD at the physical pion mass

We report on lattice QCD calculations of the nucleon isovector axial, scalar, and tensor charges. Our calculations are performed on two 2+1-flavor ensembles generated using a 2-HEX-smeared Wilson-clover action at the physical pion mass and lattice spacings $a\approx$ 0.116 and 0.093 fm. We use a wide range of source-sink separations - eight values ranging from roughly 0.4 to 1.4 fm on the coarse ensemble and three values from 0.9 to 1.5 fm on the fine ensemble - which allows us to perform an extensive study of excited-state effects using different analysis and fit strategies. To determine the renormalization factors, we use the nonperturbative Rome-Southampton approach and compare RI'-MOM and RI-SMOM intermediate schemes to estimate the systematic uncertainties. Our final results are computed in the MS-bar scheme at scale 2 GeV. The tensor and axial charges have uncertainties of roughly 4%, $g_T=0.972(41)$ and $g_A=1.265(49)$. The resulting scalar charge, $g_S=0.927(303)$, has a much larger uncertainty due to a stronger dependence on the choice of intermediate renormalization scheme and on the lattice spacing.


I. INTRODUCTION
Nucleon charges quantify the coupling of nucleons to quark-level interactions and play an important role in the analysis of the Standard Model and Beyond the Standard Model (BSM) physics. The isovector charges, g X , are associated with the β-decay of the neutron into a proton and are defined via the transition matrix elements, p(P, s)|ūΓ X d|n(P, s) = g Xūp (P, s)Γ X u n (P, s) (1) where the Dirac matrix Γ X is 1, γ µ γ 5 and σ µν for the scalar (S), the axial (A) and the tensor (T) operators, respectively. They are straightforward to calculate in lattice QCD since they receive only connected contributions arising from the coupling of the operator to the valence quarks, i.e. there are no contributions from disconnected diagrams. Lattice calculations of these charges were recently reviewed by FLAG [1], and we note some calculations of them in the last few years in Refs. [2][3][4][5][6][7][8][9][10][11][12][13][14][15]. The nucleon axial charge is experimentally well determined; the latest PDG value is g A = 1.2724 (23) [16]. In addition to its role in beta decay, the axial charge gives the intrinsic quark spin in the nucleon, and its deviation from unity is a sign of chiral symmetry breaking. Since the axial charge is so well measured, it is considered to be a benchmark quantity for lattice calculations, and it is essential for lattice QCD to reproduce its experimental value.
Unlike the axial charge, the nucleon scalar and tensor charges are difficult to directly measure in experiments. Thus, computations of those observables within lattice QCD will provide useful input for ongoing experimental searches for BSM physics. The generic BSM contributions to neutron beta decay were studied in Ref. [17], where it was shown that the leading effects are proportional to these two couplings; thus, calculations of g S and g T are required in order to find constraints on BSM physics from beta-decay experiments. The tensor charge is also equal to the isovector first moment of the proton's transversity parton distribution function (PDF), 1 δu−δd . Constraining the experimental data with lattice estimates of the tensor charge reduces the uncertainty of the transversity PDF significantly [18]. In experiment, there are multiple observables that could be used to constrain g T [19], and the overall precision will be greatly improved by the SoLID experiment at Jefferson Lab [20], providing a test of predictions from lattice QCD. In addition, the tensor charge controls the contribution of the quark electric dipole moments (EDM) to the neutron EDM, which is an important observable in the search for new sources of CP violation. The scalar charge is related via the conserved vector current relation to the contribution from the difference in u and d quark masses to the neutron-proton mass splitting in the absence of electromagnetism, g S = δM QCD N /δm q [21]. In this paper, we present a lattice QCD calculation of the isovector axial, scalar, and tensor charges of the nucleon using two ensembles at the physical pion mass with different lattice spacings. This paper is organized as follows. In Sec. II, we describe the parameters of the gauge ensembles analyzed, the lattice methodology, and the fits to the two-point functions used to extract the energy gaps to the first excited state on each ensemble. We discuss different analysis methods for estimating the three bare charges and eliminating the excited-state contaminations and present a procedure for combining the multiple results in Sec. III. The procedure we follow for determining the renormalization factors for the different observables using both RI -MOM and RI-SMOM schemes is described in Sec. IV. In Sec. V, we give the final estimates of the renormalized charges and discuss the continuum and infinite volume effects. Finally, we give our conclusions in Sec. VI. In Appendix A, we show analysis results of the bare charges using the many-state fit, which is an alternative model for excited-state contributions based on the contributions of noninteracting N π states with relative momentum ( p) 2 < ( p max ) 2 . In Appendix B, we list the bare charges determined on the two ensembles studied in this work, along with data used in previous publications [2,3].

A. Correlation functions
To determine the nucleon matrix elements in lattice QCD, we compute the nucleon two-point and three-point functions at zero momentum, Here, we place the source at timeslice 0, the sink at timeslice T , and insert the operator O X at the intermediate timeslice τ . The latter is the isovector current O X =qΓ X τ 3 2 q, where q is the quark doublet q = (u, d) T , and χ = abc (ũ T a Cγ 5 1+γ4 2d b )ũ c is a proton interpolating operator constructed using smeared quark fieldsq. We use Wuppertal smearing [22],q ∝ (1 + αH) N q, where H is the nearest-neighbor gauge-covariant hopping matrix constructed using the same smeared links used in the fermion action; the parameters are chosen to be α = 3.0 on both ensembles, N = 60 on the coarse ensemble, and N = 100 on the fine ensemble. The spin and parity projection matrices are defined 1 as Γ pol = 1 2 (1 + γ 4 )(1 − iγ 3 γ 5 ). In order to compute C X 3 , we use sequential propagators through the sink [23]. This has the advantage of allowing for any operator to be inserted at any time using a fixed set of quark propagators, but new backward propagators must be computed for each source-sink separation T . The three-point correlators have contributions from both connected and disconnected quark contractions, but we compute only the connected part since for the isovector flavor combination the disconnected contributions cancel out.

B. Simulation details
We perform our lattice QCD calculations using a tree-level Symanzik-improved gauge action and 2 + 1 flavors of tree-level improved Wilson-clover quarks, which couple to the gauge links via two levels of HEX smearing [24]. We use two ensembles at the physical pion mass: one with size 48 4 and lattice spacing a ≈ 0.116 fm which we call coarse, and another with 64 4 and a ≈ 0.093 fm which we call fine. Both volumes satisfy m π L ≈ 4. On the coarse ensemble, we perform measurements on 212 gauge configurations using source-sink separations T /a ∈ {3, 4, 5, 6, 7, 8, 10, 12} ranging roughly from 0.4 to 1.4 fm. In addition, we make use of all-mode-averaging (AMA) [25,26] to reduce the computational cost through inexpensive approximate quark propagators. For T /a ∈ {3, 4, 5}, we use approximate samples from 96 source positions per gauge configurations and high-precision samples from one source position for bias correction, and for T /a ∈ {6, 7, 8, 10, 12} we use double those numbers. On the fine ensemble, we perform the calculations on 442 gauge configurations using source-sink separations T /a ∈ {10, 13, 16} ranging roughly from 0.9 to 1.5 fm. AMA is applied with 64 sources with approximate propagators and one source for bias correction per gauge configuration. Table I   being the number of source positions used on each gauge configuration. The factor of 2 in Nmeas accounts for the use of forward-and backward-propagating states. Finally, N HP meas refers to the number of measurements made with high-precision.

C. Fitting two-point functions
On each ensemble, we perform two-state fits to the two-point correlation functions at zero momentum, where a i and E i denote the amplitudes and the energies of the ground-and first-excited states, respectively. For comparison, we also perform one-state fits with C 2 (t) = a 0 e −E0t only. The blue and red points in Fig. 1 show the dependence of aE 0 and a∆E 1 = a(E 1 − E 0 ) on the start time slice t min /a for the coarse (left) and fine (right) ensembles. The values for aE 0 were obtained using both the one-and two-state fits. The shaded blue and red bands indicate our preferred estimates of aE 0 and a∆E 1 , respectively. Those correspond to the two-state fits with t min /a = 4 for the coarse ensemble and t min /a = 5 for the fine ensemble. The quality of the resulting fits are shown in Fig. 2 by plotting the two-point function divided by its fitted ground-state contribution Table II gives a summary of the estimated fit parameters on both the coarse and fine ensembles.

III. ESTIMATION OF BARE CHARGES
Performing a spectral decomposition of the two-and three-point functions by inserting a complete set of states I = n |n n| into Eq. (2) and Eq. (3) leads to where we use the shorthand χ β = χ β (0) and similarly for O X . When the time separations τ and T − τ are large, excited states are exponentially suppressed and the ground-state denoted by n, n = 0 dominates. In this limit, the ratio of C X 3 (τ, T ) and C 2 (T ) yields the bare charge: where ∆E n ≡ E n − E 0 is the energy gap between nth excited state and ground state. Increasing T suppresses excitedstate contamination, but it also increases the noise; the signal-to-noise ratio is expected to decay asymptotically as e −(E− 3 2 mπ)T [27]. The ratio R X (τ, T ) produces at large T a plateau with "tails" at both ends caused by excited states. In practice, for each fixed T , we average over the central two or three points near τ = T /2, which allows for matrix elements to be computed with errors that decay asymptotically as e −∆E1T /2 .
Excited-state contamination is a source of significant systematic uncertainties in the calculation of nucleon structure observables. These contributions to different nucleon structure observables have been studied recently using baryon chiral perturbation theory (ChPT) [28][29][30][31]. Contamination from two-particle N π states in the plateau estimates of various nucleon charges, which becomes more pronounced in physical-point simulations, has been studied in Refs. [29,31]. It was found that this particular contamination leads to an overestimation at the 5-10% level for sourcesink separations of about 2 fm. This suggests that the source-sink separations of ∼ 1.5 fm reached in present-day calculations may not be sufficient to isolate the contribution of the ground-state matrix element with the desired accuracy. On the other hand, in Ref. [30] a model was used to study corrections to the LO ChPT result for the axial charge; it was found that high-momentum N π states with energies larger than about 1.5M N can be the cause for the underestimating of the axial charge observed in Lattice QCD calculations. These contributions, however, cannot be estimated in chiral perturbation theory. Refs. [28][29][30][31] find that multiple low-lying nucleon-pion states give important contributions to R X (τ, T ), which is in stark contrast to the commonly-used fit model based on a single excited state.
In the remainder of this section, we discuss the analysis methods we employ to study and suppress excited-state contributions to the axial, scalar, and tensor charges. We start with estimating the bare charges using the standard 'ratio method' in Sec. III A. In Sec. III B, we discuss the use of the summation method in addition to presenting a two-state fit model to the summations, which was inspired by the calculation in Ref. [10] that quotes a percent-level uncertainty for g A . Furthermore, we employ a two-state fit to the ratios R X (τ, T ), which is presented in Sec. III C. Finally, in Sec. III D, we explain the procedure we follow to combine the estimates from the different fit strategies and extract final values for the bare charges.

A. Ratio method
The ratio method is a simple approach that allows for excited-state effects to be clearly seen. We know that the excited-state contributions to R X (τ, T ) decay as e −∆E1T /2 which results eventually in a plateau when the source-sink separation is large enough. We observe on both the coarse and fine ensembles that the axial and scalar charges plateau as expected with increasing T /2, however, this does not happen in the case of the tensor charge, indicating that this method fails to reliably control excited states for g T .

B. Summation method
For studying the excited-state contributions, we use in addition to the aforementioned ratio method, the summation method [32,33]. The summation method allows improving the asymptotic behavior of excited-state contributions through summing ratios at each source-sink separation T . The summed ratios can be shown to be asymptotically linear in the source-sink separation, We choose τ 0 /a = 1. The matrix element can then be extracted from the slope of a linear fit to S X (T ) at several values of T . The leading excited-state contaminations decay as T e −∆E1T . For performing the fits of the summation method on the coarse ensemble, we vary the fit range by fixing the maximum source-sink separation included in the fit to T max /a = 12 and changing the minimal source-sink separation, T min /a. The obtained results for the three charges on the coarse ensemble are displayed as red squares in the upper right panels of Figs. 3 to 5 which demonstrate the dependences of g bare X on T min /a. Here, we see that the obtained g bare A shows a slight increase when increasing from the shortest T min and g bare T shows a somewhat larger decrease, whereas g bare S is flat. We eventually reach a plateau in all cases. The fit quality is measured by computing the p-value and the open symbols refer to fits with p-value < 0.02. The red squares in the lower right panels of Figs. 3 to 5 show the results for the summation method on the fine ensemble including all three available source-sink separations, which leads to one summation point at T min /a = 10.
The numerous source-sink separations used for calculations on the coarse ensemble allow us to perform the fit to the summations in Eq. (9) including contributions from the first excited state. This leads to two additional fit parameters c 1 and c 2 where ∆E 1 is estimated from two-state fit to the two-point correlation function. The fit function becomes In Fig. 6, we show the results of such a fit for all three charges. As before, we fix T max /a = 12 and vary T min . The results in Fig. 6 show that the fits for the different charges are stable although with relatively large statistical errors. The endcaps of the error bars refer to the resulting statistical uncertainties when fixing ∆E 1 in Eq. (10) to its central value whereas the vertical lines of the error bars result from taking the uncertainties in ∆E 1 into consideration when evaluating the fit in Eq. (10). We see that fixing ∆E 1 to its central value has little to no effect on the final results.
C. Two-state fit of the ratio In this section, we study including the contribution from a single excited state when fitting the ratio, R X (τ, T ). This is performed using the fit function Here, ∆E 1 is estimated from two-state fit to the two point function. We perform the stability analysis for this method by fitting to all points with τ ∈ [τ 0 , T − τ 0 ] and varying τ 0 . The circles with the outer statistical uncertainties in the plots of Fig. 7 show the resulting unrenormalized isovector charges as we vary τ 0 for the coarse (left column) and fine (right column) ensembles. The fit range includes source-sink separations satisfying T ≥ 2τ 0 ; this means that for the coarse ensemble, as τ 0 is increased the shorter source-sink separations (which have the most precise data) will be excluded from the fit. We notice that for g bare A , there is no significant dependence on τ 0 . The estimates for g bare S show a noisier signal on the fine ensemble. The signal for g bare T on the fine ensemble shows an upward trend in the central value for increased τ 0 while the statistical uncertainties are decreasing; this is normally not expected, whereas the signal on the coarse ensemble shows no to little dependence on τ 0 . The inner error bars of Fig. 7 show the uncertainties when ∆E 1 is fixed to its central value. The difference between the inner and outer statistical uncertainties for the axial charge shows that the uncertainty on the energy gap makes a large contribution to the uncertainties of the final results, particularly when including small time separations in the fit. This may be because the small time separations are more sensitive to the model parameters used to remove excited-state contributions. This observation applies also to the tensor charge but less to the scalar charge.

D. Combining different analyses
We have so far applied four methods for analyzing the excited-state contributions to the different observables on each ensemble, namely For each method, we have plotted the estimated charges as functions of a Euclidean time separation T /2, T min , or τ 0 , which we will generically call δt. For each bare charge, we want to choose a preferred δt for each method and then combine the results from all methods to obtain a final result. In order to reduce the number of case-by-case decisions, we have devised a procedure that we will follow to accomplish this. Our procedure is designed to fulfill the following requirements ordered in decreasing importance: • Fit ranges with poor fit quality are excluded, since that indicates the data are not compatible with the fit model and therefore the fit result is not trustworthy.
• Estimations should be taken from the asymptotic plateau regime, where there is no significant dependence on δt.
• Smaller statistical uncertainties are preferred.
• Larger time separations are preferred so that we reduce the residual excited-state contamination. T min /a T min /a In the following, we outline the first part of the procedure which aims to find a preferred δt from each analysis method.    On the fine ensemble, we do not have small values of δt for the ratio and summation methods. When δt 1 = δt min , this suggests that the plateau could start earlier than our available data. In this case, we choose to take δt 2 determined on the coarse ensemble for the same method and same charge, and use it (scaled to account for the different lattice spacings) as δt 2 on the fine ensemble.
The above procedure gives multiple estimates for each observable: at most one from each method. The obtained estimates of the charges for the coarse and fine ensembles are listed in Tab. III and Tab. IV, respectively. In those tables, we also outline for each case the obtained δt min which is the smallest available δt, δt 0 resulting from the first step in the above procedure and δt 1 from the second step. For cases where δt min = δt 1 , this indicates that there is no significant residual excited-state contamination. This is always the case for two-state fits, indicating that the data are compatible with the single-excited-state model. There are cases in the two tables where we have no remaining data after the second or the third step of the above procedure to define a δt f and therefore we leave those fields empty, as no reliable result could be obtained. We notice that we obtain similar δt 1 for the ratio and summation methods which indicates that it is appropriate to compare ratios at separation T with summation points at T min = T /2. In this case (and it can be seen in Fig. 3-5), the summation method provides more precise results than the ratios; this is in contrast to the usual comparison of ratio at separation T and summation at T min = T , which finds that the summation method has larger uncertainties. The values for the axial, scalar, and tensor charges in both tables show consistency within error bars between the different methods. The statistical uncertainties differ between the different fit strategies; in particular we obtain relatively large error bars for the scalar charge on both ensembles.
For obtaining a final estimate of the charges, we combine the different analysis methods by performing a weighted average to determine the central value. The statistical uncertainty is determined using bootstrap resampling. We test the compatibility of the central value with the set of analysis methods using a correlated χ 2 . If the reduced χ 2 is greater than one, then this indicates the different analysis methods are not in agreement, and the corresponding systematic uncertainty can be accounted for by scaling the statistical uncertainty by χ 2 /dof. We list our final estimates of the charges on both ensembles in Tab. V. In this table, the given uncertainties are obtained from bootstrap resampling and all the χ 2 values are acceptable. We obtain the largest χ 2 /dof = 1.04 for g bare S from the fine ensemble.

IV. NONPERTURBATIVE RENORMALIZATION
We determine renormalization factors for isovector axial, scalar, and tensor bilinears using the nonperturbative Rome-Southampton approach [34], in both RI -MOM [34,35] and RI-SMOM [36] schemes, and (for the scalar and tensor bilinears) convert and evolve to the MS scheme at scale 2 GeV using perturbation theory. Our primary data are the Landau-gauge quark propagator where ψ is a u or d field, and the Landau-gauge Green's functions for operator O, In our case, O is an isovector quark bilinear and there is only one Wick contraction, which is a connected diagram. We evaluate both of these objects using four-dimensional volume plane-wave sources, yielding an average over all translations in the lattice volume. From these, we construct our main objects, the amputated Green's functions, These renormalize as Λ R O = (Z O /Z ψ )Λ O . We will not determine Z ψ directly; instead, we will take ratios to determine Z O /Z V and compute Z V from pion three-point functions.

A. Conditions and matching
The RI -MOM scheme uses kinematics p = p, whereas RI-SMOM uses p 2 = (p ) 2 = q 2 , where q = p − p. In both cases the scale is defined as µ 2 = p 2 . Note that a comparison of RI-MOM and RI-SMOM renormalization was previously done using chiral fermions in Ref. [37].
For the vector current, the operator is V µ =ψγ µ ψ. In RI -MOM, the renormalization condition is where P µν = δ µν − p µ p ν /p 2 is the projector transverse to p, and for RI-SMOM the condition is 1 12q 2 Tr q µ Λ R Vµ (p , p) / q = 1.
Imposing the vector Ward identity, both of these imply that the quark field renormalization condition must be although we do not evaluate this explicitly. For the axial current, the operator is A µ =ψγ µ γ 5 ψ. In RI -MOM, the condition is and for RI-SMOM, it is 1 12q 2 Tr q µ Λ R Aµ (p , p)γ 5/ q = 1.
Each of these is related by a chiral rotation to the corresponding condition on the vector current. This implies that in the chiral limit, the renormalized axial current will satisfy the axial Ward identity, and therefore no matching to MS is needed. For the scalar bilinear, the operator is S =ψψ. In RI -MOM, the condition is and for RI-SMOM, it has the same form, For RI -MOM, the matching to MS is known to three loops [38,39], and for RI-SMOM it is known to two loops [40]. The anomalous dimension is obtained from the quark mass anomalous dimension via γ S = −γ m ; we use the four-loop MS result [41,42].
We write the tensor operator as T µν =ψ 1 2 [γ µ , γ ν ]ψ. In RI -MOM, Gracey [39] starts from the decomposition and then imposes the condition Σ (1),R T (p 2 ) = 1. Note that chiral symmetry breaking allows more terms to appear, but they won't contribute to any relevant trace. As Gracey notes, this term can be isolated via This can be rewritten to obtain the renormalization condition in a simple form: For RI-SMOM, the condition is For RI -MOM, the matching to MS is known to three loops [39], and for RI-SMOM it is known to two loops [40]. We use the four-loop MS anomalous dimension [43] 2 .

B. Vector current
Following e.g. Refs. [24,46], we determine Z V by computing the zero-momentum pion two-point function C 2 (t) and three-point function C 3 (t), where the latter has source-sink separation T = L t /2 and an operator insertion of the time component of the local vector current at source-operator separation t. The charge of the interpolating operator gives the renormalization condition for 0 < t 1 < T < t 2 < L t , where R(t) = C 3 (t)/C 2 (T ). We choose t 2 = t 1 + T ; the difference results in a large cancellation of correlated statistical uncertainties, so that precise results can be obtained with relatively low statistics; see Fig. 8. Results on the coarse ensemble are much noisier than on the fine one, although the statistical errors are still below 1%. We take the unweighted average across the plateau, excluding the first and last three points. This yields Z V = 0.9094(36) on the coarse ensemble and Z V = 0.94378(10) on the fine one.

C. Axial, scalar, and tensor bilinears
We use partially twisted boundary conditions, namely periodic in time for the valence quarks rather than the antiperiodic condition used for sea quarks. The plane-wave sources are given momenta p = 2π L (k, k, k, ±k), k = 2, 3, . . . , L 4a . By contracting them in different combinations, we get data for both RI -MOM kinematics, p −p = 0, and RI-SMOM kinematics, p − p = 2π L (0, 0, 0, ±2k). We used 54 gauge configurations from each ensemble. However, the modified boundary condition rendered one configuration on the coarse ensemble exceptional and the multigrid solver 2 Note that the sign of the three-loop term proportional to N 2 f disagrees between the proceedings of Baikov and Chetyrkin [43] and the first three-loop calculation, done by Gracey [44]. However, in an appendix of a later publication by Chetyrkin and Maier [45], the sign agrees with Gracey, and therefore we use that sign.  was unable to converge; therefore, we omitted this configuration and used only 53 on the coarse ensemble. In addition, on the coarse ensemble we also performed a cross-check using different kinematics, p, p ∈ { 2π L (k, k, 0, 0), 2π L (k, 0, k, 0)}, which ensure that in the RI-SMOM setup the components of p − p are not larger than those of p and p . Since the primary kinematics have p and p oriented along a four-dimensional diagonal and the alternative kinematics have them oriented along a two-dimensional diagonal, these setups will sometimes be referred to as 4d and 2d, respectively.
After perturbatively matching the RI -MOM or RI-SMOM data to the MS scheme and evolving to the scale 2 GeV, there will still be residual dependence on the nonperturbative scale µ 2 due to lattice artifacts and truncation of the perturbative series. To control these artifacts, we perform fits including terms polynomial in µ 2 and also, following Ref. [47], a pole term. Our fit function has the form A + Bµ 2 + Cµ 4 + D/µ 2 ; the constant term A serves as our estimate of the relevant ratio of renormalization factors Z O /Z V . We also consider fits without the pole term, i.e. with D = 0. We use two different fit ranges: 4 to 20 GeV 2 and 10 to 30 GeV 2 .
The main results on the two ensembles are shown in Fig. 9. The RI-SMOM data are generally very precise (more so than the RI -MOM data), which makes the fit quality very poor in many cases. If the covariance matrix from the RI -MOM data is used when fitting to the RI-SMOM data, then the fit qualities are good except for some of the fits without a pole term for the axial and tensor bilinears. For the RI -MOM data, the fit quality is good when using a pole term and also good for the scalar bilinear when omitting the pole term. Therefore, we elect to always include the pole term in our fits for Z A /Z V and Z T /Z V . For Z S /Z V we use fits both with and without the pole term, however the fit with a pole term to the RI -MOM data is very noisy and therefore we exclude it.
To account for the poor fit quality for some of the RI-SMOM fits, we scale the statistical uncertainty of the estimated ratio of renormalization factors by χ 2 /dof whenever this is greater than one. For each intermediate scheme, we take the unweighted average of all fit results as the central value, the maximum of the statistical uncertainties, and the root-mean-square deviation of the fit results as the systematic uncertainty. We combine results from both schemes in the same way to produce our final estimates, with the constraint that both schemes are given equal weight. These estimates are also shown in Fig. 9. For Z S /Z V there is a large discrepancy between the two intermediate schemes, which leads to a large systematic uncertainty. This discrepancy is smaller on the fine ensemble, suggesting that it is caused by lattice artifacts. Figure 10 shows the second set of kinematics on the coarse ensemble. These data do not reach as high in µ 2 ; therefore, we choose to fit to a single range of 4 to 15 GeV 2 . We use the same fit types as for the first set of kinematics, and the results (which can seen from the values of the curves at µ 2 = 0) are consistent with the final estimates from the first set of kinematics.
Our final estimates of the renormalization factors, after adding errors in quadrature, are given in Table VI. The uncertainty on Z S is more than 10% and we obtain percent-level uncertainties on Z A and Z T . In our previous publications using this lattice action [2,3,50], we used different values for these renormalization factors, which are listed in Table VII. These previous values were all obtained using an RI ( ) -MOM type scheme. Because of our large  For the fits that include a pole term, the fit curve without the pole term is also plotted, in the range 0 < µ 2 < 6 GeV 2 . The fits for Z S /Z V without a pole term are shown using desaturated colors.  uncertainty, Z S is in agreement with the previous value. The latter is also in agreement with our result from only the RI -MOM scheme. Our result for Z T is also consistent with the previous value. However, we find that Z A is 5-7% higher than the values that we previously used, a discrepancy of three standard deviations on the coarse lattice spacing and more than six on the fine one. The previous values would imply that Z A /Z V is within about one percent of unity for both lattice spacings, which is very difficult to reconcile with Fig. 9.

V. RENORMALIZED CHARGES
Multiplying the bare charges in Table V by the renormalization factors in Table VI and adding the uncertainties in quadrature, we obtain the renormalized charges on the two ensembles, shown in Table VIII. The final values should be obtained at the physical pion mass, in the continuum and in infinite volume. Since both ensembles have pion masses very close to the physical pion mass and have large volumes, we neglect these effects as their contribution to the overall uncertainty is relatively small. With two lattice spacings, we are unable to fully control the continuum limit; instead, we choose to account for discretization effects by taking the central value from the fine ensemble and quoting an uncertainty that covers the spread of uncertainties on both ensembles, i.e. δg X = max(δg f X , |g c X − g f X | + δg c X ), where g c X and g f X denote the charge computed on the coarse and fine ensembles, respectively. We obtain g T = 0.972 (41).
The overall uncertainties for the axial and tensor charges are roughly 4%. The scalar charge has a much larger uncertainty, due to the large uncertainty in the renormalization factor and the large difference in central values between the two ensembles. Results on these two ensembles can be compared with our earlier calculations using the same lattice action and heavier pion masses [2,3], reevaluating those earlier works based on the more extensive study of excited-state effects in Section III and using the renormalization factors from Section IV. For g A , the summation method with T min ≈ 0.7 fm was found to be acceptable; therefore, we reuse the summation-method results from Ref. [3], which had T min ≈ 0.9 fm. For g T , we found that the ratio method with the middle separation (T ≈ 1.2 fm), as used in Ref. [2] was inadequate; instead we will use the summation method. Finally, for g S the large statistical uncertainty means that the sourcesink separation used in Ref. [2] with the ratio method was larger than necessary, and here we will take the shortest separation (T ≈ 0.9 fm) rather than the middle one.
The comparison with our earlier results is shown in Fig. 11. In these plots, the ensembles used for a study of short time-extent effects are excluded and for two ensembles at m π ≈ 250 MeV of size 32 3 ×48 and 24 3 ×48, we have increased statistics. The data show no significant dependence on the pion mass, which justifies our neglect of this effect in the final values of the charges. If we assume that finite-volume effects scale as m 2 π e −mπL / √ m π L as for the axial charge in chiral perturbation theory at large m π L [51], then the finite-volume correction can be obtained by multiplying the difference between the two volumes at m π ≈ 250 MeV by 0.28 and 0.23 for the coarse and fine physical-pion-mass ensembles, respectively. One can see that this effect is also small compared with the final uncertainties.
This comparison provides the opportunity to revisit our earlier result for g A [3], which was unusually low. This was partly caused by the lower value of Z A , but the value obtained for m π = 149 MeV is still two standard deviations below the physical-point coarse ensemble. It appears that this is a statistical fluctuation, since the methodology has not been significantly changed.

VI. SUMMARY AND OUTLOOK
We have computed the nucleon isovector axial, scalar, and tensor charges using two 2+1-flavor ensembles with a 2-HEX-smeared Wilson-clover action. Both ensembles are at the physical pion mass and have lattice spacings of 0.116 and 0.093 fm. We have demonstrated control over excited-state contamination by using eight source-sink separations in the range from roughly 0.4 to 1.4 fm on the coarse ensemble and three source-sink separations in the range from 0.9 to 1.5 fm on the fine ensemble. The shorter source-sink separations are useful for the summation method but larger ones are needed for the ratio method. In addition, the choice of T is observable-dependent: if excited-state effects are drowned out by noise, then shorter separations are more useful. We have studied a range of different fitting strategies to extract the different charges of the nucleon from ratios of correlation functions, namely the ratio, two-state fit to the ratios, summation method, two-state fit to the summations (only on the coarse ensemble). We have studied the stability of the different analysis methods and designed a procedure for combining the multiple estimates obtained for each observable and giving an estimate of its final value. We have observed consistency between the different analysis methods, although within larger error bars for the scalar charge. We have determined the renormalization factors for the different observables using the nonperturbative Rome-Southampton approach and compared between the RI -MOM and RI-SMOM intermediate schemes to estimate the systematic uncertainties.
Our final results are given in Eqs. (27)(28)(29). The axial and tensor charges show overall uncertainties of roughly 4%. The obtained scalar charge, however, shows a much larger uncertainty, due to the large uncertainty in the renormalization factor and the large difference in the central values we observe between the the coarse and fine ensembles. In this study, since both ensembles have pion masses very close to the physical pion mass and have large volumes, we neglect the pion-mass dependence and finite volume effects. We have shown that this is justified when comparing our results to earlier calculations using the same lattice action and heavier pion masses. This calculation supersedes the earlier ones since it improves on them by working directly at the physical pion mass, using much higher statistics, and performing a more extensive study of excited-state effects.
In our calculation, we have found a large discrepancy for Z S between the two intermediate renormalization schemes; it would be therefore useful to verify whether this goes away at finer lattice spacings, and to compare against other approaches such as the Schrödinger functional [52] or position-space [53] methods.

ACKNOWLEDGMENTS
We thank the Budapest-Marseille-Wuppertal collaboration for making their configurations available to us. Calculations for this project were done using the Qlua software suite [54], and some of them made use of the QOPQDP adaptive multigrid solver [55,56]. Fixing to Landau gauge was done using the Fourier-accelerated conjugate gradient algorithm [57].
This research used resources on the supercomputers JUQUEEN [58], JURECA [59], and JUWELS at Jülich Supercomputing Centre (JSC) and Hazel Hen at the High Performance Computing Centre Stuttgart (HLRS). We In the two-state fit presented in Sec. II C, the obtained E 1 is much higher than the lowest expected N π or N ππ state. For this reason, in addition to using the two-state fit model, we also implement a many-state model for the excited-state contributions.
Inspired by [60], the many-state fit models the contributions from the first few N π noninteracting states with relative momentum ( p ) 2 < ( p max ) 2 . The noninteracting levels in a finite cubic volume with periodic boundary conditions are determined by where L denotes the spatial extent of the lattice and n is a three-vector of integers. We are interested in states with quantum numbers equal to that of a proton i.e. I(J P ) = 1/2(1/2 + ). However, the state of a pion and nucleon both at rest does not contribute since its parity is opposite that of the nucleon. The shift between free and interacting energy levels is small relative to the gap to the single nucleon state, as shown in [30]. This justifies the use of noninteracting finite-volume spectrum for the values of ∆E n = E n − m N .
The obvious difficulty in performing such a fit comes from the many fit parameters needed to parametrize the matrix elements and the overlap of the nucleon interpolating operator onto each of the N π states. In order to reduce the number of fit parameters involved in the many-state fit, we assume that the coefficient for ground-to-excited transitions is the same for all states, and the off-diagonal transition matrix elements between different excited states are small but that excited states in the two-point function in the denominator of the ratio are important. This yields a formula like the following 3 R X (τ, T ) = g bare X + b X n =0 | n| 2 ≤| nmax| 2 e −∆E n τ + e −∆E n (T −τ ) + c X n =0 | n| 2 ≤| nmax| 2 e −∆E n T .
(A2) where the three parameters are g bare X , b X , and c X . Also, the above n fulfills | n| 2 ≤ | n max | 2 and we exclude the state where the nucleon and pion are at rest, n = 0. We perform the many-state fit using four different values of n 2 max , n 2 max ∈ {5, 10, 20, 50}. Figure 12 shows a summary of the estimated unrenormalized isovector axial, scalar, and tensor charges as functions of τ 0 using this approach applied on both the coarse (left column) and fine (right column) ensembles. We notice that the estimated charges at short τ 0 depend significantly on n 2 max and that increasing n 2 max results in decreasing the statistical uncertainties of the estimated charges. In addition, Fig 12 shows that the obtained charges for different n 2 max values tend to be consistent at the largest τ 0 . When comparing to our final estimates of the bare charges in Tab. V (gray bands), we notice that the estimates from the many-state fit approach are consistent within error bars and that the many-state fit leads to larger statistical uncertainties for g bare A and g bare S compared to other analysis methods. The strong dependence on n 2 max at short τ 0 suggests that this method may be relatively unreliable and that a more sophisticated model such as the one in Ref. [30] is needed to extend into the resonance regime.  Appendix B: Table of bare charges   Table IX contains the bare charges used (after renormalizing) in Fig. 11, based on this work (physical-point ensembles), data from previous publications [2,3], and increased statistics at m π ≈ 250 MeV.