CalLat Collaboration ††thanks: Not all in California

# Detailed analysis of excited state systematics in a lattice QCD calculation
of $g_{A}$

Jinchen He (何晋琛)  Department of Physics, University of Chinese Academy of
Sciences, Beijing, 100049, P. R. China  Department of Physics, University of
California, Berkeley, CA 94720, USA  David A. Brantley  Physics Division,
Lawrence Livermore National Laboratory, Livermore, CA 94550, USA  Chia Cheng
Chang (張家丞)  Interdisciplinary Theoretical and Mathematical Sciences Program
(iTHEMS), RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Department of Physics, University of California, Berkeley, CA 94720, USA  Ivan
Chernyshev  Department of Physics, University of California, Berkeley, CA
94720, USA  Department of Physics, University of Washington, Seattle, WA
98195, USA  Evan Berkowitz  Department of Physics, University of Maryland,
College Park, MD 20742, USA  Dean Howarth  Physics Division, Lawrence
Livermore National Laboratory, Livermore, CA 94550, USA  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Christopher Körber  Institut für Theoretische Physik II, Ruhr-Universität
Bochum, D-44780 Bochum, Germany  Nuclear Science Division, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  Aaron S. Meyer  Department of
Physics, University of California, Berkeley, CA 94720, USA  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Henry Monge-Camacho  Department of Physics and Astronomy, University of North
Carolina, Chapel Hill, NC 27516-3255, USA  Escuela de Física , Universidad de
Costa Rica, San José, San Pedro, 11501, Costa Rica  Enrico Rinaldi  Physics
Department, University of Michigan, Ann Arbor, MI 48109, USA  Theoretical
Quantum Physics Laboratory, RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan
Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS),
RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Chris Bouchard  School of
Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, UK  M.A. Clark
NVIDIA Corporation, 2701 San Tomas Expressway, Santa Clara, CA 95050, USA
Arjun Singh Gambhir  Design Physics Division, Lawrence Livermore National
Laboratory, Livermore, CA 94550, USA  Nuclear Science Division, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, USA  Christopher J. Monahan
Department of Physics, The College of William & Mary, Williamsburg, VA 23187,
USA  Theory Center, Thomas Jefferson National Accelerator Facility, Newport
News, VA 23606, USA  Amy Nicholson  Department of Physics and Astronomy,
University of North Carolina, Chapel Hill, NC 27516-3255, USA  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Pavlos Vranas  Physics Division, Lawrence Livermore National Laboratory,
Livermore, CA 94550, USA  Nuclear Science Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  André Walker-Loud  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Physics Division, Lawrence Livermore National Laboratory, Livermore, CA 94550,
USA  Department of Physics, University of California, Berkeley, CA 94720, USA
https://callat-qcd.github.io

###### Abstract

Excited state contamination remains one of the most challenging sources of
systematic uncertainty to control in lattice QCD calculations of nucleon
matrix elements and form factors. Most lattice QCD collaborations advocate for
the use of high-statistics calculations at large time separations ($t_{\rm
sep}\gtrsim 1$ fm) to combat the signal-to-noise degradation. In this work we
demonstrate that, for the nucleon axial charge, $g_{A}$, the alternative
strategy of utilizing a large number of relatively low-statistics calculations
at short to medium time separations ($0.2\lesssim t_{\rm sep}\lesssim 1$ fm),
combined with a multi-state analysis, provides a more robust and economical
method of quantifying and controlling the excited state systematic
uncertainty, including correlated late-time fluctuations that may bias the
ground state. We show that two classes of excited states largely cancel in the
ratio of the three-point to two-point functions, leaving the third class, the
transition matrix elements, as the dominant source of contamination. On an
$m_{\pi}\approx 310$ MeV ensemble, we observe the expected exponential
suppression of excited state contamination in the Feynman-Hellmann correlation
function relative to the standard three-point function; the excited states of
the regular three-point function reduce to the 1% level for $t_{\rm sep}>2$ fm
while, for the Feynman-Hellmann correlation function, they are suppressed to
1% at $t_{\rm sep}\approx 1$ fm. Independent analyses of the three-point and
Feynman-Hellmann correlators yield consistent results for the ground state.
However, a combined analysis allows for a more detailed and robust
understanding of the excited state contamination, improving the demonstration
that the ground state parameters are stable against variations in the excited
state model, the number of excited states, and the truncation of early-time or
late-time numerical data.

††preprint: JLAB-THY-21-3350††preprint: RIKEN-iTHEMS-Report-21

##  I  Introduction

Lattice QCD (LQCD) calculations of nucleon matrix elements have reached a
level of maturity for inclusion in the most recent Flavour Lattice Averaging
Group (FLAG) review Aoki  _et al._ (2020). It is now common that results are
obtained with multiple lattice spacings, multiple volumes and pion masses at
or near the physical pion mass. Control over the continuum, infinite volume
and physical pion mass extrapolations are necessary to compare LQCD results
amongst themselves as well as with experiments.

However, there is an additional source of systematic uncertainty in the
calculations which must be brought under control before the extrapolations can
be relied upon, and that is the excited state contamination of the correlation
functions. The source of the issue is tied to the well-known signal-to-noise
(S/N) problem Lepage (1989). At early time, where the stochastic noise is
under control, the correlation functions have significant contamination from
excited states, while at large time, where there is ground state saturation,
the noise overwhelms the signal. To date, there are no calculations of nucleon
three-point functions performed at light pion masses with sufficient
statistics that the ground state matrix element can be determined without
addressing the contributions from excited states.

The FLAG review summarizes challenges in controlling the excited state
contamination in these calculations and the strategies various groups use to
do so. In this article, we will focus on two important points that are not
discussed in the FLAG review and which are generally lacking in the
literature:

  1. 1.

Stability of the ground state matrix element under the truncation of data
and/or the variation in the number of excited states used in the analysis;

  2. 2.

Quantification of the excited state contamination.

The first issue is generally not discussed because most calculations utilize
too few values of fixed source-sink separation times ($t_{\rm sep}$) for data
truncation to be feasible. A larger number of matrix elements are required to
perform multi-state fits, so a limited number of $t_{\rm sep}$ values used
limits the number of excited states that can be removed. In the flavor physics
LQCD community, stability of the extracted ground state observable over an
appreciable time range (sufficiently large to rule out correlated
fluctuations) has long been recognized as crucial for ruling out excited state
contamination in two-point and three-point calculations Lepage  _et al._
(2002); Wingate  _et al._ (2003); Gamiz  _et al._ (2009); Bazavov  _et al._
(2016, 2019); any residual contamination to the ground state observable from
excited states will either be observed as a trend in time of its value, or is
smaller than the precision with which the observable has been extracted. For
nucleon correlation functions, it is even more crucial that this stability be
demonstrated as they are more susceptible to correlated late-time fluctuations
through their degrading S/N ratios. Care must be taken, however, to
demonstrate this stability before the exponential growth of the noise erases
the ability to detect time dependence in the observable.

The second issue is typically addressed qualitatively as most calculations
rely upon numerical results with larger values of the source-sink separation
time, where the size of the excited-state contributions are relatively smaller
and the stochastic noise is larger, thus limiting the ability to obtain a
controlled, quantitative understanding of the excited state contributions to
the correlation function.

In order to address these issues, we have generated results with a large
number of short to intermediate values of the source-sink separation time (13
values for $t_{\rm sep}\approx 0.18-1.22$ fm) on an a09m310 ($a\approx 0.09$
fm, $m_{\pi}\approx 310$ MeV) ensemble. The large number of $t_{\rm sep}$
values with precise numerical results allows us to include up to five states
in the correlation function analysis while performing a variety of data cuts,
as discussed in detail in Sec. III.

We focus our analysis and discussion on $g_{A}$, the nucleon matrix element of
the axial current in the forward limit, as this matrix element has proved to
be one of the most challenging regarding control of excited state
contamination. The first LQCD calculation of $g_{A}$ with relatively light
dynamical quarks ($m_{\pi}\gtrsim 350$ MeV) appeared in 2005 Edwards  _et al._
(2006), resulting in a value that had a 7% statistical uncertainty and agreed
with the experimental value after extrapolation to the physical pion mass.
This led the community to anticipate $g_{A}$ would soon become a precision
benchmark quantity for LQCD.

However, subsequent calculations confounded these expectations with the
results remaining roughly independent of the pion mass (and below the physical
value) or, worse, trending away from the physical value as the pion mass was
reduced Yamazaki  _et al._ (2009); Bratt  _et al._ (2010); Alexandrou  _et
al._ (2011). It was speculated that the issue might be due to finite volume
corrections which were much larger than predicted by chiral perturbation
theory ($\chi$PT) Jaffe (2002); Cohen (2002); Yamazaki  _et al._ (2008). As
groups investigated the sensitivity of the extracted matrix elements as a
function of $t_{\rm sep}$, and added an excited state in the fit model for the
correlation function, it became clear that the dominant unresolved issue was
contamination from excited states Green  _et al._ (2011); Dinter  _et al._
(2011); Capitani  _et al._ (2012); Green  _et al._ (2014); Bhattacharya  _et
al._ (2014); Bali  _et al._ (2014). After this, a number of calculations were
performed Capitani  _et al._ (2012); Horsley  _et al._ (2014); Bali  _et al._
(2015); Abdel-Rehim  _et al._ (2015); Bhattacharya  _et al._ (2016);
Alexandrou  _et al._ (2017); Capitani  _et al._ (2019); Berkowitz  _et al._
(2017a); Chang  _et al._ (2018); Berkowitz  _et al._ (2018a); Walker-Loud  _et
al._ (2020); Gupta  _et al._ (2018); Hasan  _et al._ (2019); Alexandrou  _et
al._ (2020); Park  _et al._ (2021) that were in agreement with the physical
value of $g_{A}$ Zyla  _et al._ (2020); Bopp  _et al._ (1986); Yerozolimsky
_et al._ (1997); Liaud  _et al._ (1997); Mostovoi  _et al._ (2001); Schumann
_et al._ (2008); Mund  _et al._ (2013); Mendenhall  _et al._ (2013); Märkisch
_et al._ (2019); Beck  _et al._ (2020). Following our computation with the
Feynman-Hellmann method described in Ref. Bouchard  _et al._ (2017), other
groups have also utilized a larger number of $t_{\rm sep}$ values at the
physical pion mass Hasan  _et al._ (2019); Alexandrou  _et al._ (2020) and
found an improved understanding of excited states.111To the best of our
knowledge, the use a large number of $t_{\rm sep}$ values was first advocated
for in Refs. Detmold  _et al._ (2012a, b), which demonstrated the benefit with
LQCD calculations of heavy-hadron axial matrix elements using five $t_{\rm
sep}$ values.

Despite this progress, there remains some tension in the literature, in
particular between our results Berkowitz  _et al._ (2017a); Chang  _et al._
(2018); Berkowitz  _et al._ (2018a); Walker-Loud  _et al._ (2020) and those
from the PNDME Collaboration Bhattacharya  _et al._ (2016); Gupta  _et al._
(2018), both of which are the only results to utilize three (CalLat) or four
(PNDME) lattice spacings and physical pion masses. Both sets of results were
generated with mixed actions that use the $N_{F}=2+1+1$ Highly Improved
Staggered Quarks (HISQ) Follana  _et al._ (2007) in the sea-quark sector
generated by the MILC collaboration Bazavov  _et al._ (2013) and, in the
former case, also by the CalLat Collaboration Miller  _et al._ (2020, 2021).
The CalLat results are generated with a Möbius Brower  _et al._ (2017) Domain-
Wall Kaplan (1992); Shamir (1993); Furman and Shamir (1995) Fermion (MDWF)
valence action Berkowitz  _et al._ (2017b) and are computed with
$a\approx\\{0.09,0.12,0.15\\}$ fm lattice spacings while the PNDME results are
generated with a tadpole-improved Lepage and Mackenzie (1993) clover-Wilson
valence action with $a\approx\\{0.06,0.09,0.12,0.15\\}$ fm.

In Ref. Gupta  _et al._ (2018), it was shown that the tension between the
CalLat and PNDME results is driven by the PNDME results on the $a\approx 0.06$
fm ensembles, which tend to pull the final result to a smaller value,
suggestive that the discrepancy may be a discretization effect. However, it
was pointed out that the lever-arm in values of $t_{\rm sep}$ between the
smallest and largest source-sink separation times was the smallest on the
$a\approx 0.06$ fm ensembles, and that the high-correlation between
neighboring time-slices on these fine ensembles makes them more susceptible to
correlated fluctuations Walker-Loud . Further, subsequent analysis by PNDME
demonstrated an under-reported excited-state fitting systematic in their
results which seems to alleviate the tension Jang  _et al._ (2020).

Another difference between the CalLat and PNDME results is that the PNDME
results rely on the more common fixed source-sink separation method while the
CalLat results utilize a variant of the summation method Maiani  _et al._
(1987) which can be derived with the Feynman-Hellmann Theorem Bouchard  _et
al._ (2017). As was shown in Ref. Capitani  _et al._ (2012), the summation
method suppresses excited states more than the standard fixed source-sink
separation method, and as we will discuss in some detail in Sec. II, the
Feynman-Hellmann derived correlator suppresses excited states even more than
the summation method. The excited state systematic uncertainty, therefore,
deserves more scrutiny. This has been recognized by the community, illustrated
by the focus on excited state contamination in the most recent review on
nucleon structure from the annual lattice field theory symposium Ottnad
(2021).

The challenge of controlling calculations of $g_{A}$ and other matrix elements
has inspired a series of papers aimed at understanding the excited state
contamination by utilizing chiral perturbation theory Tiburzi (2009); Bär and
Golterman (2013); Bar (2015); Tiburzi (2015); Bär (2016); Bar (2017). This
work has led to further ideas to try to improve the calculation of the nucleon
axial form-factor Bali  _et al._ (2020); Bar (2020). It is worth noting,
however, that $\chi$PT predicts that the excited state contributions should
shift the correlation function above its asymptotic value, while numerical
results from all calculations show that the ground state limit is approached
from below, and thus, there is a significant discrepancy between this
theoretical prediction and the numerical data. While there are some
significant indications that $SU(2)$ baryon chiral perturbation theory,
without explicit delta degrees of freedom, is not a converging expansion even
at the physical pion mass Walker-Loud  _et al._ (2009); Walker-Loud (2008,
2013, 2014); Chang  _et al._ (2018); Walker-Loud  _et al._ (2020); Drischler
_et al._ (2019), one might anticipate that the predictions from chiral
perturbation theory should be at least qualitatively correct.

In this work, we take a data-driven approach and ask, given a large dataset,
what can we learn about the excited state contamination of the nucleon axial-
vector three-point function? We begin with a summary of the spectral
representation of the three point functions in Sec. II, then we turn to our
numerical results and analysis in Sec. III. We offer some observations and
conclusions in Sec. IV, and we present extensive details of our results and
analysis in the appendices.

##  II  Spectral Decomposition

Lattice QCD calculations are performed in Euclidean space in a mixed time-
momentum basis. In this paper, we focus on the forward matrix element at zero
momentum. Most LQCD calculations are performed with a local creation operator
and a momentum-space annihilation operator. With such a setup, the relevant
two-point correlation function at zero momentum and time separation $t_{\rm
sep}$ is given by

$\displaystyle C_{2}(t_{\rm sep})$
$\displaystyle=\sum_{\mathbf{x}}\langle\Omega|N(t_{\rm
sep},\mathbf{x})N^{\dagger}(0,\mathbf{0})|\Omega\rangle$
$\displaystyle=\sum_{n=0}^{\infty}|z_{n}|^{2}e^{-E_{n}t_{\rm sep}}$
$\displaystyle=|z_{0}|^{2}e^{-E_{0}t_{\rm
sep}}\left[1+\sum_{n=1}^{\infty}|r_{n}|^{2}e^{-{\Delta}_{n0}t_{\rm
sep}}\right]\,.$ (1)

In this expression, we assume that the overlap factors,
$z_{n}=\langle{\Omega}|N|n\rangle$, used to create ($N^{\dagger}$) and
annihilate ($N$) the states with quantum numbers of the nucleon from the
vacuum ($|{\Omega}\rangle$) are conjugate to each other.222We are using the
non-relativistic normalization $\langle n|m\rangle={\delta}_{nm}$ and
$\mathds{1}\equiv{|{\Omega}\rangle}{\langle{\Omega}|}+\sum_{n=0}^{\infty}{|n\rangle}{\langle
n|}$. We assume that the contributions from the thermal bath or alternate
temporal boundary conditions are sufficiently suppressed that they can be
ignored. In the last line, we have defined the energy splitting and ratio of
overlap factors:

$\displaystyle{\Delta}_{mn}=E_{m}-E_{n}\,,$ $\displaystyle
r_{n}=\frac{z_{n}}{z_{0}}\,.$ (2)

The parameterization of Eq. (7a) recasts all excited-state parameters with
respect to the ground-state, and yields a more universal set of excited-state
distributions, simplifying the estimation of their starting values in a
frequentist minimization or prior distribution in a Bayesian minimization Hörz
_et al._ (2021).

