The pion-nucleon sigma term from lattice QCD

We present an analysis of the pion-nucleon $\sigma$-term, $\sigma_{\pi N}$, using six ensembles with 2+1+1-flavor highly improved staggered quark action generated by the MILC collaboration. The most serious systematic effect in lattice calculations of nucleon correlation functions is the contribution of excited states. We estimate these using chiral perturbation theory ($\chi$PT), and show that the leading contribution to the isoscalar scalar charge comes from N$\pi$ and N$\pi\pi$ states. Therefore, we carry out two analyses of lattice data to remove excited-state contamination, the standard one and a new one including N$\pi$ and N$\pi\pi$ states. We find that the standard analysis gives $\sigma_{\pi N} = 41.9(4.9)$ MeV, consistent with previous lattice calculations, while our preferred $\chi$PT-motivated analysis gives $\sigma_{\pi N} = 59.6(7.4)$ MeV, which is consistent with phenomenological values obtained using $\pi$N scattering data. Our data on one physical pion mass ensemble was crucial for exposing this difference, therefore, calculations on additional physical mass ensembles are needed to confirm our result and resolve the tension between lattice QCD and phenomenology.


I. INTRODUCTION
This Letter presents results for the pion-nucleon σterm, σ πN ≡ m ud g u+d S ≡ m ud N (k, s)|ūu +dd|N (k, s) calculated in the isospin symmetric limit with m ud = (m u + m d )/2 the average of the light quark masses. It is a fundamental parameter of QCD that quantifies the amount of the nucleon mass generated by the u-and dquarks. The scalar charge g q S is determined from the forward matrix element of the scalar densityqq between the nucleon state: g q S = N (k = 0, s)|Z Sq q|N (k = 0, s) , where Z S is the renormalization constant and the nucleon spinor has unit normalization. The connection between g q S and the rate of variation of the nucleon mass, M N , with the mass of quark with flavor q is given by the Feynman-Hellmann (FH) relation [1][2][3] ∂M N ∂m q = N (k, s)|qq|N (k, s) = g q S /Z S .
The analytic continuation can be further improved in the framework of Roy-Steiner equations [37][38][39][40][41][42][43][44][45], whose constraints on σ πN become most powerful when combined with pionic-atom data on threshold πN scattering [46][47][48][49][50]. Slightly updating the result from Refs. [39,41] to account for the latest data on the pionic hydrogen width [48], one finds σ πN = 59.0(3.5) MeV. In particular, this determination includes isospin-breaking corrections [51][52][53][54] to ensure that σ πN coincides with its definition in lattice QCD calculations [42]. The difference from Refs. [19][20][21] traces back to the scattering lengths implied by Refs. [22,23], which are incompatible with the modern pionic-atom data. Independent constraints from experiment are provided by low-energy πN crosssections, including more recent data on both the elastic reactions [55][56][57] and the charge exchange [58][59][60][61], and a global analysis of low-energy data in the Roy-Steiner framework leads to σ πN = 58 (5) MeV [45], in perfect agreement with the pionic-atom result. In contrast, so far lattice QCD calculations [62][63][64][65][66][67][68][69][70] have favored low values σ πN ≈ 40 MeV (with the exception of Ref. [71]), and it is this persistent tension with phenomenology that we aim to address in this Letter. There are two ways to calculate σ πN using lattice QCD, which are called the FH and the direct methods [72]. In the FH method, the nucleon mass is obtained as a function of the bare quark mass m ud (equivalently M 2 π ) from the nucleon 2-point correlation function, and its numerical derivative multiplied by m ud gives σ πN . In the direct method, the matrix element ofūu +dd is calculated within the ground-state nucleon. Both methods have their challenges. In the FH method, one needs to calculate the derivative about the physical m ud , which is computationally very demanding. Most calculations extrapolate from heavier masses or fit the data for M N versus M 2 π to an ansatz motivated by χPT and evaluate its derivative at m ud . On the other hand, the signal in the matrix element is noisier since it is obtained from a 3-point function with the insertion of the scalar density. In both methods, one has to ensure that all excited-states contamination (ESC) has been removed. Both methods give σ πN ≈ 40 MeV-see Fig. 4, review by the Flavour Lattice Averaging Group (FLAG) in 2019 [72], and the two subsequent works [69,70].
Here, we present a new direct-method calculation. Our main message is that N π and N ππ excited states, which have not been included in previous lattice calculations, can make a significant contribution. We provide motivation for this effect from heavy-baryon χPT [73,74], and show that including the excited states in fits to the spec-tral decomposition of the 3-point function increases the result by about 50%. Such a change brings the lattice result in agreement with the phenomenological value.

