On the effect of excited states in lattice calculations of the nucleon axial charge

Excited-state contamination is one of the dominant uncertainties in lattice calculations of the nucleon axial-charge, $g_A$. Recently published results in leading-order chiral perturbation theory (ChPT) predict the excited-state contamination to be independent of the nucleon interpolator and positive. However, empirical results from numerical lattice calculations show negative contamination (downward curvature), indicating that present-day calculations are not in the regime where the leading-order ChPT predictions apply. In this paper we show that, under plausible assumptions, one can reproduce the behavior of lattice correlators by taking into account final-state $N \pi$ interactions, in particular the effect of the Roper resonance, and by postulating a sign change in the infinite-volume $N \to N \pi$ axial-vector transition amplitude.

Excited-state contamination is one of the dominant uncertainties in lattice calculations of the nucleon axial-charge, gA. Recently published results in leading-order chiral perturbation theory (ChPT) predict the excited-state contamination to be independent of the nucleon interpolator and positive [1][2][3]. However, empirical results from numerical lattice calculations show negative contamination (downward curvature), indicating that present-day calculations are not in the regime where the leading-order ChPT predictions apply. In this paper we show that, under plausible assumptions, one can reproduce the behavior of lattice correlators by taking into account final-state N π interactions, in particular the effect of the Roper resonance, and by postulating a sign change in the infinite-volume N → N π axial-vector transition amplitude.

I. INTRODUCTION
In the past decades, outstanding progress has been made in reproducing the properties of the strong interaction by numerically calculating QCD correlators on a Euclidean spacetime lattice. One goal of such calculations is to extract various aspects of nuclear structure from the underlying theory, and a target quantity here is the nucleon axial charge, g A , defined via where A a µ ≡ QΓ a A,µ Q and Γ a A,µ = T a γ µ γ 5 . Here T a = τ a /2 are the generators of SU (2) isospin and Q is a doublet containing the up and down quarks. We have also introduced single nucleon states with momentum p and spin σ, σ as well as their corresponding spinors u σ (p) and u σ (p). These are also isospin doublets built from the proton and neutron. In this work we use Euclidean conventions for the gamma matrices, {γ µ , γ ν } = 2δ µν .
The axial charge is in many ways an ideal quantity for lattice QCD (LQCD). In particular, it can be directly accessed from plateaus in Euclidean correlators and does not contain the noisy quark-disconnected diagrams. However, as a nuclear quantity, it suffers from the signal-to-noise problem and this is only made worse in the three-point function required to create a nucleon, couple it to the axial current and then annihilate it. For some time now, lattice calculations of g A have been prone to underestimate the quantity. 1 Possible explanations for this include underestimated systematic uncertainties from extrapolation to the physical pion mass, from finite-volume effects and from excited-state contamination. This work is concerned with the latter.
two-point correlators, constructed to satisfy b n e −∆En(T −t) + e −∆Ent + c n e −∆EnT , (2) where we have dropped subleading exponentials as we explain in detail in the following section. Here we introduce T as the nucleon source-sink separation, t as the current insertion time, and ∆E n as the gap between the nucleon mass and the (n − 1)th excited state. The coefficients b n and c n are related to finite-volume matrix elements as given in Eqs. (10) and (11) below. The excited-state contribution to R(T, t) has been recently studied in both non-relativistic [1] and relativistic [2,3] baryon chiral perturbation theory (ChPT). In both cases the authors find that the leading-order (LO) ChPT predictions are independent of the form of the nucleon interpolators. 2 This leads to the universal prediction that b n > 0, and thus that the excited-state contamination is positive. Since the predictions for b n and c n depend only on g A , the pion decay constant, f π , and known kinematic quantities, the ChPT expressions could in principle be used to remove the leading excited-state contribution in order to more accurately extract g A .
To make use of the LO ChPT results, however, one must ensure that these describe present-day numerical LQCD data. As g A is often extracted from the central value, R(T, T /2), or by fitting a constant to a range of central values, determining the T values needed for R(T, T /2) to enter the regime of LO ChPT is particularly useful. If the source-sink separation is too small, then the set of finite-volume states needed to estimate R(T, T /2) goes beyond the region described by the leading-order prediction. Indeed, the curvature of nearly all available numerical LQCD data for R(T, T /2) as a function of T is negative, indicating negative excited-state contamination, in contradiction with the LO ChPT prediction. 3 Similarly, at fixed T , R(T, t) is consistently observed to have negative curvature as a function of the currentinsertion time, t. We take this as strong evidence that, in present day LQCD calculations, the values of T are too small for R(T, t) to be well described by the LO ChPT results.
In this paper we show that, under plausible assumptions, one can reproduce the qualitative behavior of numerical LQCD correlators by including the contributions of higher-energy states, taking into account N π final-state interactions, and postulating a sign change in the infinite-volume axial-vector transition amplitude, N π, out|A µ |N . Using experimentally-determined N π scattering data in a generalization of Lüscher's quantization condition [11,12], we predict the energies of the finite-volume excited states entering ∆E n . We then use a generalization of the Lellouch-Lüscher formalism, again with experimental scattering data, to relate the finitevolume matrix elements in b n and c n to infinite-volume matrix elements involving N π asymptotic states [13]. To complete the construction, we estimate the remaining infinite-volume matrix elements in a model based on LO ChPT, supplemented by the scattering data.
Within this set-up we find that a large number of excited states give an important contribution to R(T, T /2) for realistic T values, and that a sign flip in the axialvector transition can readily accommodate the empirically observed negative excited-state contamination [see Figs. 9 and 10 below]. We find that, for physical pion masses, T 2 fm is needed to enter the regime where LO ChPT describes the lattice correlators.
This analysis suffers from various limitations that prevent us from offering reliable quantitative predictions. The most important limitation is the neglect of N ππ states. Here we only study the energies and matrix elements of finite-volume N π states. Both the Lüscher quantization condition and the Lellouch-Lüscher formalism hold only for energies below three-particle production threshold, but in this work we also include energies above N ππ threshold where the relations develop uncontrolled systematic uncertainties. There is evidence that the breakdown of the formalism turns on slowly as one crosses multi-particle thresholds, 4 but in the vicinity of the Roper resonance, the neglected three-particle states could have a significant contribution. [See also the discussion in the paragraph following Eq. (16) below.] Other limitations of this study include the modeling of the infinite-volume matrix elements, explained in detail in Sec. IV C, as well as the restriction to physical pion masses. The latter is a natural limitation given our approach of working with experimental scattering data. As a result, the predictions for R(T, t) discussed in Sec. V are most directly applicable to ensembles near the physical point.
As an aside we comment that, in order to have a solid theoretical foundation for this work, it was necessary to make contact with the LO ChPT results derived in Refs. [2,3]. In these earlier publications, ∆E n is approximated using non-interacting N π states in a finite volume, so that the work is concerned only with predicting the coefficients, b n and c n . Since we are using the Lellouch-Lüscher formalism to predict the coefficients in this study, it was necessary to first understand how this formalism can be used to reproduce the LO ChPT results. We were able to make this connection in detail, rederiving some of the expressions reported in Refs. [2,3]. This is interesting in its own right as it shows how the Lellouch-Lüscher formalism provides a shortcut for extracting ChPT predictions of these and related quantities. In particular, the numerous one-loop diagrams needed to determine b n in Ref. [3] are replaced in the present approach by five tree-level diagrams. Details are given in the appendix.
The remainder of this article is organized as follows. In the following section we define the correlators, the ratio R and the parameters ∆E n , b n and c n that describe the excited states. In Sec. III we use experiemental partial wave data to estimate the interacting energy gaps ∆E n associated with N π states. Then in Sec. IV we give estimates for the coefficients b n and c n . This leads to estimates of the excited state contamination for typical present-day lattice set-ups, presented in Sec. V. In the appendix we detail the derivation of various ChPT expressions used in the main text.

