Nucleon axial coupling from Lattice QCD

We present state-of-the-art results from a lattice QCD calculation of the nucleon axial coupling, $g_A$, using M\"obius Domain-Wall fermions solved on the dynamical $N_f = 2 + 1 + 1$ HISQ ensembles after they are smeared using the gradient-flow algorithm. Relevant three-point correlation functions are calculated using a method inspired by the Feynman-Hellmann theorem, and demonstrate significant improvement in signal for fixed stochastic samples. The calculation is performed at five pion masses of $m_\pi\sim \{400, 350, 310, 220, 130\}$~MeV, three lattice spacings of $a\sim\{0.15, 0.12, 0.09\}$~fm, and we do a dedicated volume study with $m_\pi L\sim\{3.22, 4.29, 5.36\}$. Control over all relevant sources of systematic uncertainty are demonstrated and quantified. We achieve a preliminary value of $g_A = 1.285(17)$, with a relative uncertainty of 1.33\%.


Introduction
The nucleon axial coupling, g A , parameterizes the interaction strength of the nucleon with the weak axial current in the Standard Model, and is one of the fundamental properties which governs nuclear physics. For example, the β-decay rate of a neutron to a proton is governed by the strength of g A , and confonting the sub-percent determination of g A obtained from experiments [1] with precise determination of this coupling from Lattice QCD serves as a test for the weak structure of the Standard Model. Interpretation of potential signals from observation of neutrinoless double-beta decay in nuclear targets also depends on the strength of the axial coupling to the nucleon, and in general, Speaker, e-mail: chiachang@lbl.gov multi-nucleon systems. Demonstration of control over all sources of systematic uncertainty related to single nucleon g A is needed before calculating precision multi-body corrections to g A . Precision tests of neutrino properties also depend on g A as an input, and in general F A (Q 2 ), the nucleon axial form factor, which parameterizes the nucleon axial coupling as a function of momentum transfer. In particular, future experimental results from DUNE and T2HK will be very interesting in light of recent results from T2K which exclude no CP violation in the leptonic sector with 90% confidence level [2], lending evidence to scenarios involving baryogenesis through leptogenesis. In the regime of quasielastic scattering, only coupling to a single nucleon is required, however determination of the shape of the form factor from lattice QCD at high precision is desired when interpreting next generation experimental results. Robust determination of the form factor at zero momentum transfer provides a cross check with experimental measurements.
The calculation of g A from lattice QCD has, however, proven to be notoriously challenging and results have historically been in tension with the value determined from experiment. Insight of previous works has let us identify two main challenges that needed to be overcome in order to obtain a precise and accurate value of g A : 1) the nucleon suffers from an exponentially problematic signal-tonoise problem in comparison to pions, leading to large statistical uncertainties and poor estimation of systematics at moderate values of source-sink separation, and 2) the excited state nucleon structure is dense compared to pion quantities and objectively demonstrating control over excited state contributions is required to remove this leading systematic uncertainty.
In this work we present a calculation of the nucleon axial coupling, g A = 1.285 (17), based on a new computational method inspired by the Feynman-Hellmann theorem which alleviates the above concerns. This result is commensurate with the value determined from experiment, g PDG A = 1.2723 (23), to a precision of 1.33%. We summarize our computational method in Sec. 2, the lattice setup in Sec. 3, the construction of correlation functions in Sec. 4, the analysis of the correlation functions in Sec. 5, our renormalization procedure in Sec. 6, the physical point extrapolation in Sec. 7, sources of systematic error in Sec. 8, and conclude in Sec. 9.