II. LATTICE METHODOLOGY AND EXCITED STATES
The construction of all nucleon 2-and 3-point correlations functions is carried out using Wilson-clover fermions on six 2+1+1-flavor ensembles generated using the highly improved staggered quark (HISQ) action [75] by the MILC collaboration [76]. In each of these ensembles, the u-and d-quark masses are degenerate, and the s-and c-quark masses have been tuned to their physical values. Details of the six ensembles at lattice spacings, a ≈ 0.12, 0.09, and 0.06 fm, and M π ≈ 315, 230, and 138 MeV are given in Table I and in Table II, and of the analysis in App. A. To obtain flavor-diagonal charges g q S , two kinds of diagrams, called connected and disconnected and illustrated in Fig. 1, are calculated. The details of the methodology for the calculation of the connected contributions (isovector charges) using this clover-on-HISQ formulation are given in Refs. [77,78] and of the disconnected ones in Ref. [77].
The main focus of the analysis is on controlling the ESC. To this end, we estimate σ πN using two possible sets of excited-state masses, M 1 and M 2 , given in Table I. These M i are obtained from simultaneous fits to the zero momentum nucleon 2-point, C 2pt , and 3-point, C 3pt , functions using their spectral decomposition truncated to four and three states respectively: Here A i are the amplitudes for the creation or annihilation of states by the nucleon interpolating operator used on the lattice, N = abc u aT Cγ 5 (1 + γ 4 )d b u c , with color indices {a, b, c} and charge conjugation matrix C. The nucleon source-sink separation is labeled by τ and the operator insertion time by t. The issue of ESC arises because N couples not only to the ground-state nucleon but to all its excitations including multihadron states with the same quantum numbers. In the current data, the signal in C 3pt S extends to τ ≈ 1.5 fm, at which source-sink separation the contribution of excited states is significant as evident from the dependence on {τ, t} in the ratio R S (τ, t) = C 3pt S (t, τ )/C 2pt (τ ) shown in Fig. 2. In the limits t → ∞ and (τ − t) → ∞, the ratio R S (τ, t) → g S . Fits to C 3pt using Eq. (3) with the key parameters M i left as free parameters have large fluctuations. We, therefore, remove ESC and extract the ground-state matrix element, 0|S|0 , using simultaneous   fits to C 2pt and C 3pt with common M i . Statistical precision of the data allowed, without overparameterization, four states in C 2pt (labeled {4} or {4 N π }), and three states in C 3pt (labeled {3 * }). We also dropped the unresolved 2|S|2 term in Eq. (3). Keeping it increases the errors slightly but does not change the values. Using empirical Bayesian priors for M i and A i given in Table V, we calculate σ πN for two plausible but significantly different values of M 1 and M 2 in Table I that give fits with similar χ 2 . A similar strategy has been used in the analysis of axial-vector form factors, where also the N π state gives a large contribution as discussed in Refs. [79,80].
Data for C 3pt , by Eq. (3), should be (i) symmetric about τ /2, and (ii) converge monotonically in τ for sufficiently large τ , especially when a single excited state dominates. These two conditions are, within errors, satisfied by the data shown in Fig. 2. In the simultaneous fits, M 1 and M 2 are mainly controlled by the 4-state fits to C 2pt , however, as discussed in Refs. [78,80], there is a large region in M i>0 in which the augmented χ 2 of fits with different priors for M i is essentially the same, i.e., many M i are plausible. This region covers the towers of positive parity N π, N ππ, . . ., multihadron states, labeled by increasing relative momentum k, that can contribute and whose energies start below those of radial excitations. To obtain guidance on which excited states give large contributions to C 3pt , we carried out a χPT analysis.
We study two well-motivated values of M 1 and M 2 for the analysis of C 3pt . The "standard" strategy (called the {4, 3 * } fit) imposes wide priors on M i>0 , mostly to stabilize the fits, while the {4 N π , 3 * } fits use narrow-width priors for M 1 centered about the noninteracting energy of the almost degenerate lowest positive parity multihadron states, N (1)π(−1) or N (0)π(0)π(0). Thus, the label N π implies that the contribution of both states is included. Details of extracting the M i from these two four-state fits, {4} and {4 N π } to just C 2pt , can be found in Ref. [80].  Table I. The fits (see Fig. 2) and the χ 2 with respect to C 3pt data are equally good, however, the results for the isoscalar charge g u+d S differ significantly. The {4, 3 * } fit leads to a result consistent with σ πN ≈ 40 MeV, whereas the {4 N π , 3 * } fit gives ≈ 60 MeV. The major difference comes from the disconnected quark loop diagram shown in Fig. 1, and is strongly M π dependentthe effect of the N π states is hard to resolve in the M π ≈ 315 MeV data, debatable in the 230 MeV data, and clear in the M π = 138 MeV data.
It is important to point out that the values of M 1 and M 2 used in both fit strategies are an effective bundling of the many excited states that contribute into two. In fact, as mentioned above, many combinations of M 1 and M 2 between {4, 3 * } and {4 N π , 3 * } (see Table I) give fits with equally good χ 2 values. Ref. [80] showed that for τ 1.0 fm and for both fit strategies, the dominant ESC in C 3pt comes from the first excited state. Thus, operationally, our two results for σ πN should be regarded as: what happens if the first "effective" excited state has M 1 ≈ 1220 MeV (motivated by χPT and corresponding to the lowest theoretically possible states N π or N ππ) versus 1600 MeV obtained from the standard fit to the 2-point function. To further resolve all the excited states that contribute significantly and their energies in a finite   box requires much higher precision data on additional M π ≈ 135 MeV ensembles. In short, while our {4 N π , 3 * } analysis reconciles the lattice and the phenomenological values, it also calls for validation in future calculations.