II. EXTRACTING gA FROM THE LATTICE
Various methods exist for using numerical LQCD to determine g A . Common to all approaches is the determination of two-and three-point correlators of the form where O α , O β are proton interpolating fields, A 3 µ is the third isospin component of the axial vector current, and Γ and Γ are projectors. In this work we restrict attention to states that have zero three-momentum in the finitevolume frame. Defining , and performing a spectral de-composition, we reach where we have assumed T > t > 0 and have used the shorthand O β ≡ O β (0) and similar for A 3 µ . To treat the fields equivalently we have Fourier transformed O α over spatial volume but have also divided by volume to preserve the definitions. Throughout this work all finitevolume states are normalized as n|n = 1.
We next observe that the lowest state in the sum, denoted by n, m = 1, is the single nucleon state. From this follows that the ratio of the n, m = 1 terms in C 3 (T, t) and This relies on the definitions of Γ and Γ . These are constructed to ensure that the result holds. It follows that g A can be accessed by identifying a plateau in the ratio Substituting the spectral decompositions, Eqs. (5) and (6), taking T t 0 and expanding the denominator, we find , with E n the energy of the (n − 1)th excited state, m N the nucleon mass and M π the pion mass. Here we have introduced L as the linear spatial extent of the volume and have used the fact that finite-volume corrections to the nucleon mass are exponentially suppressed. We neglect such corrections throughout.
In Eq. (9) we have also introduced where Note that the definition for b n , Eq. (10), directly arises from the coefficient on the first exponential, e −∆En(T −t) , whereas the factor multiplying the second exponential, e −∆Ent , has a different definition. However, as long as Γ µ is anti-hermitian and Γ is hermitian, then Euclidean definitions of charge-conjugation and time-reversal invariance imply R(T, t) = R(T, T − t) = R * (T, t). Thus the two coefficients are identically equal and we take b n as the coefficient for both source-to-current and current-to-sink time dependence. As can be seen by comparing the definitions, Eqs. (10) and (11), the matrix elements required to access the source-to-sink coefficient, c n , are more complicated than those needed for b n . The first term in the definition of c n , proportional to c 2,n defined in Eq. (12), arises from expanding the excited-state contamination of C 2 (T ) in the denominator. This term depends on the same matrix elements that appear in the definition of b n and can be studied using the same approach. The second term in c n , c 3,n defined in Eq. (13), arises from source-to-sink contributions in C 3 (T, t) and is thus more complicated. This term turns out to be numerically suppressed in LO ChPT and thus unimportant in our qualitative study. With this in mind, in this work we simply use the LO ChPT result for c 3,n and only apply the Lellouch-Lüscher like analysis to b n and c 2,n .
The ellipsis in Eq. (9) stands for terms suppressed by additional factors of e −∆Emt , e −∆Em(T −t) or e −∆EmT . These neglected terms arise for two reasons. One contribution is from higher orders in the expansion of C 2 (T ) in the denominator. This expansion is not required but is a good approximation and simplifies the resulting expressions. The second neglected contribution is from terms in Eq. (5) with n = m, with both indices corresponding to excited states. Such terms involve two-to-two (rather than one-to-two) axial-vector matrix elements and are expected to be suppressed relative to those we keep. We caution that these two-to-two transitions are not necessarily volume suppressed. For example, LO ChPT predicts the same volume dependence for c 2,n and c 3,n . This is the case because, at leading order, the current mediating the two-to-two transition couples only to one of the two particles. When the current couples to both particles an extra factor of volume suppression does arise. 5 The aim of this work is to estimate the value of the sum in Eq. (9) for given T and t. In the following section we study ∆E n and in Sec. IV we turn to b n and c n .