A Feyman-Hellmann inspired method
For previous calculations of g A [3][4][5][6][7][8], and in general, for processes in which an external current mediates the interaction (this excludes, for example, neutral meson mixing [9,10]), the three-point correlation function is constructed using the fixed-sink method. The fixed-point correlation function is constructed from a sequential propagator, in which a regular one-to-all propagator is inverted again at the sink, at fixed time location. As a result, multiple fixed-sink correlation functions need to be calculated in order to reveal the source-sink separation dependence of the three-point correlation function, adding linearly to the computational cost per sink location.
For this work we use a novel method which, using the same computational cost as a single fixedsink time, gives access to the complete sink time dependence. Having the full time dependence of the three-point correlation function availed us to the use of both exponentially more precise data at small source-sink separation times, and the full functional dependence needed to control excited state effects, directly addressing the two main challenges encountered while calculating nucleon related quantities. Furthermore, while the fixed-sink method allows for reuse of computational resources in the calculation of different currents within the same state, our method instead will allow us to reuse the resources for the calculation of the same current within different states, such as g A for multi-nucleon states or the axial charge of hyperons.

Feynman-Hellmann correlation function
Here we summarize the Feynman-Hellmann theorem inspired method used to perform our calculation, which is discussed in detail in Ref. [11]. The Feynman-Hellmann theorem relates matrix elements to variations in the spectrum [12][13][14][15], where the Hamiltonian is given by H = H 0 + λH λ . We apply the Feynman-Hellmann theorem directly on to the analogous lattice energy, the effective mass, where C(t) is the two-point correlation function. In the literature, there exists similar methods which adopt the name of the Feynman-Hellmann method [16,17] or a background field method [18]. The Feynman-Hellmann method that we employ for this calculation differs from the cited methods in that similar methods numerically calculate the derivative of the two-point correlation function of Eq. 3, while we algebraically resolve this derivative. Numerical evaluation of the two-point correlation function requires the introduction of a background field at different values of the coupling λ, and modifies the derivative in Eq. 3 as a finite difference at additional computational cost for each value of the coupling. The algebraic derivative of the two-point correlation function may be derived from the QCD generating functional with an external source J(t) ≡ d 3 x j(t, x) with coupling λ, and nucleon interpolating operators N such that The first term contains the vacuum expectation value of the current and vanishes for the axial current. The second term is a time-ordered three-point function integrated over the current insertion time t and yields four contributions: 1) current insertion between the source and sink which contains the signal of interest, 2) current insertion outside the source and sink, and 3) and 4) contact terms when the current is on the same time slice as the source or sink. Note that a direct application of the Feynman-Hellmann theorem on the effective mass yields contribution 1), which reproduces the summation method [4] if the summation method is applied to the fixed-sink correlator for all sink locations.

Spectral decomposition
The spectral decomposition of the derivative of the two-point correlator given in Eq. (4), including the undesired contributions is derived in detail in Sec. II A of Ref. [11]. We reproduce here the relevant conclusions of that paper. The expression for the spectral decomposition of the three-point correlation function generated from the Feynman-Hellmann method (Eq. (4)), may be simplified to the following form, As a result of the sum over the current insertion, the t in the above expression is the source-sink time separation. The first term in the square brackets that receives the (t−1) temporal enhancement contains g nm = n|J|m / √ 4E n E m , where n = m = 0 is the bare nucleon axial coupling, and is the dominant signal in the entire correlator. The overlap factor has a factor of energy absorbed into the definition as well, z n = Ω|N|n / √ 2E n . The rest of the sum over n for the (t − 1) term, and the sum over n m contains the excited state contamination when the current insertion is between the source and sink. We see that the additional artifacts introduced by the Feynman-Hellmann method (when compared to the fixed-sink method), may be cleanly parameterized by one additional tower of parameters d n . This is illuminated by introducing the substitutions performed to obtain Eq. (5), The first two terms of Eq. (6) are from contact terms, and give rise to an analogous overlap factor Z J:n , where the interpolating operator, which still overlaps with the nucleon, consists of three quarks and a quark-antiquark pair. The last two contributions come from when the current is inserted after the nucleon annihilation operator, resulting in baryon-meson transition matrix elements, Z jn , and overlap factors on mesonic states, J j . In principle, these extra contributions are difficult to determine. However, the entirety of d n is time independent, and can be simply reparameterized. In addition, note that, As a result, inspecting the correlator at t = 1 yields an extremely good order of magnitude estimate for the values of d n .