###  II.1 Three-point correlators

The matrix elements of interest (with Dirac and isospin structure ${\Gamma}$)
are determined through an analysis of three-point correlation functions which
are also computed in a mixed time-momentum basis. The most common strategy is
to use a sink with fixed definite spatial momentum at $t=t_{\rm sep}$ with a
current insertion ($j_{\Gamma}$) at ${\tau}$, between the source
($N^{\dagger}$) at $t=0$ and the sink ($N$). In the limit of zero momentum and
zero-momentum transfer, the three-point function is given by

$\displaystyle C_{\Gamma}(t_{\rm sep},{\tau})$
$\displaystyle=\sum_{\mathbf{y,x}}\langle\Omega|N(t_{\rm
sep},\mathbf{y})\,j_{\Gamma}({\tau},\mathbf{x})\,N^{\dagger}(0,\mathbf{0})|\Omega\rangle$
$\displaystyle=\sum_{n}|z_{n}|^{2}\ g^{\Gamma}_{nn}\ e^{-E_{n}t_{\rm sep}}$
$\displaystyle\phantom{=}+2\sum_{n<m}z_{n}z^{\dagger}_{m}\ g^{\Gamma}_{nm}\
e^{-(E_{n}+\frac{{\Delta}_{mn}}{2})t_{\rm sep}}$
$\displaystyle\phantom{=}\qquad\times\cosh\left[{\Delta}_{mn}\left({\tau}-\frac{t_{\rm
sep}}{2}\right)\right]\,,$ (3)

where the matrix elements of interest are given by

$\displaystyle g^{\Gamma}_{mn}=\langle m|j_{\Gamma}|n\rangle\,.$ (4)

In this limit, the correlation functions are real, so $z^{\dagger}_{n}=z_{n}$
and $g^{\Gamma}_{nm}=g^{\Gamma}_{mn}$, which are used to simplify Eq. (II.1).

Typically one determines the ground state matrix element, $g^{\Gamma}_{00}$,
by performing the computation for 1, 2, 3 or sometimes 4 values of $t_{\rm
sep}$ in the range $t_{\rm sep}\approx 1-1.5$ fm while utilizing local three-
quark interpolating fields for the creation and annihilation operators. The
correlation functions are then analyzed using models with zero, one, or more
excited states Capitani  _et al._ (2012); Horsley  _et al._ (2014); Bali  _et
al._ (2015); Abdel-Rehim  _et al._ (2015); Bhattacharya  _et al._ (2016);
Alexandrou  _et al._ (2017); Capitani  _et al._ (2019); Berkowitz  _et al._
(2017a); Chang  _et al._ (2018); Berkowitz  _et al._ (2018a); Walker-Loud  _et
al._ (2020); Gupta  _et al._ (2018); Hasan  _et al._ (2019); Alexandrou  _et
al._ (2020); Park  _et al._ (2021), with some analysis using up to four states
Yoon  _et al._ (2017); Gupta  _et al._ (2017); Jang  _et al._ (2020). There
are several challenges and shortcomings with this strategy, which we will
summarize. However, one should note that this method has been used, as so far,
it has seemed the most economical. Improvements upon this strategy would
require a substantial increase in computational complexity, effort, and
computing resources.

The ground state (gs) matrix element $g^{\Gamma}_{00}$ is typically determined
by constructing the ratio correlation function

$R_{\Gamma}(t_{\rm sep},{\tau})=\frac{C_{\Gamma}(t_{\rm
sep},{\tau})}{C_{2}(t_{\rm sep})}\,,$ (5)

which gives the ground state matrix element in the limit

$\lim_{t_{\rm sep}\rightarrow\infty}R_{\Gamma}(t_{\rm sep},{\tau}\approx
t_{\rm sep}/2)=g^{\Gamma}_{00}\,.$ (6)

The energies and overlap factors can be constrained from the two-point
function, leaving the three-point function to constrain the matrix elements
$g^{\Gamma}_{nn}$ and $g^{\Gamma}_{nm}$.333Some groups fix the overlap factors
and energies from the two-point correlation function and then pass the central
or correlated values into a three-point function analysis. Instead, we perform
a simultaneous fit to the two- and three-point functions.  From Eq. (II.1),
one observes that only the transition terms ($n\neq m$) are sensitive to the
current insertion time ${\tau}$. To isolate the ground state matrix element,
the minimum number of values of $t_{\rm sep}$ required is at least one greater
than the number of states used in the analysis in order to have a 1-degree-of-
freedom fit. For example, a two-state fit requires the determination of
$g^{\Gamma}_{00}$, $g^{\Gamma}_{11}$ and $g^{\Gamma}_{01}$: the latter can be
constrained from the ${\tau}$ dependence leaving the remaining $t_{\rm sep}$
dependence to constrain the two former matrix elements. In this minimal
scenario, computations which utilize three values of $t_{\rm sep}$ (or less)
are not able to perform a systematic study on omitting values of $t_{\rm
sep}$, or changing the number of states used in the analysis. This prohibits a
verification that the full uncertainty on $g^{\Gamma}_{00}$ associated with
excited state contamination has been correctly captured. In other words, we
would have to assign an unquantified systematic uncertainty to the ground
state matrix element.

Another significant challenge is that lattice computations of three-point
functions are typically performed for $t_{\rm sep}\gtrsim 1$ fm, which is
roughly the time separation when the stochastic noise becomes significant.
Again, because of the limited values of $t_{\rm sep}$ typically used, one
cannot determine if the data at this time is susceptible to a correlated
fluctuation or not, which, if present, would cause a bias in the results. We
will return to this point in Sec. III.

To understand the various sources of excited state (es) contamination, we
reorder Eq. (5) in a way which deliberately disentangles the different types
of excited state contributions. For the two-point correlation function, it is
straightforward to separate the ground state from the excited states

$\displaystyle C_{2}(t_{\rm sep})$ $\displaystyle=C_{2}^{\rm gs}(t_{\rm
sep})+C_{2}^{\rm es}(t_{\rm sep})\,,$ (7a) $\displaystyle C_{2}^{\rm
gs}(t_{\rm sep})$ $\displaystyle=|z_{0}|^{2}e^{-E_{0}t_{\rm sep}}\,,$ (7b)
$\displaystyle C_{2}^{\rm es}(t_{\rm sep})$ $\displaystyle=\sum_{n\geq
1}|z_{n}|^{2}e^{-E_{n}t_{\rm sep}}\,.$ (7c)

For the three-point functions, we define the $n$-to-$n$ as scattering (sc)
states and the $n$-to-$m$ as transition (tr) states such that

$\displaystyle C_{\Gamma}(t_{\rm sep},{\tau})$ $\displaystyle=C_{\Gamma}^{\rm
gs}(t_{\rm sep})+C_{\Gamma}^{\rm sc}(t_{\rm sep})+C_{\Gamma}^{\rm tr}(t_{\rm
sep},{\tau})\,,$ (8a) $\displaystyle C_{\Gamma}^{\rm gs}(t_{\rm sep})$
$\displaystyle=|z_{0}|^{2}g^{\Gamma}_{00}e^{-E_{0}t_{\rm
sep}}=g^{\Gamma}_{00}C_{2}^{\rm gs}(t_{\rm sep})\,,$ (8b) $\displaystyle
C_{\Gamma}^{\rm sc}(t_{\rm sep})$ $\displaystyle=\sum_{n\geq
1}|z_{n}|^{2}g^{\Gamma}_{nn}e^{-E_{n}t_{\rm sep}}\,,$ (8c) $\displaystyle
C_{\Gamma}^{\rm tr}(t_{\rm sep},{\tau})$
$\displaystyle=\sum_{n<m}z_{n}z^{\dagger}_{m}\ 2g^{\Gamma}_{nm}\
e^{-(E_{n}+\frac{{\Delta}_{mn}}{2})t_{\rm sep}}$
$\displaystyle\qquad\qquad\times\cosh\left[{\Delta}_{mn}\left({\tau}-\frac{t_{\rm
sep}}{2}\right)\right]\,.$ (8d)

The ratio correlation function can then be expressed as

$\displaystyle R_{\Gamma}(t_{\rm sep},{\tau})$
$\displaystyle=\frac{C_{\Gamma}^{gs}(t_{\rm sep})+C_{\Gamma}^{\rm sc}(t_{\rm
sep})+C_{\Gamma}^{\rm tr}(t_{\rm sep},{\tau})}{C_{2}(t_{\rm sep})}$
$\displaystyle=g^{\Gamma}_{00}+\frac{C_{\Gamma}^{\rm sc}(t_{\rm
sep})-g^{\Gamma}_{00}C_{2}^{\rm es}(t_{\rm sep})}{C_{2}(t_{\rm sep})}$
$\displaystyle\phantom{=}+\frac{C_{\Gamma}^{\rm tr}(t_{\rm
sep},{\tau})}{C_{2}(t_{\rm sep})}$
$\displaystyle=g^{\Gamma}_{00}+\frac{\sum_{n\geq
1}(g^{\Gamma}_{nn}-g^{\Gamma}_{00})|z_{n}|^{2}e^{-E_{n}t_{\rm
sep}}}{C_{2}(t_{\rm sep})}$ $\displaystyle\phantom{=}+\frac{C_{\Gamma}^{\rm
tr}(t_{\rm sep},{\tau})}{C_{2}(t_{\rm sep})}\,.$ (9)

Consider the leading excited state contamination to $g^{\Gamma}_{00}$ arising
from $g^{\Gamma}_{11}$, and $g^{\Gamma}_{01}$

$\displaystyle R_{\Gamma}(t_{\rm sep},{\tau})$ $\displaystyle\approx
g^{\Gamma}_{00}+|r_{1}|^{2}(g^{\Gamma}_{11}-g^{\Gamma}_{00})e^{-{\Delta}_{10}t_{\rm
sep}}$
$\displaystyle\phantom{=}+2r_{1}^{\dagger}g^{\Gamma}_{01}e^{-{\Delta}_{10}\frac{t_{\rm
sep}}{2}}\cosh\left[{\Delta}_{10}\left({\tau}-\frac{t_{\rm
sep}}{2}\right)\right]$ $\displaystyle\phantom{=}+\cdots\,,$ (10)

where the $\cdots$ includes terms from higher excited states as well as from
the first excited state, but further suppressed by extra powers of
$\exp(-{\Delta}_{10}t_{\rm sep})$. The scattering and two-point excited state
contributions, in addition to cancelling against each other for same-sign
values of $g^{\Gamma}_{11}$ and $g^{\Gamma}_{00}$, are suppressed by the full
$t_{\rm sep}$, $\exp(-{\Delta}_{10}t_{\rm sep})$. In contrast, the transition
excited states are only suppressed by half the time separation,
$\exp(-{\Delta}_{10}\frac{t_{\rm sep}}{2})$, and so they are expected to be
the dominant source of excited state contamination.

###  II.2 Feynman-Hellmann and Summed Correlators

Rather than analyzing the three-point functions, one can construct a
correlation function in which the current insertion time is summed over,
cutting $\tau_{c}$ timeslices from either end of the correlation function

$\displaystyle{\Sigma}_{\Gamma}(t_{\rm sep},{\tau}_{c})$
$\displaystyle=\sum_{{\tau}={\tau}_{c}}^{t_{\rm
sep}-{\tau}_{c}}C_{\Gamma}(t_{\rm sep},{\tau})$
$\displaystyle=\sum_{n}|z_{n}|^{2}(t_{\rm
sep}+1-2{\tau}_{c})g^{\Gamma}_{nn}e^{-E_{n}t_{\rm sep}}$
$\displaystyle\phantom{=}+2\sum_{n<m}z_{n}z^{\dagger}_{m}g^{\Gamma}_{nm}e^{-(E_{n}+\frac{{\Delta}_{mn}}{2})t_{\rm
sep}}$
$\displaystyle\phantom{=}\quad\times\frac{\sinh\left[\frac{{\Delta}_{mn}}{2}(t_{\rm
sep}+1-2{\tau}_{c})\right]}{\sinh\left(\frac{{\Delta}_{mn}}{2}\right)}\,.$
(11)

The original implementation of this idea included a sum over all timeslices,
including the contact terms (${\tau}=0$ and ${\tau}=t_{\rm sep}$) as well as
“out-of-time-order” regions (${\tau}<0$ and ${\tau}>t_{\rm sep}$) Maiani  _et
al._ (1987). Given a set of fixed source-sink separation datasets, one is of
course free to perform the sum in a variety of ways, e.g. excluding data from
the out-of-time-order regions, excluding the source and sink time (setting
${\tau}_{c}=1$ in Eq. (II.2)) or also cutting time near the source and sink
(${\tau}_{c}>1$).

The summed correlation function at large $t_{\rm sep}$

${\Sigma}_{\Gamma}(t_{\rm sep},{\tau}_{c})=|z_{0}|^{2}e^{-E_{0}t_{\rm
sep}}\bigg{\\{}\\\ (t_{\rm
sep}+1-2{\tau}_{c})\Big{[}g^{\Gamma}_{00}+|r_{1}|^{2}g^{\Gamma}_{11}e^{-{\Delta}_{10}t_{\rm
sep}}\Big{]}\\\
+r_{1}^{\dagger}g^{\Gamma}_{01}\frac{e^{\frac{{\Delta}_{10}}{2}(1-2{\tau}_{c})}-e^{-{\Delta}_{10}(t_{\rm
sep}+\frac{1}{2}-{\tau}_{c})}}{\sinh(\frac{{\Delta}_{10}}{2})}\bigg{\\}}+\cdots\,,$
(12)

can be used to determine the leading excited state contamination. Note that
the ground state and the scattering ($n$-to-$n$) states are relatively
enhanced by $t_{\rm sep}+1-2{\tau}_{c}$. The transition matrix elements
($m$-to-$n$) lead to a $t_{\rm sep}$-independent term and those that depend
upon $t_{\rm sep}$ become exponentially suppressed by the full excited state
gap ($e^{-{\Delta}_{10}t_{\rm sep}}$) rather than half the gap, as with the
three-point function, Eq. (II.1), as noted in Ref. Capitani  _et al._ (2012).
Thus, one expects that the excited state contamination of the summed
correlation function is smaller than for the standard fixed source-sink
separation time three-point correlation function, up to this $t_{\rm
sep}$-independent term.

This summed correlation function has received some attention in the literature
Gusken  _et al._ (1989); Bulava  _et al._ (2012); Capitani  _et al._ (2012);
de Divitiis  _et al._ (2012). Recently, there have been calculations which
utilize a Feynman-Hellmann approach, by performing a computation in the
presence of background fields and extracting the matrix elements through the
linear response of the spectrum to the background field Chambers  _et al._
(2014, 2015); Savage  _et al._ (2017). In Ref. Bouchard  _et al._ (2017), it
was shown that the application of the Feynman-Hellmann theorem to the
effective mass directly leads to a derivative of the summed correlation
function Maiani  _et al._ (1987), relating the matrix element to the spectrum
without the need for an explicit background field. We call this the Feynman-
Hellmann (FH) correlation function444In the original implementation Bouchard
_et al._ (2017), the sum over the current time is over all time-slices, which
we denote ${\tau}_{c}={\rm none}$, as in the original summation method Maiani
_et al._ (1987).

$\displaystyle{\rm FH}_{\Gamma}(t_{\rm sep},{\tau}_{c},dt)$
$\displaystyle=\sum_{{\tau}={\tau}_{c}}^{t_{\rm
sep}-{\tau}_{c}}\frac{R_{\Gamma}(t_{\rm sep}+dt,{\tau})-R_{\Gamma}(t_{\rm
sep},{\tau})}{dt}$
$\displaystyle=\frac{1}{dt}\left[\frac{{\Sigma}_{\Gamma}(t_{\rm
sep}+dt,{\tau}_{c})}{C_{2}(t_{\rm sep}+dt)}-\frac{{\Sigma}_{\Gamma}(t_{\rm
sep},{\tau}_{c})}{C_{2}(t_{\rm sep})}\right],$ $\displaystyle{\rm
FH}_{\Gamma}(t_{\rm sep},{\tau}_{c})$ $\displaystyle\equiv{\rm
FH}_{\Gamma}(t_{\rm sep},{\tau}_{c},dt=1)\,.$ (13)

Since the FH correlation function is constructed from
${\Sigma}_{\Gamma}(t_{\rm sep},{\tau}_{c})$, it enjoys the larger suppression
of excited states, with the leading excited state contamination scaling as
$\exp(-{\Delta}_{10}t_{\rm sep})$ rather than with ${\Delta}_{10}/2$ (the
$t_{\rm sep}$ independent pieces exactly cancel in the numerical derivative).
Additionally, the numerical derivative serves to both isolate the ground
state, whose contribution grows linearly in $t_{\rm sep}$, as well as to
further suppress the scattering ($n$-to-$n$) and transition ($m$-to-$n$)
excited states which do not strongly differ from one timeslice to the next.
This stronger suppression of excited states is what allowed us to utilize
earlier Euclidean time data Berkowitz  _et al._ (2017a); Chang  _et al._
(2018); Berkowitz  _et al._ (2018a) than is common in the three-point
correlation function analysis and to enjoy the benefits of the lower
stochastic noise. We will show this in some detail in Sec. III.1.