III. ESTIMATING THE EXCITED-STATE ENERGIES
The finite-volume quantization condition derived by Lüscher [11,12] has since been extended to include moving frames, non-identical and non-degenerate particles, coupled two-particle channels, and particles with spin [16][17][18][19][20][21][22][23][24][25][26][27]. These extensions can be used to estimate the finite-volume energies that appear in R(T, t). In particular, in the range m N + M π < E n < m N + 2M π , the finite-volume energies can be determined using the Lüscher quantization condition by inputting the experimentally determined phase shift for N π scattering.
It is useful to consider these energies relative to the energies of non-interacting particles in a finite-volume. The non-interacting levels are determined by constraining the momentum to satisfy p = 2πn/L, where L is the linear extent of the volume and n is a three-vector of integers. This constraint is appropriate to a cubic finite spatial volume with periodic boundary conditions, and in this work we restrict ourselves to this simplest set-up.
In Fig. 1 we display the non-interacting energies as a function of M π L, given by E n ∈ {ω π,n +ω N,n }, {ω π,n +ω π,m +ω N,n+m }, · · · , where ω π,n ≡ M 2 π + (2π/L) 2 n 2 , and where the ellipsis indicates four-(or more) particle states. As described in the figure caption, we are interested in states that have the quantum numbers of a single nucleon, I(J P ) = 1 /2( 1 /2 + ). For this reason the state with a pion and nucleon both at rest does not contribute. This state only couples to the s-wave and thus has negative parity due to the intrinsic parity of the pion. Also apparent from Fig. 1 is that, for physical pion masses and realistic values of M π L, the Lüscher formalism only rigorously describes, at most, the first excited state. For E n > m N + 2M π , an extended three-particle formalism is required. This has recently been developed by one of us for three-pion states, and the extension to general three-particle systems is underway [28,29]. Because the three-particle formalism is not yet directly applicable to N π → N ππ, in this work we restrict attention to the two-particle formalism, but also apply it above threshold where the predictions suffer from systematic uncertainties. As we explain more in Sec. V, an important conclusion of our analysis is that finite-volume states in the vicinity of the Roper resonance can contribute significantly to R(T, t). Given that the Roper has a ∼ 40% branching fraction to N ππ [30], three-particle states certainly need to be included to offer reliable quantitative predictions in this region. However, barring delicate cancellations between two-and three-particle states, the qualitative conclusions presented here are expected to hold. Three-particle contributions may also be enhanced when the energy of an N π pair within N ππ is close to the delta resonance. Indeed, in the ChPT analysis of Ref.
[1] the delta resonance was also included and was found to reduce the value of R(T, T /2). Finally we stress that one can only use the Lüscher quantization condition with scattering amplitudes that take the unitary form of elastic two-particle scattering, Eq. (21). Thus, in this approximation we must also neglect the inelasticity in the two-particle scattering amplitude. [See also Footnote 6.] In the following subsections we present our prediction for the finite-volume energy gaps ∆E n . First, in Sec. III A, we give the quantization condition for general two-particle systems and show how it can be reduced to describe the N π states of interest. Then, in Sec. III B, we use the experimental phase-shift data to predict the finite-volume spectrum.

A. Reducing the quantization condition
The quantization condition for particles with spin is a straightforward generalization of Lüscher's original result. Indeed a wide class of generalizations are all described by the same basic form [11,12,[16][17][18][19][20][21][22][23][24][25][26][27]31] det Here M is the two-to-two scattering amplitude and F is a known geometric function. This result describes any twoparticle system, with any number of two-particle channels with identical or non-identical particles, degenerate or non-degenerate masses and with arbitrary spin. To describe a specific system one need only specify the exact definitions, and in particular the index space, for M and F . As a preliminary example we consider a system with one channel of two non-identical scalars with masses M π and m N . In this case both M and F have two sets of spherical harmonic indices. The scattering amplitude is a diagonal matrix in this space, whose entries are related in the standard way to scattering phase shifts. F , by contrast, has on-and off-diagonal entries. This encodes the mixing of partial waves due to the reduced symmetry of the box. F can be written as a sum-integraldifference [11,12,[16][17][18] where k = |k|,k = k/k, and the sum runs over all k = (2π/L)n, n ∈ Z 3 . In Eq. (18) an ultraviolet regulator is needed to make the quantity well defined. Since the sum and integral have the same ultraviolet divergence, a universal result is recovered as the regulator is removed. Here p is the magnitude of CM frame momentum for particles with energy E and masses M π and m N To incorporate spin in this system it is most straightforward to first work in the basis where the nucleon is polarized along some fixed direction in its CM frame. This new degree of freedom, denoted by σ, can be accommodated with two simple modifications. First, the amplitude gains an additional index, M = M ,m ,σ ; ,m,σ . Second, the kinematic matrix F is multiplied with a Kronecker delta, δ σ σ . This completely defines the scalarnucleon quantization condition. Indeed, the arbitraryspin quantization condition is given by simply multiplying the F matrices with Kronecker deltas [26,27].
Next, to connect with experimental phase shifts, it is convenient to change to the basis of total angular momentum, J, orbital angular momentum, , and azimuthal component of total angular momentum, µ. The basis change is effected by contracting both sets of indices with standard Clebsch-Gordan coefficients. The amplitude in the new basis can be written 6 Note that the conservation of orbital angular momentum is special to the meson-baryon system. Generally orbital angular momenta will mix, but in this case conservation of total angular momentum implies that could at most couple with ± 1. Since changing by one unit flips parity, this coupling vanishes and is conserved. F in the new basis is given by [25,26] We make one final simplification before introducing approximations. One can show that the imaginary parts of M −1 and F perfectly cancel in Eq. (17), giving We now reduce the quantization condition to a determinant of a finite-dimensional matrix by ignoring high partial waves. It turns out that, in the even-parity sector, we reach the simplest possible truncation by neglecting δ J, for ≥ 3. Then the system is truncated to the = 1 space. In this space F J µ ;J µ is a 6 × 6 matrix: a 2 × 2 block for = 1, J = 1/2, a 4 × 4 block for = 1, J = 3/2. To determine its specific form we first note that where Z 00 is the Lüscher zeta-function described in Ref. [12] and q ≡ pL/(2π). The fact that F is proportional to the identity matrix in the = 1 subspace is preserved when we change to the J basis. Thus, both matrices in the quantization condition are diagonal and the final result is two independent, one dimensional equations for J = 1/2 or 3/2. These can be reexpressed as where n is an arbitrary integer and We comment that this has the same form as the s-wave, scalar quantization condition. The quantity φ is often referred to as the pseudophase. The fact that the J = 1/2 and J = 3/2 sectors decouple can be explained by examining the symmetry group of the finite-volume system. For the case of one scalar and one spin-half particle in a finite cubic box with zero total momentum, the symmetry group is 2 O ⊗ S 2 and the irreps are G ± 1 , G ± 2 and H ± [31]. If we neglect ≥ 3 and thus also neglect J ≥ 5/2, then we find a perfect correspondence between finite-and infinite-volume irreps G − 1 ≡ (J = 1/2) and H − ≡ (J = 3/2). This implies that, within this approximation, the two partial-waves cannot mix, as we have seen by explicit calculation.