Domain-wall fermions on gradient-flowed HISQ
The calculation is performed on the MILC collaboration's 2+1+1 flavor Highly-Improved Staggered Quark (HISQ) ensemble [19,20]. The MILC ensembles are the only set of publicly available gauge configurations that allow for control over the continuum limit, infinite volume, and physical pion mass extrapolation. The MILC configurations used in this project span three lattice spacings, a ∼ {0.15, 0.12, 0.09} fm allowing for control over a 2 discretization corrections. Three pion masses are generated by the MILC collaboration, m π ∼ {130, 220, 310} MeV, while additional heavier pion masses of m π ∼ {350, 400} MeV are generated by the CalLat collaboration to control interpolation to the physical pion mass. A volume study is performed at a ∼ 0.12 fm, m π ∼ 220 MeV, where all input parameters are held fixed at three spatial volumes m π L ∼ {3.2, 4.3, 5.4}. Formally the HISQ action has leading discretization errors starting at O(α S a 2 , a 4 ), however improved link-smearing greatly suppresses taste-changing interactions, leading to numerically smaller coefficients. The gluons are described by the tadpole-improved [21], one-loop Symanzik gauge action [22] with leading discretization errors starting at O(α 2 S a 2 , a 4 ). The valence quarks are calculated using the Möbius domain-wall action [23][24][25]. The mixed action set up is chosen because of the good chiral properties of the domain-wall action, which significantly suppresses chiral symmetry breaking and discretization effects from using a staggered sea. The Möbius kernel allows for moderate values of the fifth dimensional extent L 5 to be chosen while still keeping the residual chiral symmetry breaking to less than 10% of the input light quark mass. Due to the very small residual chiral symmetry breaking, numerically, the domain-wall action has  Table 1: HISQ gauge configurations and valence sector parameters. The HISQ ensembles used in this work (e.g. a15m310 stands for the ensemble with a ∼ 0.15 fm and m π ∼ 310 MeV) along with the number of configurations N cfg , lattice volume, approximate lattice spacing a, approximate HISQ taste-5 pion mass, and approximate value of m π,5 L. The taste-5 pion is the one protected by γ 5 -symmetry, and does not receive additive renormalizations. Values are obtained from Table I of Ref. [19] with increased number of configurations. Mobius domain-wall propagators are generated at a number of sources per configuration N src , with the fifth dimensional extent L 5 /a, such that m res is minimized at aM 5 , with the Mobius kernel defined by b 5 and c 5 , and valence light-quark masses am val. l . We also list the width σ smr and iteration count N smr of the SHELL_SOURCE and the GAUGE_INV_GAUSSIAN smearing algorithm in Chroma.
discretization errors starting at O(a 2 , α S a 2 ). The valence pseudoscalar masses are tuned to within 2% of the HISQ masses, yielding a unitary theory in the continuum limit.
We add an additional layer of gauge field smearing with the 4-dimensional Wilson flow procedure with a dimensionless flow time of t g f = 1.0 [26,27]. The new scale introduced by the gradient flow is approximately l g f ∼ 8t g f a. We observe that the gradient flowed ensembles are more continuumlike, and are stochastically less noisy, because the new scale decouples UV contamination from the correlator. In addition, we have evidence supporting the fact that our action is not over smeared. For the axial coupling, we observe that the ratio of g A /g V is flow time independent, as shown in Fig. 1, demonstrating that the continuum extrapolation is flow time independent. Additionally, we performed a complete flow time study on F K /F π calculated at t g f = {0.2, 0.4, 0, 6, 0.8, 1.0}. The ratio of pseudoscalar decay constants after chiral-continuum extrapolation is demonstrably flow time independent, and consistent with the FLAG average [28]. Table 1 lists parameters pertinent to the gauge configurations, and valence sector. A complete discussion of flow time dependence, and an in-depth review for this mixed action set up may be found at Ref. [29].