##  III  Lattice Calculation

Figure 1:  Left: (Psychedelic Moose/Water Buffalo plot) We plot the numerical
results of $R_{A_{3}}(t_{\rm sep},{\tau})$ for $t_{\rm sep}=2$ (the single
gray point in the middle) to $t_{\rm sep}=14$, the top (red) dataset. In
addition, we plot the resulting posterior description of the correlation
function from our 5-state fit as the (correspondingly colored) fit bands. The
(gray) squares in the upper left/right of the “moose antlers” are not included
in the analysis, as indicated by the break in the fit band from the inner
region. The horizontal (gray) band is the ground state matrix element,
$\mathring{g}_{A}$. Right: The ${\rm FH}_{A_{3}}(t_{\rm sep},{\tau}_{c}=1)$
numerical data from the global analysis is plotted along with the posterior
description of the correlation function. As a comparison, for even values of
$t_{\rm sep}$, we plot $R_{A_{3}}(t_{\rm sep},{\tau}=t_{\rm sep}/2)$ with open
(colored) symbols (the numcerical data in the middle of the “moose”). The
time-axis is converted from lattice units (top) to fm (bottom) using our
scale-setting Miller  _et al._ (2021).

For the present study, we use results from our MDWF on gradient-flowed HISQ
action Berkowitz  _et al._ (2017b) on the a09m310 ensemble, which has a
lattice spacing of $a\approx 0.09$ fm and a pion mass of $m_{\pi}\approx 310$
MeV. We use the same parameters as in Refs. Berkowitz  _et al._ (2017b, a);
Chang  _et al._ (2018) except for the quark-smearing, we use the parameters
${\sigma}_{\rm smr}=3.5$ and $N_{\rm smr}=45$ that go in the
GAUGE_COVARIANT_GAUSSIAN smearing routine in Chroma Edwards and Joo (2005).

We generate 16 sources per gauge configuration on 784 configurations. We
generate the three-point functions using a sequential-propagator through the
sink at 13 values of $t_{\rm sep}$,

$t_{\rm sep}/a_{09}\in\\{2,3,4,5,6,7,8,9,10,11,12,13,14\\}\,.$ (14)

Our sources and sinks are generated with a local three-quark interpolating
field using only the upper-spin components of the quark field in the Dirac-
Pauli basis (lower components for the negative parity states), which gives the
largest overlap onto the ground state of the nucleon at rest Basak  _et al._
(2005a, b). In App. A, we present further details of our computation,
including the cost-benefit analysis of improving the stochastic sampling by
combining 8 coherent sinks Bratt  _et al._ (2010) for each sequential
propagator, the use of spin up and spin down sources and sinks (versus
utilizing a spin-projector that isolates one of the spin states) and the use
of time-reversed negative parity correlators.

###  III.1  Full results

We begin with a presentation of our final results which come from a fully-
correlated Bayesian constrained curve-fit Lepage  _et al._ (2002) with a five-
state model to describe

$\displaystyle C_{2}(t_{\rm sep})\,,$ $\displaystyle R_{A_{3}}(t_{\rm
sep},{\tau})\,,$ $\displaystyle R_{V_{4}}(t_{\rm sep},{\tau})\,,$
$\displaystyle{\rm FH}_{A_{3}}(t_{\rm sep},{\tau}_{c}=1)\,,$
$\displaystyle{\rm FH}_{V_{4}}(t_{\rm sep},{\tau}_{c}=1)\,.$

The final result is obtained with all values of ${\tau}$ between the source
and sink time, ${\tau}=[1,t_{\rm sep}-1]$. For $R_{\Gamma}(t_{\rm
sep},{\tau})$, the results are symmetrized about ${\tau}=t_{\rm sep}/2$ and
half the data (plus ${\tau}=t_{\rm sep}/2$ point for even values of $t_{\rm
sep}$) is used in the analysis.

In the left panel of Fig. 1 (the “Psychedelic Moose”), we plot the numerical
results for the ratio of the three-point function generated with the
$A_{3}=\bar{q}{\gamma}_{3}{\gamma}_{5}{\tau}_{3}q$ current, divided by the
two-point function at the given value of $t_{\rm sep}$ (cfr. Eq. (5)). We also
plot the resulting posterior description with our 5-state model. The fit
quality is good and visually one can see that the model accurately describes
the numerical results over the full range of $t_{\rm sep}$ and ${\tau}$ used
in the analysis. The horizontal (gray) band is the value of the ground state
matrix element $\mathring{g}_{A}$. In the right panel we plot the ${\rm
FH}_{A_{3}}(t_{\rm sep},{\tau}_{c}=1)$ data that are used in the global fit as
well as the resulting posterior distribution of this correlation function.

Figure 2:  We plot the numerical ${\rm FH}_{A_{3}}(t_{\rm sep},{\tau}_{c})$
data for various values of ${\tau}_{c}$ as well as the posterior
reconstruction of these correlation functions from the global analysis that
uses ${\tau}_{c}=1$.

In Fig. 2, we explore the ${\rm FH}_{A_{3}}(t_{\rm sep},{\tau}_{c})$ data as
the number of data near the source and sink time are cut from the sum over
current insertion time denoted by increasing ${\tau}_{c}$. The posterior fit
bands are from the global analysis that uses ${\tau}_{c}=1$. As ${\tau}_{c}$
is increased, for a fixed $t_{\rm sep}$, one observes that the excited state
contribution becomes larger. This can be understood by looking at the leading
excited state contribution to ${\rm FH}_{\Gamma}(t_{\rm sep},{\tau}_{c})$ that
depends upon ${\tau}_{c}$

${\rm FH}_{\Gamma}(t_{\rm sep},{\tau}_{c})\ni
e^{{\Delta}_{10}({\tau}_{c}+\frac{1}{2})}r_{1}^{\dagger}g_{01}^{\Gamma}\frac{e^{-{\Delta}_{10}t_{\rm
sep}}-e^{-{\Delta}_{10}(t_{\rm sep}+1)}}{\sinh\frac{{\Delta}_{10}}{2}}.$

As is evident from Fig. 1, the transition matrix element $g_{01}^{A_{3}}<0$;
an observation common to all LQCD calculations of $g_{A}$. Therefore, the
leading excited state contamination that depends upon ${\tau}_{c}$ is negative
and grows exponentially with increasing ${\tau}_{c}$, consistent with the
results.

In order to have confidence in our final results for $m_{N}=E_{0}$,
$\mathring{g}_{A}$ and $\mathring{g}_{V}$,555With our action, the axial-vector
and vector renormalization factors are equal to $10^{-4}$ so the renormalized
axial charge is given by the ratio of the bare matrix elements,
$g_{A}=\mathring{g}_{A}/\mathring{g}_{V}$ Chang  _et al._ (2018). as well as
the ability to quantify and control the excited state contamination, we
discuss our analysis strategy and the stability of our results under model
variation and data truncation.

###  III.2  Analysis strategy

Figure 3:  Sensitivity of the extracted spectrum and $\mathring{g}_{A}$ on the
model of excited states used as a function of $t_{\rm sep}^{\rm min}$ in the
$C_{2}(t_{\rm sep})$ correlation function. For a given $t_{\rm sep}^{\rm
min}$, from left to right, offset for visual clarity, we plot the prior
(vertical (colored) box) and posterior values (filled markers) of
$\mathring{g}_{A}$, $E_{0}$ and $E_{n}$ for the harmonic-oscillator (HO),
inverse $n$ ($1/n$) and inverse $n^{2}$ ($1/n^{2}$) excited state models, Eq.
(III.2). The prior width on $\mathring{g}_{A}$ and $E_{0}$ are larger than the
displayed $y$-limits, and thus not shown. The horizontal light (red) bands (to
guide the eye) are from the chosen fit denoted with the black border on the
markers. The horizontal dashed lines denote the non-interacting energy levels
of the $P$-wave $N\pi$ states. See the text for more details.

We have two goals with this work:

  1. 1.

Identify the ground state mass and matrix elements, $\mathring{g}_{A}$ (and
$\mathring{g}_{V}$), with a complete systematic uncertainty;

  2. 2.

Obtain a quantitative understanding of the excited state contamination of the
correlation functions used in the analysis.

The former is possible without the latter through a demonstration that the
extracted ground state mass and matrix elements are invariant under
modifications of the fit model (i.e. the number of excited states and the
model of their mass-gap) and truncations in the time range used in the
analysis.

A robust identification of the excited state spectrum and matrix elements
would require the use of a variational basis that includes multi-hadron
operators Dudek  _et al._ (2013); Lang and Verduci (2013); Wilson  _et al._
(2015). For this reason, our second goal is to obtain a quantitative
understanding of the sum of all excited state contributions, and in
particular, to quantify how they contaminate the ground state values. To
achieve this goal it is essential that the fit model accurately describes the
correlation function over the full range of time separations without over-
fitting the data.

In App. B, we provide a detailed description of the analysis and the
sensitivity of the ground state posteriors under variations in the number of
states used, the model of excited states and the fit ranges in $t_{\rm sep}$
and ${\tau}$. Here, we summarize the findings of this study.

We begin with a discussion of the model for the excited states. The lowest-
lying excited states consist roughly of either a nucleon-pion in a relative
P-wave or a nucleon with a two-pion excitation. With $m_{\pi}L\approx 4$ the
non-interacting energy levels of these two states are nearly identical and so
a calculation without multi-hadron operators can not distinguish them.
Therefore, they are treated as a single excited state. We then use three
different models for the spectrum of excitations and parametrize the energy
gaps ${\Delta}E_{n}$, where $n$ indicates the $n^{th}$ state, as:

$n^{th}$ energy level $\displaystyle
E_{n}=E_{0}+\sum_{l=1}^{n}{\Delta}E_{l}\,,$ harmonic oscillator (HO)
$\displaystyle{\Delta}E_{n}=2m_{\pi}\,,$ inverse $n$ ($1/n$)
$\displaystyle{\Delta}E_{n}=2m_{\pi}/n\,,$ inverse $n^{2}$ ($1/n^{2}$)
$\displaystyle{\Delta}E_{n}=2m_{\pi}/n^{2}\,.$ (15)

For each model, for each level $n<N_{\rm max}$, we set the prior as
${\Delta}\tilde{E_{n}}=({\Delta}E_{n},m_{\pi})$ where the first entry is the
prior mean and the second is the prior width. For the highest state, which we
expect to be a “garbage can”, we set the prior width to be $5m_{\pi}$. We
observe two important facts when we vary the excited state model:

  1. 1.

The ground state posteriors are insensitive to the model used;

  2. 2.

The excited state posteriors are mostly insensitive to the model used, even
when the posterior is in significant conflict with the prior;

both strong indicators that the extracted energy levels are dictated by the
numerical data and not the priors used. In Fig. 3, we show the sensitivity of
$\mathring{g}_{A}$ and the spectrum on the excited state model versus the
$t_{\rm sep}^{\rm min}$ of the two-point correlation function (in all cases
for the full five-state, five-correlation function analysis). The lowest panel
is the relative weight determined from the three models at a given $t_{\rm
sep}^{\rm min}$

$w_{\rm i}=\frac{e^{{\rm logGBF}_{i}}}{\sum_{j\in{\rm models}}e^{{\rm
logGBF}_{j}}}\,,$ (16)

where logGBF is the log of the Gaussian Bayes Factor. The next panel gives the
$Q\in[0,1]$ value which is a measure of the fit-quality. For the values of
$\mathring{g}_{A}$ and $E_{n}$, the light (red) horizontal band is the value
of these quantities from the chosen fit from the HO model at $t_{\rm sep}^{\rm
min}=3$, to guide the eye, which is also denoted by the markers with a black
border. The vertical bars represent the prior width for the given quantity and
are aligned in the same vertical column as their corresponding posterior
values (for $\mathring{g}_{A}$ and $E_{0}$, the prior widths are larger than
the displayed $y$-limits and so they are not shown).666The energy spectrum is
priored with a series of ordered energy splittings. Therefore, to construct
the priors shown in Fig. 3, we plot
$\tilde{E}_{n}=\hat{E}_{0}+\sum_{l=1}^{n-1}{\Delta}\hat{E}_{l}+{\Delta}\tilde{E}_{n}$
where the priors are denoted with a tilde and the posteriors are denoted with
a hat. Finally, the horizontal dashed lines give the non-interacting $P$-wave
$N(\mathbf{q})\pi(-\mathbf{q})$ energy levels, which are quite dense. The
lowest non-interacting $N(\mathbf{q+p})\pi(\mathbf{-q})\pi(\mathbf{-p})$ level
where all hadrons are at rest is nearly degenerate with the lowest
$N(\mathbf{q})\pi(-\mathbf{q})$ level. In particular,

  1. 1.

The extracted posteriors are consistent between models which all have a high
fit-quality, even when the posterior value is in tension with the prior. This
is a strong indicator the extracted spectrum is highly constrained by the
numerical data and not the priors, even for the high-lying energy levels.

  2. 2.

The first excited state is consistent with the lowest lying $N\pi$ state
(which is nearly degenerate with the $N\pi\pi$ state). The higher lying states
have an uncertainty that spans several anticipated energy levels, indicating
that they are likely a linear combination of eigenstates.

|  |
---|---|---
($t_{\rm sep}^{\rm min}/t_{\rm sep}^{\rm max}$ sensitivity) | (even $t_{\rm sep}$ only) | (odd $t_{\rm sep}$ only)
Figure 4:  Stability of the determination of $\mathring{g}_{A}$ under data
truncation. The horizontal band is the value of $\mathring{g}_{A}$ determined
from the full analysis, to guide the eye. The bottom part of each panel
represents the fit quality $Q\in[0,1]$. The left panel shows the stability as
we increase $t_{\rm sep}^{\rm min}$ while holding $t_{\rm sep}^{\rm max}=14$
to the left of the vertical dashed line and, similarly, the result as $t_{\rm
sep}^{\rm max}$ is increased with $t_{\rm sep}^{\rm min}=2$ to the right of
the dashed line. The middle and right panels show the $t_{\rm sep}^{\rm max}$
stability when only even or odd values of $t_{\rm sep}$ are used in
$R_{\Gamma}(t_{\rm sep},{\tau})$ and ${\rm FH}_{\Gamma}(t_{\rm
sep},{\tau}_{c}=1,dt=2)$.

Next, we examine the stability of $\mathring{g}_{A}$ as we remove data from
the fit (all fits still use the set of five correlation functions). We examine
the sensitivity as we increase $t_{\rm sep}^{\rm min}$ in the analysis, as we
reduce $t_{\rm sep}^{\rm max}$ and also through using only even or odd values
of $t_{\rm sep}$ in the three-point functions. For each choice of which values
of $t_{\rm sep}$ to use, FH${}_{\Gamma}(t_{\rm sep},{\tau}_{c}=1,dt)$ is
constructed from the set of $R_{\Gamma}(t_{\rm sep},{\tau})$. When we use
consecutive values of $t_{\rm sep}$, we take $dt=1$ in FH${}_{\Gamma}(t_{\rm
sep},{\tau}_{c}=1,dt)$. When we use only even or odd values of $t_{\rm sep}$,
we take $dt=2$.

The left panel of Fig. 4 shows $\mathring{g}_{A}$ as a function of the minimum
and maximum value of $t_{\rm sep}$ used in the analysis. When $t_{\rm
sep}^{\rm min}$ is varied, $t_{\rm sep}^{\rm max}$ is held fixed at 14. When
$t_{\rm sep}^{\rm max}$ is varied, $t_{\rm sep}^{\rm min}$ is fixed at 2. In
the middle panel of Fig. 4, we show the value of $\mathring{g}_{A}$ as a
function of $t_{\rm sep}^{\rm max}$ when only even values of $t_{\rm sep}$ are
used in the analysis of $R_{\Gamma}(t_{\rm sep},{\tau})$ and ${\rm
FH}_{\Gamma}(t_{\rm sep})$. The right panel is the same as the middle panel
but we only use the odd, rather than the even, values of $t_{\rm sep}$ in the
analysis.

When the three-point functions are computed with a sufficiently large number
of $t_{\rm sep}$ values, the ground state parameters are very stable under
data truncation and the excited state model. Of note, it is sufficient to use
every other value of $t_{\rm sep}$ at this $a\approx 0.09$ fm lattice spacing.

###  III.3  Excited State Breakdown

Given a model that is demonstrated to accurately describe the correlation
functions over the full range of $t_{\rm sep}$ and ${\tau}$ used in the
analysis, we can separate the various sources of excited state contamination
into the “scattering” ($n$-to-$n$) and “transition” ($m$-to-$n$) sources as
well as those arising from the excited states of the two-point function (see
Sec. II.1). While we can accurately describe these various sources of excited
state contamination, we do not claim to have a rigorous determination of the
spectrum, since the creation and annihilation operators we have used are
purely local three-quark operators which are known to have poor overlap with
the nucleon-pion scattering states Lang and Verduci (2013). However, since
only the transition excited states depend upon the current insertion time
${\tau}$, we can confidently separate the excited states into these various
classes of excited state contamination.