B. Predicting the spectrum from the experimental N π phase shift
To predict the finite-volume spectrum of N π states we use experimental data made available by the George Washington University Institute for Nuclear Studies. Their data analysis center is available online at http://gwdac.phys.gwu.edu/analysis/pin_ analysis.html. In this study we use their partial wave analysis WI08 solution. The relevant phase shift data are plotted in Fig. 2. For detailed information about the experimental data set and the WI08 fit solution see Refs. [32,33].
Substituting this phase shift into Eq. (26) we reach the prediction for the two-particle energies, shown in Fig. 3. Note that, relative to the gap to the single nucleon state, the shift is relatively small between free-and interacting levels. This means that it makes little difference whether one uses the free or interacting finite-volume spectrum for the values of ∆E n that enter R(T, t). 7 Also apparent from Fig. 3 is that no avoided level crossing is visible. This is because the Roper resonance is too broad to generate such an effect. It follows that, near the physical point, no direct association between LQCD energies and the resonance can be made and a careful Lüscher based analysis will be needed to extract resonance properties from LQCD. 7 In this work we do not plot an explicit comparison but, as we comment in Sec. V below, if one uses LO ChPT for the infinitevolume matrix elements, then the effect of interactions in the energies and Lellouch-Lüscher factors affects the prediction for R(T, T /2) at the percent level. To better understand these results consider the form of the pseudophase curves, plotted together with the experimental phase shift for M π L = 4 in Fig. 4. The interacting energies, at this L value, are given by the intersections of the curves. This shows that there are universal features for the levels predicted by certain types of phase shifts. In particular, for any phase shift that slowly rises from 0 to π, the spectrum is given by a smooth deformation of the free levels. When δ(p) is near 0 or π the energies coincide with free values. As one follows a given interacting level from high energies to low (by increasing M π L) it rises by one free level. This implies that, for any slowly rising phase shift, interacting levels tend to be separated FIG. 4. Experimental scattering phase together with Lüscher pseudophase curves, (nπ − φ(q)), for MπL = 4. Each intersection gives an interacting level in terms of q 2 , which can be converted to energy via from their neighbors on each side by levels of the free theory. Also, the rise of the phase shift from 0 to π results in exactly one additional energy level relative to the free theory.
Finally, as we have already stressed above, in this prediction of the energy levels we neglect the effects of crossing three-particle production threshold. Roughly speaking crossing this threshold has two effects. First, threeparticle states appear on top of the two-particle states shown in the figure. Second, the positions of all energies are modified relative to those predicted by the twoparticle Lüscher formula. Strictly it does not make sense to distinguish between two-and three-particle states. All finite-volume states will have both two-and threeparticle components once the energy exceeds 2M π + m N . However, for sufficiently weak two-to-three couplings, the levels are well described by being two-particle or threeparticle like in certain regions, with avoided level crossings occurring whenever a given level changes from one type to the other. The overlap of the interpolator on a given state is also expected to be suppressed when the state has a large three-particle component, possibility with the exception of energies near the Roper. These observations, and the limitation of the formalism available, motivate us to use the effective spectrum plotted in Fig. 3 in our study of excited-state contamination.

