Ab-initio calculation of the proton and the neutron's scalar couplings for new physics searches
Sz. Borsanyi, Z. Fodor, C. Hoelbling, L. Lellouch, K.K. Szabo, C. Torrero, L. Varnhorst
AAb-initio calculation of the proton and the neutron’sscalar couplings for new physics searches
Sz. Borsanyi , Z. Fodor , , , C. Hoelbling , L. Lellouch , ∗ ,K. K. Szabo , , C. Torrero , † , L. Varnhorst , Department of Physics, University of Wuppertal, D-42119 Wuppertal, Germany, Institute for Theoretical Physics, Eötvös University, H-1117 Budapest, Hungary Jülich Supercomputing Centre, Forschungszentrum Jülich, D-52428 Jülich, Germany Aix Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France ∗ To whom correspondence should be addressed; e-mail: [email protected] † Present address: U-Hopper Srl, 38122 Trento, Italy
Many low-energy, particle-physics experiments seek to reveal new fundamen-tal physics by searching for very rare scattering events on atomic nuclei. Theinterpretation of their results requires quantifying the non-linear effects ofthe strong interaction on the spin-independent couplings of this new physicsto protons and neutrons. Here we present a fully-controlled, ab-initio cal-culation of these couplings to the quarks within those constituents of nuclei.We use lattice quantum chromodynamics computations for the four lightestspecies of quarks and heavy-quark expansions for the remaining two. We de-termine each of the six quark contributions with an accuracy better than 15%.Our results are especially important for guiding and interpreting experimentalsearches for our universe’s dark matter. a r X i v : . [ h e p - l a t ] J u l any low-energy experiments, that search for new elementary particles or interactions,are based on detecting very rare scattering events on atomic nuclei. These include searchesfor weakly interacting massive particles (WIMPs), that could constitute the elusive dark mat-ter (DM) of our universe. They also comprise the search for charged-lepton, flavor-violating(cLFV) processes, whose observation would be a clear sign of new, fundamental physics. Am-bitious experiments are under construction in both these areas. These may very well lead tothe discovery of new phenomena that are not described by the standard model (SM) of particlephysics. For a recent review of DM direct detection experiments, see ( ), and ( ) for cLFVsearch experiments.On the theoretical side, the difficulty resides in making accurate predictions for the ratesexpected in those experiments and in interpreting their results. Atomic nuclei are composed ofprotons, p , and neutrons, n , collectively known as nucleons, N . To predict the expected ratesand interpret the measurements, one must quantify the interactions of the hypothetical new par-ticles with nucleons. The challenge is that nucleons are themselves complex, nonlinear boundstates of quarks, antiquarks and gluons. One must be able to accurately describe these boundstates to relate, what is observed in experiment, to the fundamental parameters that describe theinteractions of the new particles with quarks and antiquarks. These couplings are particularlyimportant for the many direct DM search experiments that attempt to detect WIMPs in the spin-independent channel. They are also related to the low-energy coupling of the Higgs boson, andof the ambient Higgs field, to nucleons. In that sense, they are also associated with contributionsto the mass of nucleons. These couplings are known as nucleon σ -terms.We compute the u , d , s and c -quark σ -terms of the proton and neutron, using ab initiocalculations based on quantum chromodynamics (QCD), the fundamental theory of the stronginteraction. The Euclidean Lagrangian of QCD is L = 1 / (2 g ) Tr G µν G µν + (cid:80) q ¯ q [ γ µ ( ∂ µ + iG µ ) + m q ] q , where γ µ are the Dirac matrices, q , the quark fields and index that runs over quarkflavors, and g is the coupling constant. The m q are the quark masses and G µν = ∂ µ G ν − ∂ ν G ν + i [ G µ , G ν ] , with G µ the gluon field, a × matrix in QCD. At the energies typical of quarksand gluons inside hadrons, this theory exhibits highly nonlinear behavior. Thus we introduce ahypercubic spacetime lattice on which the above Lagrangian is discretized ( ), keeping quarksup to and including the charm. The discretization puts quark variables on the lattice sites andgauge variables on the links between neighboring sites. The discretized theory is equivalent to afour-dimensional statistical physics system. The Feynman path integral, which is used to definethe quantum theory, can thus be evaluated numerically using powerful, importance-samplingmethods ( ). The calculation is performed for a number of lattice spacings a and spatial sizes L . The final results are obtained after taking the limits L → ∞ and a → .For the b and t quarks, we use a different strategy, based on a sequence of heavy-quarkeffective field theories (HQET) ( ). These are obtained from six-flavor QCD by sequentiallyintegrating out the most massive quark left in the theory, starting from the top. As long as themasses of the heavy quarks, Q = t, b, . . . , integrated out are much larger than the typical QCDscale, Λ QCD , this can be done systematically in perturbation theory, up to power correctionswhich begin at order (Λ QCD /m Q ) . The systematic error made in this approach is determined2y the powers of α s ( m Q ) and (Λ QCD /m Q ) appropriate for the order at which the calculation isperformed, with Q corresponding to the lightest quark integrated out.The nucleon σ -terms are conveniently parametrized by the dimensionless ratios, f Nud = m ud (cid:104) N | ¯ uu + ¯ dd | N (cid:105) / (2 M N ) = σ πN /M N and f Nq = m q (cid:104) N | ¯ qq | N (cid:105) / (2 M N ) = σ qN /M N ( ),where N can be either a p or an n , the initial and final nucleon states have identical mo-menta and are normalized relativistically, q can be any one of the six quark flavors, and m ud =( m u + m d ) / . Note that in the isospin limit, where m u = m d , f nud = f pud and f ns = f ps and wewill call these quantities f Nud and f Ns , respectively. σ -terms corresponding to the light u , d and s quarks have a long history, as they concernone of the earliest low-energy theorems established with current algebra ( ). This has led todeterminations of these σ -terms, based on πN -scattering data and chiral perturbation theory( – ). There also exist direct lattice calculations of the combined u and d σ -term ( – ), aswell as some of σ sN ( – ), and two of the charm ( , , ). However, these calculationsdo not exhibit full control over systematic errors or have too low statistical precision. Herewe perform high-statistics lattice calculations of the σ -terms for the four lightest quarks, inwhich all relevant limits are taken with full control over uncertainties. We use the Feynman-Hellmann theorem, which relates the sigma terms to the partial derivative of M N with respectto the corresponding quark mass. For q = ud, s , we proceed in two steps:1. We calculate the logarithmic derivatives of M N with respect to the π + and K χ mesonmasses, F NP ≡ ( ∂ ln M N /∂ ln M P ) , with P = π + , K χ . According to leading-order SU(3)chiral perturbation theory ( χ PT), M π + ∝ m ud and M K χ ≡ ( M K + + M K − − M π + ) / ∝ m s , so we expect these derivatives to be close to the desired σ -terms.2. We compute the Jacobian, J , for the coordinate transformation (ln M π + , ln M K χ ) → (ln m ud , ln m s ) and recover the σ -terms via ( f Nud , f Ns ) T = J · ( F Nπ + , F NK χ ) T .Thus, to obtain the σ -terms, we study M N as a function of M π + , M K χ and m c , as well as M π + and M K χ as functions of m ud and m s . For that purpose, we have used 33, N f = 1 +1 + 1 + 1 , 3HEX-smeared, clover-improved Wilson ensembles, with a total of , gaugeconfigurations and staggered ensembles, with a total of configurations, each separatedby 10 rational-hybrid-Monte-Carlo trajectories. The staggered ensembles have u, d, s and c masses straddling their physical values at three values of the bare gauge coupling, β = 6 /g ,corresponding to lattice spacings in the range a = 0 . − .
10 fm , and spatial extents around . The clover ensembles have pion masses in the range of −
420 MeV , the s and c masses straddling their physical values, four gauge coupling corresponding to lattice spacingsin the range a = 0 . − .
10 fm and spatial extents up to . Details are given in ( ). Inaddition, as M N is well-known from experiment ( ), we use it to fix the lattice spacing. Inthe u, d, s sector, the physical mass point is defined by setting M π + and M K χ to their physicalvalues ( ).In Fig. 1, we plot our clover determinations of M N versus the values of M π and M K χ ,together with the experimentally measured values for these quantities that define the physical3oint. We fit these results for M N to various polynomial, Padé and χ PT-motivated functions,of M π + and M K χ , that include discretization and finite-volume corrections ( ). From each fitwe determine the corresponding pair ( F Nπ + , F NK χ ) at the physical point. A similar study of M π + and M K χ versus m ud and m s determines the four matrix elements of J at the physical point ( ).Combining these results, as described above, yields the desired f Nud and f Ns . The statistical errorson these results are calculated using 2000 bootstrap samples. To obtain systematic errors, weconsider 6144 different analyses each leading to a result for f Nud and f Ns and to a total goodnessof fit ( ). These variations are chosen to probe the main sources of systematic error. The resultsare then combined into distributions whose means give our central values and whose widthsdetermine our systematic errors ( , , ) (see Table 1). Note that our analysis of the Jacobianalso provides a very precise determination of the mass ratio, m s /m ud = 27 . , wherethe errors are as in Table 1. From f Nud we further determine the individual u and d contributionsto the proton and neutron as in ( ) (see Table 1).To determine f Nc , we use 9 staggered ensembles at three β , corresponding to lattice spac-ings a = 0 . − .
12 fm . For each β , the ensembles feature three values of m c equal to (0 . , , . times the physical value, m ( φ ) c , with other quark masses held at their physicalvalues. We find that the most reliable way to obtain f Nc is to consider the ratios of the nucleoncorrelator at different m c . The exponential fall-off of this ratio gives the difference of M N atthese two values of m c . From these ratios we get the desired derivative either by a simple Taylorexpansion or by a heavy-quark-motivated expansion as detailed in ( ) (see Table 1).The b and t contributions are determined using the HQ expansions discussed above. Herewe use the next-to-next-to-next to leading order result of ( ) to reduce the correction termsto O ( α s ( m b ) , (Λ QCD /m b ) ) , leading to an uncertainty from the HQ expansion of around . ( ). The results are given in Table 1. With the same approach, we can also com-pute f Nc . But then the HQ expansion uncertainties on the c , b and t contributions becomes O ( α s ( m c ) , (Λ QCD /m c ) ) , i.e. of order . We obtain f Nc | HQ = 0 . , where the er-rors are as in Table 1. The difference to the full lattice result is f Nc − f Nc | HQ = 0 . ,in good agreement with the dimensional estimate of the possible HQ corrections, i.e. the HQexpansion works as expected here.In Fig. 2 we show our results for the σ -terms of the proton, as fractions of the total protonmass. This representation is renormalization scheme and scale independent, as described in ( ).The σ -terms of the neutron differ from those of the proton only in the small u and d contri-butions, at the present level of precision. In addition, we provide a computer code, based onour results, that allows to make predictions and to interpret results of low-energy, experimentalsearches for new fundamental physics. In particular, using it to sum all contributions, we obtainthe low-energy coupling of the Higgs to the nucleon, f hN = 0 . , in agreement with,but significantly more precise than the result in ( ).Beyond its importance for direct dark matter searches and charged-lepton flavor violation,our calculation also allows us to complete the quantitative picture of how protons and neutronsacquire mass, described, for instance, in ( , ). Instants after the Big Bang, the universe is ahot gas of elementary particles. Quarks and gluons are massless and interact through the strong4nteraction with a strength that has been measured at the Large Hadron Collider (LHC). Asthe universe expands and cools, it undergoes a transition ( , ) during which the Higgs fieldacquires a non-zero expectation value. At that point, elementary fermions get a mass throughtheir interactions with the background Higgs field. As the universe cools down further, in turntop, bottom and charm quarks and antiquarks vanish through decay and annihilation, and sub-sequently appear only as fleeting quantum fluctuations. The universe can then be describedby a theory of the strong interaction in which these particles are completely absent, as long astheir fluctuations are subsumed into an increase in the strength of the strong interaction. Thisis the context in which the calculations of ( , ) are performed. As the universe continuesto expand, it undergoes another transition, the QCD crossover ( ): the strongly interacting up,down, strange quarks and antiquarks and gluons become confined within bound states. In par-ticular, protons and neutrons form out of up and down quarks, with strange quarks, antiquarksand gluons contributing through fluctuations. As shown in ( ), roughly of the mass ofthese two bound states comes about due to the energy stored in the quantum fluctuations withinthem, while less than are induced by the up and down quark masses. The calculation per-formed here reminds us that, within this , (cid:80) q = t,b,c f Nq = 21 . of M N is actually dueto quantum fluctuations of the massive c , b , t quarks, with the remainder being due to gluon andstrange quark-antiquark fluctuations. And, as ( ) confirmed, the permil difference betweenneutron and proton mass arises from a subtle cancellation of electromagnetic effects and effectsdue to the difference of up and down quark masses.5 cknowledgments We are indebted to S. Collins, S. Dürr, J. Lavalle, H. Leutwyler and E. Nezri for informativediscussions and correspondence. Computations were performed on JUQUEEN and JURECA atForschungszentrum Jülich, on Turing at the Institute for Development and Resources in Inten-sive Scientific Computing (IDRIS) in Orsay, on SuperMUC at Leibniz Supercomputing Centrein München, on HazelHen at the High Performance Computing Center in Stuttgart. This projectwas supported, in part by the Excellence Initiative of Aix-Marseille University - A*MIDEX(ANR-11-IDEX-0001-02), a French “Investissements d’Avenir” program, through the Chaired’Excellence program and the OCEVU Laboratoire d’Excellence (ANR-11-LABX-0060), bythe DFG Grant SFB/TR55, by the Gauss Centre for Supercomputing e.V and by the GENCI-IDRIS supercomputing Grant No. 52275. 6ucleon Individual p and nf Nud f pu f Ns f pd f Nc f nu f Nb f nd f Nt σ -terms of the nucleons, in units of the corresponding nucleonmass. The statistical (SEM) and systematic uncertainties on the last digits are given in the firstand second set of parentheses, respectively. M π [MeV ] M N [ M e V ] a~0.1022 fma~0.1022 fm QCD+QEDa~0.0888 fma~0.0769 fma~0.0637 fm M K2 -M π /2 [MeV ] M N [ M e V ] a~0.1022 fma~0.1022 fm QCD+QEDa~0.0888 fma~0.0769 fma~0.0637 fm Figure 1: Example of the dependence of the nucleon mass on pion mass (left) and reduced-kaon mass (right) squared. The circles correspond to our simulation results for these massesat different values of the lattice spacing, indicated in the caption. The black cross in each plotcorresponds to the physical value for these masses, as given in ( , ). The black curves withgray bands represent a typical fit to our results, with error bands. The values of ( F Nπ + , F NK χ ) obtained from the fit are given by the slope of the curves at the black cross. Note that allsimulation points have been corrected, using the result of the fit, for the M K χ or M π and latticespacing and volume dependencies that are not shown. All data points represent the mean ± SEM. 7 q q proton uds cbt q u a r k s o n l y p r e s e n t a s v i r t u a l q u a r k - a n t i q u a r k p a i r s statistical uncertaintysystematic uncertainty Figure 2: σ -terms of the proton shown as fractions of the total proton mass. Small blackdots indicate the number of valence quarks per quark flavor. If the proton were an elementaryfermion, the sum of these terms would be the mass of the proton. However, our calculationshows that this contribution is only (offset part) of the total, with less than comingfrom valence flavors. 8 upplementary Material Let | N ( (cid:126)p, r ) (cid:105) be a relativistically-normalized, free nucleon state with three-momentum (cid:126)p , he-licity r and mass M N . It is an eigenstate of the renormalized QCD Hamiltonian, H QCD = (cid:82) d x H QCD ( x ) with eigenvalue E (cid:126)p = (cid:112) M N + (cid:126)p . Then, applied to the momentum-independentvariation of H QCD with respect to the renormalizaed mass of quark q , m Rq , the Feynman-Hellmann theorem ( – ) implies ∂M N ∂m Rq = 12 (cid:88) r =1 , (cid:104) N ( (cid:126)p, r ) | ∂ H QCD ∂m Rq ( x ) | N ( (cid:126)p, r ) (cid:105) c , (S1)where the subscript “ c ” indicates that we are referring to the connected part of the matrix ele-ment. Now, according to ( ): m Rq ∂ H QCD ∂m Rq ( x ) = m Rq (¯ qq ) R ( x ) . so that the σ -terms are given by σ Nq = ∂M N ∂ ln m Rq = 14 M N (cid:88) r =1 , (cid:104) N ( (cid:126)p, r ) | m Rq (¯ qq ) R ( x ) | N ( (cid:126)p, r ) (cid:105) c (S2) = (cid:104) N | m q ¯ qq | N (cid:105) / (2 M N ) , where, in the last line, we introduce a short-hand notation for the previous expression, whichtakes into account the fact that m q ¯ qq is renormalization scale and scheme dependent. From thisone can define scalar, quark contents: f Nq ≡ ∂ ln M N ∂ ln m Rq = σ Nq M N (S3)Note that both the σ − terms and the quark contents are renormalization scheme and scale inde-pendent. We thank Heiri Leutwyler for correspondence on the formulation of the Feynman-Hellmann theorem in quan-tum field theory. .2 Nucleon-Higgs coupling in the standard model In the standard model, the Higgs field φ obtains a vacuum expectation value φ = v √ |(cid:104) | φ | (cid:105)| . The quark masses, like all fundamental fermion masses, originate from a Yukawa term g f ¯ f φf in the fundamental Lagrangian, where g f is that fermion’s Yukawa coupling and the fermionmass is given by m f = g f φ . (S4)Expanding φ = φ + h around its vacuum expectation value, it is obvious that the Yukawaterm is also responsible for the fermion-Higgs interaction, g f ¯ f hf which, in principle, can bemeasured independently. We can thus form the ratio f f = g f φ m f , (S5)which the standard model predicts to be one for fundamental fermions. A similar ratio can beformed for fundamental gauge bosons and it has been verified experimentally to be one for the b and t quarks, the τ , W and Z ( ). We can slightly rewrite this expression by noting that,according to (S4), ∂m f (cid:48) ∂ ln g f = g f δ ff (cid:48) φ , so that f f (cid:48) f = ∂ ln m f (cid:48) ∂ ln g f (S6)is the fraction of the mass of fermion f (cid:48) that couples to the Higgs field via the Yukawa coupling g f .The advantage of the definition (S6) is that it generalizes naturally. Returning to (S3), wenote that the scalar quark content is, in terms of fundamental parameters of the standard model: f Nq = ∂ ln M N ∂ ln g q . Thus, it can be interpreted as the mass fraction of the nucleon that couples to the Higgs field viathe Yukawa coupling g q of quark q . Within the framework of the standard model, one may alsointerpret f Nq as the fraction of the nucleon mass originating from the Higgs field via the Yukawacoupling g q . This definition is scale and scheme independent, but it is not unique and we willtherefore not pursue it any further.A large part of the nucleon mass does not couple to the Higgs field via the g q . The bulk ofthis “rest” originates from QCD dynamics and there are standard techniques for decomposingthis contribution further ( , , ). Quantum electrodynamics contributes at the permil level2 ) and there are even smaller contributions due to the weak interaction and the lepton sector,including its coupling to the Higgs. Since the latter two do not take part in the strong interaction,these contributions are suppressed by factors α ∼ − for the lepton sector and G F Λ ∼ − for the weak dynamics respectively. They can safely be ignored at the precision of ourcurrent results.Note that the sum over the scalar quark contents of the nucleon, f Nh = (cid:88) q f Nq , (S7)is an observable quantity. It denotes the strength of the coupling of the Higgs to nucleons, inthe limit of vanishing momentum transfer. However, measuring it experimentally is a hugechallenge, since this requires isolating the highly-suppressed, Higgs exchange contribution tothe low-momentum-transfer scattering of a standard model particle off a nucleon. Of course, inHiggs portal dark-matter models, these interactions are the primary detection channel. In otherextensions of the standard model, scalar exchange interactions with nucleons may require linearcombinations of scalar quark contents that differ from the one given in (S7). For this reason, inSec. 8 we provide a tool to compute arbitrary linear combinations of quark contents from ourresults, taking the full correlation matrix into account. Different parts of the analysis rely on different gauge configurations. Some were obtained witha Wilson fermion action and others with a staggered fermion one. These actions were chosenbecause each has properties that are better suited to different aspects of the calculation.
In order to compute the dependence of the nucleon mass on the meson masses used to interpolateto the physical u , d and s quark mass point, we use N f = 4 × , 3HEX-smeared, clover-improved Wilson fermions and a tree-level improved Symanzik gauge action (for details see( )). Compared to ( ), two new ensembles were generated at the finest lattice spacing, withstrange quark masses significantly differing from the physical value, so as to give us a largerlever arm in the strange quark mass direction. The full list of these ensembles is provided inTable S1. We mainly use the pure QCD ( e = 0 ) ensembles. Pseudoscalar meson masses, M P are deter-mined from a fit to multiple, zero-three-momentum correlators with a common value of M P : C P P ( t ) = (cid:104) P ( t ) P † (0) (cid:105) and C A P ( t ) = (cid:104) A ( t ) P † (0) (cid:105) , (S8)3 /g e am u am d am s L × T m π [MeV] m π L × trajectories3.2 0 -0.0686 -0.0674 -0.068 ×
413 6.9 13.2 0 -0.0737 -0.0723 -0.058 ×
353 5.9 43.2 0 -0.0733 -0.0727 -0.058 ×
356 5.8 13.2 0 -0.0776 -0.0764 -0.05 ×
294 4.9 43.2 0 -0.0805 -0.0795 -0.044 ×
238 4.0 123.2 0 -0.0806 -0.0794 -0.033 ×
266 4.4 123.2 0 -0.0686 -0.0674 -0.02 ×
488 8.1 43.2 0 -0.0737 -0.0723 -0.025 ×
411 6.8 43.2 0 -0.0776 -0.0764 -0.029 ×
336 5.6 43.2 0 -0.077 -0.0643 -0.0297 ×
438 7.3 43.2 0 -0.073 -0.0629 -0.0351 ×
469 7.8 43.2 0 -0.077 -0.0669 -0.0391 ×
405 6.7 43.2 1.00 -0.0859 -0.0792 -0.0522 ×
298 3.7 53.2 1.00 -0.0859 -0.0792 -0.0522 ×
295 4.9 43.2 1.00 -0.0859 -0.0792 -0.0522 ×
295 7.3 43.2 1.00 -0.0859 -0.0792 -0.0522 ×
295 12.2 13.3 0 -0.0486 -0.0474 -0.048 ×
422 6.1 13.3 0 -0.0537 -0.0523 -0.038 ×
348 5.1 23.3 0 -0.0535 -0.0525 -0.038 ×
349 5.0 23.3 0 -0.0576 -0.0564 -0.03 ×
275 4.0 123.3 0 -0.0576 -0.0564 -0.019 ×
293 4.2 123.3 0 -0.0606 -0.0594 -0.024 ×
200 4.3 203.4 0 -0.034 -0.033 -0.0335 ×
403 5.0 43.4 0 -0.0385 -0.0375 -0.0245 ×
321 4.0 43.4 0 -0.0423 -0.0417 -0.0165 ×
219 4.1 43.5 0 -0.0218 -0.0212 -0.0215 ×
426 4.4 43.5 0 -0.0254 -0.0246 -0.0145 ×
348 5.4 43.5 0 -0.0268 -0.0262 -0.0115 ×
310 4.8 83.5 0 -0.0269 -0.0261 -0.0031 ×
317 4.9 83.5 0 -0.0285 -0.0275 -0.0085 ×
266 4.1 83.5 0 -0.0302 -0.0294 -0.0049 ×
199 4.1 43.5 0 -0.027 -0.027 -0.027 ×
280 4.4 33.5 0 -0.028 -0.028 +0.009 ×
282 4.4 3.5Table S1: List of 3HEX clover ensembles used to determine the meson-mass dependence ofthe proton mass. 4ith P = ¯ qγ q (cid:48) , A = ¯ qγ γ q (cid:48) and q and q (cid:48) are distinct quark flavors chosen amongst u , d and s . We fit these correlation functions to C P P ( t ) = d cosh( M P ( t − T / and C A P ( t ) = f sinh( M P ( t − T / , (S9)where T is the lattice temporal extent. Nucleon correlation functions are those of ( ). Thecorresponding masses were obtained by fits to single, decaying exponentials.Fit ranges are five lattice spacings long and the initial time slices, t i , are matched, on differ-ent ensembles, to have a constant ratio between estimated excited-state effects and the relative,statistical error on the nucleon mass. Specifically, we take the smallest t i > t − ln( (cid:15) i ) / ∆ M ,where (cid:15) i is the statistical error of the nucleon mass extracted with initial time slice t i , ∆ M =500 MeV is an estimate of the mass difference to the first excited nucleon state and t is aparameter which we vary in our analysis from − . to − .
25 fm . t KS prob. M p KS prob. M n -1.40 0.59 0.86-1.35 0.87 0.99-1.30 0.54 0.85-1.25 0.32 0.60Table S2: The Kolmogorov-Smirnov probability of the CDF of the proton and neutron mass fitqualities compared to a uniform distribution for various values of the offset t (see text).Table S2 lists the Kolmogorov-Smirnov (KS) probabilities from a comparison of the CDFsof the nucleon mass fit qualities to a uniform distribution as detailed in ( ). KS probabilitiesfor the multi-channel meson fits are around . because of their highly correlated nature. Fittingchannels individually however gives fully compatible results and KS probabilities that are in thesame range as those for M p and M n . Therefore, we are confident that, by considering the fullrange of offset values, t , from Table S2 we obtain a conservative estimate of remnant excitedstate effects.For the extrapolation to the physical point, we used the ensembles reported in Table S1. Toget a better handle on the finite-volume dependence of the nucleon mass, we added the four e = 1 , β = 3 . ensembles at M π ∼
290 MeV that differ only in their volume L × T , with L ∈ { , , , } , compared to our set of pure QCD ensembles. Since these ensembleswere tuned to the isospin symmetric point, we can use the neutron mass instead of the nucleonmass, the connected-pseudoscalar-meson mass average, ( M uu + M dd ) / , instead of M π and M K − M dd instead of M K χ . We checked that the exact offset value does not substantially change our results. /g am ud am s am c L × T m π L × trajectories3.84 0.00151556 0.0431935 0.511843 × × × × × × × × × × × We use staggered ensembles near the physical point to extract the dependence of the pseu-doscalar meson masses, M P , on the quark masses, m q . The full list of these ensembles used inthe analysis is provided in Table S3.The Goldstone-pion mass is close to the physical one in all of these ensembles. The volumesare matched and M π L > , which implies that the finite-volume corrections to the pion massand decay constant are below the permil level ( ). The kaon mass and decay constant receiveeven smaller corrections. Since we do not determine any other observables from these config-urations, an infinite-volume extrapolation is not necessary in this part of the analysis. Furtherdetails on the action used can be found in ( , ).The dependence of the nucleon mass on the charm quark mass is obtained directly from a setof nine 4-stout-smeared staggered ensembles (“charm ensembles”), at three values of β = 3 . , . , . . At each β one ensemble, which we call the central ensemble, was tuned to thephysical point with a deviation of less then 4 % in M π /f π , (2 M K − M π ) /f π and the bare charmquark mass was set to . times the bare strange-quark mass ( ). Two other ensembles wereobtained from each central ensemble by varying the bare charm-quark mass to . and . times the value on the central ensemble and leaving all other parameters fixed. In each of these9 ensembles, we have generated 64 configurations separated by trajectories each.The sources and corresponding propagators are the standard ones provided, for instance, bythe MILC code ( ). 6 .4 Determination of hadron masses on the staggered ensembles For the ensembles in Table S3, we obtain pseudoscalar masses and decay constants from fitsto the pseudoscalar propagators starting at a fixed (in physical units) t min = 1 . or t min =2 . . Fit ranges are ten lattice spacings long. In Table S4 we list the KS probabilities thatresult from comparing the CDFs of the pion and kaon fit qualities to a uniform random distri-bution. t min [fm] KS prob. M π KS prob. M K t min .To extract masses from a staggered propagator, c t , we either use a standard staggered twostate fit or, alternatively, we construct a time-shifted propagator: d t ( M ) ≡ c t + e M c t +1 (S10)The time-shifted propagator is useful if there is a region in t with negligible backward contribu-tions, m ( T / − t ) (cid:29) , and negligible excited states except the parity partner of the nucleon.Here m is the ground state mass and T the lattice temporal extent. The staggered propagator inthis region behaves as c t = e − mt (cid:0) c + ( − t c e − ∆ t (cid:1) , (S11)with the mass difference to the staggered parity partner, ∆ , and the matrix elements c and c of the ground state and staggered parity partner, respectively.The contribution of the staggered parity partner to the time-shifted propagator (S10) can, inprinciple, be cancelled out by setting M = m + ∆ : d t ( m + ∆) = c e − mt (1 + e ∆ ) . (S12)We can determine M self-consistently by defining a local effective mass l t ( M ) ≡ ln( d t ( M ) /d t +1 ( M )) , (S13)as well as its average, ¯ l ( M ) = 1 n t + n − (cid:88) t =0 l t ( m ) (S14)in the signal region, t ≤ t ≤ t + n , and by minimizing t + n − (cid:88) t ,t =0 ( l t ( M ) − ¯ l ( M ))( C − ) t t ( l t ( M ) − ¯ l ( M )) (S15)7ith respect to M , where C is the correlation matrix corresponding to l t ( M ) − ¯ l ( M ) . Theresulting parameter is M opt and we call the corresponding propagator, d t ( M opt ) , the optimaltime-shifted propagator.It is important to note that M opt is just one number per ensemble that fixes the relativecontributions of c t and c t +1 in the time-shifted propagator. In particular, it does not directly enterany further stages of analysis. After having determined M opt , we proceed to extract the groundstate mass from d t ( M opt ) in a standard fashion. In fact, the time-shifted propagator d t ( M ) does have the correct asymptotic time behavior for any constant M and thus an inaccuratedetermination of M opt does not invalidate the ground-state-mass extraction from d t ( M opt ) . Ifcancellation of the staggered parity partner is not achieved to within the statistical accuracy ofthe correlator, or additional excited states or backward contributions are present, the consequentfit to determine the ground state mass from d t ( M opt ) will simply fail with a bad fit quality.On our charm ensembles, we determine mass differences from ratios of optimal time-shiftedpropagators. The plateau length is always 8 lattice spacings and the plateau start is either t min =0 . or t min = 1 . . σ -terms We determine the dependence of the nucleon mass on the pseudoscalar meson masses fromour 3HEX ensembles (Table S1). We define the isospin symmetric physical point of QCD by M N = ( M p + M n ) / .
919 MeV ( ) and M π = 134 . , M K χ = 685 . ( ). Inorder to estimate the dependence of our result on the range of the chiral expansion, we imposedtwo different cuts on the maximal pion mass M π ≤ /
420 MeV entering our analysis.
The nucleon mass is extrapolated to the physical point with an ansatz M N = M φN × (cid:89) i (cid:16) c i ( v i − v ( φ ) i ) (cid:17) t i , (S16)where the c i are the fit parameters, v i the fit variables, t i = ± (Taylor/Padé) and generically x ( φ ) denotes the physical value of the observable x . We perform fits to functions that all containthe pseudoscalar mass dependencies v = M π , with t = 1 , and v = M K χ , with either t = ± .They also include the finite-volume term v = M / π L − / e − M π L , with t = 1 . In addition, theycontain either of the two next-to-leading-order terms in M π , namely v = M π with t = 1 , or v = M π with t = 1 , corresponding to either a chiral or generic Taylor expansion in M π ( ).They further include either no discretization term or the formally leading discretization terms v = α s a ( M π − M ( φ )2 π ) with t = 1 and v = α s a ( M K χ − M ( φ )2 K χ ) with t = 1 ; or the formallysubleading v = a ( M π − M ( φ )2 π ) , t = 1 , and v = a ( M K χ − M ( φ )2 K χ ) , t = 1 ( ). Note that8e set the scale with M N , so that discretization terms proportional to just α s a or a would beredundant. The total number of distinct analyses is 4(plateau ranges) × M π cut) × M π ) × M K χ ) × = 96 . The way in which we combinethese analyses to give a final central value and systematic error is described in Sec. 3.4. InFig. S1 we present the variation of our final observables, f Nud and f Ns , that result from applyingthese different fit procedures. Results from all different fit procedures are in good agreement.The leading sources of systematic error come from the continuum limit and the pion mass cuts. From the nucleon fits we can extract the lattice spacings, which are reported in Table S5. Acomparison to ( ), where M Ω was used to set the scale on a set of ensembles which has a largeoverlap with the current one, reveals perfect agreement. β a [fm] c =35(13)(5)GeV − is compatible with the numerical predictions of ( ). To further investigatepossible finite-volume effects beyond the leading order, we performed two complete auxiliaryanalyses with explicit pion finite-volume effects according to ( ), one with fixed and one withfitted prefactors. The influence on central values and errors was found to be insignificant and,in the case of the fitted prefactor, the additional term was found to be compatible with zero. The nucleon fits have fit qualities in the range Q = 0 . − . , with an average fit quality ¯ Q = 0 . . In Figs. 1 and S2 we display, for one sample fit, the dependence of the nucleon masson M π , M K χ and the inverse of the spatial lattice extent L . Defining pion and reduced kaon9 .03 0.04 0.05 0.05 0.06 0.07 f s t =-1.25 Q=0.40t =-1.3 Q=0.36t =-1.35 Q=0.42t =-1.4 Q=0.24no a Q=0.38 α a Q=0.34a Q=0.34M π <360 MeV Q=0.37M π <420 MeV Q=0.34M π Q=0.36M π Q=0.35M K χ Pade Q=0.33M K χ Taylor Q=0.38t min =1.9fm Q=0.26t min =2.3fm Q=0.45no a m s Q=0.23a M π m s Q=0.46a M K χ m s Q=0.31a (M π ,M K χ ) m s Q=0.43a m ud Q=0.38a (1,M π ) m ud Q=0.33a (1,M K χ ) m ud Q=0.37a (1,M π ,M K χ ) m ud Q=0.35no M π M K χ mix Q=0.37M π M K χ mix Q=0.31TOTAL Q=0.36 f ud π ,K χ N N
Figure S1: Variation of f Nud and f Ns from different restrictions on the full analysis procedure.The top part of the panel corresponds to variations in the determination of the meson massdependence of the nucleon mass. From top to bottom these are: restriction to a single plateaurange, to a single scaling behavior, to a single pion cut, to either Taylor or chiral expansion in M π and to Taylor or Padé expansion in M K χ . The middle part displays the effect of variationsin the computation of the elements of the mixing matrix J . From top to bottom these are:restriction to a single plateau range, to a single continuum fit form for m s , to a single continuumfit form for m ud and to the inclusion/exclusion of crossterms in the m ud , m s and scale settingfits. The last row shows the final result including all fits. The average fit quality of the analyses ¯ Q is given in each case. The percentage of the total systematic error due to each variation aloneis also displayed (percentages add up to in quadrature).10 /10 1/6 1/4 1/3 L -1 [fm -1 ] M N [ M e V ] a~0.1022 fma~0.1022 fm QCD+QEDa~0.0888 fma~0.0769 fma~0.0637 fm Figure S2: Lattice-spatial-extent dependence of the nucleon mass in one sample fit. The fit hasbeen used to shift the simulation results to the physical point in all of the other variables onwhich these results depend. σ -terms as, σ πN = M π ∂M N ∂M π (cid:12)(cid:12)(cid:12)(cid:12) M Kχ and σ K χ N = M K χ ∂M N ∂M K χ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) M π , (S17)we obtain σ πN = 42 . . .
4) MeV and σ K χ N = 50 . . .
8) MeV . We also define thelogarithmic derivatives, F Nπ = σ πN M N = ∂ ln M N ∂ ln M π (cid:12)(cid:12)(cid:12)(cid:12) M Kχ and F NK χ = σ K χ N M N = ∂ ln M N ∂ ln M K χ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) M π , (S18)for which our final results are F Nπ = 0 . and F NK χ = 0 . . Fromthese fits we also obtain an estimate of the nucleon mass at vanishing M π keeping M K χ fixed,11 N | M π =0 = 896(1)(2) MeV , and, conversely, an estimate of the nucleon mass at vanishing M K χ keeping M π fixed, M N | M Kχ =0 = 889(3)(3) MeV . These results, as well as our value forthe nucleon mass in the SU (3) chiral limit M N | M π = M Kχ =0 = 848(3)(3) MeV , should be inter-preted with great care, since it is not clear that our interpolation and extrapolation functions arevalid down to vanishing quark masses. To leading order in chiral perturbation theory, M π and M K χ are proportional to m ud and m s .The logarithmic derivatives F Nπ and F NK χ defined in (S18) are therefore equal to the nucleonquark content, f ud = ∂ ln M N ∂ ln m ud (cid:12)(cid:12)(cid:12)(cid:12) m s and f s = ∂ ln M N ∂ ln m s (cid:12)(cid:12)(cid:12)(cid:12) m ud (S19)in leading order. To go beyond leading order, we need to determine the relation between pseu-doscalar meson masses and quark masses around the physical point. For this purpose we usethe N f = 2 + 1 + 1 staggered ensembles presented in Sec. 2.3, which bracket the physical point. Transforming the F NP , P ∈ { π, K χ } into the f q , q ∈ { ud, s } is achieved via the Jacobian matrix, J , with elements J P,q = ∂ ln M P ∂ ln m q , (S20)evaluated at the physical point so that f q = ∂ ln M N ∂ ln m q = (cid:88) P ∂ ln M N ∂ ln M P J P,q = (cid:88) P F NP J P,q . (S21)While the F NP are defined at the physical point, including the charm quark mass, the N f =2+1+1 ensembles, on which we calculate the Jacobian, are at a variety of different charm quarkmasses. Thus, when discussing the Jacobian, we are careful about recording m c dependence.Since the computation of this Jacobian is best done with staggered fermions, whose massesare multiplicatively renormalized, we parametrize this dependence by the ratio r cs = m c /m s ,where renormalization factors cancel in mass-independent schemes. Then we fit the light andstrange quark mass as functions of M π , M K χ and r cs : m ud ( M π , M K χ , r cs ) and m s ( M π , M K χ , r cs ) . (S22)12rom these fit functions, we determine the following logarithmic derivatives at the physicalpoint: N q,π = ∂ ln m q ∂ ln M π (cid:12)(cid:12)(cid:12)(cid:12) M Kχ ,r cs , N q,K χ = ∂ ln m q ∂ ln M K χ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) M π ,r cs , N q,r = ∂ ln m q ∂ ln r cs (cid:12)(cid:12)(cid:12)(cid:12) M π ,M Kχ . (S23)In principle we need to compute the Jacobian between the two sets of parameters ( M π , M K χ , r cs ) and ( m ud , m s , m c ) . However, the elements of the Jacobian matrix that we are ultimately inter-ested in are those at fixed m c , i.e. J P,ud = ∂ ln M P ∂ ln m ud (cid:12)(cid:12)(cid:12)(cid:12) m s ,m c and J P,s = ∂ ln M P ∂ ln m s (cid:12)(cid:12)(cid:12)(cid:12) m ud ,m c . (S24)They are obtained from (S23) as J π,ud = N s,K χ N ud,π N s,K χ − N ud,K χ N s,π ,J π,s = − N ud,K χ (1 − N s,r ) − N s,K χ N ud,r N ud,π N s,K χ − N ud,K χ N s,π ,J K χ ,ud = − N s,π N ud,π N s,K χ − N ud,K χ N s,π ,J K χ ,s = N ud,π (1 − N s,r ) − N s,π N ud,r N ud,π N s,K χ − N ud,K χ N s,π . (S25) The matrix J is a scheme and scale independent quantity. Each element has the form of m q /M P × ∂M P /∂m q , where m q is a quark mass and M P , either M π + or M K χ . Any multi-plicative renormalization factor cancels in these ratios provided that it is independent of thequark masses. However, to allow for a global fit of m ud and m s to the functions of (S22), in-volving ensembles with different lattice spacings, a specific scheme has to be introduced. Wechose the simplest possible scheme: fixing the value of the renormalized strange quark massat the physical point to a constant independent of quark masses and lattice spacing. The exactprocedure is detailed in Sec. 4.3 In order to set the scale, we interpolate the pion decay constant, f π , to the physical mass pointdefined by the physical values, M ( φ )2 π and M ( φ )2 K χ , of the squared pion and reduced kaon masses,and by r ( φ ) cs , which we set to . . Since all of our results are dimensionless ratios, scale setting13s not critical and especially a variation in r ( φ ) cs has negligible effect on the result. Varying theinterpolation function has no effect on our results as is evident from Fig. S1.To determine J we study the behavior of the functions m ud ( M π , M K χ , r cs ) and m s ( M π ,M K χ , r cs ) . We expand these around the physical mass point, again defined by M ( φ )2 π , M ( φ )2 K χ and r ( φ ) cs . The fit functions we employ have the form m Rud ( M π , M K χ , r cs ) = c ud + c ud ∆ π + c ud ∆ K χ + c ud ∆ π + c ud ∆ K χ + c ud ∆ π ∆ K χ + c udc ∆ r cs (S26)and m Rs ( M π , M K χ , r cs ) = c s + c s ∆ π + c s ∆ K χ + c s ∆ π + c s ∆ K χ + c s ∆ π ∆ K χ + c sc ∆ r cs , (S27)where ∆ π = M π − M ( φ )2 π , ∆ K χ = M K χ − M ( φ )2 K χ , ∆ r cs = r cs − r ( φ ) cs and m Rq are renormalizedquark masses. They are related to the bare quark masses m q , which are parameters of theaction, by a multiplicative renormalization factor, Z , so that m Rud = Zm ud ( M π , M K χ , r cs ) and m Rs = Zm s ( M π , M K χ , r cs ) . We choose a scheme in which the renormalization factors do notdepend on the quark masses, and consequently the squared pion and reduced kaon masses,but only on the gauge coupling β . To avoid an explicit determination of the renormalizationfactors, we divide the above expansions by the value of the renormalized strange quark mass atthe physical mass point for each β , yielding: m ud ( M π , M K χ , r cs ) m s ( M ( φ )2 π , M ( φ )2 K χ , r ( φ ) cs ) = r + d ud ∆ π + d ud ∆ K χ + d ud ∆ π + d ud ∆ K χ + d ud ∆ π ∆ K χ + d udc ∆ r cs (S28)and m s ( M π , M K χ , r cs ) m s ( M ( φ )2 π , M ( φ )2 K χ , r ( φ ) cs ) = 1 + d s ∆ π + d s ∆ K χ + d s ∆ π + d s ∆ K χ + d s ∆ π ∆ K χ + d sc ∆ r cs , (S29)where r is the ratio of the strange and to light quark mass at the physical point and d qij = c qij /m s ( M ( φ )2 π , M ( φ )2 K χ , r ( φ ) cs ) , for q = ud, s and for all values of ij appearing in Eqs. (S28) and(S29). Note that all renormalization factors cancel out. The value m s ( M ( φ )2 π , M ( φ )2 K χ , r ( φ ) cs ) isdifferent for each gauge coupling β . This is a manifestation of the need for renormalization. Inpractice, we introduce fit parameters m ( φ ) s [ β ] , where β in rectangular brackets indicates that oneindependent parameter per gauge coupling is considered.The above expansions are subject to discretization artefacts. We consider these to contributeto the ratio r and to the linear terms in ∆ π and ∆ s . Discretization artefacts on the quadraticterms are subleading. The final fit functions are m ud ( M π , M K χ , r cs ) m ( φ ) s [ β ] = r + a r + ( d ud + e ud a )∆ π + ( d ud + e ud a )∆ K χ + d ud ∆ π + d ud ∆ K χ + d ud ∆ π ∆ K χ + d udc ∆ r cs (S30)14nd m s ( M π , M K χ , r cs ) m ( φ ) s [ β ] = 1 + ( d s + e s a )∆ π + ( d s + e s a )∆ K χ + d s ∆ π + d s ∆ K χ + d s ∆ π ∆ K χ + d sc ∆ r cs . (S31)We perform various fits, all of which contain the fit parameters r , r , d q , d q and d qc . Thediscretization terms e q , e q and the higher order interpolation terms d q are optionally presentand we perform fits with all possible variations, with each term either present or absent. Theadditional, higher-order terms, d q and d q , turn out to be irrelevant and can be omitted entirely.Finite-volume effects on all relevant quantities are safely below 1 permil on all of our en-sembles ( ) and cannot be detected with the statistical accuracy of our data.The M π and M K χ dependence of the quark masses for a single representative analysis aredisplayed in Fig. S3. The total number of distinct analyses for the Jacobi matrix is 2(plateau ranges) × m s ) × m ud ) × f Nud and f Ns ,resulting from these different fit procedures. Results from all fit procedures are in good agree-ment, with the leading contribution towards the systematic error coming from the variation ofthe continuum extrapolation terms in m s . One of the fit parameters in (S30), r , is the ratio of light to strange quark masses at the physicalpoint. This is a phenomenologically interesting parameter that we can determine. It provides anadditional crosscheck of our fit procedure. For its inverse, the strange to light quark mass ratio,we obtain m s m ud = 27 . (S32)which is in good agreement with the current PDG world average . from ( ) and with themost precise lattice values . +10 − ) from ( ) and . from ( ). The continuumextrapolation of this quantity for a single representative analysis is shown in Fig. S4. The quark mass fits have fit qualities in the range Q = 0 . − . , with an average fit qualityof Q = 0 . . The final results on the elements of the mixing matrix are given in Fig. S5.15 .0370.0380.039 m ud / m s φ M π [MeV ] m s / m s φ a~0.093 fm M K χ [MeV ] a~0.063 fma~0.077 fm Figure S3: Dependence of the quark masses, m ud and m s , on the squared pseudoscalar-mesonmasses, M π and M K χ , for one of our analyses. Data points are corrected using the fit to tuneall, except the plotted variables, to their physical value.16 .08 a [fm ] m ud / m s Figure S4: Example of a continuum extrapolation of the quark mass ratio m s /m ud .17 .0 0.2 0.4 0.6 0.8 1.0 log M π log m ud = 0 . log M K χ log m ud = − . log M π log m s = 0 . log M K χ log m s = 1 . LO χ PTstatistical errorsystematic error
Figure S5: Final results for the elements of the Jacobian matrix (S20). There are only smallcorrections to the leading order χ PT predictions.Evidently, the mixing matrix provides only a small correction to the leading order prediction ofchiral perturbation theory, where the mixing matrix is the identity.Applying the Jacobian matrix to the vector, ( F Nπ + , F NK χ ) T , we obtain the final results listedin Table 1. We use a weighted average and standard deviation of the results from the individualanalyses to determine central value and systematic errors ( , ). For our main result, we usethe Akaike information criterion (AIC) ( ) to determine the relative weight of the analyses. Inorder to check that the choice of weight does not significantly alter the result, we have plottedthe cumulative distribution function of f Nud and f Ns in Figs. S6 and S7, with a flat, unit weight,with a weight equal to the quality of fit, Q , and with the AIC weight. The choice of weight doesnot substantially influence the final result.Our results are largely compatible with other recent lattice determinations ( , , , – ) as well as with the seminal calculations by Gasser, Leutwyler and Sainio ( ), while recentphenomenological determinations that obtain f Nud from πN scattering data ( , , – ) givesomewhat higher values. Scalar quark contents, the nucleon mass and the QCD trace anomaly are related by a sum rule( ). This sum rule allows one to compute the heavy-quark contents from the light-quark ones,in the heavy-quark limit. This is achieved by considering a succession of effective field theoriesof QCD in which the heavy top, bottom and charm quarks are integrated out in turn. Thus,18 .03 0.04 0.05 f udN r e l a t i v e w e i gh t AIC weightQ weightflat weight
Figure S6: The cumulative distribution function of f Nud over all analyses, obtained using threedifferent weight functions. In the top part of the figure, the resulting central values and totalerrors, corresponding to each weight function, are plotted.19 .05 0.06 0.07 f sN r e l a t i v e w e i gh t AIC weightQ weightflat weight
Figure S7: The cumulative distribution function of f Ns over all analyses, obtained using threedifferent weight functions. In the top part of the figure, the resulting central values and totalerrors, corresponding to each weight function, are plotted.20onsider QCD with N f flavors treated as light and with one heavy quark, Q , to be integratedout. Then, matching the N f + 1 and N f theories at scale m Q yields( ) yields f NQ = 23 β (1 − ¯ f N f ) (cid:34) O (cid:32)(cid:18) Λ QCD m Q (cid:19) , α s ( m Q ) (cid:33)(cid:35) (S33)where β = 11 − (2 / N f is the leading-order coefficient of the QCD β -function, in N f -flavorQCD, and ¯ f N f = N f (cid:88) i =1 f Nq i (S34)is the sum of the scalar quark contents over the N f quark flavors, q i , not integrated out. Recently,higher-order QCD corrections to (S33) have been computed to O ( α s ) ( ), yielding f Q = (cid:88) n =0 ( b N f n − c N f n ¯ f N f ) α ns (cid:34) O (cid:32)(cid:18) Λ QCD m Q (cid:19) , α s ( m Q ) (cid:33)(cid:35) ) (S35)where N f = 5 for Q = t , N f = 4 for Q = b and N f = 3 for Q = c . The numerical values ofthe b n and c n are given in Table S6 for N f = 3 , and . n b n c n b n c n b n c n ).We leverage these results in various ways. First, we can obviously use our results for thelight quark contents to determine various ¯ f N f and thus the heavy-quark contents. This methodis used to compute f Nb and f Nt in Sec. 7.Now, if only the top quark is integrated out, the neglected corrections in (S35) are neglible.If also the bottom quark is removed from the dynamics, these corrections are still very smallsince (Λ QCD /m b ) ∼ . ( α s ( m b ) is smaller). The heavy-quark expansion is possibly less-well behaved for the charm. Indeed, (Λ QCD /m c ) ∼ ( α s ( m c ) < ∼ is smaller), whichis no longer negligible and is comparable to our lattice uncertainties. Thus, we also compute f Nc directly on the lattice, using the Feynman-Hellmann theorem. As shown in Sec. 6, theresult obtained is compatible with the one given by the heavy-quark expansion, within the naiveestimate of O ((Λ QCD /m c ) ) corrections. This indicates that, even for the charm, the heavy-quark expansion behaves as expected. 21he heavy-quark expansion also helps in the direct determination of the scalar, charm-quarkcontent on the lattice. Since the nucleon mass depends on the charm-quark mass only verymildly, we vary this mass by ± around its physical value and measure the correspondingvariations of the nucleon mass: ∆ + = M N ( m c = m φc ) − M N ( m c = m φc ) , ∆ − = M N ( m c = m φc ) − M N ( m c = m φc ) . (S36)These are then combined to determine the quantity ˆ f Nc = 2 ∆ + + ∆ − M N ( m c = m φc ) (S37)which, using a second order Taylor series expansion in m c around the physical point, can beshown to equal the scalar, charm-quark content f Nc , up to terms of order ( δm c /m c ) = 1 / .Insight from the heavy-quark expansion provides an alternative expansion. First we notethat α s changes roughly between . and . when changing the scale from / m φc to / m φc ( ), implying a relative variation of f Nc in this range by ∼ according to (S35). To a goodapproximation f Nc ( m c ) , in this range, is therefore constant, i.e. f Nc ( km c ) = ∂∂ ln m ln M ( m ) (cid:12)(cid:12)(cid:12)(cid:12) m = km c (cid:39) cst (S38)for k = 0 . or k = 1 . . This implies M ( m c ) = M ( m φc ) (cid:18) m c m φc (cid:19) f Nc . (S39)Taylor expanding this expression to second order around f c = 0 and plugging the result into(S36), we find that the quantity, ˜ f Nc = 1 M ( m φc ) ln ln ln (cid:32)(cid:18) ln 43 (cid:19) ∆ + + (cid:18) ln 54 (cid:19) ∆ − (cid:33) , (S40)is equal to the charm-quark content f Nc , up to terms of order ( f Nc ) ∼ × − , assuming f Nc tobe constant in the ± region around the physical charm quark mass. Since the heavy-quarkexpansion tells us that this assumption is valid to O (3%) , we conclude that ˜ f Nc = f Nc + O (3%) .The difference between ˆ f Nc extracted with the Taylor series (S37) and ˜ f Nc extracted with theheavy quark expansion (S40) provides an estimate of the systematic error due to replacing thederivative of the nucleon mass, with respect to m c , by a finite difference.22 Lattice computation of the c sigma term The direct computation of the charm-quark content from our lattice ensembles poses a differentset of challenges. On the one hand, locating the physical point precisely is not critical, as de-tailed in the previous section. On the other hand, one needs to vary the charm-quark mass over asignificant range to obtain a signal. The strategy that we employ takes these two considerationsinto account. Instead of performing a combined fit to the dependence of the nucleon mass onlattice spacing, quark masses and volume and using the result to compute its derivative withrespect to m c , we directly determing the charm-quark content from finite differences (S37,S40)at each lattice spacing. The central ensembles, at each lattice spacing, are tuned to the phys-ical mass point to within less than 4 % in the light and strange quark masses leading to tinycorrections covered by the variation of finite-difference forms (S37,S40). Furthermore, sincethe ensembles at one lattice spacing share all parameters except for the charm-quark mass, and M π L > , finite-volume effects are irrelevant too. We compute the mass differences ∆ ± of (S36) directly from the ratio of nucleon correlatorsfrom the two relevant ensembles. Staggered excited states in the two nucleon correlators aresuppressed by using the time-shifted propagator described in Sec. 2.4, for each of the ensem-bles separately. Fig. S8 shows an effective mass plateau for ∆ M = ∆ + + ∆ − , for our β = 3 . ensembles. We identify the onset of the plateau to be slightly below t min ∼ . and corre-spondingly consider two values, t min (cid:39) . and t min (cid:39) .
95 fm , for determining the massdifference. The variation of the result with respect to t min enters into the systematic error.We determine the lattice spacing either by interpolating the pseudoscalar decay constant f π as describes in Sec. 4.3 or by directly using the nucleon mass from the central, physicalensemble at each β . The difference between these two procedures also enters into the systematicerror estimate.In order to estimate cutoff uncertainties, we perform three different continuum extrapola-tions of f Nc . Using values from all three lattice spacings, we perform either a constant or linearextrapolation in a . In addition, we also perform a constant extrapolation using the two finestlattice spacings only. The result from two of these extrapolations is plotted in fig. S9. Thespread between the results of these methods again enters the systematic error. We perform a total of 2( ˆ f Nc , ˜ f Nc ) × × × =24 different analysis procedures. The way in which we combine these analyses to give a finalcentral value and systematic error is described in Sec. 6.4.Fig. S10 gives a breakup of the systematic error into its different components. As one cansee, all restrictions of the fit procedure are in agreement with the final result and the main23 ta -0.0200.020.04 D M e ff Figure S8: Effective mass plateau and extracted mass from the more aggressive fit with t min ∼ . at β = 3 . . 24 .08 a [fm ] f c N Figure S9: Continuum extrapolation of the scalar, charm-quark content of the nucleon, f Nc . Theblue curve and point correspond to a linear extrapolation in a using all three lattice spacingswhile the grey curve and black point corresponds to a constant extrapolation to the results fromour two finest lattices. 25 .05 0.06 0.07 0.08 f c N t min =0.8fm Q=0.56t min =0.95fm Q=0.632 β , no a Q=0.703 β , a Q=0.543 β , no a Q=0.55HQ Q=0.49Taylor Q=0.70f π scal Q=0.56nucl scal Q=0.62TOTAL Q=0.60 Figure S10: Variation of f Nc from different restrictions to the full analysis procedure. Fromtop to bottom these are: restriction to a single plateau range, to a single scaling behavior, to asingle definition for f Nc and to a single scale-setting method. The last row shows the final resultincluding all fits. The average fit quality of the analyses ¯ Q is given in each case. The percentageof the total systematic error due to each variation alone is also displayed (percentages add up inquadrature). 26ontribution to the systematic uncertainty originates from varying the continuum fit. As discussed in Sec. 6, heavy-quark effective theory provides us with a crosscheck of the charmsigma term up to a precision of O ((Λ QCD /m c ) ) ∼ . We enter, into the heavy-quark ex-pansion (S35), α s ( m c ) = 0 . originating from m c ( m c ) = 1 . +25 − )GeV ( ) and thenumerically integrated 5-loop beta function ( ) with α s (91 . . ( ). Thisresults in f Nc | HQ = 0 . − . f (S41)where ¯ f denotes the sum of quark contents of the lighter quarks ¯ f = (cid:88) q = ud,s f Nq . (S42)Our light quark results give ¯ f = 0 . , (S43)yielding f Nc | HQ = 0 . (S44)This result is in excellent agreement with the one from the direct lattice determination. The continuum extrapolations have fit qualities in the range Q = 0 . − . , with an averagefit quality ¯ Q = 0 . .As in the case of the light and strange quark contents, we plot the cumulative distributionfunction of f Nc from all analyses with different weight functions (see Fig. S11). The choice ofthe weighting function does not significantly affect our result and we take the AIC weight forproducing our final result.The lattice calculation of f Nc ( , , ) is compatible with our number. So is the determi-nation of ( ), in which systematic errors are not estimated. b , and t σ -terms Having checked the validity of the heavy-quark expansion (S35) in the case of the charm, weuse it to compute f Nb and f Nt . We use as inputs α s ( m b/t ) = 0 . / . originating from m b ( m b ) = 4 . +4 − )GeV and m t ( m t ) = 160 . +4 . − . )GeV ( ) and the numerically integrated5-loop beta function ( ) with α s (91 . . ( ). Then, expression (S35) for f Nb and f Nt becomes: f Nb = 0 . − . f ,f Nt = 0 . − . f (S45)27 .06 0.07 0.08 f cN r e l a t i v e w e i gh t AIC weightQ weightflat weight
Figure S11: The cumulative distribution function of f Nc over all analyses, obtained using threedifferent weight functions. In the top part of the figure, the resulting central values and totalerrors, corresponding to each weight function, are plotted.28here ¯ f N f denotes the sum of scalar, quark contents defined in (S34).We input our lattice results into these equations to obtain the final numbers reported in Ta-ble 1. The statistical error originates from the lattice input while systematic error estimates fromthe heavy-quark expansion (Λ QCD /m b ) = 0 . , from the lattice and from α s are combined inquadrature to give the systematic error. The errors on f Nb and f Nt are both dominated entirelyby lattice errors on the light, strange and charm quark contents used as input quantities. σ -terms with full corre-lations We provide a C routine that computes arbitrary linear combinations of the scalar, quark contentsfor the proton, the neutron and the nucleon, while retaining the full correlation between thosequantities. When called without arguments, it returns the individual quark contents as wellas the Higgs-nucleon coupling, f Nh , and brief instructions on how to obtain any other linearcombination.The code is available as an ancillary file or via download from http://particle.uni-wuppertal.de/hch/lincomb.c References . M. Schumann, J. Phys.
G46 , 103003 (2019). . S. Mihara, Springer Proc. Phys. , 149–156 (2019). . K. G. Wilson, Phys. Rev.
D10 , [,45(1974)], 2445–2459 (1974). . S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth, Phys. Lett.
B195 , 216–222 (1987). . M. A. Shifman, A. I. Vainshtein, V. I. Zakharov, Phys. Lett.
B78 , 443–446 (1978). . Please see Supplementary Online Material. . S. Weinberg, Phys. Rev. Lett. , 616–621 (1966). . J. Gasser, H. Leutwyler, M. E. Sainio, Phys. Lett.
B253 , 252–259 (1991). . M. M. Pavan, I. I. Strakovsky, R. L. Workman, R. A. Arndt, PiN Newslett. , 110–115(2002). . M. Hoferichter, J. Ruiz de Elvira, B. Kubis, U.-G. Meißner, Phys. Rev. Lett. , 092301(2015). . R. Horsley et al. , Phys. Rev.
D85 , 034506 (2012). . S. Durr et al. , Phys. Rev. Lett. , 172001 (Oct. 2015). . G. S. Bali et al. , Phys. Rev.
D93 , 094504 (2016).29 . Y.-B. Yang, A. Alexandru, T. Draper, J. Liang, K.-F. Liu, Phys. Rev.
D94 , 054503 (2016). . A. Abdel-Rehim et al. , Phys. Rev. Lett. , 252001 (2016). . C. Alexandrou et al. , 1909.00485, arXiv: (Sept. 2019). . H. Ohki et al. , Phys. Rev.
D87 , 034509 (2013). . W. Freeman, D. Toussaint, Phys. Rev.
D88 , 054503 (2013). . P. Junnarkar, A. Walker-Loud, Phys. Rev.
D87 , 114510 (2013). . M. Gong et al. , Phys. Rev.
D88 , 014503 (2013). . C. Patrignani et al. , Chin. Phys.
C40 , 100001 (2016). . S. Aoki et al. , Eur. Phys. J.
C77 , 112 (2017). . S. Dürr et al. , Science , 1224–1227 (2008). . S. Borsanyi et al. , Science , 1452–1455 (2015). . R. J. Hill, M. P. Solon, Phys. Rev.
D91 , 043505 (2015). . M. Hoferichter, P. Klos, J. Menéndez, A. Schwenk, Phys. Rev. Lett. , 181803 (2017). . K. Kajantie, M. Laine, K. Rummukainen, M. E. Shaposhnikov, Phys. Rev. Lett. , 2887–2890 (1996). . F. Csikor, Z. Fodor, J. Heitger, Phys. Rev. Lett. , 21–24 (1999). . Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, K. K. Szabo, Nature , 675–678 (2006). . P. Güttinger, Zeitschrift für Physik , 169–184, ISSN : 1434-601X (Mar. 1932). . W. Pauli, in Handbuch der Physik, vol. 24 (Springer, Berlin, 1933), p. 162. . H. Hellmann, Zeitschrift für Physik , 180–190, ISSN : 0044-3328 (Mar. 1933). . R. P. Feynman, Phys. Rev. , 340–343 (1939). . X.-D. Ji, Phys. Rev.
D52 , 271–281 (1995). . “Combined measurements of Higgs boson production and decay using up to 80 fb − ofproton–proton collision data at √ s =
13 TeV collected with the ATLAS experiment”, tech.rep. ATLAS-CONF-2018-031 (CERN, Geneva, July 2018), ( https : / / cds . cern .ch/record/2629412 ). . X.-D. Ji, Phys. Rev. Lett. , 1071–1074 (1995). . Y.-B. Yang et al. , arXiv: (2018). . G. Colangelo, S. Durr, C. Haefeli, Nucl. Phys.
B721 , 136–174 (2005). . S. Borsanyi et al. , Phys. Lett.
B730 , 99–104 (2014). . R. Bellwied et al. , Phys. Rev.
D92 , 114505 (2015). . C. T. H. Davies et al. , Phys.Rev.Lett.104:132003,2010 , 132003 (2010).30 . MILC, MILC code Version 7 , 2013, ( ). . G. Colangelo, A. Fuhrer, S. Lanz, Phys. Rev.
D82 , 034506 (2010). . M. Tanabashi et al. , Phys. Rev.
D98 , 030001 (2018). . A. Bazavov et al. , Phys. Rev.
D90 , 074509 (2014). . S. Durr et al. , Phys. Lett.
B701 , 265–268 (2011). . H. Akaike, IEEE Transactions on Automatic Control , 716 (1974). . C. Alexandrou, C. Kallidonis, Phys. Rev.
D96 , 034511 (2017). . C. Alexandrou et al. , Phys. Rev.
D95 , [erratum: Phys. Rev.D96,no.9,099906(2017)], 114514(2017). . X.-L. Ren, X.-Z. Ling, L.-S. Geng, Phys. Lett.
B783 , 7–12 (2018). . J. M. Alarcon, J. Martin Camalich, J. A. Oller, Phys. Rev.
D85 , 051503 (2012). . Y.-H. Chen, D.-L. Yao, H. Q. Zheng, Phys. Rev.
D87 , 054019 (2013). . M. Hoferichter, J. Ruiz de Elvira, B. Kubis, U.-G. Meißner, Phys. Rept. , 1–88 (2016). . P. A. Baikov, K. G. Chetyrkin, J. H. Kühn, Phys. Rev. Lett.118