In Fig. 5, we plot the percent contamination of various sources of excited
state contamination to the ground state for $R_{A_{3}}(t_{\rm
sep},{\tau}=t_{\rm sep}/2)/\mathring{g}_{A}$ (left) and FH${}_{A_{3}}(t_{\rm
sep},{\tau}_{c}=1)/\mathring{g}_{A}$ (right). The legend keys for
$R_{A_{3}}^{\rm es}(t_{\rm sep},{\tau}=t_{\rm sep}/2)/\mathring{g}_{A}$
correspond to

2pt es: $\displaystyle\quad-\frac{C_{2}^{\rm es}(t_{\rm sep})}{C_{2}(t_{\rm
sep})}\,,$ 3pt sc: $\displaystyle\quad\phantom{-}\frac{C_{A_{3}}^{\rm
sc}(t_{\rm sep},{\tau}=t_{\rm sep}/2)}{\mathring{g}_{A}C_{2}(t_{\rm sep})}\,,$
3pt tr: $\displaystyle\quad\phantom{-}\frac{C_{A_{3}}^{\rm tr}(t_{\rm
sep},{\tau}=t_{\rm sep}/2)}{\mathring{g}_{A}C_{2}(t_{\rm sep})}\,.$ (17)

The 2pt + sc is the sum of the first two terms and the 3pt is the sum of all
excited states. The FH${}_{A_{3}}(t_{\rm sep},{\tau}_{c}=1)/\mathring{g}_{A}$
keys correspond to the same breakdown of various excited state contributions
after they are passed through Eq. (II.2). From these plots, there are several
observations and conclusions one can make:

  1. 1.

There is a significant cancellation between the scattering excited states and
two-point excited states;

  2. 2.

The transition excited states are the dominant source of excited states and
they are relatively suppressed in ${\rm FH}_{A_{3}}(t_{\rm sep},{\tau}_{c})$
as compared to $R_{A_{3}}(t_{\rm sep},{\tau}=t_{\rm sep}/2)$;

  3. 3.

The total excited state contamination at $t_{\rm sep}\approx 1$ fm is $\approx
8\%$ for $R_{A_{3}}(t_{\rm sep},{\tau}=t_{\rm sep}/2)$ and $\approx 1\%$ for
${\rm FH}_{A_{3}}(t_{\rm sep},{\tau}_{c})$.

Figure 5:  Ratio of excited state contributions to the ground state, $g_{A}$,
for $R_{A_{3}}(t_{\rm sep},{\tau})$ (left, 3pt) and ${\rm FH}_{A_{3}}(t_{\rm
sep},{\tau}_{c}=1)$ (right, fh). Each panel contains three types of data
points: the $n$-to-$n$ scattering (sc) excited states, $C_{\Gamma}^{\rm
sc}/(g_{A}C_{2})$, the $n$-to-$m$ transition (tr) excited states,
$C_{\Gamma}^{\rm tr}/(g_{A}C_{2})$, and the excited states from the two-point
(2pt es) function, $-C_{2}^{\rm es}/C_{2}$ (and the corresponding Feynman-
Hellmann versions of these ratios in the right figure). The decomposition of
the numerical data into these three classes of excited state contamination
depends upon our posterior reconstruction of the correlation functions (2pt,
3pt and FH). Because we are looking down the middle of $R_{A_{3}}$ (along
$\tau=t_{\rm sep}/2$), we only show data at even $t_{\rm sep}$s in the left
panel. In addition to plotting the corresponding posterior contribution from
each class of excited state, we plot two hatched bands. The white hatched band
is the sum of the 3pt sc and 2pt es contributions, which are observed to
largely cancel, as indicated in Eq. (II.1). The gray hatched band is the sum
of all excited state contributions (the full posterior distribution, minus the
ground state, normalized by the ground state). The 3pt sc contribution is
positive for $R_{A_{3}}(t_{\rm sep},{\tau})$ and it becomes negative for ${\rm
FH}_{A_{3}}(t_{\rm sep},{\tau}_{c}=1)$ at $t_{\rm sep}>0.5$ fm. The opposite
is true for the 2pt es contribution. The sum of the scattering and two-point
excited state contributions is opposite in sign to the transition excited
state contributions, unlike for $R_{A_{3}}(t_{\rm sep},{\tau})$ in which they
are the same sign, leading to an even stronger suppression of the total
excited state contributions beyond the expected suppression, as described in
the text. The horizontal red dashed lines represent the threshold of a 1%
contribution of excited state contamination to the ground state value.

In Eq. (II.1) we see, from the proportionality of the $n^{\rm th}$ excited
state to $g_{nn}^{\Gamma}-g_{00}^{\Gamma}$, an explicit cancellation between
the scattering (3pt sc) and two-point excited state (2pt es) contamination.
For this cancellation to be significant, as observed in the posterior
determination of the two classes of excited states, one explanation is that
$g_{nn}^{\Gamma}\approx g_{00}^{\Gamma}$ for all $n$. For the vector current,
we know that in infinite volume, the vector operator measures the iso-vector
charge of the system, which in the isospin limit, is exactly equal to
$g_{00}^{V_{4}}$. It has further been demonstrated that in the forward
scattering limit, the finite volume interaction amplitudes conspire with the
finite volume Lellouch-Lüscher factors to ensure that
$g_{nn}^{\Gamma}=g_{00}^{\Gamma}$ for two-particle states Briceño  _et al._
(2019).

This result, while not surprising, demonstrates charge conservation in finite
volume by allowing one to relate $n$-point Greens functions to $(n-1)$-point
Greens functions through the Ward-Takahashi equations. In the case of the
axial-vector matrix element, it is plausible that in the forward scattering
limit, the corrections to such a relation arising from the partially-conserved
nature of the axial-vector current would also vanish, such that in this limit,
the Lellouch-Lüscher factors are similarly removed. For example, for the
nucleon, we know the induced pseudo-scalar form factor vanishes in the forward
limit. We leave a detailed study of this question to future work.

In the left panel of Fig. 5 we see that the scattering excited state
contribution (the upper (red) band/data) and the contribution from the excited
states coming from the two-point function (the middle (purple) band/data) are
roughly equal and opposite in sign. The hatched curve, which is mostly white
with gray hash lines, is the sum of these two contributions and it lies
between them. We observe that the sum of these two sources of excited state
contributions decay to a 1% correction (the horizontal dashed (red) lines) at
$t_{\rm sep}\approx 1$ fm. In contrast, the transition (3pt tr) excited states
depicted by the lowest (blue) band/data do not decay to the 1% level until
$t_{\rm sep}\approx 2.2$ fm.

In the right panel of Fig. 5 we show the same excited state breakdown for the
Feynman-Hellmann correlation function. The sign of the scattering and two-
point excited state contributions both change compared to $R_{A_{3}}(t_{\rm
sep},{\tau}=t_{\rm sep}/2)$. This sign change can be understood from Eq.
(II.2). It is also interesting to note that the magnitude of the transition
excited state contributions at $t_{\rm sep}\approx 1$ fm go from $\approx
7.5\%$ for $R_{A_{3}}(t_{\rm sep},{\tau}=t_{\rm sep}/2)$ to $\approx 2.5\%$
for FH${}_{A_{3}}(t_{\rm sep},{\tau}_{c}=1)$. Finally, for $R_{A_{3}}(t_{\rm
sep},{\tau}=t_{\rm sep}/2)$, the sum of the scattering and two-point excited
state contributions are equal in sign to the transition excited states, while
for FH${}_{A_{3}}(t_{\rm sep},{\tau}_{c}=1)$ they are opposite in sign at
intermediate and large $t_{\rm sep}$. Thus, the total excited state
contamination decays to the 1% level for the Feynman-Hellmann correlation
function at $t_{\rm sep}\approx 1$ fm while this does not happen until $t_{\rm
sep}>2$ fm for $R_{A_{3}}(t_{\rm sep},{\tau}=t_{\rm sep}/2)$.

There is no proof that this fortunate cancellation must happen. However, an
examination of our results in Ref. Chang  _et al._ (2018) shows a consistent
picture that the excited state contamination of the FH${}_{A_{3}}(t_{\rm
sep},{\tau}_{c}={\rm none})$ correlation function decays to the 1% level by
roughly 1 fm over a broad range of pion masses, $130\lesssim m_{\pi}\lesssim
400$ MeV.

###  III.4  Comparison with late-time only results

In this section, we perform the correlator analysis on late time separation
data, $t_{\rm sep}=[10,12,14]$, which is directly comparable to the results in
Refs. Bhattacharya  _et al._ (2016); Gupta  _et al._ (2018) that were
performed on the same ensembles, and otherwise mimicking the more common
strategy used for example in Refs. Capitani  _et al._ (2012); Horsley  _et
al._ (2014); Bali  _et al._ (2015); Abdel-Rehim  _et al._ (2015); Alexandrou
_et al._ (2017); Capitani  _et al._ (2019); Park  _et al._ (2021).

The results are obtained with a simultaneous fit to the two- and three-point
functions using the same priors as our main fit, given in App. B.1. We also
demonstrate the sensitivity of the ground-state parameters posterior
distributions to changes in the input priors as in App. B.2. Because we have
only three values of $t_{\rm sep}$, we use a minimal number of excited states:
with one excited state, there is only one degree of freedom in describing the
$t_{\rm sep}$ dependence of $R_{\Gamma}(t_{\rm sep},{\tau})$, after which, the
analysis relies upon the priors for the excited state matrix elements.

Figure 6:  We plot the numerical data for $R_{A_{3}}(t_{\rm sep},{\tau})$ and
resulting posterior fit bands obtained from a combined two-point and three-
point analysis when only $t_{\rm sep}=[10,12,14]$ are used in the analysis.
The inner (colored) numerical data are used in the analysis. There is a break
in the fit band and the outer (gray) band and data indicate regions of
${\tau}$ not used in the analysis. The horizontal gray band is the resulting
posterior value of the ground state, $\mathring{g}_{A}$.

With the quark smearing we have used we find that we must apply a relatively
aggressive truncation on the current insertion time in order for the model to
describe the numerical data, restricting the analysis to the “center” of the
current insertion time ${\tau}$. We illustrate the region of ${\tau}$ used for
each $t_{\rm sep}$ in Fig. 6 and depict the stability of the ground state
axial matrix element in Fig. 7. The optimal result is chosen as the fit with
the largest amount of numerical data while maintaining a good fit quality
(cfr. the lower panel of Fig. 7). We refer to the optimal choice of data
included for this 2-state model fit as $\tau_{\rm inc}=\tau_{\rm inc}^{\rm
opt}$, but we also show results for the ground state posteriors when one more
($\tau_{\rm inc}^{\rm opt}+1$) or one less ($\tau_{\rm inc}^{\rm opt}-1$)
value of ${\tau}$ is included in the fit. While not depicted, the ground state
posteriors are also stable under variation of the choice of $t_{\rm sep}^{\rm
min}$ in the two-point function.

Figure 7:  The stability of $\mathring{g}_{A}$ versus the number of states and
current insertion times included in the analysis. $\tau_{\rm inc}^{\rm opt}$
denotes our optimal choice of included data. The bottom panel shows the
corresponding fit quality $Q$. The dark horizontal band indicates the optimal
fit from the late-$t_{\rm sep}$ only analysis while the light horizontal band
above it represent the optimal fit from the full analysis presented in Sec.
III.1.

We observe that for this specific dataset, the late-$t_{\rm sep}$ result is in
$\approx 2{\sigma}$ tension with the result from the full analysis that
includes more values of $t_{\rm sep}$, and is clearly stable over a variety of
data truncations and model variations. This is indicative that late-$t_{\rm
sep}$ data can be subject to correlated fluctuations which are difficult to
identify without having results at many values of $t_{\rm sep}$.

The same systematic effect can be seen in a simpler case using a fit to the
two-point correlation function. In Fig. 8, we show the result of the ground
state energy extracted from the two-point correlation function under varying
$t_{\rm sep}^{\rm min}$. We highlight the best fit ground-state energy, which
is supported by a robust plateau in the dimension of $t_{\rm sep}^{\rm min}$
with high $Q$-value. However, we see that at late time beyond 1 fm, a second
plateau develops approximately two standard deviations above our best fit.
This second plateau cannot be physical and simply reached later in time due to
latent excited state contamination, as the 2-pt correlator is positive
definite and therefore required to approach the ground state from above. This
second “stable” plateau is also the logical result for a one-state fit.
Analogous to the three-point analysis, if the model fails to describe excited-
state contaminations, the model must extract ground-state parameters at late
time. These values are then sensitive to uncontrolled statistical
fluctuations. Also note that we observe the fluctuation in the effective mass
occurs approximately 0.2 fm later than the matrix elements, see the top panel
of Fig. 11 in App. B. Similar to the three-point analysis, in the absence of a
more holistic view as granted by analyzing more source-sink separation data,
there is no measure on the size of how large the potential underestimation of
errors are.

Figure 8:  The extracted ground state energy as a function of $t_{\rm
sep}^{\rm min}$ and the number of states ($n_{s}$) used in the analysis. The
late-$t_{\rm sep}$ data is subject to a correlated fluctuation, as observed in
the $2{\sigma}$ increase in $E_{0}$ around $t_{\rm sep}^{\rm min}=12$.

Without using more than three values of $t_{\rm sep}$, it is not possible to
identify if one is susceptible to such a fluctuation, and therefore, it is not
possible to fully quantify the uncertainty on the posterior distribution of
the ground state parameters.

###  III.5 Comparison with the 2-state FH analysis

In this final section, we compare our result with that from Ref. Chang  _et
al._ (2018) on the same a09m310 ensemble. Those results were obtained using
the Feynman-Hellmann method as described in Ref. Bouchard  _et al._ (2017),
which is the same as Eq. (II.2) with the sum over current insertion time
running over the full time extent, including contributions from contact
operators at ${\tau}\in\\{0,t_{\rm sep}\\}$ and from out-of-time region,
${\tau}<0$ and/or ${\tau}>t_{\rm sep}$. We applied a two-state model and a
frequentist analysis on the data from the two-point and Feynman-Hellmann
correlation functions.

Figure 9:  Comparison of $\mathring{g}_{A}$ determined from the three analysis
strategies described in the text with our result from Ref. Chang  _et al._
(2018) which utilized the Feynman-Hellmann correlator with a sum over all
current insertion times Bouchard  _et al._ (2017) (note that the quark
smearing in the present analysis is different from that used in Ref. Chang
_et al._ (2018)).

In Fig. 9 we show the extracted value of $\mathring{g}_{A}$ with three
different fully-correlated analysis strategies:

  1. 23s

: Fit to $C_{2}(t_{\rm sep})$, $R_{\Gamma}(t_{\rm sep},{\tau})$ and
FH${}_{\Gamma}(t_{\rm sep},{\tau}_{c}=1)$;

  2. 2s

: Fit to $C_{2}(t_{\rm sep})$ and FH${}_{\Gamma}(t_{\rm sep},{\tau}_{c}=1)$;

  3. 23

: Fit to $C_{2}(t_{\rm sep})$ and $R_{\Gamma}(t_{\rm sep},{\tau})$.

The horizontal gray band is the result from Ref. Chang  _et al._ (2018). These
results show the consistency between the four computational and analysis
strategies. Even though the FH${}_{\Gamma}(t_{\rm sep},{\tau}_{c})$
correlation functions are constructed from the $R_{\Gamma}(t_{\rm
sep},{\tau})$ correlation functions, they are subject to a different pattern
of excited state contamination (cfr. Sec. II). Therefore it is useful to
include them in the fully-correlated analysis. The combined fit yields a
robust and consistent fitting strategy and precise extractions of hadronic
matrix elements at zero-momentum transfer. Whether an analogous strategy can
be pursued for non-zero momentum transfer is left for future work.

In Fig. 10, we present the posterior correlation function on top of the
numerical data for both the results in Ref. Chang  _et al._ (2018) as well as
from the FH${}_{A_{3}}(t_{\rm sep},{\tau}_{c}=1)$ dataset in the present work.
The FH${}_{A_{3}}(t_{\rm sep},{\tau}_{c}={\rm none})$ dataset from Ref. Chang
_et al._ (2018) is observed to have significantly less excited state
contamination at early $t_{\rm sep}$ as compared to FH${}_{A_{3}}(t_{\rm
sep},{\tau}_{c}=1)$. This is what enabled a two-state analysis in Ref. Chang
_et al._ (2018). The combined fit (23s) in the present work enables a
determination of significantly more excited state parameters through the
precise early-$t_{\rm sep}$ data. Moreover, as shown in this section, a
confident extraction of the ground state posterior distributions is attained:
as depicted in the figure, both strategies yield precise and consistent
extractions of the large-$t_{\rm sep}$ extrapolation of the results. The
former result is more economical for obtaining a precise value of
$\mathring{g}_{A}$, while the method in this work also enables a determination
of non-zero momentum transfer results which can be used to determine the form
factors.