IV. ESTIMATING THE FINITE-VOLUME MATRIX ELEMENTS
In this section we use experimental scattering data, together with LO ChPT and a model to describe N π final-state interactions, in order to estimate the finitevolume matrix elements entering b n and c n . The finitevolume two-particle states, denoted by |n , arise from insertions of the identity in Eqs. (5) and (6) and always appear as an outer product of the form |n n|. This is exactly the structure that is readily accommodated by the generalized Lellouch-Lüscher formalism, as we now describe.
The original work by Lellouch and Lüscher gives the relation between a finite-volume matrix element and the K → ππ decay rate [13] ππ, E n , L|H(0)|K, where M K is the kaon mass, H is the weak hamiltonian density in position space, H out ≡ ππ, out|H(0)|K is the corresponding infinite-volume matrix element and where φ is defined in Eq. (27) above and δ ππ is the s-wave phase-shift for elastic pion scattering. The finite-volume states on the right-hand side of Eq. (28) are unit normalized, whereas the infinite-volume states within H out satisfy the standard relativistic normalization. In the derivation of Lellouch and Lüscher the box size must be tuned so that the two-pion and kaon states are degenerate, E n = M K . As has become increasingly clear through various subsequent studies, see for example Refs. [13,17,18,22,23,27,[34][35][36][37], the conversion factor, |R|, can be understood as a relation between the two-particle states defined in finite and infinite volume. This means that essentially the same relation holds even after relaxing a number of assumptions going into the original derivation. In particular one can define a Lellouch-Lüscher like relation for operators with generic quantum numbers and with nonzero energy and momentum. 8 For products of matrix elements that involve an outerproduct of finite-volume states, the relation can be derived without taking magnitudes of matrix elements [37]. Thus, information about the relative sign of transition processes can also be obtained.
In the context of this study, the relevant relation, given in Ref. [27], takes the form L) is a matrix generalization of |R|, defined in the following subsection. This is the conversion needed for b n . The analog for c 2,n is given by The key limitation of Eqs. (30) and (31) is that these only hold for E n < 2M π + m N . For such two-particle energies the relation is valid up to exponentially suppressed corrections of the form e −MπL , but above three-particle threshold a generalized form with three-particle matrix elements is required. As in the previous section, here we again apply the two-particle formalism outside of its region of validity. We expect this to give a qualitative indication of the nature of excited-state contamination, but only by applying a rigorous three-particle formalism can one reach a reliable quantitative prediction, especially in the vicinity of the Roper.
Applying Eqs. (30) and (31), with R(E n , L) determined using experimental scattering data, it remains only to analyze the matrix elements of the nucleon interpolating operator, O, and of the axial-vector current, A µ , in infinite volume. In this way, the details of the finite-volume set-up are factored out. To explain this in detail we find it convenient to assume specific forms for the projectors entering Eqs. (3) and (4). In particular we take Γ = u + (0)ū + (0)/m N and Γ µ = δ 3µ 2iΓ, were u and u are standard nucleon spinors, already used in Eq. (1). One then finds where Here the first factor, B(E n ), is understood as a rowvector on the J, , µ index space labeling the two-particle state. It depends on the spin-projected interpolator as well as the kinematic factors ω π and ω N , evaluated at momentum p as defined in Eq. (20). The middle factor, C(E n , L), is a matrix on the J, , µ space. We discuss its definition in detail in the following subsection. Finally A(E n ), understood as a column on the same index space, is the infinite-volume axial-vector transition amplitude. We comment that all three of the factors entering Eq. (32) are dimensionless, real functions. The latter claim holds due to the diagonal matrices e −iδ included in the definitions. Here δ is a diagonal matrix of N π scattering phase shifts. For example, the J = 1/2, = 1 entry is plotted in Fig. 2. Watson's theorem states that the complex phase of a two-particle matrix element (below threeparticle production threshold) is equal to the elastic N π scattering phase in the same channel [38]. Thus the phase matrices in the definitions cancel those in the infinitevolume matrix elements. Above three-particle threshold this no longer holds, but in this work we model the matrix elements with a form satisfying this two-particle unitarity constraint. In other words we build in the approximation that Watson's theorem persists above threshold.
[See Sec. IV C for details.] Similarly the factors of e iδ in Eq. (35) cancel the intrinsic phase in R(E n , L) as we show in the next section. This ensures that C(E n , L) is also a real function.
In the following subsection we give the matrix definition of R(E n , L) and C(E n , L) and explain that, in the present case, one can truncate these to single entries by applying the same truncation used for the scattering amplitude in the previous section. In Sec. IV B we then use experimental scattering data to calculate the interacting values of C(E n , L). Finally in Sec. IV C we use a model, based in LO ChPT supplemented by the experimental scattering data, to estimate both B(E) and A(E). We then apply these results in Sec. V, to give predictions for the excited-state contamination to g A .

A. Reducing the Lellouch-Lüscher-like relation
We begin this subsection by defining R(E n , L), introduced in Eq. (30) above. In this equation the right-hand side should be understood as the product of a column vector 0|O β (0)|N π, in , followed by the matrix R(E n , L), followed by a row vector N π, out|A 3 µ (0)|N . Each of these quantities is defined on the J, , µ index space, where the three labels correspond to total angular momentum, orbital angular momentum, and azimuthal total angular momentum respectively. The matrix R(E n , L) is defined by with M and F defined in Eqs. (21) and (22) respectively. R has both on-and off-diagonal elements and, in the context of Eq. (30), gives a linear combination of infinitevolume matrix elements that equals a particular finite-volume matrix element. The same matrix structure holds in Eq. (32).
To truncate R we first observe that the operator A 3 3 (0) acting on the infinite-volume single-nucleon state generates a state which couples to both J = 1/2 and J = 3/2. In the corresponding finite-volume matrix element this state couples to two-particle finite-volume states in the G − 1 = 1/2 ⊕ · · · and H − = 3/2 ⊕ · · · representations. Thus, if we choose the two-particle state to transform in the G − 1 irrep, then the right-hand sides of Eqs. (30) and (32) will contain one term, depending on the J = 1/2 two-particle scattering state Given this truncation we are left only to determine the on-diagonal J = 1/2, = 1 entry of R. In principle this single entry depends on the full matrix structure of F −1 and M, since they enter via a matrix inverse. However, if we apply the p-wave truncation on M, as in the previous section, then M and F both truncate to single entry matrices. We find [13] R where δ(p) = δ J=1/2, =1 (p), is the N π phase shift, shown in Fig. 2.
To understand the phase in R, we recall from Watson's theorem that, at energies where only two-particle elastic scattering can occur, the complex phase of zero-to-two and one-to-two transition amplitudes is given by the twoto-two strong scattering phase [38]. Thus the phase in R perfectly cancels the phase in the matrix element We conclude by discussing the rescaled quantity C(E n , L), defined in Eq. (35). Substituting Eq. (42) into the definition and simplifying, we reach In next section we will also be interested in the noninteracting limit, and thus define where q ≡ pL/(2π) was already introduced above. Note that in Eqs. (44) and (45) we have implicitly extended the definition of C(E, L) to all energies. As is clear from Eqs. (30) and (32), the quantity only has physical meaning when evaluated at the energies of the finitevolume spectrum. However, understanding the continuous form of the function is useful for predicting how C(E n , L) will vary with the strength of the particle interaction.