Lattice correlation functions
Standard nucleon interpolating operators are used to when constructing the correlation functions [30,31],N where u(x) and d(x) are up-and down-quark field operators at x, while the Γ and P are spin projection operators. We double statistics by constructing the spin averaged correlation functions. The local axial and vector operator is used to construct the derivative correlator in Eq. (4), from which we obtain the bare values of g A and g V , The Pauli τ 3 matrix acts on isospin space, and projects the isovector combination, yielding g A relevant to nuclear beta-decay. A detailed discussion on the construction of the lattice correlation function, and their associated Wick contractions are presented in Sec. II B of Ref. [11].
Multiple sources per configuration are created in order to increase our statistical sampling, as reported in the first column under valence parameters of Table 1. For each configuration, a random origin is chosen, (x 0 , y 0 , z 0 , t 0 ), along with its antipode in space, (x 0 , y 0 , z 0 , t 0 ) + L/2(1, 1, 1, 0), on which evenly spaced time locations are chosen for the multiple sources. Once generated, the correlators are shifted to t 0 = 0 and averaged over different sources. We observe a reduction in statistical uncertainty inversely proportional to √ N src , indicating no correlation between different sources. We further double the statistics by combining the time-reversed negative-parity nucleon correlator with the forward-propagating positive-parity nucleon correlator.
Correlators are constructed using a gauge invariant Gaussian smeared source [32] to increase the overlap with the ground state, and increase statistics by constructing and analyzing correlators with smeared and point sinks. The last two columns of Table 1 under valence parameters provide the smearing width and iteration count of the procedure.
Finally, the dataset is free of autocorrelations as a function of Monte Carlo time as demonstrated by Fig. 2. The error of the mean is stable as a function of bin size at all time slices for the a09m310 and a15m310 ensembles. The independence of bin size is representative of all the data. As a result, we do not bin any of our data.

Correlation function analysis
In the following section, we discuss our correlator fit procedure and provide evidence for full control over contamination of the ground state from contributions from excited states. A short discussion on the analysis of the pion mass and decay constants is also provided because a correlated analysis is performed in a later step when the chiral-continuum limit is taken. All correlator fits are performed using the Python library lsqfit [33].