III. EXCITED STATES IN χPT
The contributions of low-momentum N π and N ππ states to C 2pt and C 3pt can be studied in χPT [81][82][83][84][85][86][87][88], a low-energy effective field theory (EFT) of QCD that provides a systematic expansion of R S (τ, t) in powers of Q/Λ χ , where Q denotes a low-energy scale of order of the pion mass, GeV is the typical scale of QCD. In contrast to the isovector scalar charge considered in Ref. [83], we find large contributions from the N (k)π(−k) and N (0)π(k)π(−k) states, which can give up to 30% corrections to R S and thus affect the extraction of g u+d S and σ πN in a significant way.
The diagrams contributing to R S are shown in Fig. 5, where we assume N (x) to be a local nucleon source with well defined transformation properties under chiral symmetry. The chiral representation of this class of sources has been derived in Refs. [81,82,89]. Details of the calculation at next-to-next-to-leading order (N 2 LO) in χPT and the expansion of N in terms of heavy nucleon and pion fields are summarized in App. B. The crucial observation is that the isoscalar scalar source couples strongly to two pions, so that loop diagrams with the scalar source emitting two pions, which are consequently absorbed by the nucleon, are suppressed by only one chiral order, Q/Λ χ . These diagrams have both N π and N ππ cuts, which give rise to ESC to Euclidean Green's functions. A second important effect is that the next-to-leadingorder (NLO) couplings of the nucleon to two pions, parameterized in χPT by the low-energy constants (LECs) c 1,2,3 , are sizable, reflecting the enhancement by degrees of freedom related to the ∆(1232). When the pions couple to the isoscalar source, these couplings give rise to large N 2 LO corrections that are dominated by N ππ excited states and have the same sign as the NLO correction. Since, in the isospin-symmetric limit, the isovector scalar source does not couple to two pions, the NLO diagrams and the N 2 LO diagrams proportional to c 1,2,3 do not contribute to the isovector 3-point function, whose leading ESC arises at O(Q 2 /Λ 2 χ ). A detailed analysis showing that the functional form of the ESC predicted by χPT matches the lattice data and fits for sufficiently large time separations τ is given in App. B. In particular, the NLO and N 2 LO ESC can each reduce σ πN at a level of 10 MeV, thus explaining the {4 N π , 3 * } fits, i.e., a larger value when ESC is taken into account.