Figure 10:  Comparison of our FH${}_{A_{3}}(t_{\rm sep},{\tau}_{c}={\rm
none})$ numerical results and analysis from Ref. Chang  _et al._ (2018) using
a two-state frequentist analysis (top data/band) to the FH${}_{A_{3}}(t_{\rm
sep},{\tau}_{c}=1)$ numerical data and posterior description of the
correlation function from the fully-correlated five-state Bayesian analysis to
$C_{2}(t_{\rm sep})$, FH${}_{{\Gamma}}(t_{\rm sep},{\tau}_{c}=1)$ and
$R_{\Gamma}(t_{\rm sep},{\tau})$ in the present work (bottom data/band).

##  IV  Observations and Conclusions

In this work, we have computed the three-point correlation functions that are
used to determine the nucleon axial charge for 13 values of the source-sink
separation time in the range $t_{\rm sep}\approx 0.17-1.22$ fm on an ensemble
with $a\approx 0.09$ fm and $m_{\pi}\approx 310$ MeV. This large numerical
dataset (Fig. 1) has enabled us to robustly determine the ground state mass
and matrix elements with a fully-quantified systematic uncertainty arising
from excited states. It further enabled us to quantify the excited state
contribution to the correlation function with greater detail than has been
previously achieved in the literature. We were able to demonstrate that the
ground state parameters are stable against the model of excited states (Eq.
(III.2) and Fig. 3), the number of excited states used in the analysis, and
truncations of either small $t_{\rm sep}$ or large $t_{\rm sep}$ results as
well as through the use of only even or odd values of $t_{\rm sep}$ results
(Fig. 4).

A re-writing of the spectral decomposition of the three-point correlation
functions revealed a prospective cancellation between two classes of excited
states arising from the $n$-to-$n$ scattering states and those arising from
the two-point correlation function used to construct $R_{\Gamma}(t_{\rm
sep},{\tau})$, Eq. (II.1) (there is an exact cancellation for the vector
matrix element for which $g_{nn}^{V_{4}}=g_{00}^{V_{4}}=1$ up to
renormalization). For our calculation, this cancellation seems to materialize
for $g_{A}$ (Fig. 5) leaving the dominant excited state contributions to be
the $n$-to-$m$ transition terms. In the standard fixed source-sink separation
method, the excited state contributions do not decay down to the 1%
contamination level until $t_{\rm sep}\gtrsim 2$ fm while in the Feynman-
Hellmann correlation function, Eq. (II.2), they decay to the 1% level at
$t_{\rm sep}\approx 1$ fm. We further observed that the Feynman-Hellmann
correlation function constructed with a sum over the entire time extent of the
lattice Bouchard  _et al._ (2017) (FH${}_{\Gamma}(t_{\rm sep},{\tau}_{c}={\rm
none})$), rather than just between the source and the sink
(FH${}_{\Gamma}(t_{\rm sep},{\tau}_{c}=1)$), leads to an even further
suppression of excited state contamination (Fig. 10). This stronger
cancellation of excited states enabled us to compute $g_{A}$ with a $\approx
1\%$ uncertainty by utilizing early-time data Berkowitz  _et al._ (2017a);
Chang  _et al._ (2018); Berkowitz  _et al._ (2018a); Walker-Loud  _et al._
(2020), which demonstrated this higher suppression of excited states over a
broad range of pion masses, $130\lesssim m_{\pi}\lesssim 400$ MeV.

By utilizing a large number of early to mid-time data, one is able to detect
if the late-time data is subject to a correlated fluctuation which might
otherwise bias the ground state matrix elements (Figs. 7 and 8). With only
three values of $t_{\rm sep}$ at late time, it is not possible to perform such
a data truncation study that could identify this issue. We have found that
including the Feynman-Hellmann correlation functions in a global fully
correlated analysis, even though they are constructed from the three-point
functions as in Eq. (II.2) and are highly correlated with $R_{\Gamma}(t_{\rm
sep},{\tau})$, improves the stability and precision of the extracted ground
state parameters (Fig. 9). This improvement arises because the excited states
present themselves differently in the two sets of matrix-element correlation
functions.

Even though the Feynman-Hellmann correlation function has greater noise than
the standard three-point correlator at equal $t_{\rm sep}$ and statistics
(compare, for example, $t_{\rm sep}=1$ fm. in Fig. 1), the more rapid decay of
excited states demonstrated in Sec. III.3 means that this strategy mitigates
the exponential degradation of the signal-to-noise of the nucleon’s two- and
three-point correlation functions as $t_{\rm sep}$ is increased, which
requires exponentially more computational resources to control the stochastic
precision. The strategy we present in this work offers a more economical
method of obtaining the ground state matrix elements than that which is more
commonly advocated for in the literature, which is to use high-statistics
calculations at $t_{\rm sep}\approx 1-2$ fm or larger. In future work, we will
investigate the same strategy for non-zero momentum transfer correlation
functions which are used to determine the nucleon form factors.

Data availability: The computations were performed utilizing LaLiBe lal  which
utilizes the Chroma software suite Edwards and Joo (2005) with QUDA solvers
Clark  _et al._ (2010); Babich  _et al._ (2011) and HDF5 The HDF Group
(1997-NNNN) for I/O Kurth  _et al._ (2015). They were efficiently managed with
METAQ Berkowitz (2017); Berkowitz  _et al._ (2018b) and status of tasks logged
with EspressoDB Chang  _et al._ (2020). The final extrapolation analysis
utilized gvar Lepage (2020) and lsqfit Lepage (2021). The analysis and data
for this work can be found at https://github.com/callat-qcd/project_fh_vs_3pt.

###### Acknowledgements.

We thank O. Bär, R. Briceño, M. Hansen and A. Jackura for enlightening
conversations and correspondence. We thank B. Hörz for help cross-checking the
three-point code in LaLiBe and generating some of the results. We thank the
MILC Collaboration for use of their HISQ gauge ensembles. Computing time for
this work was provided through the Innovative and Novel Computational Impact
on Theory and Experiment (INCITE) program and the LLNL Multiprogrammatic and
Institutional Computing program for Grand Challenge allocations on the LLNL
supercomputers. This research utilized the NVIDIA GPU-accelerated Summit
supercomputer at Oak Ridge Leadership Computing Facility at the Oak Ridge
National Laboratory, which is supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC05-00OR22725 as well as the
Lassen supercomputer at Lawrence Livermore National Laboratory. This work was
supported in part by the Berkeley Physics International Education Program
(JH), the Berkeley Laboratory Undergraduate Research Program (IC); the NVIDIA
Corporation (MAC), the Alexander von Humboldt Foundation through a Feodor
Lynen Research Fellowship (CK), the U.S. Department of Energy, Office of
Science, Office of Nuclear Physics under Award Numbers DE-AC02-05CH11231 (CCC,
CK, AWL), DE-AC52-07NA27344 (DAB, DH, ASG, PV), DE-FG02-93ER-40762 (EB), DE-
AC05-06OR23177 (CJM); DE-SC00046548 (ASM); the DOE Early Career Award Program
(AWL), and the U.K. Science and Technology Facilities Council grants
ST/S005781/1 and ST/T000945/1 (CB);

##  Appendix A Discrete symmetry systematics

“Coherent sink technique” Bratt  _et al._ (2010): We reduce the numerical cost
of the computations by solving for a single sequential-propagator from many
sequential sinks simultaneously. We found that we can combine eight sinks into
a single coherent sink generated with two sources per timeslice, with a 10-20%
loss in statistical precision as compared to solving a single sequential-
propagator for each of the eight sources separately. For each $t_{0}$, a
random origin (O) is chosen and then the antipode (A) location is also chosen
Berkowitz  _et al._ (2018c)

$\displaystyle s_{\rm O}(t_{0})$ $\displaystyle=(x_{0},y_{0},z_{0})$
$\displaystyle s_{\rm A}(t_{0})$ $\displaystyle=\left[s_{\rm
O}(t_{0})+\frac{L}{2}(1,1,1)\right]\ {\rm mod}\ L\,.$

We repeat this for four values of $t_{0}$ spaced by $T/4$. We generate all 16
sources by running a second set of eight sources shifted by $T/8$ from the
first set of sources.

Spin-averaging: We find that combining the spin-up–to–spin-up and spin-
down–to–spin-down correlation functions (with a $+$ sign for $V_{4}$ and a $-$
sign for $A_{3}$) leads to a near perfect $\sqrt{2}$ reduction in the
stochasitc uncertainty of the numerical data. We further observe that the non-
symmetric behavior of $R_{\Gamma}(t_{\rm sep},{\tau})$ about ${\tau}=t_{\rm
sep}/2$ for larger values of $t_{\rm sep}$, which must vanish in the infinite
statistics limit, is less pronounced when we perform this spin-averaging.

Time-reversal symmetry: We find that combining the backwards temporal
propagation of the negative-parity two- and three-point functions – the time-
reversed correlation functions (negative parity three-point functions with
negative values of the source-sink separation time $t_{\rm sep}=t_{\rm
snk}-t_{\rm src}<0$) with the positive parity three-point functions generated
with $t_{\rm sep}>0$ leads to a near-perfect $\sqrt{2}$ reduction in
stochastic uncertainty, allowing us to make use of both the positive and
negative parity components of the quark propagators.

##  Appendix B Analysis Details

In this appendix, we discuss in detail the analysis of the various correlation
functions and what led us to our final strategy presented in Sec. III.

To analyze the two-point (2pt), Eq. (7a), three-point (3pt), Eq. (5) and
Feynman-Hellmann (FH), Eq. (II.2) correlation functions and determine the
parameters of the fit model, we perform a maximum-likelihood Bayesian
analysis. We explore analyzing three combinations of correlation functions in
a global (simultaneous) analysis:

  1. 23s

: $C_{2}(t_{\rm sep})$, $R_{\Gamma}(t_{\rm sep},{\tau})$ and
FH${}_{\Gamma}(t_{\rm sep},{\tau}_{c}=1)$;

  2. 2s

: $C_{2}(t_{\rm sep})$ and FH${}_{\Gamma}(t_{\rm sep},{\tau}_{c}=1)$;

  3. 23

: $C_{2}(t_{\rm sep})$ and $R_{\Gamma}(t_{\rm sep},{\tau})$.

###  B.1 Prior selection

The first step in the analysis is to choose prior distributions for the
parameters. In order to estimate the ground state priors, we use the effective
mass, the effective overlap and an effective $g_{A}$ plot:777We do not show
the effective $g_{V}$ plot since $g_{V}$ is very close to 1 with our lattice
action.

$\displaystyle m_{\rm eff}(t_{\rm sep})$
$\displaystyle=\ln\left(\frac{C_{2}(t_{\rm sep})}{C_{2}(t_{\rm
sep}+1)}\right)\,,$ $\displaystyle z^{2}_{\rm eff}(t_{\rm sep})$
$\displaystyle=e^{m_{\rm eff}(t_{\rm sep})t_{\rm sep}}C_{2}(t_{\rm sep})\,,$
$\displaystyle g_{A}^{\rm eff}$ $\displaystyle={\rm FH}_{A_{3}}(t_{\rm
sep},{\tau}_{c}=1)\,.$ (18)

In Fig. 11, we plot these effective quantities which all asymptote to their
ground state values in the large $t_{\rm sep}$ limit. We choose conservative
ground state priors to be

$\displaystyle\tilde{E}_{0}=0.50(5)\,,$ $\displaystyle\tilde{g}_{A}=1.2(2)\,,$
$\displaystyle\tilde{z}_{0}=0.00034(34)\,,$
$\displaystyle\tilde{g}_{V}=1.0(2)\,,$ (19)

which are plotted as the wide gray horizontal bands. We also plot the
resulting posterior distribution of the effective quantities resulting from
the parameters of our final analysis.

Figure 11:  The effective mass (Top), overlap factor $z$ (Middle) and $g_{A}$
correlation functions (Bottom). The wide horizontal gray bands are the chosen
ground state priors, while the blue overlay bands are the reconstructed
effective quantities using the posterior distributions from the final
analysis.

For the excited state energies we explore three models of excited states, Eq.
(III.2). As shown in Fig. 3, the posterior energies are largely insensitive to
the model. We therefore focus the discussion on the spectrum of our chosen
model, ${\Delta}E_{n}=2m_{\pi}$ for all $n$. The lowest-lying excitation is a
nucleon-pion P-wave or a nucleon with two-pions at rest, up to interaction
energies which are a small fraction of the total energy. For our $m_{\pi}L$,
these two energy levels are practically degenerate and therefore modeled as a
single excitation. We prior all the ${\Delta}E_{n}$ with a log-normal
distribution, $\ln({\Delta}E_{n})=(\ln(2m_{\pi}),0.5)$ such that the resulting
energies, $E_{n}=E_{0}+\sum_{l=1}^{n}{\Delta}E_{l}$ are ordered.

While the creation and annihilation operators are conjugate to each other,
this does not fix the absolute sign of $z_{n}$. Further, there is a redundancy
in sign of the transition matrix elements: if all the overlap factors are
taken to be positive, $z_{n}>0$, a negative contribution will manifest in a
negative value of $g_{nm}$. Only the combination $z_{n}z^{\dagger}_{m}g_{nm}$
has a well defined sign.

To be conservative, we prior the central values of the excited state overlap
factors with a central value of $0$. For the first excited state, we choose a
slightly smaller prior width with respect to the ground state and for the
higher excited states we use again a slightly smaller width:

$\displaystyle\tilde{z}_{1}=0(0.00025)\,,$ $\displaystyle\tilde{z}_{n\geq
2}=0(0.00015)\,.$ (20)

These slight reductions are motivated by the use of a smeared quark source,
which suppresses excited state overlap factors compared to the ground state.

For the vector matrix elements, the conserved charge protects the charge of
all states to be $g^{V}_{nn}=1$, even in finite volume Briceño  _et al._
(2019). For the transition matrix elements, we postulate that these are the
same order of magnitude, but with an unknown sign. For the axial-vector matrix
elements, we postulate the excited state matrix elements and transition matrix
elements are of the same order of magnitude as $g_{A}$. This leads us to the
prior values of

$\displaystyle g^{V}_{nn}=1.0(2)\,,$ $\displaystyle g^{A}_{nn}=1.2(2)\,,$
$\displaystyle g^{V}_{nm}=0(1)\,,$ $\displaystyle g^{A}_{nm}=0(1)\,.$ (21)

A complete list of the prior and posterior distributions of all fit parameters
is provided in App. C.

When performing a multi-exponential fit, it is expected that the highest state
used in the analysis serves as a “garbage” can that is contaminated by the
tower of more highly-excited states not included in the analysis. The 3pt and
FH correlation functions have different parametric dependence upon the excited
states. Therefore, when exploring the parameter space of fits, if the number
of states used for example in the FH correlation functions differs from the
2pt function, we allow the highest lying state in each correlation function to
have different priors. Specifically, if the 2pt uses five states and the FH
uses three states, then the energy and overlap priors of the third FH state
are decoupled from the third state of the 2pt function. When the number of
states used is the same, we find we are able to describe the correlation
functions well when keeping the highest garbage-can state the same in all
correlation functions.

###  B.2 Sensitivity analysis

For each set of correlation functions, the best fit is chosen after a careful
study of the posterior distribution sensitivity on input fit parameters
including fit ranges, number of states in the fit model, prior widths, and
model dependence of the excited state spectrum. In the following sections, we
discuss the best fit under the context of three different fit strategies, 23,
2s and 23s, and then the costs and benefits of these different strategies.

####  B.2.1 Fit region and $n_{\textrm{state}}$ stability analysis

Due to the structure of excited state contamination, the posterior
distributions are most sensitive to changes in $t_{\rm sep}^{\rm min}$, the
minimum source-sink time, and the number of excited states
$n_{\textrm{states}}$ ($n_{s}$) included in the model. In this section we
discuss the stability of the best fit under changing these two dimensions.

Two-point with Feynman-Hellmann analysis (2s)

Figure 12:  Stability analysis for the axial coupling for the combined two-
point and Feynman-Hellmann fit for (Top) varying $t_{\rm sep}^{\rm min}$ of
the two-point function only and (Bottom) Feynman-Hellmann correlation function
only. In both cases, results of $n_{s}=[1,2,3]$ state fit ansatzes are shown.
The best fit for a given $n_{s}$ is highlighted by a solid marker and
corresponding horizontal band. The filled horizontal red band for the
$n_{s}=2$ result highlights the best fit out of the entire explored parameter
space. The corresponding $Q$-values are also provided.