B. Predicting the Lellouch-Lüscher factors
In this section we give numerical predictions for the values of C(E n , L) and in doing so also build some intuition about the meaning of this quantity.
We begin with the non-interacting version, C NI . This is plotted in Fig. 5 as a function of the dimensionless squared momentum, q 2 . The energies for which this curve has physical meaning correspond to q 2 = n 2 with n ∈ Z 3 . At these values our rescaled Lellouch-Lüscher factor takes on particularly simple values where ν n is the degeneracy of the nth state, equivalently the number of integer vectors that satisfy n 2 = n. The first few values of ν n are given in Table I. These degeneracies are also indicated by the horizontal tick marks crossing each vertical line in Fig. 5. We next turn to the interacting values of C(E n , L), plotted in Fig. 6. In contrast to C NI , the factor in the interacting case depends independently on E and L. Here We also indicate the two-particle energies corresponding with certain q 2 values. The interacting curve (blue) is lower than the non-interacting (gray) due to a shift proportional to dδ(p)/dp in the inverse. For ease of comparison we have plotted the phase shift over the same range for each MπL value. However, the most important effect for the interacting value of C is the shift in energies. Since the curve is rapidly oscillating, these shifts result in large discrepancies between the interacting and non-interacting C values.  we plot C(E, L) as a function of q 2 at fixed M π L = 4 [ Fig. 6(a)] and M π L = 5 [ Fig. 6(b)]. Note that two effects are important in the shift from non-interacting to interacting systems. First, the curve characterizing the Lellouch-Lüscher factors is reduced, due to the addition of a term proportional to dδ/dp in the inverse. Comparing the light gray and blue curves, we see that this has a relatively small numerical effect. Second, the physically relevant values on the curve (the finite-volume energies) are shifted to new locations. Interestingly, this second shift has a large effect for both M π L = 4 and M π L = 5 (assuming physical pion masses). In particular, the contribution to R(T, t) from the seventh excited state is significantly enhanced due to the large value of C(E n , L).

C. Modeling the matrix elements
We now turn to the remaining building blocks entering the definitions of b n and c n : the overlap factor B(E n ) and the axial-vector transition amplitude A(E n ).
As already mentioned above, the matrix elements within B(E) depend on the details of the interpolator, O + , and cannot be constrained using experimental data. The precise definition of O + depends on the set-up of a particular lattice calculation, and the only clear defining property is that it has the quantum numbers to annihilate a single nucleon. In a variational approach, this operator is designed to reduce the size of the two-particle matrix elements shown in Eq. (30). Thus by choosing a sufficiently large basis one can make this contribution arbitrarily small. 9 By contrast, A(E) is a physical matrix element that can in principle be accessed in a scattering experiment. In fact, since we are focusing on the case where both the nucleon and the N π state are at rest, parity guarantees that the corresponding matrix element with a vector current must vanish. Thus we can equally well think of A(E) as the matrix element of the left-handed (vector-minusaxial) current. The kinematics of A(E) correspond to a process such as p π − → p W − → p e −νe , in which the current is evaluated at time-like four-momentum. Such a transition is difficult to extract in experiment and the present data is insufficient to constrain the amplitude.
We also note that the value of this matrix element at space-like momenta is of great experimental relevance for determining QCD effects in neutrino-nucleon scattering, ν p → pπ + − . 10 In this work we rely on LO ChPT, together with a model that incorporates the experimental N π scattering data, in order to gain insight into the behavior of both the interpolator overlap B(E) and the axial-vector transition amplitude A(E). Beginning with B(E), we first give the expression predicted by LO covariant baryon ChPT, derived in the appendix

whereḡ
is a convenient shorthand introduced in Ref. [2], and f π = 93 MeV is the pion decay constant. This result predicts the part of b n that depends on O to be negative and have magnitude ∼ 10 −3 . In addition, as was already pointed out in Refs. [1][2][3], the leading-order prediction is independent of the details of the interpolator used. Again, we stress that this only holds for a three-quark interpolating field and, in particular, does not apply to any interpolator built from multi-particle operators, for example N π-or N ππ-like operators. Eq. (47) is expected to break down at higher energies, when next-to-leading-order ChPT corrections become important. For instance, if O is optimized to couple to a single-nucleon state then, depending on the nature of the Roper, the overlap of the operator with two-particle states may also be enhanced in the vicinity of the resonance. This enhancement is not visible in LO ChPT.
To model the effect of the Roper we first consider the Bethe-Salpeter equation for the two-particle matrix element within B(E) where the subscript 2PI refers to the sum of all diagrams that are two-particle irreducible in the s-channel. To reach the full matrix element, the 2PI quantity must be attached, via a two-particle loop, to the two-to-two scattering amplitude, M. In the loop both the pion [∆(k)] and nucleon [S(P −k)] propagators should be fully dressed and evaluated at the off-shell momenta sampled by the integral. The scattering amplitude is also sampled at off-shell momenta. 10 See for example Ref. [39].
We cannot evaluate this quantity without making a number of approximations. First we evaluate the k 0 integral, approximating the matrix element and M to have no important analytic structure, so that we need only encircle the poles in the two-particle loop. This gives where S is a smooth function below three-particle production threshold. If we assume the dominant contribution comes form the first term, and drop the smooth part then we are left with a contribution in which both the matrix element and the scattering amplitude are projected on shell. We find where and where PV indicates a principal-value pole prescription. The subscript R indicates that this loop integral requires a regulator, an artifact that has been introduced by our approximations. In this work we choose to regulate by subtracting the integrand evaluated at threshold This subtraction is motivated by the observation that the second term in Eq. (49) should not play a role at low energies. Note also that the on-shell restriction projects the scattering amplitude down to its I(J P ) = 1 /2( 1 /2 + ) component.
To complete the construction of our model we use the fact that the diagrams in the LO ChPT calculation of B(E) are two-particle irreducible, and thus also give the leading-order contribution to the 2PI restriction of this quantity. This leads us to define where we have introduced the free parameter γ to partially compensate the ad hoc procedure that lead us to this expression. Here we have also included the phase factor e −iδ that is needed to cancel the phase that ap-pears in the two-particle matrix element. 11 To evaluate both the e −iδ factor and the scattering amplitude M we use the experimentally determined N π scattering phase, plotted in Fig. 2. Note also that we must discard a small imaginary part that arises in this model.
In the case of γ = 0 we omit the e −iδ phase factor and thus recover the leading order ChPT result. For γ > 0 the rising phase shift causes the matrix element to flip sign roughly in the region of the Roper, with the energy value at the node dependent on the specific value chosen for γ. For γ < 0 the sign of the ChPT prediction is preserved and a peak is observed in the vicinity of the Roper. Past this energy range the matrix element can flip sign, but for γ < −3 the crossing is well outside the relevant energy range. In Fig. 7 we plot the energy dependence predicted by various values of γ.
We now turn to the matrix element of the axial-vector current, A(E). As we derive in the appendix, the LO ChPT prediction for this quantity is To estimate this quantity beyond ChPT we apply the same model used for the overlap factor (56) As with B(E, γ), α = 0 gives the LO ChPT prediction, α < 0 preserves the sign of the matrix element and enhances the magnitude near the resonance, and α > 0 gives a zero-crossing roughly in the range of the resonance. Also as with B(E, γ), we include the phase factor needed to cancel the phase in the matrix element, and discard a small imaginary contribution that arises as an artifact of our model. For α = 0 the phase factor is not included.
In Fig. 8 we plot the energy dependence of the axial transition amplitude for various choices of α. Note that we restrict attention to a range of α that is smaller than that considered for γ. The LO ChPT prediction for B(E) has a magnitude that decreases with energy whereas A(E) is nearly constant at higher energies. This has the consequence that varying α over a given range has a larger effect than varying γ. We choose the parameters such that the models have a maximum magnitude roughly within a factor of two of the maximum predicted by LO ChPT. 11 This simple phase structure is only strictly valid below threeparticle production threshold. In the present model we are using the elastic form for the two-particle scattering amplitude and this has the consequence that the phase is preserved also above production threshold. This is consistent with the neglect of threeparticle states in the Lüscher quantization condition. We stress again that, unlike with the overlap factor, the functional form of A(E) is independent of the lattice set-up and is of direct experimental relevance. In this study we are particularly interested in whether A(E) has a node (zero crossing) at some energy, as predicted by the α > 0 models. Interestingly, such a node is observed in the CLAS scattering data for ep → e π + π − p in their analysis of the electromagnetic transition amplitude as a function of photon virtuality, Q 2 [40]. This node is not directly relevant for A(E) because (1) it concerns the electromagnetic transition and (2) it is for space-like momenta. It is nonetheless interesting to note that such crossings are observed.
Given that LO ChPT predicts the same sign for B(E) and A(E), and thus a positive value for b n = B(E n )C(E n , L)A(E n ), and given also that the curvature in LQCD correlator data indicates important contributions from states with b n < 0, we postulate that a node in A(E) might provide a reasonable explanation for the apparent discrepancy. More generally one can identify four basic scenarios: (i) neither B(E) nor A(E) cross zero in the relevant energy range, (ii) both cross zero, (iii) only the overlap factor B(E) has a node, or (iv) only the transition amplitude A(E) has a node. The first two cases lead to positive excited-state contamination and thus fail to describe present day numerical LQCD correlator data. The third and fourth scenarios can both explain the empirically observed excited state-contamination, as we explain in the next section.