IV. ANALYSIS OF LATTICE DATA
Examples of fits with strategies {4, 3 * } and {4 N π , 3 * } to remove ESC and obtain g u+d,bare S are shown in Fig. 2 and the results summarized in Table I. The final results are obtained from fits to the sum of the connected and disconnected contributions. These values overlap in all cases with the sum of estimates from separate fits to g u+d, conn S and g 2l S . From the separate fits, we infer that most of the difference between the two ESC strategies comes from the disconnected diagrams, which we interpret as due to the N π/N ππ contributions through quark-level diagrams such as shown in Fig. 1. Figure 3 shows data for σ πN = m bare ud g u+d,bare S as a function of a and M 2 π . The chiral-continuum (CC) extrapolation is carried out using the N 2 LO χPT expression [41]: Neglected finite-volume corrections can also be estimated in χPT, see App. B and Refs. [78,90], indicating a correction of less than 1 MeV for the a09m130 ensemble. Figure 3 shows two chiral fits based on the N 2 LO χPT expression for σ πN . The {2, 2a, 3, 4} fit keeps terms proportional to {d 2 , d a 2 , d 3 , d 4 } with all coefficients free. In the {2, 2a, 3 χ , 4, 4L} fit we use the χPT value for d 3 = d χ 3 , which does not involve any LECs, and include the d 4L term. Each panel also shows the six data points obtained with the {4, 3 * } and {4 N π , 3 * } strategies and the fits to them. In each fit d a 2 comes out consistent with zero. The results for σ πN at the physical point M π = 135 MeV from the various fits are essentially given by the a09m130 point. We have neglected a correction due to flavor mixing inherent in Wilson-clover fermions since it is small as shown in App. D. Our final result, σ πN = 59.6(7.4) MeV, is the average of results from the {2, 2a, 3, 4} and {2, 2a, 3 χ , 4, 4L} fits to the {4 N π , 3 * } data given in Fig. 3, which overlap. In App. B, we consider more constrained fit variants, which show that the fit coefficients of the M 2 π and M 4 π log M 2 π terms are also broadly consistent with their χPT prediction.

V. CONCLUSIONS
Results for σ πN were reviewed by FLAG in 2019 [72], and there have been two new calculations since as summarized in App. E, and shown in Fig. 4. The ETM collaboration [69], using the direct method on one physical mass 2+1+1-flavor twisted mass clover-improved ensemble, obtained σ πN = 41.6(3.8) MeV; the BMW collaboration using the FH method and 33 ensembles of 1+1+1+1-flavor Wilson-clover fermions [70], but all with M π > 199 MeV, find σ πN = 37.4(5.1) MeV. The χPT analysis of the impact of low-lying excited N π states in the FH and direct methods is the same, and as shown in Fig. 3, it mainly affects the behavior for M π 135 MeV. Our work indicates that previous lattice calculations give the lower value σ πN ≈ 40 MeV because in the FH analysis [70] the fit ansatz (Taylor or Padé) parameters are determined using M π ≥ 199 MeV data, and in the direct method, the N π/N ππ states have not been included when extracting the ground-state matrix element [69].  Gasser   The BMW 20 result from 1+1+1+1-flavor lattices is listed along with the other 2+1+1-flavor calculations for brevity. Following the FLAG conventions, determinations via the direct approach are indicated by squares and the FH method by triangles. Also, the symbols used for lattice estimates that satisfy the FLAG criteria for inclusion in averages are filled green, and those not included are open red. The references from which lattice results have been taken are: JLQCD 18 [68], χQCD 15A [65], BMW 15 [64], ETM 14A [71], ETM 19 [69], and BMW 20 [70]. Phenomenological estimates using πN scattering data (blue filled circles) are from Gasser 91 [20], Pavan 02 [26], Alarcon 11 [28], Hoferichter 15 [39], and Ruiz de Elvira 17 [45].