Bare g A and g V fits
The bare axial and vector couplings of each ensemble are extracted from a simultaneous and correlated fit to the two point correlation function and the effective mass derivatives given by Eq. (3) for the axial and vector current insertions, each with smeared-and point-sink operators. The six correlators are fit together because the overlap factors and the spectrum are shared. We first perform a two-state constrained fit in order to survey a wide range of source-sink separation times for all six correlators.
The central values of the resulting posterior distributions are then used as initial guesses for the final unconstrained correlator fit. The initial step of constrained curve fitting does not affect the final result, but simply provides an efficient way to explore large parameter spaces. The initial guesses motivated by the mean of the posterior distributions also only decreases the number of steps required for convergence. The final unconstrained fit is performed using a two-state fit for all six correlators. The preferred correlator fits are the ones which pass a set of stringent requirements, resulting in an objective selection process. We first assess the quality of all fits by considering only results with p-values greater than 0.05 in order to discriminate against fits of poor quality. The p-value in practice does an adequate job eliminating fits which are in great tension with data. This is usually seen when the fit is performed too aggressively at very small source-sink separations, where "very small" is of course, a relative statement. A good p-value, defined to be greater than 0.05, however, is an inadequate condition for a good fit. To demonstrate control over excited state contamination, we vary the fit regions over different time separations, and demand that the preferred fit lie in the region of stability. The reasoning is that if the data is adequately described by our fit ansatz, which in the case of this calculation, we limit to the inclusion of only 2 states. This is because an unconstrained fit to more than two states leads to numerical instability. For the fit ansatz for g A and g V , we directly implement the ratio given by Eq. (3) by constructing the difference in the ratio of Eq. (5) and the standard two-point spectral decomposition without performing any additional manipulation to the fit ansatz.
In Fig. 3, we show example fits for the a12m220 and a09m220 ensembles, which pass this test. The figure demonstrates that the preferred fit (in black) lies in the region of stability while the fit 1.25 1.23 1.25 regions are varied over the two point C(t), g A correlator G A (t), and g V correlator G V (t). The horizontal bands highlights the preferred fit, while the green (blue) points to the right lie consistently within the preferred fit. We also observe that by going to larger values of G A (t) t min , the uncertainty on the ratio g A /g V increases. This is because dropping more precise data leads to a more uncertain result. In general the uncertainty of the ratio is more sensitive to changes in the fit region to the g A correlator because the uncertainty ong V is approximately an order of magnitude smaller thang A , while the twopoint correlator only indirectly affects the ratio. However, dropping too much data from the two-point correlators (not shown) will eventually lead to catastrophic numerical instabilities because the overlap factors can only be disentangled from information provided by the two-point correlators. On the other hand, stability plots also reveal fits that are too aggressive, as shown by the fit at G A (t) t min = 3 for the a12m220 ensemble. We observe that the more aggressive fit is more than one standard deviation lower than the preferred fit, lending evidence to residual excited state artifacts contaminating in the ground state parameter. The different preferred fits also all pass the minimum p-value (indicated by the dashed red line), as demonstrated by the bottom portion of Fig. 3. Analysis of the a15m310 ensemble was performed in previous work [11] under the Bayesian framework with up to 8 states and is consistent with the result presented here.
We propagate the resulting uncertainty of the correlator fits though bootstrap resampling. The simultaneous fit already accounts for correlations between the bare couplings, however we also account for correlations between the couplings and the pion mass and pion decay constants, which become relevant when the extrapolation to the physical point is performed. We take the central value of our preferred fit as the initial guess for the bootstrap routine, and sample the final distribution with 5000 bootstrap draws. We save a master list of random integer values for every ensemble, in which we coordinate all bootstrap routines for all calculations performed on these gradient-flowed HISQ ensembles. The bootstrap list and complete bootstrapped correlator results are saved in the collaboration Post-greSQL database hosted at NERSC. From the bootstrap distribution, we plot the histogram of the ratio of couplings, and the resulting fit on top of the correlator data. Fig. 4 provides an example of this check for the a12m220 ensemble. We observe that the correlator fit is numerically stable under bootstrap, showing no outliers and is clearly Gaussian distributed as expected. Plotting the fit result on top of the correlator data also reveals that the fit ansatz describes the data very well, capturing the curvature at early source-sink separation times, indicating that excited states are well described by a  2 state fit. Here we note that the nucleon suffers from a much larger signal-to-noise problem when compared to mesons, including heavy-light systems. This is reflected in the observation of correlated fluctuations which are observed at moderate values of source-sink separation time. On the other hand, when one studies heavy-light systems, while the signal damps away faster than for nucleon systems, the noise is also damped by a mass dominated by the bottom quark scale, which yields a noisy, but well-behaved uncertainty in the correlator. For the example shown in Fig. 4, the correlated fluctuation can be seen in the bottom right panel. We observe that around 1 fm, the g V correlator exhibits a downward fluctuation of approximately one standard deviation. Therefore, aside from generating exponentially more data which will with brute force push the fluctuation to larger values of t, a controlled correlator fit to nucleon quantities can only be performed if one has access to small values of t. This is regularly performed for the two-point correlator, but will be more expensive if a fixed-sink method is used to construct three-point correlators.