Fig. 12 shows the dependence on the axial-vector matrix element fit parameter,
$\mathring{g}_{A}$, under varying two-point and FH $t_{\rm sep}^{\rm min}$,
and their corresponding $n_{\textrm{states}}$. Typically, stability plots
demonstrate that the best fit lies at a locally optimal point in the parameter
space. The simplest strategy is to simply fix all but one parameter, such as
$t_{\rm sep}^{\rm min}$ or $n_{\textrm{states}}$ in this section. If results
surrounding the best fit are absent of spurious correlations, then the
stability plots provide evidence that systematic errors arising from the fit
procedure are all accounted for. However, in the case of the FH form of the
correlation function, the time dependence observed is analogous to a two-point
correlation function. And similar to the two-point function, the best fit for
a given number of excited states, requires a more careful choice of $t_{\rm
sep}^{\rm min}$. For example, a two-point fit with many excited states may
successfully describe the data at small source-sink separation times, but will
surely fail when the model is simplified. Therefore, for the stability plots
shown in Fig. 12, the analysis varies $t_{\rm sep}^{\rm min}$ around the best
fit for a given number of states in the model. For simplicity, we vary the
$t_{\rm sep}^{\rm min}$ of the vector and axial-vector matrix elements in
tandem.

After the summation, the three-point correlation function only has one degree
of freedom per time slice. Coupled with the fact that the signal-to-noise is
exponentially worse compared to the two-point function (see Fig. 11), this
means that in order to prevent over-fitting the data, simpler fit models with
less excited states should be chosen. Therefore, the stability shown in Fig.
12 is checked for fits with 1, 2, and 3 states. Beyond three states, the model
will have comparable or more parameters than data. While this can in principle
be alleviated with careful choices of priors, strategies that rely more
heavily on prior information may also inadvertently introduce possible
systematic errors.

For a single state, we observe the fit to be stable at approximately a source-
sink separation time of 1 fm and identifies the length scale in which all
excited-states have decayed to below the noise given a percent-level
determination of the matrix element. This observation is consistent with what
is observed in our previous work Bouchard  _et al._ (2017); Chang  _et al._
(2018). With two- and three-state fits, we observe that the best fit $t_{\rm
sep}^{\rm min}$ can encompass progressively shorter source-sink separation
times, in agreement with expectation. Furthermore, the best fit for the three
models, and neighboring results, are all consistent within one standard
deviation with no appreciable systematic trend.

Figure 13:  Stability analysis for the axial coupling for the combined two-
point and Feynman-Hellmann fit for varying $t_{\rm sep}^{\rm max}$ of the two-
point function. The best fit highlighted by a solid marker and corresponding
horizontal band is identical to the one shown in Fig. 12. The corresponding
$Q$-values are provided.

A two-state fit allows us to capture some of the excited state dependence in
order to not rely entirely on the correlator reaching the plateau region.
Simultaneously, fitting to two states avoids the possibility of over-fitting
since an introduction of three additional fit parameters (ground state and
first excited state matrix elements) can be extracted from 10 data points in
this specific analysis ($t_{\rm sep}^{\rm min}$ of 5 to 14). We assume that
the overlap and energies are well constrained by the two-point correlation
function in this counting. The three-state fit introduces 6 new matrix element
parameters and is at the edge of what is naively allowed by data. Fitting to
three states allows $t_{\rm sep}^{\rm min}$ to encompass down to $t_{\rm
sep}=4$, resulting in 12 data points. However, since lattice correlation
functions have an exponential signal-to-noise problem, data points do not
carry equal weight in determining the posterior distributions. In particular,
the weighted loss function penalizes larger time separation data with the
inverse of the variance. Referring back to Fig. 11, we see for this specific
example, data beyond $t_{\rm sep}=9$ have little impact on the outcome of the
result. This can be seen for example in Fig. 13 where the results are given
for varying $t_{\rm sep}^{\rm max}$ of the two-point function. As a result,
the three-state fit with 6 extra parameters is effectively constrained by
approximately 6 data points. It follows that more complicated fit functions
run into the danger of over-fitting. We conclude that fits to the FH
correlation function are best performed with two states.

Two-point with three-point analysis (23)

Figure 14:  Stability analysis for $\mathring{g}_{A}$ for the combined two-
point and three-point fit for (Top) varying $t_{\rm sep}^{\rm min}$ of the
two-point function and (Bottom) three-point correlation function. In both
cases, results of $n_{s}=[4,5,6,7]$ state fit ansatzes are shown. The filled
marker and horizontal green band for the $n_{s}=5$ result highlights the best
fit. The corresponding $Q$-values are provided. The Bayes Factors $w$ for a
given $t_{\rm sep}^{\rm min}$ are shown in the bottom panel of each plot.

Rather than constructing the FH combination, we explore the possibility of
fitting directly to the three-point correlator as a function of both source-
sink and current-sink insertion time. Due to the complexity involved in
choosing a two-dimensional fit region, we simplify the decision making process
by fitting to all valid insertion-sink time separation times (the contact
terms are dropped) as shown in Fig. 15. This choice reduces the $t_{\rm
sep}^{\rm min}$ stability study to again a single dimension. To capture this
curvature, a large number of excited states are required. However, unlike the
FH correlator, the three-point correlator supplies many data points as shown
in Fig. 1, and can be used to fully capture the complicated excited-state
structure. Additionally, conservative smearing on the interpolating operator
yields a two-point correlator that provides corroborating support to a large
number of excited states when fit to small source-sink separation times. Due
to the large number of states, stability with respect to $t_{\rm sep}^{\rm
min}$ are less sensitive to $n_{\textrm{states}}$. As a result, we further
simplify the study by showing only stability with respect to the final result,
instead of first identifying $n_{\textrm{states}}$-dependent best fits shown
in Fig. 12.

Figure 15:  This figure depicts the regions in $t_{\rm sep}$ and ${\tau}$ used
in the 23 analysis. The numerical results have been averaged about
${\tau}=t_{\rm sep}/2$ denoted by the light-shaded squares on the right half
of the figure. The outer-most (red) squares for each $t_{\rm sep}$ are the
points ${\tau}=0$ or ${\tau}=t_{\rm sep}$ which are excluded from the
analysis. For $t_{\rm sep}>10$, the values of ${\tau}=1$ and ${\tau}=t_{\rm
sep}-1$ are also excluded from the optimal fit. Moving in, the (black) squares
represent the values of ${\tau}$ closest to the source/sink times included in
the optimal analysis (with all “inner” values of ${\tau}$ also included),
${\tau}_{\rm inc}^{\rm opt}$. Cutting an additional value of ${\tau}$ closest
to the source/sink times is ${\tau}_{\rm inc}^{\rm opt}-1$ denoted by the
lighter (blue) squares. The lighter (green) squares denote ${\tau}_{\rm
inc}^{\rm opt}-2$ and the innermost (yellow) squares denote ${\tau}_{\rm
inc}^{\rm opt}-3$. When cutting values of ${\tau}$, we always keep one (for
even $t_{\rm sep}$) or two (for odd $t_{\rm sep}$) values of ${\tau}$ in the
middle of the source-sink separation time.

Fig. 14 shows the stability of the best fit for the two-point and three-point
simultaneous fit under varying $t_{\rm sep}^{\rm min}$ for the two-point and
three-point function. Similar to the FH fit in the previous section, we
simplify the analysis by varying the $t_{\rm sep}^{\rm min}$ for the vector
and axial-vector matrix elements simultaneously. The best fit presumably can
be slightly improved if this condition is relaxed. Due to the choice of
fitting all current insertion times (aside from the source and sink time), the
best fit is observed to require a 5 state model, and includes source-sink
separation times that are commensurate to the inclusion of small current-sink
times. The best fit lies in the region of stability, and is chosen to
incorporate the most data possible given a fit model. Alternatively, given a
fit region, the simplest model is chosen (fewest number of states). For
example, the best fit with a three-point $t_{\rm sep}^{\rm min}=3$ includes 5
states even though 6 or 7 states yields similar ground state posterior
distributions. This decision is corroborated with the set of Bayes Factors
normalized to a fixed fit region, Eq. (16). For example, Fig. 14 suggests that
5 state and 7 state fit have comparable probability for reproducing the
underlying data.

Figure 16:  Stability analysis for the combined $C_{2}(t_{\rm sep})$ and
$R_{\Gamma}(t_{\rm sep},{\tau})$ analysis for varying regions of current-
insertion time. The best fit given by $\tau^{\textrm{opt}}_{\textrm{inc}}$ is
defined in Table 3 (see also Fig. 15 for a graphical depiction of
$\tau^{\textrm{opt}}_{\textrm{inc}}-n$). The filled marker and horizontal
green band highlights the best fit and is identical to the best fit in Fig.
14. The $Q$-values are provided.

While we make the choice of fitting data over all current insertion times, the
ground-state posterior distributions are stable when extracted from a subset
of the data. In Fig. 16, the fit region with respect to current insertion time
is symmetrically truncated, keeping however, at least one data point per
source-sink separation time (e.g. $t_{\rm sep}=3$ includes ${\tau}=[1,2]$, so
a ${\tau}_{c}=1$ will not eliminate any data since otherwise the entire
dataset for $t=3$ would be eliminated). It is also observed that while
curvature in $\tau$ is dependent predominantly on excited-state behavior,
aggressively truncating data still leads to larger statistical uncertainty in
ground state parameters since less data is being included in the analysis.

We conclude that a simultaneous fit to the two- and three-point correlators is
best performed by fitting to the maximum amount of data while choosing the
simplest model which can describe the data. In particular, we observe that the
under this strategy, the accompany two-point correlator also provides
sufficient constraints on the excited state overlap and energy parameters.

Two-point, three-point and FH Analysis (23s)

Figure 17:  Stability analysis for the combined two-point, three-point, and
Feynman-Hellmann analysis. These plots follow what is shown in Fig. 14.

Finally, we perform a simultaneous fit to the two-point correlator along with
all FH and three-point correlators. Since the FH correlator exposes different
excited-state dependence when compared to the three-point correlator, it may
be possible to extract a more robust calculation of the ground state
parameters. Unfortunately, the overall strategy of the previous two studies
are incompatible with one another. In the case of the FH analysis, the
strategy is to fit simpler models in order to avoid overfitting, while for the
three-point analysis overfitting is much less of a concern and instead a
majority of clean data is fitted with more complex models. A successful
combined fit will need to reconcile these differences.

We take a simple approach by recognizing that the excited state information
extracted from the three-point analysis can be used to constrain a more
complex FH fit function. Following this logic, the best combined fit follows
the same 5 state model as the two-point with three-point correlator fit
discussed previously, while the FH fit is now modelled by 5 states with the
intention of relying on the three-point correlator to constrain high excited
state parameters. Fig. 17 demonstrates the stability of the combined best fit
under changes in $t_{\rm sep}^{\rm min}$ of the two-point, three-point and FH
correlators. We observe that in the combined fit the ground-state parameters
are insensitive to changes in $t_{\rm sep}^{\rm min}$ for all datasets,
including the FH correlation functions. This observation corroborates the
hypothesis that the three-point correlator lends support to high excited-state
contributions which is consistent with the predicted spectral decomposition
for both correlation functions. Additionally, while the best fit $t_{\rm
sep}^{\rm min}$ for the three-point is kept the same as the three-point with
two-point fit, the best-fit FH $t_{\rm sep}^{\rm min}$ now extends down to
$t_{\rm sep}=3$ (previously, in the 2s best fit $t_{\rm sep}^{\rm min}$ is 5)
due to the inclusion of more states in the model.

Figure 18:  Stability analysis for the combined two-point, three-point, and
Feynman-Hellmann analysis for varying regions of current-insertion time. This
plot follows what is shown in Fig. 16.

Finally, sensitivity under varying current insertion time for the three-point
correlator is studied. Unlike the simpler three-point with two-point
correlator analysis, we have to drop one additional data point away from the
contact interaction at $t=11$ to $t=14$ for the ground-state posterior
distributions to be insensitive to changes in fit region. Fig. 18 shows the
varying fit region with respect to the best fit. In particular, the
${\tau}^{\rm opt}_{\rm inc}+1$ fit is to all current insertion dependence
between the source and sink, ${\tau}=[1,t_{\rm sep}-1]$ for all $t_{\rm sep}$.
The colored regions in Fig. 15 highlight the various ${\tau}_{\rm inc}$
regions with the colored boxes. The black region denotes the best fit region.
We observe that the best fit lies in a region that is insensitive to varying
subsets of the three-point correlator. Conversely, the combined fit suggests
that at large values of $t_{\rm sep}$, where the distribution of the
correlator become under-sampled, the three-point data shows signs of being
inconsistent with the FH data, leading to the instability seen in ${\tau}^{\rm
opt}_{\rm inc}+1$ fit.

We conclude that the three-point correlator provides enough information to
constrain the excited states of the FH correlator. A combined fit is
therefore, best performed by rooting the calculation in a two-point with
three-point strategy and extending the analysis to encompass as much of the FH
correlator as is describable by the truncated spectral decomposition.

####  B.2.2 Prior width analysis

Figure 19:  Sensitivity analysis of the best fit result from the combined two-
point, three-point, and Feynman-Hellmann fit under varying prior widths for
the ground-state (red) and together with the first-excited state (blue). The
$x$-axis indicates that the study performed variations from 0.1 to 10 times
the best fit prior width. The filled marker and horizontal green band
highlights the best fit and is identical to the result shown in Fig. 17.

Our overall strategy is to extract ground-state parameters from a multi-
exponential fit, subject to prior constraints. Since the objective function
now depends on prior distributions in addition to data, we check that the
posterior distributions of interest are insensitive to our prior knowledge.
The purpose of introducing priors in this context is not to supplement
additional information, but to constrain the search space of the numerical
minimization for faster convergence. Lower computation costs allow us to more
thoroughly investigate the sensitivity of ground state posterior
distributions, which in turn lends to more robust results. We enforce the
expectation that choices of priors should not yield changes in the posterior
distributions of ground-state parameters.

Fig. 19 demonstrates the robustness of our quoted results under variations of
the ground-state and first-excited state prior widths for the combined fit.
The study indicates that the extracted matrix elements are unconstrained by
prior distributions until the widths are reduced by a factor of 10, while
broadening the prior distribution by a factor of 10 leaves the matrix elements
unchanged. Similar conclusions hold for the vector matrix element for the 2s,
and 23 strategies.

##  Appendix C Prior and posterior distributions

In Table 1, we list the prior and posterior distributions of the energies and
overlap factors used in our final analysis. The matrix element prior and
posterior distributions are given in Table 2. We use the $(\mu,{\sigma})$
format to report the central value $\mu$ and the width ${\sigma}$ of the
distributions. In Table 3, we list the range of $t_{\rm sep}$ values used in
the three sets of correlation functions that are analyzed.