MeV
To conclude, a χPT analysis shows that the low-lying N π and N ππ states can make a significant contribution to g u+d S . Including these states in our analysis (the {4 N π , 3 * } strategy) gives σ πN = 59.6(7.4) MeV, whereas the standard analysis ({4, 3 * } strategy) gives σ πN = 41.9(4.9) MeV consistent with previous analyses [72]. These chiral fits are shown in Fig. 3. Since the {4, 3 * } and {4 N π , 3 * } strategies to remove ESC are not distinguished by the χ 2 of the fits, we provide a detailed N 2 LO χPT analysis of ESC, which reveals sizable corrections consistent with the {4 N π , 3 * } analysis, restoring agreement with phenomenology. Since the effect of the N π and N ππ states becomes significant near M π = 135 MeV, further work on physical mass ensembles is needed to validate our result and to increase the precision in the extraction of the nucleon isoscalar scalar charge. The parameters of the seven 2+1+1-flavor ensembles generated using the HISQ action [75] by the MILC collaboration [76] are given in Table II. Note that the seventh ensemble, a09m310, listed has been used only for the analysis of the nucleon mass in App. C. On each of these ensembles, the calculations of the 2-and 3-point functions was carried out using tadpole improved Wilsonclover fermions as described in Ref. [78].
To reduce ESC, smeared sources using the Wuppertal method [92] were used to generate quark propagators . Columns 6-9 give the number of configurations analyzed, and the number of low-(high-) precision measurements made per configuration for the connected contributions. Columns 10-12 give the number of configurations, low-precision random sources used, and the ratio NLP/NHP for the disconnected contributions. The analysis of the isovector contributions has been presented in Ref. [78]. The ensemble a09m310 has been used only in the analysis of the nucleon mass in App. C.
with parameters given in Ref. [78]. The same smearing was used at the source and sink points.
All correlation functions are constructed using the truncated solver method with bias correction [78,93,94]. In this method, high statistics are obtained by using a low-precision (LP) stopping criterion in the inversion of the quark propagators, which was taken to be between r LP ≡ |residue| LP /|source| = 10 −3 and 5 × 10 −4 . These estimates are corrected for possible bias using high-precision (HP) measurements with r HP taken to be between 10 −7 and 10 −8 [77,78]. The number of configurations analyzed, and the number of LP and HP measurements made for the connected and disconnected contributions, are given in Table II. In our data, the bias correction term was found to be a fraction of the 1σ error in all quantities and for all six ensembles.
For the statistical analysis of the data, we first constructed bias-corrected values for the 2-and 3-point correlation functions, then averaged these over the multiple measurements made on each configuration, and finally binned these. These binned data, 250-320 depending on the ensemble, were analyzed using the single elimination jackknife process. The analysis was repeated to quantify model variation of results by choosing data with different set of source-sink separations τ and different number of points t skip , next to the source and the sink for each τ , skipped in the excited-state fits. The final result was taken to be the average over the model values, weighting each by its Akaike information criteria weight.
The bare quark mass is defined to be m bare ud = 1/2κ l − 1/2κ c , with the critical value of the hopping parameter, κ c , determined using a linear fit to (aM π ) 2 versus 1/2κ at each of the three values of a. Results for M π are given in Table II and fits to the nucleon mass are discussed in App. C. A subtle point in the renormalization of σ πN for Wilson-clover fermions is presented in App. D.
The final quoted errors are from the chiral fits shown in Fig. 3 and given in the labels. The error in each data point, σ πN = m bare ud × g u+d,bare S , combines in quadrature those in m bare ud and g u+d,bare S (see Table I), with the latter given by the appropriate fit used to remove the ESC as illustrated in Fig. 2.
In addition to the simultaneous fits to C 2pt and C 3pt to remove ESC, we have also carried out the full analysis by first calculating the M i from 4-state fits to C 2pt and using these as input in 3-state fits to C 3pt as described in Ref. [80]. The priors used for the excited state masses M i and the amplitudes A i are given in Table V in App. C. The two sets of results for M i from the two approaches (simultaneous versus individual fits) are consistent and the ground state matrix elements agree within 1σ. This agreement occurs for both strategies, {4, 3 * } and {4 N π , 3 * }. The error estimates from the simultaneous fits used to get the final results are slightly larger.