Pseudoscalar mass and decay constant fits
The chiral-continuum extrapolation can be reparameterized to depend on the dimensionless quantity circumventing the necessity of performing a scale-setting analysis. Calculation of the pion mass, m π , and the pion decay constant, F π , are performed on the same lattice actions, gauge configurations, and sources as the main analysis of this paper.
A Bayesian constrained fit with a four-state fit ansatz is performed on the pion two-point correlation function in order to extract m π and its overlap factors; in Ref. [29], we show that the oscillating state present in the Domain-Wall action are highly suppressed when the gauge fields are smeared with gradient-flow, and therefore are neglected in this analysis. A simultaneous fit to both the point-sink and smeared-sink correlators is performed, and statistics is further doubled by "folding" the meson correlators. Similarly, a Bayesian constrained fit to a constant is performed to extract m res from the m res correlator, as defined by Eq. (5) of Ref. [29]. The 5D Ward Identity is used to obtain F π from m π and m res as given in Eq. (6) of Ref. [29].
The ground state priors for the pion mass and overlap factors are determined by the long-time limit of the effective mass and scaled correlators. The scaled correlator is the raw correlator with the leading exponential scaled away, and in the long time limit is proportional to the overlap factor. The ground state prior widths are set to 10% of the prior central value, approximately two orders of magnitude larger than the width of the posterior distribution, thus leaving the ground state effectively unconstrained. The prior for the excited state energy splitting is log-normal and is set by approximately the mass splitting between the three-pion and one-pion state, with a width encompassing the two-pion to one-pion splitting within one standard deviation. The prior for m res is set by plotting the m res correlator, with a width set to 10% of the central value, and is approximately one order of magnitude larger than its posterior distribution.
The fit regions for the pion two-point correlation functions and the m res correlators are chosen in the region of stability under varying choices of source-sink separation time to ensure full control over excited state systematic uncertainty.
Uncertainties are propagated by bootstrap resampling. For the preferred fits, 5000 bootstrapped configurations are analyzed, allowing for a correlated analysis with the bootstrapped samples of g A /g V . For each bootstrap, the prior central values for all parameters are set to a value randomly drawn from their corresponding initial prior distributions. The bootstrap histogram for π for the a12m220 and a09m220 ensembles is provided in Fig. 5 as an example. We observe that the distributions are Gaussian without outliers, indicating that the constrained fits are numerically stable.

Non-perturbative renormalization
The discretization of the Dirac action leads to differences between the local current, as defined by Eq. (9,10), and the conserved currents. We correct for these differences using the non-perturbative Rome-Southampton renormalization procedure [34], with non-exceptional kinematics [35,36], and implement momentum source quark propagators [37], achieving high statistical accuracy.
Defining the incoming and outgoing momentum-space propagators as G(p) andḠ(p), and the quark bilinear matrix element V Γ (p 2 , p 1 ), the amputated vertex function follows, We match the quark bilinear matrix element to its tree-level value, indicated by the subscript '0', at scale µ, where the quark bilinear operators J Γ relevant to g A are defined by Eq. (9,10), and Z Γ is the renormalization coefficient. In particular, we enforce the SMOM condition such that µ 2 = p 2 2 = p 2 1 = (p 2 − p 1 ) 2 . The matching to the tree-level value is performed in color and Dirac spaces, as a result, taking the trace yields the renormalization condition, where F (s) Γ is the trace of the tree-level matrix element, and Λ (s) Γ is the amputated vertex function projected to the s = {γ µ , / q} intermediate schemes with q = p 2 − p 1 , We are only interested in the cases where Γ = {γ µ , γ µ γ 5 }, in which case the intermediate schemes and choice of scale changes only the definition of the wave function renormalization factor Z q because the vector and axial currents are protected by Ward identities. We circumvent the calculation of Z q by taking advantage of the fact that the vector charge is by definition, normalized and as a result, g A can be renormalized by computing the following ratio Therefore, in practice we compute the ratio of renormalization coefficients from which should be intermediate-scheme and scale independent. Fig. 6 (Left) shows a calculation of the renormalization coefficients for the a ∼ {0.15, 0.12, 0.09} ensembles with m π ∼ 310 MeV. We observe that Z A and Z V from the γ µ and / q schemes agree  well within uncertainty above 2 GeV. Additionally, the ratio of coefficients enjoys a large Rome-Southampton window. Evidence of IR contamination vanishes for momenta larger than 3.0 GeV, while there is no evidence of an onset of growing discretization uncertainty for momenta through 4.5 GeV. Ratios of renormalization coefficients for all ensembles are commensurate with unity to one part in 10,000. Therefore, we do not quote the uncertainty related to these coefficients. Additionally, we study Z A /Z V as a function of flow time and find the different to be insignificant.
Because quark-bilinear matrix elements are gauge variant, we perform these calculations under the Landau gauge. Landau gauge fixing however, is incomplete because different Gribov copies yields different local minima. We quantify this uncertainty by performing random global gauge transformations to the configurations before recalculating the renormalization coefficients. Fig. 6 (Right) show this study on the a15m310 ensemble. We observe that the Gribov uncertainty is sub-dominant when compared to the statistical uncertainty.

