Nucleon Sigma Terms with $N_f = 2 + 1$ O($a$)-improved Wilson fermions

We present a lattice-QCD based analysis of the nucleon sigma terms using gauge ensembles with $N_f = 2 + 1$ flavors of ${\cal O}(a)$-improved Wilson fermions, with a complete error budget concerning excited-state contaminations, the chiral interpolation as well as finite-size and lattice spacing effects. We compute the sigma terms determined directly from the matrix elements of the scalar currents. The chiral interpolation is based on SU(3) baryon chiral perturbation theory using the extended on-mass shell renormalization scheme. For the pion nucleon sigma term, we obtain $\sigma_{\pi N} = (43.7\pm3.6)$ MeV, where the error includes our estimate of the aforementioned systematics. The tension with extractions based on dispersion theory persists at the 2.4-$\sigma$ level. For the strange sigma term, we obtain a non-zero value, $\sigma_s=(28.6\pm9.3)$ MeV.

Introduction.The scalar matrix element of the nucleon is an important observable, and plays a crucial role in interpreting the results of dark-matter directdetection experiments.Especially appealing candidates for cold dark matter are weakly interacting massive particles (WIMP), as they naturally reproduce the observed relic abundance of dark matter through annihilation processes in the early universe.In particular for Higgs-portal models, in which the WIMP-nucleus interaction is mediated by the Higgs boson, the spin-independent crosssection for WIMP-nucleus recoil experiments is sensitive to the values of the scalar matrix element [1].The lightquark scalar matrix element 1   σ πN ≡ m l ⟨N |ūu + dd|N ⟩ = m l (∂m N /∂m l ), where m l ≡ (m u + m d )/2, also known as the pionnucleon sigma term, is of special interest.Phenomenologically, σ πN is accessible via πN -scattering amplitudes at the Cheng-Dashen point [2].Historically, the value for σ πN ∼ 45 MeV derived in [3] was prevalent for a long time, a value compatible with most lattice determinations.However, new analyses using constraints from pionic hydrogen and deuterium led to a much larger value of σ πN = 59.1(3.5)MeV [4], consistent with the EFT analysis of [5] and in agreement with [6] based on low energy πN -scattering (see Ref. [7] for a review).By contrast, lattice calculations for σ πN [8][9][10][11][12][13][14][15][16][17][18][19], discussed in detail in the FLAG report [20], have largely confirmed the lower estimate, while being in tension with the latest dispersive analysis at the level of 3−4 standard deviations. 2 Very recently, it was suggested that the discrepancy is alleviated via an explicit treatment of N π and N ππ excited 1 We take the nucleon at rest and use the state normalization ⟨N ⃗ p ′ s ′ |N ⃗ p s⟩ = (2π) 3 δ ss ′ δ(⃗ p − ⃗ p ′ ).Also, throughout this work we assume exact isospin symmetry. 2See Refs.[21][22][23][24] for further efforts to extract σ πN from collections of lattice data for the light quark mass dependence of m N .
states in the analysis [25].As a related quantity, the strangeness matrix element a pure sea-quark effect, has often been discussed together with the pion-nucleon sigma term.Their linear combination is to first order in (m l − m s ) proportional to the nucleonhyperon mass splitting.The σ 0 value inferred from this observation, assuming a negligible strangeness content σ s of the nucleon, corresponds to a small value for σ πN .In [26] however, corrections to σ 0 were calculated which bring the associated σ πN estimate back into agreement with its Cheng-Dashen-theorem based determinations without the need to invoke a large σ s value.We perform a direct determination of the nucleon sigma terms from a lattice calculation of the matrix element of the scalar current.Our final estimates are based on a simultaneous chiral, continuum and infinite volume extrapolation of the pion-nucleon and strange sigma terms.We average the individual fits with weights based on the Akaike information criterion (AIC) [27,28] to provide a full error budget accounting for variations in the treatment of excited state contaminations, discretization errors, finite-volume effects and the quark-mass dependence.
Simulation details.We employ the N f = 2 + 1 ensembles [29] generated as part of the Coordinated Lattice Simulations (CLS) initiative with non-perturbatively O(a)-improved Wilson fermions [30] and the tree-level improved Lüscher-Weisz gauge action [31], correcting for the treatment of the strange quark determinant using [32].Table I gives details of the ensembles used in this work.In particular, lattice spacings range from 0.050 fm to 0.086 fm.
The two-point and three-point functions needed to extract the scalar matrix elements of the nucleon read x,y e iqy ⟨Ψ β (x, t s )S q (y, t)Ψ α (0)⟩, (5) where S q denotes the scalar density, The interpolating operator for the proton, is built using Gaussian-smeared quark fields [33] q and spatially APE-smeared gauge links in the covariant Laplacian ∆ [34].
The pertinent Wick contractions for the three-point function lead to the connected and disconnected contributions, . For the connected part, we employ extended propagators via the "fixed-sink" method, requiring additional inversions for each chosen value of t s [35].In order to reduce the cost of the inversions, we apply the truncated solver method with bias correction [36][37][38].For the connected part, the polarization matrices Γ ′ , Γ read The disconnected three-point function is constructed from the quark loop L q and the nucleon two-point function where Note that for forward scalar matrix elements (q = 0), the vacuum expectation value of the current insertion must be subtracted, Additionally, we improve the signal by averaging over all three different polarizations and by averaging over forward and backward propagating nucleons.[44] and the lattice spacings from [45].Using Eq. ( 16) the largest source-sink separations correspond to 1.4 fm and 1.5 fm for the two finer and coarser lattices, respectively.
stochastically using four-dimensional noise vectors η.We improve the precision of the quark loops using a variation of the frequency splitting method [39] that combines the one-end-trick [40] with a generalized hopping parameter expansion [41] and hierarchical probing [42] (for more details see App.C Ref. [43]).Let G S ≡ ⟨N |S q |N ⟩ denote the nucleon scalar form factor at vanishing momentum transfer.It can be extracted from the ratio of correlation functions Indeed, let ∆ be the the energy gap between the lowest excited state and the ground state.Performing the spectral decomposition in Eq. ( 14) and taking the limit of t, (t s − t) ≫ ∆ −1 , we obtain We extract the ground-state contribution for each flavor combination of the scalar current corresponding to σ πN , σ s and σ 0 .Errors are computed using the bootstrap method on binned data with a bin size of two.For the conversion to physical units, we first express dimensionful quantities in units of t 0 using Ref. [45] (see Tab. I) and finally use the value from [20] √ t 0 = 0.14464(87) fm (16) to calibrate the scale.FIG. 1. Left: Results of linear fits to the summed correlator on ensemble D200 with the starting time slice given on the x-axis.
The blue shaded area is the weighted average using Eq. ( 18) shown as a black line in the bottom of the plot, for the particular choice of parameters from Eq. (19).Right: Fit result of an explicit two-state fit to the effective form factor.The gray band represents the result for the ground-state matrix element of that fit; it is shown together with the result of the window average (black filled square) and the result of a two-state fit to the summed correlator (black diamond).
Excited-state analysis.A major obstacle to achieving reliable and precise determinations of the ground-state matrix element is the well-known noise problem of nucleon correlation functions [46,47].For typical source-sink separations in current lattice calculations, the ratio in Eq. ( 14) will be contaminated by exponentially suppressed terms associated with resonances and multi-hadron states.Several approaches were developed to have a better control over the excited-state systematics (see [48,49] and references therein).The summation method [50][51][52] and multi-state fits are the most widely used among them.
In the summation method, the ground-state matrix element is determined from the summed ratio by fitting b 1 and G S to S(t s ).We have extended the number of source-sink separations compared to our analysis of the isovector vector form factor [53] to include smaller source-sink separations.This enables us to monitor the range of t s where the result from the linear ansatz of Eq. ( 17) stabilizes.Rather than selecting a single fit starting at a certain value t min s , we follow the procedure defined in [54] and determine G S from an average over a range of t min s values with weights with N a normalization factor.The choice of lower (t lo ) and upper (t up ) bound suppresses the excessive influence of excited states at small values of t s and the exponentially increasing noise at larger values, respectively.We find the choices t lo = 0.8 fm, t up = 1.0 fm and ∆t = 0.08 fm, (19) to give estimates for the ground state matrix element that are more robust against statistical fluctuations than choosing one particular value of t min s .Since the onset of a plateau in the extracted matrix element as a function of t min s , such as in the left panel of Fig. 1, does not entirely exclude the possibility of remnant excited-state contributions, we also apply two further analysis methods.We note that the average using Eq. ( 18) is only applied to the fit results of Eq. ( 17) for different t min s .As a cross check, we performed fits to the summed correlator including the first excited state contribution, where m10 and m11 involve matrix elements of S q from first excited to ground state and excited to excited state, respectively.The excited-state contributions are parametrically suppressed by ∆ ⋅ t s .In this case, we need priors for the energy gap ∆ in order to stabilize the fits.We choose twice the pion mass on the given ensemble as the central value and assign a total prior width of five percent.Even with a prior for the energy gap ∆, m11 is not well constrained, and we resort to a simplified fit ansatz excluding this term.
In addition to the analysis of the summed correlators, we performed fits using a two-state ansatz for the effective form factor itself.The fit function reads Similar to the analysis of the summed correlators, the gap of the first excited state is not well constrained and we are led to using priors.For the priors, we use the same setup as in the two-state fit to the summed correlator.Even though the neglected excited-state contributions in the two-state ansatz are parametrically less suppressed, we include all t s > 0.8 fm, i.e. the same value as for t lo in Eq. ( 19).Subsequently, we cut time slices at the source and sink until a good fit is achieved.For σ s the data is too noisy to perform two-state fits of the effective form factor, and we resort to plateau fits, where we fit different t s and use the value that shows convergence with t s .A two-state fit applied to data at m π = 200 MeV is illustrated in Fig. 1 (right panel), along with the results of the two other methods.Fig. 2 shows a comparison of the σ terms obtained from the different excited-state analyses, and the results are collected in Tab.V in the appendix.
While the summation method with the averaging window fixed in units of fm is adequate if the dominant excited-state contribution is only weakly dependent on the pion mass, the two other analysis methods explicitly assume the bulk of that contribution to be associated with a mass gap ∆ = O(m π ).Therefore, in terms of excited-state contamination, we essentially have two procedures, either relying on the applicability of Eq. ( 17) or, relying on assumptions about the energy gaps through priors, applying Eqs. ( 20) and (21), where the latter are both very sensitive to the prior, but give consistent results.In order to assess the systematics associated with the very different effects of excited states in the two strategies, we perform the chiral and continuum extrapolation for the window averaged summation method (fit ansatz Eq. ( 17)) and for one method using priors (fit ansatz Eq. ( 21)), and finally model average the results with equal weights, i.e. giving no preference to either strategy.
Chiral and continuum extrapolation.The calculation of the σ term in chiral perturbation theory (ChPT) proceeds via the nucleon mass using the Feynman-Hellmann theorem.The nucleon mass has been calculated in various formulations of ChPT [26,[55][56][57] up to two-loop order [58].
Since our gauge ensembles lie on a line of constant trace of the quark mass matrix (2m l + m s ), both the pion and the kaon mass change as m l is varied.Moreover, to have a handle on the quantities σ 0 and σ s , the inclusion of the strange quark into the effective theory is mandatory.We therefore use the result of SU(3) ChPT in the extended on-mass shell scheme (EOMS) of [59].The nucleon mass reads with For the η meson mass, we assume the Gell-Mann-Okubo relation 3M We fix the values of the low-energy constants (LECs) D = 0.8, F = 0.46, F ϕ = 0.108 GeV, m 0 = 938.9MeV and fit the constants b0 and b1 .For the physical point we use the isospin-limit meson masses M π = 134.8MeV and M K = 494.2MeV [60].Eq. ( 22) is derived with respect to the quark masses, yielding the quark-mass dependence of the sigma terms.For the quark mass dependence of the octet meson masses we take the leading order expression in ChPT [61].
We treat the lattice spacing dependence of the sigma term via an additional term 3 , The finite-volume dependence of the nucleon mass in SU( 2) is given in [62], wherefrom we derive We only use the finite-volume corrections due to pion loops, as terms ∼ exp(−M K L) are parametrically much more suppressed; thus we omit finite-volume corrections for σ s .Instead of using the ChPT results for the prefactors of the finite-volume corrections, we leave them as additional fit parameters, however we use as a loose prior the value obtained from SU(2) ChPT.We proceed to fit σ πN , σ s , taking into account the correlations among the sigma terms and lattice spacing.The fits are performed with variations in the upper end of the pion mass range (220, 285 or 360 MeV), and including/excluding the artifacts with respect to finite lattice spacing and to finite volume.We analyze the two data sets obtained from the excited-states analyses separately with respect to the above variations, i.e. within each data set all variations are averaged using an AIC weight w i given by where n c and n f denote the number of cut data points and number of fit parameters, respectively.The weights are normalized per data set, and finally a flat weighting is applied between the data sets.Using the procedure of [53,63] we obtain as our final estimates where the first and second errors correspond to the statistical and systematic uncertainties, respectively.More details of the averaging procedure are given in the appendix.The systematic error dominates, with the largest source of uncertainty coming from the treatment of excited states.In Fig. 3 we compare our results to those of other lattice calculations.We note a reasonable agreement among these calculations.

Conclusion.
We have calculated the nucleon sigma terms σ πN , σ 0 and σ s with a full error budget concerning excited-state contamination as well as chiral, finite-size and continuum extrapolations.Our estimate for σ πN lies close to the early estimate from N π scattering [3].It is compatible with most other lattice determinations and in excellent agreement with the σ πN determination of [19], which uses partly the same gauge ensembles but proceeds by computing the quark-mass dependence of the nucleon mass.For σ s we find a non-zero value, again compatible with most recent lattice determinations.Including the effects of different methods for the treatment of excited states into our error budget, we clearly establish this to be the largest source of systematic uncertainty.
Analyzing the data sets from the window and two-state procedure separately, see Tab.II, we observe an upwards trend for σ πN when using priors similar to [25], albeit not as pronounced.Our final central value for σ πN lies between the two values presented in [25], but is much closer to that obtained without imposing tight priors on the gap ∆ around values of order m π .A discrepancy of 2.4 σ persists with the dispersive result of [4], after applying the correction necessary to match our definition of the pion mass in the isospin-limit from Ref. FIG. 3. Comparison of our results to other lattice determinations: RQCD22 [19], NME21 [25], BMW20 [18], ETM19 [17], JLQCD18 [16], RQCD16 [15], χQCD15 [14], BMW15 [13], ETM14A [12], QCDSF12 [9], BMW11A [8].Filled circles represent results extracted from the slope of the nucleon mass with respect to the light quark mass m l , and squares represent results obtained directly from the matrix element.The gray band corresponds to the dispersive result of [4] with the correction for the isospin-limit value of the pion mass from [64] applied, i.e. σ πN = 55.9(3.5)MeV.
(the latter through the John von Neumann Institute for Computing (NIC)), as well as on the GCS Supercomputer HAZELHEN at Höchstleistungsrechenzentrum Stuttgart (www.hlrs.de)under project GCS-HQCD.
Our programs use the QDP++ library [66] and deflated SAP+GCR solver from the openQCD package [67], while the contractions have been explicitly checked using [68].We are grateful to our colleagues in the CLS initiative for sharing the gauge field configurations on which this work is based.
In the continuum, the operator m q qq is invariant under renormalization group transformations.However, Wilson fermions explicitly break chiral symmetry, and this enables mixing with other quark flavors.In the presence of chiral symmetry breaking by the regulator, flavor-nonsinglet and flavor-singlet operators represent the more adequate basis of operators to work in.Indeed, it is straightforward to show that the operators Σ (8)  = (m q,l − m q,s )(ūu + dd − 2ss), (32) are renormalized even in Wilson-action lattice QCD, where it is the bare quark masses and bare scalar operators that appear on the right-hand side.However, in order to realize O(a) improvement, O(a) counterterms are required.In particular, terms of type O(am q ) and the gluonic operator aF 2 can be included [69].
The bare quark masses are related to the hopping parameters via where κ crit is the hopping parameter at which the octet of pseudoscalar mesons becomes massless in the SU(3) symmetric theory.The values of κ crit for our action can be found in [70].Using the previously defined operators Σ (0) , Σ (8) , the light-and strange-quark operators can be reconstructed as Σ (8) (34) Σ (8) .(35) An important observation is that the ratios of quark masses appearing in these expressions can be evaluated using the PCAC (partially conserved axial current) quark masses (see e.g.[70] for their definition, including O(a) improvement).Proceeding in this way by-passes the use of the finite renormalization factor r m , which parametrizes the difference in renormalization of the SU(3) f octet and singlet quark mass combinations [69].We note that for QCD actions with an exact chiral symmetry r m = 1 [71].
O(a) improvement of Σ (0) and Σ (8)   In the following we estimate the size of the O(a) corrections for the flavor-singlet and non-singlet scalar operator in N f = 2 + 1 flavor QCD.First, we recall that the renormalization and improvement pattern of a nonsinglet combination of quark masses reads [69] ml − ms = Z m (m q,l −m q,s )[1+3b m am av q +ab m (m q,l +m q,s )], (36) while the singlet combination renormalizes as = Z m r m [(1 + ad m 3m av q )3m av q + ad m (2m 2 q,l + m 2 q,s )].
We follow the notation of [69], denoting by a hat an operator or a parameter of the theory that has been renormalized and O(a) improved.

The octet scalar operators
Let S f ′ ,f ≡ ψf ψ f ′ and λ a denote a Gell-Mann matrix.We then define S a ≡ Tr{λ a S} = ψλ a ψ to be the octet of scalar currents, and Tr S = ψψ to be the flavor-singlet current.The octet of scalar currents has no additive improvement term in the massless limit.Thus the renormalization and improvement pattern of the local discretiza-tion of the two neutral octet combinations reads [69]

Ŝ3
= Z S (1 + 3b S am av q + b S am q,l ) S 3 , (38) where Z S = 1/Z m , and the improvement coefficients are not independent [69], Taking into account these relations, we obtain The difference d m − b m is a sea-quark effect that we will neglect in the following.Furthermore, b S = −2b m has been determined on CLS ensembles in [72].We note that in perturbation theory, b m = − 1 2 + O(g 2 0 ).At β = 3.55 for instance, Ref. [72] finds b m = −0.835.Since on ensemble D200 a(m q,l − m q,s ) = 1 2κ l − 1 2κs ≃ −0.0160, we arrive at the estimates a(m q,l − m q,s ) 1  3 a(m q,l − m q,s ) The O(a) corrections are thus on the order of a few percent.They reduce slightly the weight of the light quarks and increase the magnitude of the weight of the strange quark.
The flavor-singlet scalar operator The improvement of the scalar operator Tr S ≡ ūu + dd + ss (44) in the SU(3) chiral limit is given by where the gluonic action is given by The renormalized, improved scalar operator then reads FIG. 4. Left: Results of linear fits to the summed correlator on ensemble E300 with the starting time slice given on the x-axis.
The blue shaded area is the weighted average using the weight function of Eq. ( 18) shown as a black line in the bottom of the plot, for the particular choice of parameters from Eq. ( 19).Right: Fit result of an explicit two-state fit to the effective form factor.The gray band denotes the result for the ground-state matrix element of that fit, together with the result of the window average (black filled square) and the result of a two-state fit to the summed correlator (black diamond).
Note that M S = diag(m q,l ūu, m q,l dd, m q,s ss).Again, using Z S = 1/Z m and r S = 1/r m , as well as the relations among the improvement coefficients [69] one obtains Σ(0) ≡ Tr(M ) Tr(S) = (ūu + dd)[3m av q + b m 2am av q (m q,s − m q,l ) +a(m 2 q,s − 8m q,l m q,s − 2m 2 q,l ) +b m a(m 2 q,s + m q,s m q,l − 2m 2 q,l ) −d m a(2m q,l + m q,s ) 2 ] +ss[3m av q + 4a(m q,l − m q,s )m av q b m +a(m 2 q,s − 8m q,l m q,s − 2m 2 q,l ) In order to estimate the size of the O(a) correction we take d m ≃ b m from [72] and b m ≃ 0 ≃ d m .If m q,l ≪ m q,s , then the O(a) corrections in the square brackets have the same relative size as in the octet case, Eq. (42)(43).
As for g S , we note the relation [69] g where we used the one-loop result of [73] for b g .The trace anomaly in the nucleon at rest yields (see for instance [74]; we use the non-relativistic normalization of the nucleon state at rest) The trace anomaly is related to the improvement term via with Hence the last term in Eq. ( 50) is of the order of −2 MeV for a = (3GeV) −1 .For σ πN there is a suppression by 2 3 m l 2m l +ms ≃ 0.025.Thus we expect the lattice artifacts due to the gluonic operator to be of order −0.05MeV in σ πN .This is certainly negligible, even if the perturbative estimate of g S was too small by an order of magnitude.

Correlator Analysis Details
In our calculation of the sigma terms the statistical precision of the correlator is restricted by two factors, the signal-to-noise of the disconnected contribution and the occurrence of exceptional configurations which produces outliers.The signal-to-noise problem of the disconnected part is exacerbated at small pion mass, since the absolute contribution of the connected part decreases.In the following we summarize our strategy to deal with these two problems.

Disconnected Part
We observe that the connected part for a given statistics is far more precise than the disconnected.The dis- connected part consists of the loop contribution and the two-point function.For the quark loops, we have exhausted the number of sources per configuration, for which error scaling still holds, on most of our ensembles.Consequently we concentrate on the 2-point function, i.e. improving the signal of the disconnected part by using additional sources for the 2-point functions.For the connected part we keep a matched setup between 2and 3-point functions, as the correlation plays an important role.In Fig. 4 we show a comparison of the excited state analysis between summation and explicit two-state ansätze for E300.Here the signal for the effective form factor has been improved via additional two-point functions.However, even after including additional sources statistical fluctuations are still clearly visible.The different strategies are explained in the main text, for the window average the weight function reads with the particular choice t lo = 0.8 fm, t up = 1.0 fm and ∆t = 0.08 fm. ( 56)

Outliers
In contrast to other observables in our previous analysis, we observe a small number of configurations for which the effective form factor of the scalar operator exceeds the ensemble average by a huge amount, which potentially spoils the correct estimation of the error.Removing the measurements on a (negligibly small) fraction of the configurations considerably improves the error estimate.We may assume the sampling distribution, either Jackknife or bootstrap, in the limit of a large number of measurements to be normally distributed.However on some of the ensembles the observed distribution deviates strongly from a normal distribution.The most prominent example is ensemble D200, where we identify one configuration to be the root cause of the gross overestimation of errors, i.e. an outlier in our analysis.It is well known that the mean is not a robust estimator with respect to outliers and may be highly affected by the existence of extreme values on the correlator level.We try to identify the extreme values on a per-configuration basis performing first a Jackknife analysis, where the sampling distribution should be approximately Gaussian.We essentially look for extreme deviations from the central location of the sampling distribution on each time slice and for every source-sink separation for the effective form factor corresponding to the connected and disconnected part of the sigma terms, respectively.Whenever we find a value that is more than ∼ 6σ away, we flag the Jackknife sample, i.e. configuration, and remove it from the subsequent analysis, see Fig. 5.For estimating the central location of the distribution and its standard deviation we use the median and median absolute deviation as robust replacements for the mean and standard deviation.For the latter we apply the usual correction factor Φ −1 (3/4) = 0.67449 to make contact with the standard deviation of a normal distribution.We find this procedure correctly identifies all problematic results, that either are not symmetric with respect to t s /2 and/or are very far away from the center of observations.Data from the latter category may also come from a sampling distribution that has a longer tail than the normal distribution.We therefore apply a very loose cut using ∼ 6σ, i.e. the number of flagged configurations is kept to a minimum, so as to not distort the empirical distribution.The number of flagged configurations is generally well below 1%, except for E250 with 1.5 %, which also has the smallest number of configurations amongst the ensembles analyzed.

Results for the sigma terms
The results for the three determinations described in the main text for all ensembles are collected in Tab.V.
The conversion to physical units uses the ratios t0 a 2 from [45] and √ t 0 = 0.14464(87) fm (57) at the physical point from [20].The error estimate is based on Bootstrap procedure with a sample size of 5000.

Fits and model average
In Fig. 6 we show one particular fit for the summation window averaged data based on the SU(3) formula for the nucleon mass Eq. ( 22) without any cut in the pion mass including finite size effects.The data have been corrected for finite volume effects only, while the fit is at physical kaon mass.
We derive the expression for the sigma terms from the nucleon mass with Using the lowest order ChPT expression for the quark mass dependence of the meson masses, the sigma terms read The fits are performed simultaneously to σ πN , σ s and m N , where we include the correlations among the sigma terms.We perform variations of these fits, i.e. three cuts in the pion mass, including/excluding lattice spacing, including/excluding finite volume and including both lattice spacing and finite volume corrections.The strictest pion mass cut is such that enough data points remain to perform the fit using all values of the lattice spacing.We treat the data subset selection problem using the "perfect model" method of Ref. [75].In total we thus have 12 variations on three data sets.Instead of choosing a particular fit we perform model averages over the 36 fits using their AIC weights.As described in the main text,  FIG. 6. Simultaneous fit of the window averaged data for σ πN , σs.In this variation we include corrections due to finite volume, using all available pion masses.We correct the central value for the fitted finite volume correction only.
only two data sets enter the final analysis.The weights are normalized first on each data set, and subsequently averaged using flat weights, i.e. with a factor 1/2. From these weights we build a cumulative distribution function (see Fig. 7) following Ref.[18] P We estimate the central value and the total error of the average, using the median and the difference between the 1-σ percentiles of P x .For the separation into statistical and systematic errors we assume that and that a scaling of the individual Bootstrap errors with an arbitrary constant is expected to affect σ 2 stat exclu-sively.In Fig. 7 we show the CDF for all three quantities, note that σ 0 is not fitted.The blue shaded area is the symmetric error from the percentiles centered around the median.These coincide rather well with the 1-σ percentiles of the actual distribution.The breakup into systematic and statistical error uses the fact that the systematic error does not change if all errors are inflated by some arbitrary factor (see Ref. [63] for more details).
In Tab.VI we collect the results for the different variations and the corresponding weights in the averaging procedure w i .We note that we performed the averaging procedure with all quantities expressed in units of t 0 , applying the calibration in the end to convert to physical units.We see that in general the fit quality for all variations is acceptable.The penalty terms in the AIC weights prefer variations with more data and fewer fit parameters.That is visible for the window data, where most of the  64) for all variations of the fits to the sigma-terms.The red squares and purple dots denote fits based on summation window and two-state data, respectively.
weight is on the fits using all available pion masses and including a finite volume correction.On the other hand, the two-state data prefers fits with stricter cuts in the pion mass, and again finite volume corrections.When performing the analysis separately for the summationwindow and two-state data for the sigma terms we obtain where only the total error is given.All values are within 1 − σ of our best estimate, as can be seen in Fig. 7, where the bulk of points is covered by the total errors of our best estimate.We note that the AIC averaged result is stable with respect to including models where only terms of second order in the pion-and kaon-mass are used, and a model adding polynomial fourth order terms in the chiral counting.The former turns out to have less AIC weight compared to our (third-order) estimate, while the latter needs to be stabilized using priors.In both cases, the changes in the central values are insignificant compared to our best estimate, and the error changes within a few percent, depending on the prior applied for the fourth-order term.Similarily, removing all data points with a pion mass above 285 MeV from the analysis only has very small effect on the central value and error.Moreover, we checked that the AIC average is also stable against variations in the low-energy constants F ϕ , D and F .To this end we have varied the values of the LECs given in Tab.11 of Ref. [19] within one standard deviation, and added these as additional models in the averaging.

FIG. 2 .
FIG.2.Comparison of the different extractions.The blue squares, green circles and red diamonds correspond to the extraction based on the window average of the summed correlator, the explicit two-state fit to the summed correlator, and the explicit two-state fit to the effective form factor.

FIG. 5 .
FIG.5.Jackknife distribution of the connected part for the σ-term for all source-sink-separations. Dashed are the flagged configurations by the analysis described in the main text.For D200 a total of 4 configurations are flagged out of 2000.

FIG. 7 .
FIG.7.Cumulative distribution function of Eq. (64) for all variations of the fits to the sigma-terms.The red squares and purple dots denote fits based on summation window and two-state data, respectively.

TABLE I .
Details of CLS ensembles used in this work.The pion and kaon masses are taken from

TABLE II .
Result of the model average procedure using AIC weights defined in Eq. (29) when applied exclusively to the data set denoted in the column heading.Only total errors are shown.

TABLE III .
Number of exact and sloppy sources for the calculation of the two-point function, which enters the estimate of the disconnected contribution.Ensembles in bold, number of sources is increased with respect to the statistics for the connected three-point function.

TABLE IV .
Number of flagged configurations for each ensemble.The first number refers to the connected part, while the second concerns the disconnected.

TABLE V .
Results for the sigma terms on every ensemble in MeV, where window, sum two-state and two-state, refer to the window average of the summed correlator, the two-state fit to the summed correlator and the direct two-state fit to the correlator, respectively.