Appendix B: Chiral Perturbation Theory
The corrections to the nucleon mass and the σ-term in χPT have a long history in the literature [17,18,74,[95][96][97][98][99][100][101][102][103][104][105]. The LO, NLO, and N 2 LO diagrams contributing to the isoscalar scalar charge g u+d S are shown in Fig. 5. In these diagrams, plain and dashed lines denote pions and nucleons in the interaction picture of χPT, and not the nucleon and pion eigenstates of the full theory.
These diagrams lead to the expansion where Λ χ = 4πF π ≈ 1 GeV is the breakdown scale of the chiral expansion. The NLO diagrams (a 1 ) and (b 1 ) in Fig. 5 only contribute to the isoscalar channel, implying that the isovector channel has the different expansion We will show that the same loop diagrams responsible for the NLO and N 2 LO corrections to g S also induce a sizable contribution from N π and N ππ intermediate states. We do not consider these diagrams in our calculation, as they only produce N 2 LO recoil corrections (including (g2) despite containing only LO interactions [81,82] δσπN (τ, t), defined in Eq. (B12), is evaluated for t = τ /2 = 8a (left) and t = τ /2 = 6a (right), using the parameters of the a09m130 lattice ensemble listed in Table II. |nmax| denotes the maximum momentum included in the sum in Eqs. (B5), (B8), and (B9), mn the multiplicity of the momentum state [106]. We split the NLO corrections in Eq. (B5) into the loop diagram (a1) in Fig. 5, which also contributes to the ground state, and the diagram in which the scalar source couples to pions emitted by the nucleon source N . At N 2 LO, we label by ci the contributions in Eq. (B9) and by "recoil" those in Eq. (B8). The last line is evaluated in the continuum.
To this end, we calculate the ratio R S (τ, t) using heavy-baryon χPT and expand it as including in its definition a factor m ud to make the result R = m ud R S scale independent and ensure a normalization that facilitates the comparison to σ πN . We assume N to be a local nucleon source, transforming as 1 2 , 0 ⊕ 0, 1 2 under the chiral group SU (2) L ⊗ SU (2) R . The heavy-baryon χPT realization of N was constructed in Ref. [82] N S . The NLO diagrams receive contributions from nucleon, N π, and N ππ excited states, leading to MeV is the pion decay constant, g A = 1.276 the axial charge of the nucleon [107]. The first term in Eq. (B5) is a correction to the ground-state contribution where ∆ L σ πN is the finite-volume correction to the σ-term [108]. The remaining terms reflect the excited states, with diagram (a 1 ) receiving a contribution from N π intermediate states (nucleon and pion having opposite momenta) and from an N ππ state (the nucleon at rest and the two pions carrying momenta ±k). The N π and N ππ states with zero pion momentum vanish due to the prefactor k 2 . The amplitude of the N π contribution is suppressed by (M π L) 3 compared to g (0) S , but enhanced by the large coupling of a scalar source to the pion, making it suppressed by only a single chiral order. The last line of Eq. (B5) originates from diagram (b 1 ). In this case the dominant excited state is N ππ, with the two pions at zero momentum. Diagram (g 2 ) arises from the last term in the square bracket in Eq. (B4). Though formally NLO, this diagram vanishes up to O(1/M N ) corrections, and thus the topology (g 2 ) only contributes to N 2 LO. Similarly, the diagram with the pion emitted by N and absorbed byN vanishes at NLO.
At N 2 LO there are several contributions. They come from loop corrections to the LO (diagrams (a 2 ), (b 2 ), and (c 2 )) and from diagrams with subleading interactions in the chiral Lagrangian and the scalar source coupling to two pions (diagrams (d 2 ), (e 2 ), and (f 2 )). Here, diagrams (a 2 ) and (b 2 ) are exactly canceled by the N 2 LO corrections to the 2-point function. This is in contrast to the isovector case, in which they are responsible for the leading ESC: This result agrees with the corrections to the isovector scalar charge computed in Ref. [83], once we expand in the limit M N |k|.
Diagram (c 2 ) only contributes to the ground state, and diagrams (d 2 ) and (e 2 ) are recoil corrections. A first effect of these diagrams is to shift the energy excitation of the N π state from E π to E π + E N . The remaining recoil corrections are Finally, diagram (d 2 ) receives contributions from the LECs c 1,2,3 , which contribute to πN scattering at NLO in χPT. Including diagram (c 2 ), they give The energy gap in this case is ≈ 2M π , since k = 0 is allowed. The LECs c 1,2,3 have been determined most reliably by an analysis of πN scattering using Roy-Steiner equations, matched to χPT in the subthreshold region [40,41]. We will use the N 3 LO values which correspond to N 2 LO in the scalar form factor. We neglect the N 2 LO diagrams with pions emitted by the nucleon source, represented by (g 2 ), (h 2 ), and (i 2 ) in Fig. 5, given that for a local source the NLO contribution is already small. Of these diagrams, (g 2 ) produces N 2 LO recoil corrections, (h 2 ) contains unknown LECs that appear in the expansion of the source, and (i 2 ) cancels in the ratio between the 2-and 3-point function. While we have assumed local nucleon sources, the relative importance of the diagrams in Fig. 5 depends on the details of the lattice realization of N . However, for the Gaussian smearing applied in this work with r s ∼ 0.6 fm, corrections in addition to diagram (b 1 ) scale as (r s M π ) 2 ∼ 0.2 and can therefore be neglected, see also Refs. [81,109]. In the infinite-volume limit, the ground-state pieces of Eqs. (B5), (B8), and (B9) reproduce the N 2 LO expression for the σ-term (in the form given in Ref. [41]) except for the LO contribution proportional to c 1 M 2 π , its quark-mass renormalization proportional tol 3 − 1, and the N 2 LO LEC e 1 , all of which only contribute to the ground state. All other terms in Eq. (B11) can be obtained by replacing the finite-volume sum 1/L 3 k with infinite-volume integrals in Eqs. (B5), (B8), and (B9).
To assess the importance of the N π and N ππ contributions we define  Fig. 2 (and the shape with that of the separate contributions shown in Fig. 7). We assume gS = 18 is the asymptotic value in both cases. The contributions to δσ πN from the NLO diagrams (including the formally N 2 LO correction from E N ) and from the N 2 LO diagrams are evaluated in Table III, using the parameters of the a09m130 lattice ensemble, for two choices of source-sink separation, τ = 16a and τ = 12a. We list separately the corrections from diagram (a 1 ) and (b 1 ) in Fig. 5, as the first is dominated by N π excited states, the second by N ππ. Similarly, the N 2 LO corrections from the diagrams proportional to the LECs c 1,2,3 receive contributions from N ππ states, while the recoil corrections in Eq. (B8) are dominated by N π. From Table III we see that R (1) and R (2) are of similar size and the most important contributions come from diagram (a 1 ) and from Eq. (B9). The nucleon source and recoil contributions are substantially smaller.
We also note that the ESC is amplified by the presence of several states close in energy. We can see this in the left panel of Fig. 6, where the red and orange curves denote the NLO and N 2 LO corrections, including states up to |n max | = 1, while the green and blue curves include states up to |n max | = 3. In the continuum, the effect of the entire tower of excited states can be resummed, with the result shown in the last line of Table III and by the purple line in Fig. 6. The comparison to the different |n max | values indicates the degree of convergence, which is ultimately determined by t and τ via the exponential suppression of the continuum integrals. Indeed, we see that for τ = 16a the corrections at t = 8a beyond |n max | = 3 are around 10%, but twice as large for τ = 12a, t = 6a.
Away from t ∼ τ /2, the integrals become increasingly dominated by large-momentum modes, leading to the eventual breakdown of the chiral expansion. For this reason, the expansion becomes less reliable for small and large t, which explains why the functional form predicted by χPT, see the right panel of Fig. 6, suggests a faster decrease towards the edges than observed in the lattice data, see Figs. 2 and 7. (Note that in fits to remove the excited-state effects in the lattice data, we neglect t skip points next to the source and the sink for each τ .) In the     center, however, the EFT expansion should be reliable, in particular for the τ = 16a variant, which leads to the conclusion that NLO and N 2 LO contributions can each lead to a reduction of σ πN at the level of 10 MeV. Note that Fig. 7 shows the behavior versus t and τ for the individual contributions, i.e., insertions on the u, the d, and the disconnected loop l. Table III and Fig. 6 we focused on the ensemble with the lightest pion mass, we find good agreement between the χPT expectations and the fits in Fig. 2 also for the remaining ensembles. For example, with the parameters of the a09m220 ensemble, χPT predicts the difference between the lattice data and the extrapolated value of g S to be 3 at t = τ /2 = 8a, in good agreement with the left panel in the second row of Fig. 2. Similarly, for the a06m220 and a06m310 ensembles we obtain that g S is shifted by 2.8 and 1.6, at t = τ /2 = 12a.  Fig. 3. Four of these additional fits are shown in Fig. 8, with resulting fit parameters given in Table IV. (We neglect possible discretization errors in these fits as they are not resolved in our best fits shown in Fig. 3.) We see that for the {4 N π , 3 * } strategy all fit variants lead to parameters that agree with the χPT prediction (in-cluding evidence for the nonanalytic M 3 π term from the {2, 3} and {2, 3, 4} fits), with changes in σ πN consistent with the uncertainty assigned to the least constrained fits ({2, 3, 4} and {2, 3 χ , 4, 4L}), whose average we quoted as our central result. Even in the most-constrained case {2 χ , 3 χ , 4, 4L χ } a good fit of the data is obtained, in marked contrast to the {4, 3 * } strategy, in which case the χ 2 /dof becomes unacceptable when imposing the maximum amount of chiral constraints.

While in
The expressions in Eqs. (B5), (B6), and (B9) can also be used to assess the importance of finite-volume corrections to g S . Focusing on the ground-state contributions, we can write ∆ L σ πN in Eq. (B6) as [108] (B13) Similarly, for the c i contributions (B9) we obtain in terms of the Bessel functions K 0 and K 1 . Using the parameters of the a09m130 lattice ensemble, we get  9. Chiral fits to the {4 N π , 3 * } data for the nucleon mass M0 given in Table I. Data at a seventh point, a09m310, is included from Ref. [78]. The fit in the middle panel is overparameterized. The right panel shows that the {2 χ , 3 χ , 4, 4L χ } fit, with χPT values for three of the coefficients taken from Ref. [41], is reasonable and comparable to the fit in the left panel.
implying that the finite-volume corrections are controlled at the level of about 1 MeV. We end this appendix by pointing out a subtlety. We have used the lowest-order (linear) relation between M 2 π and 1/2κ to get m bare ud as we have ensembles at only two values of M π at each a. Higher-order corrections give the LECl 3 term in Eq. (B11). Removing it changes the χPT prediction of the d 4L term in Table IV  In both types of analyses, simultaneous fits to C 2pt and C 3pt and individual fits to them, we have used the same empirical Bayesian priors for the amplitudes and masses of the three excited states determined using the procedure described in Ref. [78]. The central values of these priors and their widths for {4, 3 * } and {4 N π , 3 * } fits are given in Table V. The resulting M 0 , M 1 , and M 2 that enter in {3 * } fits to C 3pt are given in Table I, and are essentially the same from the two types of analyses. The result for the nucleon mass on the seventh ensemble, a09m310, used only for the analysis of M N in this section, is 1.09(1) GeV from both analyses.
In this section, we study the chiral behavior of the nucleon mass M N ≡ M 0 using the N 2 LO χPT result, which has a form similar to Eq. (4): with the χPT expressions for the c i given, e.g., in Ref. [41]. Even ignoring discretization and finite-volume corrections and fitting the lattice data using Eq. (C1) we face two challenges. First, as evident from the data for M 0 in Table I, there is no significant difference between the {4, 3 * } and {4 N π , 3 * } results. Therefore, we cannot comment on the impact of N π states in the analysis of M N . Second, as shown in Fig. 9, at least three parameters (c 0 and two more) are needed to fit the data. With data at essentially only three values of M π , even these fits are overparameterized. In short, data at many more values of the lattice parameters, especially in M 2 π , are needed to quantify the lattice artifacts and check the prediction of χPT. The most constrained fit, {2 χ , 3 χ , 4, 4L χ } in Fig. 9, does yield a good description of the data, but the resulting uncertainty in σ πN via the FH method is too large to compete with the direct method.

Appendix D: Renormalization
In this appendix we discuss the renormalization of the quark mass and scalar charge for fermion schemes that break chiral symmetry such as Wilson-clover fermions. We will do this using the notation and results given in Ref. [110] for an N f flavor theory with two light degenerate flavors, m u = m d ≡ m l , and N f − 2 heavier flavors denoted generically by m s . Here m i andm i denote the bare and renormalized quark masses for flavor i and we will neglect all O(a) terms as they do not effect the continuum limit. These discretization errors start at O(a) in our calculation.
We start with Eq. (26) in Ref. [110]: where m j are defined as (1/2κ i ) − (1/2κ 0 c ) with κ i the Wilson hopping parameter and κ 0 c its critical value defined to be the point at which all pseudoscalar masses vanish, i.e., the SU (N f ) symmetric limit. Using Z m to denote the flavor nonsinglet and Z m r m the isosinglet renormalization constants, one has the relation [110]  (D5) Similarly, for the scalar density S i with flavor i, and its expectation value S i in any fixed state we have from Eqs. (22)(23) in [110] with Z S and Z S r S the nonsinglet and the singlet renormalization constants.
Using Z S Z m = 1 and r S r m = 1 [110], we get the desired result showing that mixing between flavors in Wilson-like formulations gives rise to a correction to σ πN , i.e., the second term on the second line. To explain this, we need to clarify an important point about the notation. In this work we have defined κ c as the point at which M 2 π vanishes with all the m s kept at their physical values. The difference in the definition of the chiral point, κ 0 c versus κ c , leads to a difference in the definition of the bare quark masses. The connection is that the bare quark mass used in this work is the same as m l − m l (m l = 0)| fixed sm s in Eq. (D7).
We have not calculated s S s that is needed to evaluate the correction to the isoscalar scalar charge, however, it is relatively small since |1 − r m | 2% for the six ensembles and s S s < S l . At the precision at which we are working in this paper, and since the focus is on showing that there is a difference in the result depending on the excited state analyses, i.e., between {4, 3 * } and {4 N π , 3 * }, this correction is neglected. Also, note that r m = 1 for lattice QCD formulations that preserve chiral symmetry, in which case there is no correction.