Chiral-continuum extrapolation
The renormalized lattice results for the axial coupling is extrapolated to the physical point using SU(2) NNLO Heavy Baryon chiral perturbation theory (HBχPT) [38]. Since the axial coupling is a dimensionless quantity, the χPT is parameterized in terms of dimensionless parameters, such that the low-energy constants (LEC) are naturally O(1). The normalized pion decay constant is chosen such that F π ∼ 92 MeV. The lattice spacing a is made dimensionless by the gradient-flow scale w 0 , which is precisely and accurately derived from Tab. IV of Ref. [39]. The factor of 1/4π in the 2 a parameterization makes the numerical value of 2 π ∼ 2 a , allowing for a double expansion around lattice spacing and pion mass with compatible power counting. We also observe that the a parameterization results in O(1) coefficients for the leading discretization corrections.  The NNLO chiral extrapolation for g A [40] is expressed as, Additionally we include the m 4 π analytic terms and up to NNLO discretization effects in the Symanzik expansion, and NLO finite volume corrections [41] δ L ( π , m π L) ≡g A (L) − g A (∞) where F 1 , F 3 are related to Bessel functions of the second kind. The extrapolations to the physical pion mass, continuum limit, and infinite volume are performed simultaneously with the preferred fit function, Additionally, we consider the possibility of residual discretization errors of the form a 1 a/w 0 where a 1 = O(m res ), and generic one-loop discretization errors of the form s 2 α s 2 a where s 2 = O(1). We also check for model dependency by considering a Taylor expansion around the chiral point, We present our final chiral-continuum extrapolation g A in Fig. 7 (Top Left), and observe that NNLO χPT describes the data well, with a χ 2 aug. /d.o.f. = 0.43. The fits are performed under the Bayesian framework, with O(10) priors for all LECs except for terms that are quadratic in 2 a,π , which are set with O(1) priors. We use the definition of the augmented-χ 2 , and the corresponding degreesof-freedom counting introduced in Appendix B of Ref. [9]. The chiral expansion converges after including the 3 π non-analytic term, as suggested by Fig. 7 (Top Right). The curvature present in the pion mass extrapolation can be attributed to the large cancellation between the 2 π and 3 π terms of Eq. (22). A plot of the physical point extrapolation as a function of the lattice spacing is shown in the bottom left of Fig. 7. There is no evidence of lattice spacing dependence that is quadratic in 2 a . In fact, we observe minimal lattice spacing dependence, as our coarsest lattice matrix elements are within 6% of our final extrapolated result. The bottom right of Fig. 7 show the infinite volume extrapolation for the m π ∼ 220 MeV ensembles. We observe that the finite volume prediction from NLO χPT adequately describes our data, and in addition, the finite volume matrix elements are consistent with the infinite volume extrapolation, demonstrating that finite volume corrections are a small effect.

