Systematics in nucleon matrix element calculations
Jeremy Green




§ INTRODUCTION


The structure of nucleons has been studied extensively in experiments,
and nucleons also play a vital role as experimental probes. This makes
the controlled study of nucleons using lattice QCD a very attractive
goal. Some of the simplest matrix elements of interest include the
scalar and tensor charges, which control BSM contributions to neutron
beta decay, and the sigma terms, which control the sensitivity of dark
matter detection to WIMPs.

In order to make a significant impact, lattice calculations of nucleon
structure must become precision calculations, with full control over
all systematics. The traditional test observable is the isovector
axial charge g_A:

    ⟨ p(P,s')|u̅γ^μγ_5 d|n(P,s)⟩≡ g_A u̅_p(P,s') γ^μγ_5 u_n(P,s).

It is simple to compute, being an isovector forward matrix element,
and it is measured precisely in beta decay experiments[It
  should be noted, however, that over time the experimental value has
  drifted upward slightly; see the PDG's history plot.], with the
latest PDG value <cit.> being
g_A=1.2724(23). However, calculating it accurately has been
challenging: there is a long history of results being below the
experimental value. Because systematic effects in g_A are
significant and it has been studied the most extensively, the axial
charge will be the exclusive focus of this review. This review will
also focus on recent results, particularly those presented at this
conference.

This review is organized as follows. Section <ref> is a
brief overview of how matrix elements are computed. Excited-state
effects, which are a particularly challenging source of systematic
uncertainty, are discussed at length in
Section <ref>. Finite-volume effects and dependence on the
pion mass are briefly reviewed in Sections <ref> and
<ref>. Finally, a summary and outlook is given in
Section <ref>.



§ METHODOLOGY OF MATRIX ELEMENTS


The simplest approach for computing the forward hadronic matrix
element of operator 𝒪 is to use a single interpolator
χ at zero momentum. One computes two-point and three-point
functions and performs a spectral decomposition,

    C_2pt(t)    ≡⟨χ(t) χ^†(0) ⟩       
       C_3pt(τ,T)    ≡⟨χ(T) 𝒪(τ) χ^†(0) ⟩
       = ∑_n |Z_n|^2 e^-E_n t,       = ∑_n,n' Z_n' Z_n^* ⟨ n'|𝒪|n⟩
         e^-E_nτ e^-E_n'(T-τ),

where Z_n≡⟨Ω|χ|n⟩ is the overlap of the
interpolator onto state n. In the limit where all time separations
t, τ, and T-τ are large, the ground state dominates:

    C_2pt(t)    → |Z_0|^2 e^-E_0 t(1 + O(e^-Δ E t)),
    
      C_3pt(τ,T)    → |Z_0|^2 e^-E_0 t( ⟨ 0|𝒪|0⟩
      + O(e^-Δ E τ) + O(e^-Δ E(T-τ) + O(e^-Δ E T) ),

where Δ E≡ E_1-E_0 is the energy gap to the first excited
state.

In the ratio method, the prefactors are cancelled to obtain
the ground-state matrix element:

    R(τ,T) ≡C_3pt(τ,T)/C_2pt(T)→⟨ 0|𝒪|0⟩ + O(e^-Δ Eτ) + O(e^-Δ E(T-τ))
      + O(e^-Δ E T).

For forward matrix elements, the excited-state contributions are
symmetric about τ=T/2. They can be minimized by placing the
operator at the midpoint, yielding R(T2,T) = ⟨
0|𝒪|0⟩ + O(e^-Δ E T/2).

An alternative approach is the summation
  method <cit.>, which involves summing
over the operator insertion time τ. The terms that are
independent of τ grow linearly with the source-sink separation
T in the sum, whereas the terms that have an exponential dependence
on τ produce partial sums of geometric series. The derivative of
the sum yields the ground-state matrix element:

    S(T) ≡ a∑_τ R(τ,T),  d/dTS(T) = ⟨ 0|𝒪|0⟩ + O(T e^-Δ E T).

There is some flexibility about the interval over which τ is
summed. One option is to use the “interior” region where
𝒪 is between the source and the sink, possibly excluding a
fixed number of points at each end so that the sum is from τ_0 to
T-τ_0. Another option is to sum over the whole lattice, which is
the method that was first used. In two talks at the Lattice 2010
conference <cit.>, it was pointed out
that contributions from excited states decay more rapidly than for the
ratio-midpoint method; this revived interest in the summation
method. In practice, the summation method has been found to produce a
larger statistical uncertainty than the ratio method, which can negate
some of its advantage; see Fig. <ref>.



In the last few years, there has been some use of methodologies based
on the Feynman-Hellmann theorem. This states that a matrix element in
a given state can be obtained from the derivative of the state's
energy with respect to a perturbation in the Lagrangian:

    ℒ(λ)≡ℒ + λ𝒪.∂/∂λE_n(λ)|_λ=0
    = ⟨ n|𝒪|n⟩.

Discrete derivatives are sometimes used, particularly for the nucleon
sigma terms, where the theorem relates nucleon scalar matrix elements
to derivatives of the nucleon mass with respect to quark
masses. Evaluating the derivative of a two-point function exactly
leads directly to the summation method:

    -.∂/∂λlog C_2pt(t)|_λ=0
     = S(t),

where the sum is taken over all timeslices. In the large-t limit,
the time derivative of this equation yields Eq. (<ref>) for the
ground state. This result appeared in the original summation-method
paper <cit.> and has been rederived in recent
years <cit.>.



§ EXCITED-STATE CONTAMINATION


Unwanted contributions from excited states decay exponentially and
will be highly suppressed if Δ E T ≫ 1. However, the
signal-to-noise problem <cit.> prevents calculations
from being performed at large T: this ratio decays as
e^-(E_0-32m_π)T. To be more concrete about the “brute
force” approach of simply using a large source-sink separation,
assume that we are using the ratio method in the asymptotic regime,
where the statistical errors and excited-state contributions scale as

    δ_stat∝ N^-1/2e^(m_N-32m_π)T,  δ_exc∝ e^-Δ E T/2,

respectively, where N is the number of statistical
samples. Supposing that we want these two uncertainties to scale
together, i.e. δ_exc=αδ_stat≡δ for some constant α, then as δ is decreased the
source-sink separation must be increased to suppress excited state
effects. An increase in statistics is also required, both to
compensate for the reduced signal-to-noise ratio, and also to meet the
smaller target statistical error; the required statistics are given by

    N ∝δ^-(2+4m_N-6m_π/Δ E).

At the physical pion mass with Δ E=2m_π (see the next
subsection), the exponent is roughly -13, much larger than the -2
that is obtained when neglecting excited states. This situation could
be significantly improved if multilevel
methods <cit.> or other
ideas <cit.> are able to reduce the signal-to-noise
problem.

Because of the difficulty in going to large T, there has been much
effort spent on removing excited states from data at relatively small
source-sink separations: the two main approaches are improving the
interpolating operator and modeling the excited-state contributions.



 §.§ Theoretical expectations


An approximation to the finite-volume spectrum consists of the energy
levels of any number of noninteracting stable hadrons,
E=∑_i,j√(m_i^2+p_j^2), where p_j=2πLn_j and
n_j is a vector of integers. For a nucleon at rest, the leading
excitations are states with a nucleon and any number of
pions. Positive parity requires that for an odd number of pions there
must be some nonzero momenta. These noninteracting energy gaps are
shown in Fig. <ref>. At the physical pion mass with m_π
L=4, the lowest excitation is Δ E=2m_π and there are eight
Nπ and Nππ levels below 4m_π. As the box size grows, the
spectrum becomes denser.
On the other hand, at heavier pion
masses the energies rapidly increase and there are soon just a few
levels below Δ E=1 GeV. Going beyond the noninteracting
approximation, in Ref. <cit.> finite-volume
quantization conditions were applied in the Nπ sector using the
experimentally measured scattering phase shift. Significant deviations
from the noninteracting levels were found, particularly for energy
gaps between 0.4 and 0.8 GeV, however the general features of the
Nπ spectrum were not significantly changed.



Predicting excited-state contributions to two-point and three-point
functions requires knowing the spectrum E_n, the overlap factors
Z_n, and the operator matrix elements ⟨
n'|𝒪|n⟩. The key insight that allows for study using
chiral perturbation theory (ChPT) is that, if one assumes the smearing
size of the interpolator is small compared with m_π^-1, then at
leading order a single low-energy constant controls the coupling of a
local interpolator to both nucleon and nucleon-pion states, and it is
eliminated when forming ratios <cit.>.

The prediction from ChPT for the nucleon effective mass is a
percent-level excited-state contribution for T≳
1 fm <cit.>. On the other
hand, using the ratio method, an effect at the 10% level is predicted
for g_A at T=1 fm <cit.>. This effect
increases the effective lattice value of g_A, in contrast with most
numerical studies, which find that excited states decrease the
extracted value of g_A. At this conference, O. Bär presented the
first study of observables at nonzero momentum transfer, namely the
form factors G_A and G_P of the axial current; a new tree-level
diagram was found to produce a very large excited-state effect in
G_P <cit.>.

ChPT is not expected to produce accurate results for the contributions
from excited states with higher energies, particularly in the vicinity
of resonances; this means that it is only expected to be valid at
large source-sink separations where these contributions are
suppressed, namely T≳ 2 fm for g_A and similar
observables <cit.>. Deviations from ChPT for Nπ states
in the resonance regime were modeled in Ref. <cit.>;
see Fig. <ref>. Under some model scenarios, the ratio-method
result for g_A could be suppressed by excited states at short
source-sink separations and then rise to sit 1–2% above the true
value for a wide range of larger source-sink separations. Clearly,
this corresponds to multiple excited states contributing with
different signs. This sort of scenario is particularly troublesome, as
it makes systematically improving a calculation by removing
excited-state effects quite difficult.





 §.§ Numerical studies


Before discussing studies in the literature, it should be noted that
the available set of source-sink separations T and source-operator
separations τ depends on how the three-point function is
computed.


  * The most common approach uses a fixed
    sink. C_3pt(τ,T) is evaluated using a sequential
  propagator through the sink, so that T, the sink momentum, and the
  interpolating operators are fixed. All values of τ can be
  obtained and any quark bilinear operator 𝒪(τ) can be
  used. The computational cost increases with every value of T.

  * If one instead uses a fixed operator, the sequential
  propagator is evaluated through 𝒪(τ), which is thus
  fixed. All values of T can be obtained, and the sink momentum and
  interpolator can be varied. The computational cost increases with
  each operator insertion and with each value of τ. This approach
  has been used recently in some variational studies by the CSSM
  group <cit.>.

  * Rather than fixing τ, one can sum over it to obtain a
  summed operator in the sequential propagator. In this case,
  T becomes the only relevant time separation and all values can be
  obtained. The sink interpolator can be varied and the computational
  cost increases with each operator insertion. This approach has been
  used recently by CalLat <cit.> and
  NPLQCD <cit.>.

It is possible to replace the sequential propagator with a stochastic
one. This allows for increased flexibility but at the possible cost of
increased noise <cit.>.



  §.§.§ Improving the nucleon interpolator


The usual[A common alternative is to replace Cγ_5 with
  Cγ_5 P_+, where P_+≡12(1+γ_0) is a
  positive parity projector.] proton interpolator is
χ=ϵ_abc(u^T_a Cγ_5 d_b)u_c, constructed using
smeared quark fields. Standard practice is to tune the smearing width
such that the nucleon effective mass reaches a plateau as early as
possible. It is well known from spectroscopy that the variational
method <cit.>, where one finds a linear
combination of interpolating operators χ_var≡∑_i
c_i χ_i with optimized coefficients c_i, is a powerful
systematic approach for eliminating contributions from the
lowest-lying excited states to estimates of energy levels and matrix
elements <cit.>. In practice, its
effectiveness depends on the choice of interpolator basis
{χ_i}.

A simple way to produce several interpolators is to vary the smearing
width; there are two recent studies[See
  Refs. <cit.> for earlier studies.] that
used bases comprising interpolators with three different smearing
widths <cit.>. Some results from
Ref. <cit.> are shown in
Fig. <ref>. The effective mass from the
variationally optimized interpolator lies very close to that of the
standard interpolator with the widest smearing. The optimized operator
and the widest smearing also both show little sign of excited-state
effects in the plateaus for the axial charge, but narrower smearings
do show clear signs (see e.g. Fig. 15 of
<cit.>). Ref. <cit.> also found that
narrower smearings suffer from larger excited-state effects. In that
study, the variationally optimized interpolator produced smaller
excited-state effects than the largest smearing. However, a still
larger smearing might be as good as the variational interpolator. One
should take two clear lessons from these studies. The first is that
tuning the smearing width is important, since it can have a
significant effect on excited states. The second is that when studying
a computationally more expensive alternative to the standard approach,
a fair comparison should be with a well-tuned operator; it is easy to
make the standard approach appear to be worse by using a too-narrow smearing.



Beyond varying the smearing width, one can use different local
operator structures such as ϵ_abc(u_a^TCd_b)γ_5
u_c. This approach has been used to add negative-parity interpolators
to the basis, which can be important for coupling to excited states in
moving frames <cit.>. Including covariant derivatives
or the chromomagnetic field strength in the interpolator allows for a
much larger basis <cit.>. In
Ref. <cit.>[C. Egerer presented preliminary
  results at this conference.], such a basis was used, employing the
distillation method to efficiently construct the correlators. In
general, a larger basis was more effective at removing excited-state
contributions, with the largest effect seen in the tensor charge. This
reference also reports that the Laplacian-Heaviside smearing used in
distillation produces smaller excited-state effects from the standard
operator than the more commonly used Wuppertal smearing, however it is
unclear whether either smearing was tuned.



  §.§.§ Fitting excited states


A natural strategy for removing contributions from excited states is
to fit correlators [or derived quantities such as R(τ,T) or
S(T)] using a model that includes excited-state effects. By far the
most common model is based on a truncation of the spectral
decomposition to a small number of states (two or three). Generally
the two-point function provides the strongest constraints on the
energies E_n and overlaps Z_n, and the three-point function serves
to determine the matrix elements ⟨ n'|𝒪|n⟩.
Commonly, the energy gap Δ E is found to be between 0.5 and
1.0 GeV, which is usually greater than the expected lowest-lying
excitation energy shown in Fig. <ref>.



Some recent examples are given in Refs. <cit.>. In Ref. <cit.>, the extracted value of
g_A was shown to be stable when the fit range is varied; similarly,
in Ref. <cit.> stability was shown with respect to
varying fit ranges and the number of states in the fit
model. Figure <ref> shows the preferred fits for g_A on
the same ensemble (albeit with different valence quark action). In
both cases, the fits start with a minimum time separation of
3a. Given the pion mass and box size, there are more than ten
noninteracting excited energy levels with energy gap Δ E<1 GeV;
for the first points included in the fit these are only suppressed by
e^-3aΔ E>0.25. Clearly, it is difficult to associate the
“excited state” in the fit with a single actual state.

At this conference, K. Ottnad presented a different fitting
approach <cit.> that does not determine the energy gap
Δ E from the two-point function. Instead, fits are performed to
the ratios R(τ,T); fitting six different observables
simultaneously is sufficient to constrain Δ E. Interestingly,
in this case the fitted energy gap approaches the expected
lowest-lying noninteracting level as the minimum time separation
included in the fit, t_start, is increased.

In order for these fits to be trustworthy, ideally they would be
required to have good fit qualities (p-values). (This is not a
sufficient condition!)  A strong test is given in
Ref. <cit.>: when the fit is repeated on many
ensembles, the distribution of fit qualities should be uniform, which
can be checked using a Kolmogorov-Smirnov (KS) test. More than half of
the fit qualities in Ref. <cit.> are below 0.2, and
hence the KS test indicates that they are not compatible with the
uniform distribution (p<10^-3). In contrast, the fits in
Ref. <cit.> are acceptable from this point of view
(p=0.26). Some caution may be required, however, when fitting to
many variables, as the difficulty in inverting a large covariance
matrix can make it difficult to reliably estimate χ^2 and the fit
quality.



 §.§ Outlook on excited states


It is natural to ask why the picture in Fig. <ref> of
low-lying Nπ and Nππ states does not appear in the spectrum
of typical variational analyses or multi-state fits. One argument is
that the coupling of a multiparticle state to a local interpolator is
suppressed by the inverse lattice volume. However, this should be
compensated by the density of states so that in infinite volume the
interpolator will couple to continua of multiparticle states. In fact,
model predictions such as Fig. <ref> show a weak volume
dependence; this suggests that continuum spectral functions might
yield suitable fit models for excited states in large volumes.

The absence of multiparticle states is familiar from meson
spectroscopy. It has been found that a variational basis must include
nonlocal operators in order to identify the complete and correct
spectrum; see, e.g., <cit.>. Even though nucleon
structure only requires removing excited states and not obtaining
precise knowledge of them, it may still be necessary to include
nonlocal operators in order to benefit in practice from the proven
improved asymptotic approach to the ground state when using the
variational method <cit.>.

Given that state of the art nucleon structure calculations have not
reconciled their data with theoretical expectations for excited-state
effects, perhaps the safest approach is to analyze data in multiple
ways: ratio and summation methods, which don't make specific
assumptions about the spectrum of excitations, as well as fits that
can make use of shorter time separations. This was done in the
extensive excited-state study of Ref. <cit.>, which
also included a variational setup, as well as in some recent
physical-pion-mass calculations by ETMC (see e.g.<cit.>).



§ FINITE-VOLUME EFFECTS


There is a long history of attributing a low value for g_A computed
on the lattice to finite-volume effects. These effects were computed
in ChPT in Ref. <cit.>; neglecting loops with Δ
baryons, the leading contribution at large volume is

    g_A(L)-g_A/g_A∼m_π^2 g_A^2/π^2F_π^2√(π/2m_π L) e^-m_π L.

In addition to being exponentially suppressed at large m_π L, if
one fixes m_π L and decreases m_π this effect will also be
reduced. If this expression holds true, then a calculation with m_π
L=3 at the physical pion mass will have smaller finite-volume effects
than one with m_π L=4 at m_π=300 MeV.



There have been several fully-controlled studies of finite-volume
effects in g_A, i.e., the same calculation performed on ensembles
that differ only by their volume. These are
summarized[Calculations where the lattice temporal extent was
  varied together with the spatial volume are also included.] in
Fig. <ref>. For large values of m_π L and small pion
masses, no effect is observed within uncertainties at the few percent
level. As the right hand side of Eq. (<ref>) is
increased, the first significant effect is at the 5% level in the
calculation by RQCD at m_π=290 MeV and m_π L =3.4. However, this
is a negative effect rather than the positive effect predicted by
Eq. (<ref>). Encouragingly, the physical pion mass with
m_π L=3 corresponds to 0.03 on the horizontal axis, where no
effect has been detected.

Global fits to a set of ensembles — where the pion mass, lattice
spacing, and volume are all varied — provide a different approach to
study finite-volume effects. The challenge is that any failure of the
fit function to accurately describe the dependence on the other
variables is a source of systematic uncertainty in the estimate of
finite-volume effects. This is especially true because most sets of
ensembles will tend to have larger values of m_π L at larger pion
masses and on coarser lattice spacings. Using a global fit,
Ref. <cit.> assumed the volume dependence of g_A has
the form cm_π^2 e^-m_π L with c a free parameter, and found
a -0.9(5)% effect at the physical pion mass with m_π L=4. A
similar approach was used by Ref. <cit.>, and a similar
effect size was found. Finally, Ref. <cit.> assumed the
leading heavy baryon ChPT expression and allowed a higher-order term
proportional to m_π^3; this also produced a small effect at the
physical pion mass.



§ CHIRAL EXTRAPOLATION


In heavy baryon ChPT, the pion mass dependence of the axial charge is
known to take the form

    g_A(m_π) = g_0 - (g_0 + 2g_0^3)(m_π/4π F_π)^2logm_π^2/μ^2 + c_1 m_π^2 + c_2 m_π^3 + O(m_π^4),

where g_0 is the axial charge in the chiral limit and c_1,2 are
additional low-energy constants. Although the prefactor of the chiral
log is known, it is unclear how high of a pion mass can be reached
before the convergence of ChPT breaks down. As a result, recent
calculations that performed a chiral fit <cit.>, including results presented by R. Gupta and
K. Ottnad at this conference, generally preferred to use a simple
polynomial dependence on m_π^2.

The issue of chiral extrapolation can, of course, be avoided by using
only lattice ensembles with near-physical pion masses. Results using
this approach were presented at this conference by
M. Constantinou <cit.>,
Y. Kuramashi <cit.>, C. Lauer <cit.>,
Y. Lin <cit.>, and S. Ohta <cit.>.


  
    < g r a p h i c s >

  Recent calculations of the nucleon axial charge using 2
    flavours of dynamical sea quarks <cit.> (green), 2+1
    flavours <cit.> (orange), and 2+1+1
    flavours <cit.> (blue). Published results are shown with a
    filled symbol. The PDG value <cit.> is indicated
    by the vertical band.
  


Figure <ref> shows some recent determinations of the
axial charge. It is encouraging that several collaborations are now
able to reproduce the experimental value, although there is still a
tendency for results to sit below experiment and no result is even
half a standard deviation above experiment. Given that the
experimental value is known, it may be useful for future calculations
to perform a blinded analysis.



§ SUMMARY AND OUTLOOK


Excited-state contamination remains a major focus of nucleon structure
calculations. A full variational study including nonlocal
interpolators that couple well to multiparticle states could help to
determine whether current methods are adequate. In contrast, no sign
of large finite-volume effect in g_A has been observed in the
existing fully-controlled studies at low pion masses.

Discretization effects were not discussed in this review, in part
because they are not universal. In general they appear to be less
important than excited states, but they are nevertheless important for
controlling uncertainties and can have a significant impact on the
outcome such as in Ref. <cit.>.

It should be stressed that the results in Fig. <ref>
have not been filtered based on any quality criteria. Such an
evaluation is necessary for obtaining a reliable “lattice QCD
average” of any observable. However, it may be particularly difficult
to set standards for controlling excited-state effects, since analysis
strategies vary significantly. A first community attempt was made in
Ref. <cit.>, and some nucleon structure will also be
included in the next FLAG review.

Bringing simple observables like g_A under precise control over all
sources of systematic uncertainty will be an important step toward
reliable calculations of more complex observables such as the proton
charge radius and parton distribution functions. Systematics in these
observables will require further study; in particular, finite-volume
effects could be important for form factors <cit.>
and discretization effects might be significant for parton
distribution functions, which are generally not O(a) improved.



I thank everyone who sent results in advance of my talk and who
replied to my questions: C. Alexandrou, C. C. Chang, C. Egerer,
R. Gupta, J. Liang, S. Ohta, K. Ottnad, F. M. Stokes, A. Walker-Loud,
and T. Yamazaki. I also thank my colleagues at DESY for their comments
on an early version of this talk and Karl Jansen for sending comments
on a draft of these proceedings.


utphys-noitalics