V. ESTIMATING THE EXCITED-STATE CONTAMINATION
We are now ready to combine the results of the previous section to estimate the ratio R(T, t) and the excitedstate contamination to g A .
In Fig. 9 we show values for the ratio R(T, T /2)/g A , given three different scenarios for the matrix elements. In each case we show the results for both M π L = 4 (solid lines) and M π L = 6 (dashed lines). To provide comparable predictions, in both cases we sum all excited states up to an energy of 1800 MeV. In the case of M π L = 4 this corresponds to the first 8 excited states and for M π L = 6 to the first 18. In both cases we find that this number of states is both necessary and sufficient to estimate the saturated value of R(T, T /2)/g A within a few percent, for the models considered. We also see that the excited state contaminations from the two different volumes are in very good agreement.
The highest pair of curves in Fig. 9 shows the prediction from LO ChPT, but with the interacting values of the Lellouch-Lüscher factors. The excited-state contamination here is comparable to that given in Fig. 5 of Ref. [3], in particular the result of that reference that includes ten excited states. Thus we find that, if one uses LO ChPT for the overlap and the axial-vector matrix element, then the effect of interactions on the energies and Lellouch-Lüscher factors leads to a small (percent level) correction to the predicted value of R(T, T /2).
As is also stressed in Ref. [3], including excited states beyond the first few requires sampling the LO ChPT predictions outside their expected region of validity. Roper and, as the latter is not included in the ChPT prediction, it is reasonable to require a separation from this state. We thus infer that one can only predict R(T, T /2) using LO ChPT for source-sink separations large enough that the first two states dominate. For physical pion masses this may mean separations of T > 2 fm.
The middle pair of curves in Fig. 9 is the prediction from combining the LO ChPT prediction for the overlap, B(E, γ = 0), with the zero-crossing model for the axialcurrent transition amplitude, A(E, α = 1). The curve is very flat because there are large cancellations between the positive contributions from lower states and the negative contributions from higher-energy states.
The lowest pair of curves in Fig. 9 gives the scenario most consistent with observed LQCD correlators. This follows from combining A(E, α = 1) and B(E, γ = −4) [see again Figs. 7 and 8]. The negative contribution from the higher excited states overpowers the positive contribution from the first few.
In Fig. 10 we show the importance of higher excited states in the A(E, α = 1) and B(E, γ = −4) scenario. In particular, for M π L = 4, we find that the contamination predicted by summing fewer than seven states differs significantly from the saturated curve. In particular the seventh excited state has a significant contribution due to the large value of C(E n , L). The thickness of a given curve is proportional to the number of states included in the sum. The three blue curves show the result from including the first (lowest blue), first two (middle blue) and first three (highest blue) excited states. The fourth excited state is the first with bn < 0 so that the predicted contamination falls once this state is included (highest green curve). The set of green curves then indicates the sum up to the fourth state (highest green) to the sum up to the seventh state (lowest green). The effects of including additional states beyond the seventh are negligible so that the lowest green curve gives a good indication of the full excited state contamination predicted by this model.