Systematic error analysis
The robustness of the preferred chiral-continuum extrapolation is demonstrated by the stability of the extrapolated value of g A under an array of fit variations shown in Fig. 8(Left). The group of fits labeled "N (n) LO χPT" shows that while the pion mass extrapolation receives large corrections up to NNLO, the fit stabilizes onwards.
The χPT extrapolation is checked against a linear and quadratic in 2 π Taylor expansion, and shows perfect consistency. Since the extrapolation is performed under a work, we check that priors of the LECs up to the 3 π non-analytic term are not constraining. Successively doubling the prior widths of the LECs result in insignificant changes in the final result.
The robustness of the continuum extrapolation is demonstrated in the following three fits: O(10) NNNLOa, +O(α s a 2 ) disc., +O(a) disc.. The first of the three fits increase the prior width of the 4 a LEC by one order of magnitude. Because the lattice calculation is performed only on three lattice spacings, we see that the posterior distribution of this LEC is determined mainly by the prior. We observe however, that the resulting fit is still consistent with the preferred fit, and assert that χPT power counting adequately captures the uncertainty at 4 a . The next two fits in this group check for possible generic one-loop discretization errors, and the residual O(a) discretization error. We observe negligible dependence in both contributions.
The systematic uncertainty from the finite volume corrections are explored in the fit labeled "omit NLO FV", and the result is observed to be consistent with the preferred fit.
The right two panels of Fig. 8   The preferred fit lives within this group of fits and is highlighted in black. The next two fits are to a Taylor expansion, the next group labeled "2x(N n )LO width" test the dependency on priors, the next group evaluates the systematic uncertainty coming from the continuum extrapolation, and finally the systematic from truncating the finite volume correction at NLO is shown. The right two panels provides the χ 2 aug /d.o.f. and the log Bayes Factors for the various fits. The dashed vertical line in both panels highlights the value of the preferred fit. (Right) A summary of lattice results with regard to g A is compiled for LHPC05 [3], RBC/UKQCD08 [42], CLS12 [4], QCDSF13 [5], RQCD14 [6], ETMC15 [7], PNDME16 [8], LHPC17 [43],CalLat17A [44], ETMC17 [45], and CLS17 [46]. The averaged experimental determination is obtained from Ref. [1]. For simplicity, lattice results presented with statistical and systematic errors are added together in quadrature. The results with an open symbol are obtained from only one lattice spacing. The results with a † are extrapolated from the quantity g A /F π .
show that if we ignore the Taylor expansion extrapolations, the preferred fit is the most likely model which reproduces the lattice correlator data. The χPT is physically motivated, in contrast to the Taylor expansion, and is the model used in the final result that is quoted.
Our post-diction of the nucleon axial coupling is g QCD A = 1.285 ± 0.017 (27) with a final uncertainty budget of statistical 1.29% chiral extrapolation 0.21% continuum extrapolation 0.10% infinite volume 0.23% isospin breaking 0.04% total 1.33% Fig. 8 (Right) summarizes the improvement of the lattice determination of g A resulting from this work.

Conclusions and Outlook
The nucleon axial coupling has been a long standing challenge to determine using the methods of Lattice QCD. Two main challenges include overcoming the exponential noise problem, which in the case of the nucleon, leads to ill-behaved correlated fluctuations absent in meson quantities, as well as controlling the excited state contributions to the nucleon correlator. Using the Feynman-Hellmann method to generate lattice correlators, we are able to leverage data at small source-sink separation time in order to gain access to exponentially more precise data. In addition, the correlator stability plots demonstrate control over excited state contamination. As a result, we obtain a determination of g A at the percent-level commensurate with the experimental measurement. The post-diction of g A is the first step towards making a quantitative connection between nuclear physics and QCD, opening the door to calculations of nuclear quantities difficult or even inaccessible through experiment.