Table 1:  We list the prior and posterior energies and overlap factors, in lattice units, for each of the parameters that describe the two-point correlation function. The prior and posterior values are listed in $(\mu,{\sigma})$ pairs. The energy splittings are priored with a log-normal distribution with ${\Delta}E_{n}=(2m_{\pi},m_{\pi})$, except for the highest state which has a prior width of $5m_{\pi}$, with the total energy given by the sum in Eq. (III.2). For the optimal fit, all energy splittings are given a central value of $2m_{\pi}$.  Parameter | Prior | Posterior
---|---|---
$E_{0}$ | $(0.50,0.05)$ | $(0.4904,0.0016)$
${\Delta}E_{n<4}$ | $(0.29,0.14)$ | $(0.3,0.028)$
$E_{1}$ | $(0.79,0.15)$ | $(0.79,0.028)$
$E_{2}$ | $(1.07,0.21)$ | $(1.136,0.099)$
$E_{3}$ | $(1.36,0.25)$ | $(1.66,0.16)$
${\Delta}E_{4}$ | $(0.29,0.72)$ | $(1.119,0.033)$
$E_{4}$ | $(1.65,0.76)$ | $(2.78,0.16)$
$z_{0}$ | $(0.00034,0.00034)$ | $(0.0003211,0.0000039)$
$z_{1}$ | $(0,0.00025)$ | $(0.000333,0.000031)$
$z_{2}$ | $(0,0.00015)$ | $(-0.000331,0.000062)$
$z_{3}$ | $(0,0.00015)$ | $(0.000617,0.000034)$
$z_{4}$ | $(0,0.00015)$ | $(0.00032,0.00011)$
Table 2:  Priors and posteriors of the vector and axial-vector matrix elements. The values are listed in $(\mu,{\sigma})$ pairs.  $g^{A}_{nm}$ | Prior | Posterior | $g^{V}_{nm}$ | Prior | Posterior
---|---|---|---|---|---
$g^{A}_{00}$ | $(1.2,0.2)$ | $(1.253,0.019)$ | $g^{V}_{00}$ | $(1,0.2)$ | $(1.02238,0.00087)$
$g^{A}_{01}$ | $(0,1)$ | $(-0.271,0.068)$ | $g^{V}_{01}$ | $(0,1)$ | $(-0.0041,0.004)$
$g^{A}_{11}$ | $(0,1)$ | $(0.89,0.2)$ | $g^{V}_{11}$ | $(1,0.2)$ | $(0.997,0.02)$
$g^{A}_{02}$ | $(0,1)$ | $(-0.16,0.11)$ | $g^{V}_{02}$ | $(0,1)$ | $(0.002,0.013)$
$g^{A}_{12}$ | $(0,1)$ | $(-0.05,0.3)$ | $g^{V}_{12}$ | $(0,1)$ | $(0.013,0.042)$
$g^{A}_{22}$ | $(0,1)$ | $(0.71,0.53)$ | $g^{V}_{22}$ | $(1,0.2)$ | $(1.05,0.12)$
$g^{A}_{03}$ | $(0,1)$ | $(0.1,0.043)$ | $g^{V}_{03}$ | $(0,1)$ | $(0.0798,0.0075)$
$g^{A}_{13}$ | $(0,1)$ | $(0.14,0.14)$ | $g^{V}_{13}$ | $(0,1)$ | $(0.071,0.025)$
$g^{A}_{23}$ | $(0,1)$ | $(0.25,0.31)$ | $g^{V}_{23}$ | $(0,1)$ | $(-0.011,0.062)$
$g^{A}_{33}$ | $(0,1)$ | $(0.87,0.22)$ | $g^{V}_{33}$ | $(1,0.2)$ | $(0.789,0.047)$
$g^{A}_{04}$ | $(0,1)$ | $(0.61,0.21)$ | $g^{V}_{04}$ | $(0,1)$ | $(0.35,0.12)$
$g^{A}_{14}$ | $(0,1)$ | $(0.4,0.25)$ | $g^{V}_{14}$ | $(0,1)$ | $(0.33,0.13)$
$g^{A}_{24}$ | $(0,1)$ | $(-1.06,0.52)$ | $g^{V}_{24}$ | $(0,1)$ | $(-0.34,0.16)$
$g^{A}_{34}$ | $(0,1)$ | $(0.72,0.39)$ | $g^{V}_{34}$ | $(0,1)$ | $(0.63,0.24)$
$g^{A}_{44}$ | $(0,1)$ | $(-0.32,0.95)$ | $g^{V}_{44}$ | $(1,0.2)$ | $(1,0.2)$
Table 3:  We list the range of $t_{\rm sep}$ values used in our optimal fits for each of the three sets of correlation functions we analyze. We also list the values of the current insertion time, ${\tau}$, for which the optimal fit uses ${\tau}=[1,t_{\rm sep}-1]$ for the 23 analysis (2pt and 3pt). When we perform the 23s (2pt, 3pt and FH) analysis, there is mild tension when including the current times ${\tau}=1$ and ${\tau}=t_{\rm sep}-1$ for $t_{\rm sep}>10$ and so our optimal fit removes one extra current insertion time for these later $t_{\rm sep}$ values. For FH, the ${\tau}$ range indicates the values used in the summation over current time.  | (2pt, 3pt) | (2pt, FH) | (2pt, 3pt, FH)
---|---|---|---
2pt $t_{\rm sep}$ range | [3, 17] | [5, 17] | [3, 17]
3pt $t_{\rm sep}$ range | [3, 14] | – | [3, 14]
FH $t_{\rm sep}$ range | – | [5, 13] | [3, 13]
3pt ${\tau}$ range | [1,$t_{\rm sep}-1$] | [1,$t_{\rm sep}-1$] | [1,$t_{\rm sep}-1$], $t_{\rm sep}\leq 10$
|  | [1,$t_{\rm sep}-1$] | [2,$t_{\rm sep}-2$], $t_{\rm sep}>10$
FH ${\tau}$ range | – | [1,$t_{\rm sep}-1$] | [1,$t_{\rm sep}-1$]
$n_{\textrm{states}}$ | 5 | 2 | 5