VI. CONCLUSION
In this work we have studied the excited-state contamination in LQCD correlators used to extract the nucleon axial charge, g A . Combining various finite-volume formalisms with experimental scattering data, LO ChPT and a model for the infinite-volume matrix elements, we find that the excited-state behavior empirically observed in lattice correlators can be reproduced by postulating a sign change in the infinite-volume axial-vector transition amplitude, N π, out|A|N . Such nodes are observed experimentally in other transition amplitudes, but the data is insufficient to make a definitive statement about the quantity at hand.
Our findings additionally indicate that a large number of finite-volume excited states, including those at energies around the Roper resonance, give important contributions in the ratios used to access g A . This is based on mild assumptions about how the nucleon interpolators couple to states near the Roper and on the observation that the Lellouch-Lüscher factors, governing the relation between finite-and infinite-volume states, can be signif-icantly enhanced.
The results presented here serve to further emphasize the great importance of using optimal interpolators to minimize the coupling to excited states. Based on numerical LQCD calculations in the meson sector, the most promising approach seems to be the variational method, in which a large basis of operators is used to disentangle the excited states. The situation will also be improved by further advances in improving the signal-to-noise ratio in nucleon correlators.
Finally we emphasize that the nature of excited-state contamination depends heavily on the quantity under consideration. Indeed for many quantities, for example average x, the LQCD data indicates positive excitedstate contamination [6]. For this observable the same overlap factor B(E) appears, but a different transitionamplitude arises due to the differing current insertion. The observed, positive excited-state contamination can be accommodated if we suppose the matrix element has the same sign as the axial vector at low-energies and does not cross zero in the relevant energy window.
Another interesting example is the iso-singlet octet axial-vector. ChPT predictions indicate that matrix elements of this current should be highly suppressed relative to those of the iso-triplet studied here. The iso-singlet suffers from other sources of systematic uncertainty, in particular quark-disconnected diagrams. Given the potential severity of excited-state contaminations, it will be interesting to compare the systematic error budgets for these quantities as methods on both sides improve.
Working in time-momentum perturbation theory Ref. [2] found where Here we have repeated the definition ofḡ A , given in Eq. (48) above, but with ω N + ω π in place of E n . This is equivalent given the definition of the individual energies where p is defined via These definitions are needed to extend the results away from the non-interacting finite-volume energy levels.
Eq. (A4) also depends onα, a low-energy coefficient that was also introduced in Ref. [2]. This coefficient enters the relation between the lattice interpolator and the ChPT fields To cross check the result of Ref. [2] we apply the extension of the Lellouch-Lüscher approach to the relevant finite-volume matrix elements. From the discussion in the main text follows Here we have used the non-interacting value for the Lellouch-Lüscher factor as is appropriate for a leadingorder calculation. The two-particle asymptotic states on the second line are projected to definite isospin I = 1/2, m I = 1/2 definite total angular momentum J = 1/2, µ = 1/2 and definite parity P = +. This implies that orbital angular momentum is restricted to the p-wave.
To reach Eq. (A5), one must calculate the infinitevolume matrix elements in Eq. (A11) at tree level. The relevant two diagrams are shown in Fig. 11. The first of these diagrams, Fig. 11(a), does not include any interactions from the Lagrangian but instead arises from the N π contribution in the ChPT expression for the lattice interpolator [2] O ⊃α i 2f π γ 5 1 0 This is then combined with the definite isospin state 12 to reach Here we have written the right-most spinor, evaluated at the nucleon momentum p, as a boosted zero-momentum spinor. We have introduced η = sinh −1 (p/m N ) as the rapidity of the nucleon, and K i = −iγ 0 γ i /2 as the generator of the boost. In this appendix we use the Minkowski metric (including Minkowksi gamma matrices), following the conventions of Ref. [41]. As indicated by the + label on the N π state (and on the right-most spinor), this is the result in the case that the incoming nucleon is spin up.
When the incoming nucleon is spin down the spinors pick off a different combination of gamma matrices leading to a different combination of the components ofp. The result is Next performing the orbital angular-momentum projection (simply by dropping the √ 4πY * 1m factor) and forming the definite J combination we conclude the final contribution from Fig. 11(a) We now turn to Fig. 11(b). This diagram arises from the single nucleon contribution to the interpolator together with the N N π coupling in the Lagrangian Inferring the momentum-space Feynman rules and calculating the spin-up matrix element, we find We again substitute Eq. (A16), but in this case both the cosh and sinh terms contribute giving 0|O + (0)|N π, +, in ⊃ −iα √ 3m N g A f π 1 ω N + ω π − m N × cosh(η/2)p + sinh(η/2)ω π p 3 . (A24) Applying standard identities for the trigonometric functions one finds 0|O + (0)|N π, +, in ⊃ − iα √ 2f π g A (ω N − m N ) 1/2 ω N + ω π + m N ω N + ω π − m N √ 4πY * 10 (p) , In words, the contribution to the spin-up state from Fig. 11(b) is just given by the result for 11(a) multiplied by (−ḡ A ). The same relation holds for the spin-down state so that the full result for the J = 1/2 matrix element is Substituting these results into Eq. (A11) we recover Eq. (A5).
We now turn to the three-point function and in particular the ratio used to extract g A . Applying our specific choices of Γ and Γ to the definition of b n we deduce = 2i ν n 4ω π ω N L 3 0|O + (0)|N π, in N π, out|A 3 3 (0)|N 0|O + (0)|N , Here we have again applied the Lellouch-Lüscher formalism and have also substituted Eq. (A26).
To complete the calculation it remains only to substitute the result for the infinite-volume matrix element N π, out|A 3 3 (0)|N . Again this must be evaluated at tree level by calculating the three diagrams shown in Fig. 12. These diagrams include various contributions from the ChPT expression for the axial-vector current (A31) Here we are using the Euclidean convention for the axialvector current as is used throughout the main text. However, we are expressing it in terms of the Minkowski convention γ 3 to match the conventions of other quantities used in this appendix. Beginning with Fig. 12(a), we first evaluate the isospin projected result with the spin up N π state N π, +, out|A 3 The contribution from the spin down state has a similar form N π, −, out|A 3 3 (0)|N and the final result for the Fig. 12(a) contribution to the J = 1/2 matrix element is Next evaluating Fig. 12(b) we find N π, +, out|A 3 for the spin-up matrix element, for the spin-down matrix element, and 6f π (A40) for the J = 1/2 matrix element.
We conclude with Fig. 12(c), jumping straight to the projected, J = 1/2 result Combining all results, we conclude in perfect agreement with Ref. [3].