## References

  * Aoki  _et al._ (2020) S. Aoki _et al._ (Flavour Lattice Averaging Group), “FLAG Review 2019: Flavour Lattice Averaging Group (FLAG),” Eur. Phys. J. C 80, 113 (2020), arXiv:1902.08191 [hep-lat] .
  * Lepage (1989) G. Peter Lepage, “The Analysis of Algorithms for Lattice Field Theory,” in  _TASI, Boulder_ (1989) pp. 97–120, https://inspirehep.net/literature/287173.
  * Lepage  _et al._ (2002) G. P. Lepage, B. Clark, C. T. H. Davies, K. Hornbostel, P. B. Mackenzie, C. Morningstar, and H. Trottier, “Constrained curve fitting,” Nucl. Phys. B Proc. Suppl. 106, 12–20 (2002), arXiv:hep-lat/0110175 .
  * Wingate  _et al._ (2003) Matthew Wingate, Junko Shigemitsu, Christine T. H. Davies, G. Peter Lepage, and Howard D. Trottier, “Heavy light mesons with staggered light quarks,” Phys. Rev. D 67, 054505 (2003), arXiv:hep-lat/0211014 .
  * Gamiz  _et al._ (2009) Elvira Gamiz, Christine T. H. Davies, G. Peter Lepage, Junko Shigemitsu, and Matthew Wingate (HPQCD), “Neutral $B$ Meson Mixing in Unquenched Lattice QCD,” Phys. Rev. D 80, 014503 (2009), arXiv:0902.1815 [hep-lat] .
  * Bazavov  _et al._ (2016) A. Bazavov _et al._ (Fermilab Lattice, MILC), “$B^{0}_{(s)}$-mixing matrix elements from lattice QCD for the Standard Model and beyond,” Phys. Rev. D 93, 113016 (2016), arXiv:1602.03560 [hep-lat] .
  * Bazavov  _et al._ (2019) A. Bazavov _et al._ (Fermilab Lattice, MILC), “$|V_{us}|$ from $K_{\ell 3}$ decay and four-flavor lattice QCD,” Phys. Rev. D 99, 114509 (2019), arXiv:1809.02827 [hep-lat] .
  * Edwards  _et al._ (2006) R. G. Edwards, G. T. Fleming, Ph. Hagler, J. W. Negele, K. Orginos, A. V. Pochinsky, D. B. Renner, D. G. Richards, and W. Schroers (LHPC), “The Nucleon axial charge in full lattice QCD,” Phys. Rev. Lett. 96, 052001 (2006), arXiv:hep-lat/0510062 .
  * Yamazaki  _et al._ (2009) Takeshi Yamazaki, Yasumichi Aoki, Tom Blum, Huey-Wen Lin, Shigemi Ohta, Shoichi Sasaki, Robert Tweedie, and James Zanotti, “Nucleon form factors with 2+1 flavor dynamical domain-wall fermions,” Phys. Rev. D 79, 114505 (2009), arXiv:0904.2039 [hep-lat] .
  * Bratt  _et al._ (2010) J. D. Bratt _et al._ (LHPC), “Nucleon structure from mixed action calculations using 2+1 flavors of asqtad sea and domain wall valence fermions,” Phys. Rev. D 82, 094502 (2010), arXiv:1001.3620 [hep-lat] .
  * Alexandrou  _et al._ (2011) C. Alexandrou, M. Brinet, J. Carbonell, M. Constantinou, P. A. Harraud, P. Guichon, K. Jansen, T. Korzec, and M. Papinutto (ETM), “Axial Nucleon form factors from lattice QCD,” Phys. Rev. D 83, 045010 (2011), arXiv:1012.0857 [hep-lat] .
  * Jaffe (2002) R. L. Jaffe, “Delocalization of the axial charge in the chiral limit,” Phys. Lett. B 529, 105–110 (2002), arXiv:hep-ph/0108015 .
  * Cohen (2002) Thomas D. Cohen, “The Extraction of g(A) from finite volume systems: The Long and short of it,” Phys. Lett. B 529, 50–56 (2002), arXiv:hep-lat/0112014 .
  * Yamazaki  _et al._ (2008) T. Yamazaki, Y. Aoki, T. Blum, H. W. Lin, M. F. Lin, S. Ohta, S. Sasaki, R. J. Tweedie, and J. M. Zanotti (RBC+UKQCD), “Nucleon axial charge in 2+1 flavor dynamical lattice QCD with domain wall fermions,” Phys. Rev. Lett. 100, 171602 (2008), arXiv:0801.4016 [hep-lat] .
  * Green  _et al._ (2011) Jeremy Green, John Negele, Andrew Pochinsky, Stefan Krieg, and Sergey Syritsyn, “Excited state contamination in nucleon structure calculations,” PoS LATTICE2011, 157 (2011), arXiv:1111.0255 [hep-lat] .
  * Dinter  _et al._ (2011) Simon Dinter, Constantia Alexandrou, Martha Constantinou, Vincent Drach, Karl Jansen, and Dru B. Renner, “Precision Study of Excited State Effects in Nucleon Matrix Elements,” Phys. Lett. B 704, 89–93 (2011), arXiv:1108.1076 [hep-lat] .
  * Capitani  _et al._ (2012) S. Capitani, M. Della Morte, G. von Hippel, B. Jager, A. Juttner, B. Knippschild, H. B. Meyer, and H. Wittig, “The nucleon axial charge from lattice QCD with controlled errors,” Phys. Rev. D 86, 074502 (2012), arXiv:1205.0180 [hep-lat] .
  * Green  _et al._ (2014) J. R. Green, M. Engelhardt, S. Krieg, J. W. Negele, A. V. Pochinsky, and S. N. Syritsyn, “Nucleon Structure from Lattice QCD Using a Nearly Physical Pion Mass,” Phys. Lett. B 734, 290–295 (2014), arXiv:1209.1687 [hep-lat] .
  * Bhattacharya  _et al._ (2014) Tanmoy Bhattacharya, Saul D. Cohen, Rajan Gupta, Anosh Joseph, Huey-Wen Lin, and Boram Yoon, “Nucleon Charges and Electromagnetic Form Factors from 2+1+1-Flavor Lattice QCD,” Phys. Rev. D 89, 094502 (2014), arXiv:1306.5435 [hep-lat] .
  * Bali  _et al._ (2014) G. S. Bali, S. Collins, B. Gläßle, M. Göckeler, J. Najjar, R. Rödl, A. Schäfer, R. Schiel, A. Sternbeck, and W. Söldner, “Moments of structure functions for $N_{f}=2$ near the physical point,” PoS LATTICE2013, 290 (2014), arXiv:1311.7041 [hep-lat] .
  * Horsley  _et al._ (2014) R. Horsley, Y. Nakamura, A. Nobile, P. E. L. Rakow, G. Schierholz, and J. M. Zanotti, “Nucleon axial charge and pion decay constant from two-flavor lattice QCD,” Phys. Lett. B 732, 41–48 (2014), arXiv:1302.2233 [hep-lat] .
  * Bali  _et al._ (2015) Gunnar S. Bali, Sara Collins, Benjamin Glässle, Meinulf Göckeler, Johannes Najjar, Rudolf H. Rödl, Andreas Schäfer, Rainer W. Schiel, Wolfgang Söldner, and André Sternbeck, “Nucleon isovector couplings from $N_{f}=2$ lattice QCD,” Phys. Rev. D 91, 054501 (2015), arXiv:1412.7336 [hep-lat] .
  * Abdel-Rehim  _et al._ (2015) A. Abdel-Rehim _et al._ , “Nucleon and pion structure with lattice QCD simulations at physical value of the pion mass,” Phys. Rev. D 92, 114513 (2015), [Erratum: Phys.Rev.D 93, 039904 (2016)], arXiv:1507.04936 [hep-lat] .
  * Bhattacharya  _et al._ (2016) Tanmoy Bhattacharya, Vincenzo Cirigliano, Saul Cohen, Rajan Gupta, Huey-Wen Lin, and Boram Yoon, “Axial, Scalar and Tensor Charges of the Nucleon from 2+1+1-flavor Lattice QCD,” Phys. Rev. D 94, 054508 (2016), arXiv:1606.07049 [hep-lat] .
  * Alexandrou  _et al._ (2017) Constantia Alexandrou, Martha Constantinou, Kyriakos Hadjiyiannakou, Karl Jansen, Christos Kallidonis, Giannis Koutsou, and Alejandro Vaquero Aviles-Casco, “Nucleon axial form factors using $N_{f}$ = 2 twisted mass fermions with a physical value of the pion mass,” Phys. Rev. D 96, 054507 (2017), arXiv:1705.03399 [hep-lat] .
  * Capitani  _et al._ (2019) Stefano Capitani, Michele Della Morte, Dalibor Djukanovic, Georg M. von Hippel, Jiayu Hua, Benjamin Jäger, Parikshit M. Junnarkar, Harvey B. Meyer, Thomas D. Rae, and Hartmut Wittig, “Isovector axial form factors of the nucleon in two-flavor lattice QCD,” Int. J. Mod. Phys. A 34, 1950009 (2019), arXiv:1705.06186 [hep-lat] .
  * Berkowitz  _et al._ (2017a) E. Berkowitz, C. Bouchard, D. B. Brantley, C.C Chang, M.A. Clark, N. Garron, B Joó, T. Kurth, C. Monahan, H. Monge-Camacho, A. Nicholson, K. Orginos, E. Rinaldi, P. Vranas, and A. Walker-Loud, “An accurate calculation of the nucleon axial charge with lattice QCD,” (2017a), arXiv:1704.01114 [hep-lat] .
  * Chang  _et al._ (2018) C. C. Chang, A. Nicholson, E. Rinaldi, E. Berkowitz, N. Garron, D. A. Brantley, H. Monge-Camacho, C. J. Monahan, C. Bouchard, M. A. Clark, B. Joó, T. Kurth, K. Orginos, P. Vranas, and A. Walker-Loud, “A per-cent-level determination of the nucleon axial coupling from quantum chromodynamics,” Nature 558, 91–94 (2018), arXiv:1805.12130 [hep-lat] .
  * Berkowitz  _et al._ (2018a) Evan Berkowitz, M. A. Clark, A. Gambhir, K. McElvain, A. Nicholson, E. Rinalid, P. Vranas, A. Walker-Loud, C. C. Chang, B. Joó, T. Kurth, and K. Orginos, “Simulating the weak death of the neutron in a femtoscale universe with near-Exascale computing,”  (2018a), 10.1109/SC.2018.00058, arXiv:1810.01609 [hep-lat] .
  * Walker-Loud  _et al._ (2020) André Walker-Loud _et al._ , “Lattice QCD Determination of $g_{A}$,” PoS CD2018, 020 (2020), arXiv:1912.08321 [hep-lat] .
  * Gupta  _et al._ (2018) Rajan Gupta, Yong-Chull Jang, Boram Yoon, Huey-Wen Lin, Vincenzo Cirigliano, and Tanmoy Bhattacharya, “Isovector Charges of the Nucleon from 2+1+1-flavor Lattice QCD,” Phys. Rev. D 98, 034503 (2018), arXiv:1806.09006 [hep-lat] .
  * Hasan  _et al._ (2019) Nesreen Hasan, Jeremy Green, Stefan Meinel, Michael Engelhardt, Stefan Krieg, John Negele, Andrew Pochinsky, and Sergey Syritsyn, “Nucleon axial, scalar, and tensor charges using lattice QCD at the physical pion mass,” Phys. Rev. D 99, 114505 (2019), arXiv:1903.06487 [hep-lat] .
  * Alexandrou  _et al._ (2020) C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, and A. Vaquero Aviles-Casco, “Nucleon axial, tensor, and scalar charges and $\sigma$-terms in lattice QCD,” Phys. Rev. D 102, 054517 (2020), arXiv:1909.00485 [hep-lat] .
  * Park  _et al._ (2021) Sungwoo Park, Rajan Gupta, Boram Yoon, Santanu Mondal, Tanmoy Bhattacharya, Yong-Chull Jang, Bálint Joó, and Frank Winter (Nucleon Matrix Elements (NME)), “Precision Nucleon Charges and Form Factors Using 2+1-flavor Lattice QCD,” (2021), arXiv:2103.05599 [hep-lat] .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), “Review of Particle Physics,” PTEP 2020, 083C01 (2020).
  * Bopp  _et al._ (1986) P. Bopp, D. Dubbers, L. Hornig, E. Klemt, J. Last, H. Schütze, S. J. Freedman, and O. Schärpf, “Beta-decay asymmetry of the neutron and $\frac{{g}_{A}}{{g}_{V}}$,” Phys. Rev. Lett. 56, 919–922 (1986).
  * Yerozolimsky  _et al._ (1997) B. Yerozolimsky, I. Kuznetsov, Yu. Mostovoy, and I. Stepanenko, “Corrigendum: Corrected value of the beta-emission asymmetry in the decay of polarized neutrons measured in 1990,” Physics Letters B 412, 240 – 241 (1997).
  * Liaud  _et al._ (1997) P. Liaud, K. Schreckenbach, R. Kossakowski, H. Nastoll, A. Bussière, J.P. Guillaud, and L. Beck, “The measurement of the beta asymmetry in the decay of polarized neutrons,” Nuclear Physics A 612, 53 – 81 (1997).
  * Mostovoi  _et al._ (2001) Y. A. Mostovoi _et al._ , “Experimental value of G(A)/G(V) from a measurement of both P-odd correlations in free-neutron decay,” Phys. Atom. Nucl. 64, 1955–1960 (2001), [Yad. Fiz.64,2040(2001)].
  * Schumann  _et al._ (2008) M. Schumann, M. Kreuz, M. Deissenroth, F. Glück, J. Krempel, B. Märkisch, D. Mund, A. Petoukhov, T. Soldner, and H. Abele, “Measurement of the proton asymmetry parameter in neutron beta decay,” Phys. Rev. Lett. 100, 151801 (2008).
  * Mund  _et al._ (2013) D. Mund, B. Maerkisch, M. Deissenroth, J. Krempel, M. Schumann, H. Abele, A. Petoukhov, and T. Soldner, “Determination of the Weak Axial Vector Coupling from a Measurement of the Beta-Asymmetry Parameter A in Neutron Beta Decay,” Phys. Rev. Lett. 110, 172502 (2013), arXiv:1204.0013 [hep-ex] .
  * Mendenhall  _et al._ (2013) M. P. Mendenhall _et al._ (UCNA), “Precision measurement of the neutron $\beta$-decay asymmetry,” Phys. Rev. C87, 032501 (2013), arXiv:1210.7048 [nucl-ex] .
  * Märkisch  _et al._ (2019) B. Märkisch _et al._ , “Measurement of the Weak Axial-Vector Coupling Constant in the Decay of Free Neutrons Using a Pulsed Cold Neutron Beam,” Phys. Rev. Lett. 122, 242501 (2019), arXiv:1812.04666 [nucl-ex] .
  * Beck  _et al._ (2020) M. Beck _et al._ , “Improved determination of the $\beta$-$\overline{\nu}_{e}$ angular correlation coefficient $a$ in free neutron decay with the $aSPECT$ spectrometer,” Phys. Rev. C 101, 055506 (2020), arXiv:1908.04785 [nucl-ex] .
  * Bouchard  _et al._ (2017) Chris Bouchard, Chia Cheng Chang, Thorsten Kurth, Kostas Orginos, and Andre Walker-Loud, “On the Feynman-Hellmann Theorem in Quantum Field Theory and the Calculation of Matrix Elements,” Phys. Rev. D 96, 014504 (2017), arXiv:1612.06963 [hep-lat] .
  * Detmold  _et al._ (2012a) William Detmold, C. J. David Lin, and Stefan Meinel, “Axial couplings and strong decay widths of heavy hadrons,” Phys. Rev. Lett. 108, 172003 (2012a), arXiv:1109.2480 [hep-lat] .
  * Detmold  _et al._ (2012b) William Detmold, C. J. David Lin, and Stefan Meinel, “Calculation of the heavy-hadron axial couplings $g_{1}$, $g_{2}$ and $g_{3}$ using lattice QCD,” Phys. Rev. D 85, 114508 (2012b), arXiv:1203.3378 [hep-lat] .
  * Follana  _et al._ (2007) E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD, UKQCD), “Highly improved staggered quarks on the lattice, with applications to charm physics,” Phys. Rev. D 75, 054502 (2007), arXiv:hep-lat/0610092 .
  * Bazavov  _et al._ (2013) A. Bazavov _et al._ (MILC), “Lattice QCD Ensembles with Four Flavors of Highly Improved Staggered Quarks,” Phys. Rev. D 87, 054505 (2013), arXiv:1212.4768 [hep-lat] .
  * Miller  _et al._ (2020) Nolan Miller _et al._ , “$F_{K}/F_{\pi}$ from Möbius Domain-Wall fermions solved on gradient-flowed HISQ ensembles,” Phys. Rev. D 102, 034507 (2020), arXiv:2005.04795 [hep-lat] .
  * Miller  _et al._ (2021) Nolan Miller _et al._ , “Scale setting the M\”obius Domain Wall Fermion on gradient-flowed HISQ action using the omega baryon mass and the gradient-flow scales $t_{0}$ and $w_{0}$,” Phys. Rev. D 103, 054511 (2021), arXiv:2011.12166 [hep-lat] .
  * Brower  _et al._ (2017) Richard C. Brower, Harmut Neff, and Kostas Orginos, “The Möbius domain wall fermion algorithm,” Comput. Phys. Commun. 220, 1–19 (2017), arXiv:1206.5214 [hep-lat] .
  * Kaplan (1992) David B. Kaplan, “A Method for simulating chiral fermions on the lattice,” Phys. Lett. B 288, 342–347 (1992), arXiv:hep-lat/9206013 .
  * Shamir (1993) Yigal Shamir, “Chiral fermions from lattice boundaries,” Nucl. Phys. B 406, 90–106 (1993), arXiv:hep-lat/9303005 .
  * Furman and Shamir (1995) Vadim Furman and Yigal Shamir, “Axial symmetries in lattice QCD with Kaplan fermions,” Nucl. Phys. B 439, 54–78 (1995), arXiv:hep-lat/9405004 .
  * Berkowitz  _et al._ (2017b) Evan Berkowitz _et al._ , “Möbius domain-wall fermions on gradient-flowed dynamical HISQ ensembles,” Phys. Rev. D 96, 054513 (2017b), arXiv:1701.07559 [hep-lat] .
  * Lepage and Mackenzie (1993) G. Peter Lepage and Paul B. Mackenzie, “On the viability of lattice perturbation theory,” Phys. Rev. D 48, 2250–2264 (1993), arXiv:hep-lat/9209022 .
  * (58) André Walker-Loud, “CalLat elastic nucleon structure: I,” Lattice 2019, https://indico.cern.ch/event/764552/contributions/3420541/.
  * Jang  _et al._ (2020) Yong-Chull Jang, Rajan Gupta, Boram Yoon, and Tanmoy Bhattacharya, “Axial Vector Form Factors from Lattice QCD that Satisfy the PCAC Relation,” Phys. Rev. Lett. 124, 072002 (2020), arXiv:1905.06470 [hep-lat] .
  * Maiani  _et al._ (1987) L. Maiani, G. Martinelli, M. L. Paciello, and B. Taglienti, “Scalar Densities and Baryon Mass Differences in Lattice {QCD} With Wilson Fermions,” Nucl. Phys. B 293, 420 (1987).
  * Ottnad (2021) Konstantin Ottnad, “Excited states in nucleon structure calculations,” Eur. Phys. J. A 57, 50 (2021), arXiv:2011.12471 [hep-lat] .
  * Tiburzi (2009) Brian C. Tiburzi, “Time Dependence of Nucleon Correlation Functions in Chiral Perturbation Theory,” Phys. Rev. D 80, 014002 (2009), arXiv:0901.0657 [hep-lat] .
  * Bär and Golterman (2013) Oliver Bär and Maarten Golterman, “Excited-state contribution to the axial-vector and pseudoscalar correlators with two extra pions,” Phys. Rev. D 87, 014505 (2013), arXiv:1209.2258 [hep-lat] .
  * Bar (2015) Oliver Bar, “Nucleon-pion-state contribution to nucleon two-point correlation functions,” Phys. Rev. D 92, 074504 (2015), arXiv:1503.03649 [hep-lat] .
  * Tiburzi (2015) Brian C. Tiburzi, “Chiral Corrections to Nucleon Two- and Three-Point Correlation Functions,” Phys. Rev. D 91, 094510 (2015), arXiv:1503.06329 [hep-lat] .
  * Bär (2016) Oliver Bär, “Nucleon-pion-state contribution in lattice calculations of the nucleon charges $g_{A},g_{T}$ and $g_{S}$,” Phys. Rev. D 94, 054505 (2016), arXiv:1606.09385 [hep-lat] .
  * Bar (2017) Oliver Bar, “Chiral perturbation theory and nucleon–pion-state contaminations in lattice QCD,” Int. J. Mod. Phys. A 32, 1730011 (2017), arXiv:1705.02806 [hep-lat] .
  * Bali  _et al._ (2020) Gunnar S. Bali, Lorenzo Barca, Sara Collins, Michael Gruber, Marius Löffler, Andreas Schäfer, Wolfgang Söldner, Philipp Wein, Simon Weishäupl, and Thomas Wurm (RQCD), “Nucleon axial structure from lattice QCD,” JHEP 05, 126 (2020), arXiv:1911.13150 [hep-lat] .
  * Bar (2020) Oliver Bar, “$N\pi$ states and the projection method for the nucleon axial and pseudoscalar form factors,” Phys. Rev. D 101, 034515 (2020), arXiv:1912.05873 [hep-lat] .
  * Walker-Loud  _et al._ (2009) A. Walker-Loud _et al._ , “Light hadron spectroscopy using domain wall valence quarks on an Asqtad sea,” Phys. Rev. D 79, 054502 (2009), arXiv:0806.4549 [hep-lat] .
  * Walker-Loud (2008) Andre Walker-Loud, “New lessons from the nucleon mass, lattice QCD and heavy baryon chiral perturbation theory,” PoS LATTICE2008, 005 (2008), arXiv:0810.0663 [hep-lat] .
  * Walker-Loud (2013) André Walker-Loud, “Baryons in/and Lattice QCD,” PoS CD12, 017 (2013), arXiv:1304.6341 [hep-lat] .
  * Walker-Loud (2014) André Walker-Loud, “Nuclear Physics Review,” PoS LATTICE2013, 013 (2014), arXiv:1401.8259 [hep-lat] .
  * Drischler  _et al._ (2019) Christian Drischler, Wick Haxton, Kenneth McElvain, Emanuele Mereghetti, Amy Nicholson, Pavlos Vranas, and André Walker-Loud, “Towards grounding nuclear physics in QCD,” (2019), arXiv:1910.07961 [nucl-th] .
  * Hörz  _et al._ (2021) Ben Hörz _et al._ , “Two-nucleon S-wave interactions at the $SU(3)$ flavor-symmetric point with $m_{ud}\simeq m_{s}^{\rm phys}$: A first lattice QCD calculation with the stochastic Laplacian Heaviside method,” Phys. Rev. C 103, 014003 (2021), arXiv:2009.11825 [hep-lat] .
  * Yoon  _et al._ (2017) Boram Yoon _et al._ , “Isovector charges of the nucleon from 2+1-flavor QCD with clover fermions,” Phys. Rev. D 95, 074508 (2017), arXiv:1611.07452 [hep-lat] .
  * Gupta  _et al._ (2017) Rajan Gupta, Yong-Chull Jang, Huey-Wen Lin, Boram Yoon, and Tanmoy Bhattacharya, “Axial Vector Form Factors of the Nucleon from Lattice QCD,” Phys. Rev. D 96, 114503 (2017), arXiv:1705.06834 [hep-lat] .
  * Gusken  _et al._ (1989) S. Gusken, U. Low, K. H. Mutter, R. Sommer, A. Patel, and K. Schilling, “Nonsinglet Axial Vector Couplings of the Baryon Octet in Lattice QCD,” Phys. Lett. B 227, 266–269 (1989).
  * Bulava  _et al._ (2012) John Bulava, Michael Donnellan, and Rainer Sommer, “On the computation of hadron-to-hadron transition matrix elements in lattice QCD,” JHEP 01, 140 (2012), arXiv:1108.3774 [hep-lat] .
  * de Divitiis  _et al._ (2012) G. M. de Divitiis, R. Petronzio, and N. Tantalo, “On the extraction of zero momentum form factors on the lattice,” Phys. Lett. B 718, 589–596 (2012), arXiv:1208.5914 [hep-lat] .
  * Chambers  _et al._ (2014) A. J. Chambers _et al._ (CSSM, QCDSF/UKQCD), “Feynman-Hellmann approach to the spin structure of hadrons,” Phys. Rev. D 90, 014510 (2014), arXiv:1405.3019 [hep-lat] .
  * Chambers  _et al._ (2015) A. J. Chambers _et al._ , “Disconnected contributions to the spin of the nucleon,” Phys. Rev. D 92, 114517 (2015), arXiv:1508.06856 [hep-lat] .
  * Savage  _et al._ (2017) Martin J. Savage, Phiala E. Shanahan, Brian C. Tiburzi, Michael L. Wagman, Frank Winter, Silas R. Beane, Emmanuel Chang, Zohreh Davoudi, William Detmold, and Kostas Orginos, “Proton-Proton Fusion and Tritium $\beta$ Decay from Lattice Quantum Chromodynamics,” Phys. Rev. Lett. 119, 062002 (2017), arXiv:1610.04545 [hep-lat] .
  * Edwards and Joo (2005) Robert G. Edwards and Balint Joo (SciDAC, LHPC, UKQCD), “The Chroma software system for lattice QCD,” Nucl. Phys. B Proc. Suppl. 140, 832 (2005), arXiv:hep-lat/0409003 .
  * Basak  _et al._ (2005a) S. Basak, R. G. Edwards, G. T. Fleming, U. M. Heller, C. Morningstar, D. Richards, I. Sato, and S. Wallace, “Group-theoretical construction of extended baryon operators in lattice QCD,” Phys. Rev. D 72, 094506 (2005a), arXiv:hep-lat/0506029 .
  * Basak  _et al._ (2005b) Subhasish Basak, Robert Edwards, George T. Fleming, Urs M. Heller, Colin Morningstar, David Richards, Ikuro Sato, and Stephen J. Wallace (Lattice Hadron Physics (LHPC)), “Clebsch-Gordan construction of lattice interpolating fields for excited baryons,” Phys. Rev. D 72, 074501 (2005b), arXiv:hep-lat/0508018 .
  * Dudek  _et al._ (2013) Jozef J. Dudek, Robert G. Edwards, and Christopher E. Thomas (Hadron Spectrum), “Energy dependence of the $\rho$ resonance in $\pi\pi$ elastic scattering from lattice QCD,” Phys. Rev. D 87, 034505 (2013), [Erratum: Phys.Rev.D 90, 099902 (2014)], arXiv:1212.0830 [hep-ph] .
  * Lang and Verduci (2013) C. B. Lang and V. Verduci, “Scattering in the $\pi$N negative parity channel in lattice QCD,” Phys. Rev. D 87, 054502 (2013), arXiv:1212.5055 [hep-lat] .
  * Wilson  _et al._ (2015) David J. Wilson, Raul A. Briceno, Jozef J. Dudek, Robert G. Edwards, and Christopher E. Thomas, “Coupled $\pi\pi,K\bar{K}$ scattering in $P$-wave and the $\rho$ resonance from lattice QCD,” Phys. Rev. D 92, 094502 (2015), arXiv:1507.02599 [hep-ph] .
  * Briceño  _et al._ (2019) Raúl A. Briceño, Maxwell T. Hansen, and Andrew W. Jackura, “Consistency checks for two-body finite-volume matrix elements: I. Conserved currents and bound states,” Phys. Rev. D 100, 114505 (2019), arXiv:1909.10357 [hep-lat] .
  * (91) “LALIBE: CalLat Collaboration code for lattice qcd calculations,” https://github.com/callat-qcd/lalibe.
  * Clark  _et al._ (2010) M.A. Clark, R. Babich, K. Barros, R.C. Brower, and C. Rebbi, “Solving Lattice QCD systems of equations using mixed precision solvers on GPUs,” Comput.Phys.Commun. 181, 1517–1528 (2010), https://github.com/lattice/quda, arXiv:0911.3191 [hep-lat] .
  * Babich  _et al._ (2011) R. Babich, M. A. Clark, B. Joo, G. Shi, R. C. Brower, and S. Gottlieb, “Scaling Lattice QCD beyond 100 GPUs,” in  _SC11 International Conference for High Performance Computing, Networking, Storage and Analysis_ (2011) arXiv:1109.2935 [hep-lat] .
  * The HDF Group (1997-NNNN) The HDF Group, “Hierarchical Data Format, version 5,” (1997-NNNN), http://www.hdfgroup.org/HDF5/.
  * Kurth  _et al._ (2015) Thorsten Kurth, Andrew Pochinsky, Abhinav Sarje, Sergey Syritsyn, and Andre Walker-Loud, “High-Performance I/O: HDF5 for Lattice QCD,” PoS LATTICE2014, 045 (2015), arXiv:1501.06992 [hep-lat] .
  * Berkowitz (2017) Evan Berkowitz, “METAQ: Bundle Supercomputing Tasks,” (2017), https://github.com/evanberkowitz/metaq, arXiv:1702.06122 [physics.comp-ph] .
  * Berkowitz  _et al._ (2018b) Evan Berkowitz, Gustav R. Jansen, Kenneth McElvain, and André Walker-Loud, “Job Management and Task Bundling,” EPJ Web Conf. 175, 09007 (2018b), arXiv:1710.01986 [hep-lat] .
  * Chang  _et al._ (2020) Chia Cheng Chang, Christopher Körber, and André Walker-Loud, “EspressoDB: A Scientific Database for Managing High-Performance Computing Workflow,” J. Open Source Softw. 5, 2007 (2020), arXiv:1912.03580 [hep-lat] .
  * Lepage (2020) Peter Lepage, “gvar,” (2020), https://github.com/gplepage/gvar.
  * Lepage (2021) Peter Lepage, “lsqfit,” (2021), https://github.com/gplepage/lsqfit.
  * Berkowitz  _et al._ (2018c) Evan Berkowitz _et al._ , “Progress in Multibaryon Spectroscopy,” PoS LATTICE2018, 003 (2018c), arXiv:1902.09416 [hep-lat] .