Determination of quark masses from $\mathbf{n_f=4}$ lattice QCD and the RI-SMOM intermediate scheme

We determine the charm and strange quark masses in the $\overline{\text{MS}}$ scheme, using $n_f=2+1+1$ lattice QCD calculations with highly improved staggered quarks (HISQ) and the RI-SMOM intermediate scheme to connect the bare lattice quark masses to continuum renormalisation schemes. Our study covers analysis of systematic uncertainties from this method, including nonperturbative artefacts and the impact of the non-zero physical sea quark masses. We find $m_c^{\overline{\text{MS}}}(3 \text{GeV}) = 0.9896(61)$ GeV and $m_s^{\overline{\text{MS}}}(3 \text{GeV}) = 0.08536(85)$ GeV, where the uncertainties are dominated by the tuning of the bare lattice quark masses. These results are consistent with, and of similar accuracy to, those using the current-current correlator approach coupled to high-order continuum QCD perturbation theory, implemented in the same quark formalism and on the same gauge field configurations. This provides a strong test of the consistency of methods for determining the quark masses to high precision from lattice QCD. We also give updated lattice QCD world averages for $c$ and $s$ quark masses.


I. INTRODUCTION
Quark masses are fundamental parameters of the Standard Model which must be connected via theory to experimentally measured quantities. They arise in the Standard Model from interactions with the Higgs field, and precise knowledge of quark masses will be needed to test stringently the Standard Model picture of mass generation [1].
In lattice QCD simulations the bare quark masses of the theory are input parameters, and these are tuned to reproduce a set number of physical observables, typically meson masses (one for each quark mass in the simulation). These parameters are however defined at the cutoff scale of the theory and are non-universal, because they depend on the specific lattice regularisation of QCD used. To be useful, these values must then be converted to a chosen quark mass definition in a continuum regularisation of QCD at a fixed physical scale. The conversion, or mass renormalisation, factor adjusts for the different treatment of ultraviolet modes on the lattice and in the continuum and so in principle can be calculated straightforwardly by a 'matching' calculation in lattice QCD and continuum QCD perturbation theory. Lattice QCD perturbation theory [2] is very hard beyond the first order in the strong coupling constant, α s , and so this method is limited to an accuracy of several percent [3]. Higher accuracy can be achieved by methods that make use of nonperturbative calculations in lattice QCD combined * andrew.lytle@glasgow.ac.uk † christine.davies@glasgow.ac.uk ‡ URL: http://www.physics.gla.ac.uk/HPQCD with continuum QCD perturbation theory and we will compare results from two such methods here. One issue with these methods is the control of infrared nonperturbative artefacts from the lattice QCD calculation that are a source of systematic uncertainty. The conventional continuum scheme to which lattice masses are converted is the MS scheme and we will denote masses in the MS scheme by m. A scale for the mass must also be chosen and we will use 3 GeV. Having a fixed convention for quoting quark masses allows a comparison between different determinations.
One way to make the lattice QCD to continuum QCD quark mass connection is to calculate short distance physical quantities in lattice QCD that are both sensitive to the quark mass and for which continuum QCD perturbation theory (in the MS scheme) has been done to a high order. The appropriate energy scale for α s should also be large. A successful method of this type is the 'current-current correlator method' [4] that uses timemoments of heavyonium correlators, extrapolated to the continuum from lattice QCD and then matched to QCD perturbation theory accurate through O(α 3 s ) [5][6][7][8][9][10]. The advantage of this method (which we will denote the JJ method) is that nonperturbative effects (condensate contributions), that would otherwise obscure the match to perturbation theory, are suppressed by four powers of Λ/(2m h ), where m h is the heavy quark mass [11,12]. The suppression is very effective, to the point where these effects have negligible impact, because: Λ is small at around 0.3 GeV; m h is large (and can be varied to test the contribution) and 4 is a high power. Here (Λ) 4 represents the expected size of the gluon condensate 0|α s G 2 /π|0 constructed from the gluon field-strength tensor.
Uncertainties in the JJ method arise from missing higher orders in QCD perturbation theory, but these can be tested by implementing the perturbation theory at different scales [12]. This method has given 1% accurate results for charm and bottom quark masses in the MS scheme [4,[12][13][14][15][16]. The results for m c and m b can then be leveraged into an accurate result for lighter quark masses, such as the strange quark mass m s . This is done by determining fully nonperturbatively in lattice QCD the ratio of two quark masses, such as m c /m s , using the same quark formalism for both quarks [12,[17][18][19][20][21]. This ratio (in the continuum limit) is independent of the lattice quark formalism or continuum scheme and so also holds for the MS scheme at a fixed scale µ. Combining the value of the ratio m c /m s (which can now be obtained to an accuracy of better than 1% [12,21,22]) with the value for m c then yields a 1% accurate result for m s . Further ratios between strange and up/down quark masses (see, for example, [21]) can be used to cascade this accuracy down to even lighter quarks.
Since the JJ method enables the value of the quark mass in the MS scheme to be obtained for an input tuned lattice quark mass, it is equivalent to (indirectly) determining the mass renormalisation factor, Z MS m (µ), that connects the two masses [12].
Another completely different method for making the connection between lattice and MS masses is to determine ratios of appropriate matrix elements between external quark states of large virtuality, µ 2 , that can be calculated both in lattice QCD and in the MS scheme in continuum QCD perturbation theory [23]. Such calculations must be done in a fixed gauge, usually Landau gauge. The method proceeds by imposing 'momentumsubtraction' renormalisation conditions [24] on matrix elements in the lattice QCD calculation. e.g.
defines Z Γ for operator O Γ = ψΓψ, where p|O Γ |p 0 is the tree-level matrix element and p 1 | and |p 2 are external quark states. The symmetric kinematic configuration specified here (with q = p 1 − p 2 ) corresponds to the symmetric momentum-subtraction or SMOM scheme. The importance of this configuration will be discussed further below. Applying the condition of eq. (1) to a scalar operator (along with a determination of the wavefunction renormalisation factor) gives directly a mass renormalisation factor, Z SMOM m (µ), that converts the lattice quark mass to that in the SMOM scheme. Because the SMOM scheme can be implemented in the continuum it can itself then be matched to the MS scheme using continuum QCD perturbation theory (in the same gauge) [25,26]. Multiplying the lattice bare quark mass by the final Z MS m (µ) = Z MS/SMOM m (µ) × Z SMOM m (µ) gives the required m(µ). This method has been widely applied to operator renormalisation in general and not just the determination of Z m , going under the name of the 'RI-SMOM' (regularisation-independent symmetric momentum-subtraction) scheme [24]. For a review of this and the earlier RI-MOM scheme, see [27].
The RI-SMOM scheme is expected to work in a window in which Λ QCD µ π a . (2) Here the upper limit aµ 1 keeps control of lattice discretisation effects and the lower limit guards against being dominated by potentially large nonperturbative effects [28] that behave as condensates multiplied by inverse powers of µ. Nonperturbative effects were a major issue with the original RI-MOM scheme [23] which set up the kinematics for eq. (1) so that p 2 1 = p 2 2 = −µ 2 , but p 1 = p 2 so that q 2 = 0. This 'exceptional' configuration gave rise to differences, inversely proportional to µ 2 , between renormalisation factors that should be the same from chiral symmetry (such as those of the pseudoscalar and scalar operators). This was coupled in some cases to strong nonperturbative dependence of the renormalisation factors on the quark mass, see for example [29][30][31][32].
In contrast, since none of the momenta are light-like in the RI-SMOM scheme, the operators associated with it can be analysed within the Operator Product Expansion (OPE) and sensitivity to nonperturbative effects is under better control. Those associated with spontaneous chiral symmetry breaking, for example, are more benign, with behaviour as 1/µ 6 following expectations from the OPE [32,33]. The SMOM vertex functions show only small quark mass dependence. An added bonus is that the RI-SMOM to MS matching factors [25,26] for Z m are much closer to unity (through O(α 2 s )) than their RI-MOM counterparts [34,35]. This means that the RI-SMOM mass renormalisation factor can be obtained with smaller systematic uncertainty.
Nonperturbative condensate effects are still present in the RI-SMOM scheme, however, and their effects must be included in any accurate determination of the quark mass. The leading condensate contribution to Z m is chirally-symmetric and is only suppressed by 1/µ 2 . Since the associated condensate is the Landau gauge gluon condensate (also known as the gluon mass condesate) [36], 0|A 2 |0 , which is thought to be O(1)GeV 2 [37,38], this contribution could have a significant effect up to very high values of µ 2 . Such a contribution must be included in the analysis and constrained with results at multiple µ values. Here we provide a thorough analysis of systematic uncertainties in the determination of the quark mass with this method, including that of nonperturbative effects.
Using the RI-SMOM intermediate scheme we are then able to determine values for m c and m s with comparable accuracy, around 1%, to that obtained using the current-current correlator method, and using the same lattice quark formalism (highly-improved staggered quarks (HISQ)). The RI-SMOM approach has completely different systematic uncertainties, however, so that a comparison of results from the two methods is then a strong test of our understanding of systematic uncertainties, because the lattice bare quark masses are tuned to the same values in both cases.
The paper is laid out as follows: Section II describes briefly the RI-SMOM approach and Section III gives some details needed to implement it for staggered quarks; Section IV then gives results for the lattice determination of Z m in the SMOM scheme; Section V uses these results to determine the quark masses in the MS scheme. Finally Section VI compares to earlier values, giving new world averages, and concludes with prospects for future improvements.

II. THE RI-SMOM METHOD
As outlined in Section I the lattice QCD RI-SMOM approach mimics what would be done in continuum QCD in a momentum-subtraction scheme. A key part of the argument is that the calculation should be set up in a way that is regularisation-independent. Thus within the lattice QCD calculation the same answer for the quark mass in the SMOM scheme should be obtained in any quark formalism up to discretisation effects. Then the continuum limit of the lattice result also holds in the equivalent continuum SMOM scheme. The continuum SMOM to MS matching completes the conversion to the MS scheme. Within the lattice QCD calculation we must then also ensure that the tuning of quark masses and the determination of the lattice spacing are done in a regularisation-independent way. This is of course the standard practice when we determine the lattice spacing and tune lattice quark masses using physical quantities (such as hadron masses) calculated at the lattice QCD physical point (i.e. including sea quarks with physical masses) and take the value from experiment. We will return to this point below.
To determine the renormalisation factor for an operator O Γ in this framework we then need to apply renormalisation conditions to the inverse propagator (to obtain a wavefunction renormalisation factor) and to an amputated vertex function containing O Γ .
For free quarks in the continuum the inverse of the quark propagator, S(p), is The wavefunction renormalisation factor, Z q , in this scheme can be defined by [23,24] 1 so that Z q = 1 in the free theory. Vertex functions G Γ of operator O Γ (=ψΓψ) can be calculated between two external, off-shell quark lines and 'amputated' as: The renormalisation condition (eq. (1)) on Λ Γ yields Z Γ /Z q given lattice values for Λ Γ . From this we can determine Z Γ if we have Z q . Here we are interested in the mass renormalisation factor, Z m = 1/Z S obtained from the scalar quark bilinear: Again the tree-level value of Z S is 1. Here | sym indicates that p 1 and p 2 satisfy p 2 1 = p 2 2 = (p 1 − p 2 ) 2 = −µ 2 (the RI-SMOM condition), so that there is a single momentum scale. We will also be interested in the pseudoscalar operator with renormalisation condition 1 12i This method is straightforward to implement in lattice QCD. The inverse propagators and vertex functions are calculated from ensemble averages over a set of gluon fields. Note that this means that eq. (5) gives Λ Γ as the product of three ensemble averages. Z P and Z S in eqs. (6) and (7) are then defined as a ratio of ensemble averages, with uncertainties determined via a bootstrap procedure. In practice, a relatively small number of gluon field configurations are needed for good numerical precision in the renormalisation factors Z S and Z P .
Calculations can readily be done for a range of different masses for the 'valence' quarks for which propagators are calculated. We will use the same quark mass for the two sides of the vertex function but note that only quark-line connected Wick contractions appear in this calculation. It is conventional to define the RI-SMOM renormalisation constants in the limit of zero valence quark mass, and we do that here. One reason for doing this is that of consistency, since the perturbative calculations that match RI-SMOM to MS have been done for massless quarks 1 . This is discussed further below.
In practice a more important issue is that of nonperturbative quark mass dependence associated with condensate contributions. An operator product expansion (OPE) approach to the RI-SMOM scheme (where it can be rigorously applied) shows that there are contributions to the quark propagators and vertex functions used to define Z m that appear as inverse powers of µ 2 multiplied by powers of quark masses, or quark or gluon condensates or combinations of all of these [32,36]. It is important to remember that, because we are dealing with gaugenoninvariant quantities here, gauge-noninvariant condensates can also appear. These nonperturbative contributions are not part of the perturbative mass renormalisation factor, but they cannot be trivially separated from it in a lattice QCD calculation. Although the nonperturbative terms seen in the RI-SMOM scheme are wellbehaved, they are not entirely negligible at the values of µ 2 that we use here, as we will discuss in Section IV. It  [25,26] in the matching from the RI-SMOM scheme to the MS scheme. Results are given for both n f = 3 and n f = 4 with all quark masses set to zero. Note how small these coefficients are. The equivalent of c1 for the earlier RI-MOM scheme is -0.424 and for c2 with n f = 3, -0.769 [34,35].
therefore makes sense to remove them, where possible, by extrapolating in the valence quark mass to zero. This only works, of course, for cases where the effect is proportional to a power of the quark mass (and we will study these in Section IV B). The leading contribution to Z m in terms of inverse powers of µ comes from the Landau gauge gluon condensate with no powers of quark masses multiplying it and so it cannot be removed by extrapolating to zero quark mass. There are also higher order contributions of this form. This means that we have to allow for contributions of this kind in our fit ansatz for Z m and test for them by varying µ. This enables us to remove them from our determination of the MS quark mass and to allow an appropriate uncertainty in our error budget from our incomplete knowledge of these contributions. Note that the sea quark masses are not extrapolated to zero. We use calculations at physical values of the masses of the u, d, s and c quarks in the sea (with m u = m d ) to determine the lattice spacing and tune the valence masses [12]. We also calculate Z m on multiple gluon configurations with different unphysical values of the masses of the sea quarks (for a given bare coupling) to test the dependence on these parameters. As we show in Section IV A dependence of Z m on the sea quark masses is much smaller than that on the valence quark mass and barely visible. Nonperturbative contributions arising from the sea quarks, some of which depend on the sea quark masses, will be present and we have to estimate a systematic error from that effect.
We return now to the issue of the perturbative matching to MS. The renormalisation factor between the RI-SMOM scheme and the MS scheme has been worked out through O(α 2 s ) in continuum QCD perturbation theory in [24][25][26]. Writing this renormalisation factor as  Table I. These are calculated at zero (valence and sea) quark mass.
We must also account for systematic errors in the perturbative matching in the continuum from our RI-SMOM scheme with non-zero sea quark mass to the MS scheme. Sea quarks appear first at O(α 2 s ) in the matching and the largest effect present in our calculation will be for the sea c quark. We estimate the size of this effect in Appendix A. This gives an adjustment to c 2 that we will include when evaluating Z MS/SMOM m in Section V.

III. RI-SMOM WITH STAGGERED QUARKS
There are minor complications on the lattice QCD side if a staggered quark formalism is used, as here, because of the fermion doubling issue. The staggered quark action is derived from a naive transcription of the Dirac action onto the lattice in which a rotation is made to diagonalise the action in spin-space. The spin degree of freedom can then be dropped and the 16 'doublers' or tastes of the naive action become 4 tastes in the staggered action. To reconstruct the 4-taste 4-spin Dirac field then requires combining staggered quark fields, χ, over a 2 4 hypercube [40]. This has implications for the momentum-space quark field that enters into the momentum-subtraction renormalisation formalism. The full lattice Brillouin zone, in lattice units contains, for staggered quarks, both momentum and taste information [41]. To separate them we must work in a reduced Brillouin zone with an additional 4-dimensional label for each subzone. Then with B µ a 4-dimensional vector of 0s and 1s. We use the method for staggered quarks developed in [42], and here simply give an overview of that procedure. For a given momentum (in lattice units) ap in the reduced Brillouin zone, we invert the staggered Dirac operator on 16 momentum sources of the form e ipx with ap = ap + πA, where A is a 4-vector composed of 0s and 1s. Each of the resulting propagators, S(y, p) where y runs over the lattice volume, is Fourier transformed 16 times with momenta −ap + πB, with B a 4-vector of the same type as A. The results are assembled into a propagator S(ap ) ≡ S AB (ap ) = S(ap + πA, −ap + πB). (12) This is a 48 × 48 matrix, but we have kept the colour indices implicit; the matrix is diagonal in colour space on forming the ensemble average over lattice gluon fields. The propagator is also a taste-singlet [42] and so has the same properties for the purposes of the SMOM approach to those for other quark formalisms. After averaging over gluon fields the matrix is inverted for each value of p to obtain the inverse propagator.
To apply the condition in eq. (4) to determine Z q we must multiply by a representation of the matrix / p in AB space. Using the notation of [42] this is the matrix p µ (γ µ ⊗ I) that is the Fourier transform of the (tastesinglet) derivative term in the free inverse propagator. Since this derivative is improved to remove a 2 discretisation effects for our improved staggered quark action, we take ap µ = sin(ap µ ) + sin 3 (ap µ )/6 so that Z q is equal to 1 in the free case. (γ µ ⊗ I) is a matrix of 0, 1 and -1 obtained by tracing over products of gamma matrices as described in Appendix A of [42]. Then The trace is over spin, taste and colour.
The scalar operator that we use to determine the mass renormalisation factor is the local taste-singlet operator χ(x)χ(x). The vertex function for this operator is then constructed as Here (−1) x is the alternating phase factor over the lattice, (−1) x1+x2+x3+x4 . S † is the hermitian conjugate in colour space and has a permuted B index according tõ B = B + 2 (1, 1, 1, 1). To apply eq. (6) we must multiply G S,AB on both sides by the inverse propagator to give Λ S,AB and again take the trace over spin, taste and colour. For the local taste-singlet scalar this gives the simple expression For the local pseudoscalar operator the procedure is identical except that there is no (−1) x in the equivalent of eq. (14) and in the equivalent of eq. (15) multiplication by the matrix γ 5 ⊗ γ 5 is needed before taking the trace. This can be written simply as a 16 × 16 matrix with a skew-diagonal of 1s. 1/Z S = Z m is then obtained by dividing by Z q .

IV. LATTICE QCD CALCULATION
For this calculation we use ensembles of gluon field configurations generated by the MILC collaboration [43,44]. These include u, d, s and c quarks in the quark sea, with m u = m d = m l . The gluon action is fully improved to remove discretisation errors through O(α s a 2 ) [45]. The sea quarks are implemented through the Highly Improved Staggered Quark (HISQ) formalism [46,47] designed, and demonstrated, to have very small discretisation effects, at α 2 s a 2 and a 4 . We also use the HISQ formalism for our propagator and vertex function calculations. The simulation parameters for the sets (ensembles) of gluon field configurations used are given in Table II. We have sets at three different values of the bare QCD coupling, β, with finer lattice spacing as β increases. For β = 6.0, give the sea quark masses in lattice units. Sets 1-9 will be referred to in the text as 'coarse', sets 10 and 11 as 'fine' and set 12 as 'superfine'. Most of the sets we show here are used to test dependence on sea quark masses or spatial volume; our final analysis will be done using results from sets 2, 4, 9, 10, 11 and 12 (marked in bold in the  We fix the gauge field configurations to lattice Landau gauge by maximising the average trace over colour of the gluon field link. Note that this differs from the continuum Landau gauge by discretisation errors [48]. On each ensemble we then calculate quark propagators for a range of quark masses and momentum values and assemble vertex functions as described in Section II. We use the bootstrap method to determine the uncertainty in Z SMOM m from combining Z S and Z q , as well as the correlations between results obtained on a given ensemble. High precision is possible with these calculations with only a moderate number of samples of gluon field configurations. We use 20 from each set, well-spaced in Monte Carlo generation time for statistical errors below 0.1%. We have tested that the statistical errors are Gaussian by comparing the mean and median from a bootstrap distribution. We have also checked that the tolerance we use for the Dirac matrix inversion is such that tightening the tolerance has no significant effect on the results. The impact of the gauge-fixing tolerance will be discussed below.
For the RI-SMOM calculations reported here, the momenta that we use in the two propagators combined in the vertex function are given in lattice units by for integer x. Adding the additional θ/2 term through a 'momentum-twist' using phased boundary conditions [49,50] allows us to tune the value of the momenta used precisely. This means, for example, that we can tune the momenta to be the same on ensembles with different values of L s . With the definitions of eq. (16) a(p 1 − p 2 ) has the same magnitude as each of ap 1 and ap 2 which is the appropriate kinematics for the RI-SMOM scheme. We will call this magnitude aµ: We use momenta only in the spatial directions for simplicity because our lattices have a different extent in the time direction. We use a variety of x and θ values and have tested that results do not change under a change of x and θ to achieve the same ap. Note that, in keeping with eq. (10), we do not want any momentum component in lattice units to exceed π/2. This limits how high a value of µ we can reach; for example we cannot exceed a µ value of 3 GeV on the coarse lattices.
The results enable us to extract a renormalisation factor for the scalar current, Z SMOM m , in the RI-SMOM scheme for each ensemble for a variety of aµ values and HISQ quark masses used in the propagators, am. Insofar as Z m is a renormalisation factor from one QCD regularisation scheme to another, taking account of the differences in the two schemes at the cutoff, we expect Z m to behave as a power series in α s with coefficients that depend logarithmically on the ratio of the two cutoffs, i.e. on ln(aµ). Since Z m here is being determined nonperturbatively in lattice QCD, differences from this expectation arise both for small aµ and large aµ values and we will address both of these here.
For large aµ values, systematic discretisation effects can appear from the granularity of the lattice. In Z m such effects would cause systematic errors of the form (aµ) n where n is a positive power whose value depends on the quark action, a higher power corresponding to a more highly improved action. With the HISQ action we have removed tree-level a 2 terms and so we expect discretisation effects at (aµ) 2 to be suppressed by powers of α s and therefore to be relatively small [46]. The lowest order at which tree-level discretisation errors can appear is at (aµ) 4 . In fact the discretisation errors, as long as they are not too large, are benign. In the end, in order to determine a quark mass relevant to the physical world, we will perform an extrapolation in a to the continuum limit a = 0, at fixed µ, and remove discretisation errors.
Of more concern are nonperturbative effects that can have an impact at small values of µ. An operator product expansion (OPE) tells us that the vertex functions can be expanded in inverse powers of µ with coefficients that depend on condensates, vacuum expectation values of local quark and gluon operators. In the current-current correlator method, this effect was studied in [12]. There the heavy quark mass, m h , replaces µ and nonperturbative terms of the form can appear in the correlator moments. The first term contains the light quark chiral condensate ψψ and the second, the gluon condensate constructed from the gluon field-strength tensor (the heavy quark condensate being absorbed into this). Since the current-current correlator propagators and so gauge-noninvariant condensates can appear. Such condensates can be larger in magnitude than the gauge-invariant ones because powers of the Landau gauge gluon field, A µ , can appear (see, for example, [37] for the gluon case) and this can be associated with inverse powers of µ as small as µ −2 . In Section V and Appendix B we discuss how we expect such a condensate to affect Z m .
Study of the impact of the gauge-fixing tolerance provides some evidence of sensitivity to these gaugenoninvariant nonperturbative effects. For the Landau gauge-fixing, we use a tolerance of 10 −7 on the magnitude of the gradient of the gauge field on gluon configuration sets 1 through 11. This fixes the average trace of the link in Landau gauge to a few parts in 10,000. The residual effect on Z m from this gauge-fixing tolerance is at the same level as we now demonstrate. Figure 1 shows a scatter plot of bootstrap samples for Z SMOM m on coarse set 4 at µ = 2 GeV for a gauge-fixing tolerance of 10 −7 and then successive tightening of this tolerance by factors of 10 down to 10 −10 . The tighter tolerance gives a shift in the mean value of Z SMOM m by around 0.0004. Results at higher µ values show much smaller effects in a way that demonstrates their origin in nonperturbative effects. This is illustrated in Figure 2 which shows the change in Z SMOM m for a factor of 10 change in gauge-fixing tolerance as a function of µ. To cover the residual gauge-fixing effects we take an additional uncorrelated uncertainty on our Z SMOM m results of 0.0004 for µ = 2 GeV, 0.0001 for µ = 3 GeV and 0.00002 at µ = 4 GeV on sets 1 through 11. This is typically at the level of our statistical uncertainties. On set 12 we fixed to Landau gauge with a tolerance of 10 −14 and do not take any additional uncertainty from residual gauge-fixing effects. We do not consider any possible effects from Gribov copies (for studies of this in the RI-MOM scheme see, for example, [51,52]).
Another way to assess the size of nonperturbative effects is to look at differences of Z factors for operators which should have the same perturbative expansion. Since we are concentrating here on Z S it makes sense to look at the difference between Z S and Z P . This difference showed significant problems with Z P in the RI-MOM scheme because it exposed nonperturbative contributions that behaved as 1/µ 2 [31]. This behaviour can be traced to the fact that the inserted operator is carrying no momentum (q = 0) in that scheme [23]. This causes particular problems for the pseudoscalar operator and was deemed to make this vertex function of only very limited use. A related issue arises with the scalar operator in the RI-MOM scheme, however, and that is one of very strong dependence on the quark mass. These features were illustrated for the HISQ action in [53] where the vertex functions Λ P and Λ S are compared for the RI-MOM and RI-SMOM schemes as a function of momentum and quark mass, and the superior behaviour of the RI-SMOM scheme is very clear.
In the RI-SMOM scheme, as Figure 3 shows, the nonperturbative behaviour of Tr(Λ P − Λ S )/48 (proportional to (1/Z P − 1/Z S )) is quite benign, falling as 1/µ 6 with little dependence on the lattice spacing. This indicates that either Z P or Z S could be used to determine the mass renormalisation factor; we will however concentrate here on Z S . As we will see in Section IV B, the quark mass dependence of Z m derived from Z S in the RI-SMOM scheme is also much less of an issue than it was in the RI-MOM scheme.
The slope with 1/µ 6 of Tr(Λ P − Λ S )/48 in Figure 3 is O(1) GeV 6 when translated into its effect on Z P − Z S , in approximate agreement with what is seen with domainwall quarks [33]. This sets an appropriate scale to use when assessing systematic effects from nonperturbative contributions in Section V. These systematic effects will show up when results evaluated at different µ are run perturbatively in the continuum to a common scale. We will use multiple values of µ (2, 2.5, 3, 4 and 5 GeV) in our analysis and compare results for the MS mass at a reference scale of 3 GeV. Z SMOM m is dimensionless but the appropriate scale for it, µ, must be obtained in GeV by multiplying aµ by the inverse of the lattice spacing. The value of the lattice spacing is obtained from determining a dimensionful quantity that can be matched in the continuum limit at physical quark masses to an experimental value. We use the Wilson flow parameter, w 0 [54], itself fixed at the value 0.1715(9) fm using the decay constant of the π [55].
The physical quark mass limit can be approached in a number of different ways. When calculating quantities such as hadron masses, which are sensitive to low momentum scales, it is convenient to keep the bare coupling constant, α lat = g 2 /(4π), and w 0 fixed. This means that the value of a varies slightly with the sea quark masses but the variation of hadron masses is small, since they behave in a similar way to w 0 . An alternative, which is more suitable for the determination of bare parameters such as quark masses for reasons discussed in Appendix A of [12], is to fix α lat and the lattice spacing. This latter method is the one that we will implement here.
We use results from [12] where the sea quark mass dependence of w 0 /a in terms of the result at the tuned physical point is determined for different α lat . A universal linear dependence on δm sea uds (the deviation of the sum of the u/d and s sea masses from their physical values) is seen, for values of δm sea uds in units of the tuned s quark mass less than 0.5 (see Figure 10 in Appendix A). An analysis of dependence on the sea c quark mass away from the physical point is also given. We use the fits of [12] to interpolate results for w 0 /a for sets of ensembles at a given value of β to the physical quark mass point. The values that we obtain for w 0 /a in this way are given for 'coarse' (β = 6.0), 'fine' (β = 6.30) and 'superfine' (β = 6.72) sets in Table III. We will use the w 0 /a result (and the value it implies for a −1 in GeV) for all the ensembles with that value of β.
As we will see in Section IV A this approach means that Z SMOM m has little discernible sea quark mass dependence. This is expected insofar as Z SMOM m represents physics at the cut-off scale that is a function of α lat , with the light sea quark masses having only a very small effect on its perturbative expansion. We will test the impact of the c sea mass in the next section and in Appendix A.
We take a similar approach for the tuned bare quark masses for s and c, as will be discussed further in section V.
A. Sea mass dependence Table II shows the variety of ensembles on which we have calculated Z SMOM m . Note particularly how many different combinations of sea quark masses we have studied for β = 6.0 ('coarse'). This allows us to test for dependence on the sea quark masses, given the method described in Section IV for fixing the lattice spacing. No significant sea quark mass dependence is seen for the µ values that we will use for our analysis. Some of the  [54], w0, and tuned quark masses for the coarse, fine and superfine sets of ensembles as determined in [12]. These are obtained by fitting the sea quark mass dependence of these parameters and interpolating/extrapolating to physical sea quark masses. The lattice spacing is obtained from w0/a by using the value for w0 of 0.1715(9) fm determined from the pion decay constant in [55]. For the quark masses the uncertainty is split into two pieces. The first uncertainty is uncorrelated between lattice spacing values and comes from statistical/fitting errors and uncertainties in the value of w0/a. The second uncertainty is correlated between lattice spacing values because it comes from the uncertainty in w0 and from the uncertainty in the ηc or ηs meson mass as appropriate. (2)(7) ensembles have very different combinations of sea quark masses from those that would be considered suitable for a comparison to the real world. For example, set 6 has u/d and s quark masses equal at a value around 1/10th that of the physical s mass. Nevertheless even this ensemble has a Z SMOM m that agrees (for a given µ) with that from set 4 where m s is ten times larger and therefore more realistic. Note that the components of Z m , i.e. Z q and the vertex function, typically show somewhat larger changes with sea mass but the effects cancel in Z m .  Figure 4 corresponds to the u/d quark mass in the sea. This figure therefore also indicates how little variation we can expect as the u/d quark mass in the sea is varied. In our final analysis we will include results from different values of the sea quark masses to allow for small variations to be taken into account, but these results indicate that such variations are below the level of our statistical uncertainties.
A similar picture is seen for variation with the c sea mass, even though the c sea mass is much larger and 2m c is close to µ in our range of µ values, so that one could worry that an effect might be discernible. We can gauge the possible size of such an effect from the perturbative analysis of the impact of massive c quarks in the sea given in Appendix A. That shows a shift of size 0.1% for a change in m c from zero to m c at µ = 2 GeV. Figure 5 shows a comparison of results for Z SMOM m as a function of µ for sets 1 and 2 which have a substantially (30%) different c sea quark mass, along with slightly different u/d and s sea quark masses (which Figure 4 has already demonstrated have no effect) and different spatial volumes (again for which we see no effect in Section IV C). We also include results for set 3 which has a value of am c differing from that on set 2 by 1%, a size of varia-  tion which is closer to that of typical m c mistuning on our ensembles [12]. The lower plot of Figure 5 gives more detail at µ = 2 GeV and shows, as expected, no variation of Z SMOM m at the level of 0.1% for a change in am c of 30%. It also shows that there is no impact on our results at the level of our statistical errors from the slight (5%) mistuning of the c sea mass that we have on some ensembles.
Finally, in Figure 6 we compare results for two fine lattices with different sea mass values (sets 10 and 11). This plot covers three µ values we will use in our final analysis, 2, 3 and 4 GeV. Good agreement is seen between the results on sets 10 and 11 (with the largest discrepancy being 0.1% for µ = 2 GeV), testing sea-mass dependence as well as dependence on the spatial volume, to be discussed in more detail in Section IV C.   Table II). These sets have slightly different u/d and s sea quark masses, but by an amount which is much smaller than the change in the c sea mass from set 1 to sets 2 and

B. Valence mass dependence & extrapolation
In Section IV A the Z SMOM m renormalisation factors are determined for small and fixed but non-zero valence quark masses and we showed that the dependence on sea quark masses is almost negligible. Here we will show that there is a small but visible dependence for Z SMOM m on the valence quark mass. This dependence comes from the vertex function since the wave function renormalisation is almost independent of the valence quark mass. Since the impact in perturbation theory of the small valence quark masses we use should be negligible, the most likely source of valence quark mass dependence is nonperturbative, i.e. that of quark masses multiplying a condensate appearing in conjunction with inverse powers of µ as expected from the OPE.  Table II). These sets also have slightly different s and c sea quark masses and spatial volume. The plot shows results for Z SMOM m as a function of µ in GeV for the two sets, for a valence mass in lattice units of 0.0074. lence masses; that of the u/d quark mass in the sea and two and three times that value. Figure 7 shows very little visible dependence on ma but it is, however, significant (see the lower plot of Figure 7 for more detail in the superfine case). Note that the results at different values of am are correlated and we include this correlation in our fits through a covariance matrix determined by the bootstrap procedure.
As discussed in Section IV it is convenient to extrapolate in the valence quark mass to zero, so that we can convert from the SMOM scheme to the MS scheme using perturbation theory at zero valence quark mass. The extrapolation also has the advantage of removing some of the non-perturbative condensate contributions that are not part of Z SMOM m . The size of the observed mass dependence then provides an indication of the size of the remaining condensate contributions that do not depend on the quark mass.
To extrapolate to zero valence quark mass we fit Z SMOM m to a simple polynomial form in am val given below. Including only a linear term does not give a good fit over the range of am values that we use when correlations are included. We therefore add in both a quadratic and cubic term: We use a prior on Z SMOM m (the value of Z SMOM m (am) in the massless limit) of 1.0±0.5. For the coefficients d i priors of {0 ± 0.1, 0 ± 0.01, 0 ± 0.001} are used for µ = 2 GeV with the values being decreased by a factor of 2 and 4, respectively for µ = 3 GeV and 4 GeV to allow for an approximate µ −2 suppression, the smallest inverse power of µ that we expect to appear. We divide the lattice valence masses by the tuned s quark mass at that lattice spacing so that the d i are dimensionless and physical. Note that, if the linear term were set by the gauge-invariant condensate m ψψ /µ 4 , we would expect d 1 to be O(2 × 10 −4 ), which would make the mass dependence scarcely visible. Instead our priors allow for possibly larger gaugenoninvariant condensates to appear. The linear slope shows significant lattice spacing dependence and is consistent with zero on the superfine lattices, as shown in Figure 8. The coefficient of the quadratic mass dependence, d 2 , is non-zero and is shown for the superfine lattices in Figure 8 along with the results of a simple fit to the form C/µ 4 (C/µ 2 does not give a good fit, although C/µ 6 is also acceptable) with C = -0.10(3) GeV 4 , equivalent to −(0.56(4) GeV) 4 . From this (and earlier results in Section IV) we conclude that in our fits in Section V, comparing quark masses determined using Z SMOM m values at different values of µ, we should allow for condensate contributions remaining in Z SMOM m that could be as large as (1GeV) n /µ n coming from gauge-noninvariant condensates. This will allow us to include an uncertainty from such nonperturbative contributions in our determination of the mass. Since Z SMOM m is a matching factor between two different regularisations of QCD we expect it to be dominated by ultraviolet physics and not to be sensitive to the volume of the lattice. However, we have observed some infrared sensitivity in the form of nonperturbative condensate contributions. In aiming for a precise determination of Z SMOM m finite-volume effects need to be tested. This is straightforward to do on lattices that have the same β and sea quark masses and differ only by the number of lattice points in each spatial direction. Figure 9 shows such results for sets 3, 4 and 5 that have 24, 32 and 40 lattice points in each spatial direction but exactly the same parameters in the lattice QCD Lagrangian (see Table II). No significant dependence on the lattice size is seen except for very small µ (below 2.0 GeV, which is our smallest value for analysis) and for small lattices (of size L s =24 which is smaller in terms of M π L s than any of the lattices that we use for analysis). We conclude from this that finite volume effects are negligible.

V. DETERMINATION OF MASSES IN THE MS SCHEME
Our procedure here for determining the quark mass in the MS scheme has three ingredients: • A quark mass in our lattice QCD scheme tuned nonperturbatively to reproduce the mass of a given hadron from experiment; • A nonperturbative determination from lattice QCD of the mass renormalisation factor Z SMOM m that converts this mass at each value of the lattice spacing to a mass in our SMOM scheme at a given value of the scale, µ; • A perturbative calculation of the mass renormalisation factor Z MS/SMOM m (through α 2 s ) that further converts the SMOM mass at scale µ to the mass in the MS scheme at scale µ. From there we can run the mass to different scales using 4-loop running in the MS scheme [56,57]. Here m(a) is the bare lattice quark mass, in physical units, at a specific value of the lattice spacing, the first item from the list above. Z SMOM m is the mass renormalisation factor calculated nonperturbatively on lattice QCD configurations at a specific lattice spacing, allowing us to convert the lattice mass to the SMOM scheme at a scale µ. How this is calculated has been discussed in earlier sections; here we will give the results. The intermediate quark mass we obtain in the SMOM scheme, although nominally now in a continuum scheme at a physical scale, will still carry remnants of its lattice origins through discretisation errors. These must be removed by calculation at multiple values of the lattice spacing, so that an extrapolation to the continuum limit, a = 0, can be made. This could be done with the SMOM masses, but we choose to first convert to the MS scheme at scale µ by multiplying by the final factor Z MS/SMOM m . We denote the MS mass obtained this way as m(µ, a) to show that it has yet to be extrapolated to the continuum limit. for µ = 2, 2.5, 3, 4 and 5 GeV on a subset of ensembles from Table II covering coarse to superfine lattice spacings. The values are obtained by extrapolating to zero valence quark mass as discussed in the text. For each set of results we also give, in column 5, the correlation matrix between the values for different µ. Note that there is very slight mistuning of some µ values on the coarse and fine lattices (sets 2 -11) and the actual µ values are given in the column headings. Statistical errors only are given here. We include a further uncorrelated uncertainty on Zm values for sets 2 -11 as described in Section IV to account for residual gauge-fixing effects (±0.0004 at µ = 2 GeV, ±0.0001 at 3 GeV and ±0.00002 at 4 GeV). We will describe how we do this below; first we give the results that we will use for each of the ingredients of eq. (20).

A. m(a)
The bare lattice quark masses that we use are for s and c quarks and are given in [12] for the ensembles and lattice spacing values that we use here. The c quark mass, m c (a), was tuned by adjusting the lattice mass to give the physical value for the η c meson mass. The physical value for the η c mass was defined from the experimental value with a shift upwards by 2.7 MeV (less than 0.1%) to remove electromagnetic effects and to account for cc annihilation, since both of these effects are missing in our lattice QCD calculation [58,59]. The uncertainty on the physical η c mass to which we match is then increased (to 2.7 MeV) to allow for uncertainty in these corrections. The s quark mass is similarly tuned based on the physical mass of the ss pseudoscalar particle known as the η s . It is an unphysical particle since its valence quarks are artificially not allowed to annihilate, but its properties can be well determined in lattice QCD [55,60]. Its mass can be determined in terms of K and π meson masses as [55] M phys ηs = 0.6885 (22) GeV.
where the uncertainty includes a systematic error from the neglect of electromagnetism in the lattice QCD calculations.
The tuned lattice bare c and s quark masses are given in GeV in Table II of [12]. These are the values that give the physical η c or η s mass on each ensemble given a fixed value for w 0 . Since here (as explained in Section IV) we are approaching the physical mass point using a fixed lattice spacing value (since that removes sea quark mass dependence from our results, as shown in Section IV A) then we also need fixed tuned quark mass for sets of ensembles at a fixed value of β. The fits to the dependence on sea quark mass discussed in Appendix A of [12] enables us to determine the tuned c and s quark masses for physical sea quark masses at each value of β. These are the values that we will use here and they are given in Table III. The uncertainties in the tuned masses include the uncertainty from the lattice spacing. This gives a 3× smaller relative uncertainty for the c quark mass than for the s quark mass because the lattice spacing uncertainty appears with the 'binding energy' of the meson rather than its mass. For the η c the binding energy is much smaller than the mass, but for the η s it is of the same size. Table III divides the uncertainty on the tuned masses into two components: a correlated uncertainty from the value of w 0 and the value of the meson masses  Tables I and VII. The uncertainty in Z quoted here comes from the uncertainty in the value of αs (and so is 100% correlated between the values). We use α MS (n f = 4, 5.0 GeV) = 0.2128(25) [12]. Column 3 gives the multiplicative factor R(3 GeV, µ) that converts the MS mass at scale µ to the mass at the reference scale of 3 GeV. This is obtained from 4-loop running in perturbative QCD. The uncertainty is dominated by that in αs; the uncertainty from missing higher order terms in the running is negligible in comparison. The error in R is then 100% correlated or anti-correlated between the values, depending on whether R is greater than or less than 1. There is also a 100% correlation (or anticorrelation, as appropriate) with the errors in Z MS/SMOM m . Note that the values in the table are for the µ values in column 1. We calculate R and Z inside our fit function and hence allow for the fact there that the µ values are slightly mistuned on coarse and fine lattices (see Table IV). We include an additional uncertainty, as described in the text, to allow for w0/a errors feeding in to the determination of µ. This gives a (correlated) uncertainty of 0.0003 on the coarse lattices, 0.0002 on the fine lattices and 0.0008 on the superfine lattices.   The point plotted as a white circle, further offset to the left, is from the current-current correlator method [12] (along with a nonperturbative determination of mc/ms in the ms case).
Working from right to left in eq. (20) the next set of results that we need are for Z SMOM m for each ensemble that we will use in determining our continuum and chiral limit for the quark masses. We have chosen to work with multiple values of µ in order to assess the impact of nonperturbative terms on the mass renormalisation. These are µ = 2, 2.5, 3, 4 and 5 GeV. At each value of µ we determine Z SMOM m at three values of the valence quark mass, as described in Section IV B, and extrapolate to zero valence mass. Results are given in Table IV for the ensembles that we will use. The results for differ-ent µ values on a given ensemble are correlated and so we include in the Table the correlation matrix for the Z values. The correlation matrix, ρ ij , for variable x i and x j is defined by with indicating the expectation value, and σ the standard deviation. Results are given for sets 2, 4, 9, 10, 11 and 12 that we will use in our final analysis which will determine a continuum limit for the mass and allow for small residual sea quark mass effects. Note the very slight mistunings of µ from the nominal values on sets 2-11. These are allowed for in our fits.  Table VII. We use the results of Appendix A to adjust c 2 to allow for having a massive c quark in the sea. This has a very small effect for µ = 2 GeV, and even smaller one for µ = 2.5 GeV and is otherwise negligible. The resulting values for Z MS/SMOM m are given in Table V. The uncertainty in the Z values quoted there comes from the uncertainty in α s and so is 100% correlated between the values. There is also a systematic uncertainty from missing higher orders in the perturbative expansion and we will allow for that in our final fits below. By multiplying all three ingredients together as in eq. (20) we obtain values for the c or s quark mass in the MS scheme at a (nominal) scale of µ = 2, 2.5, 3, 4 or 5 GeV from each lattice ensemble. These results still contain discretisation effects from the lattice QCD component of the calculation. To remove these effects we must extrapolate to the continuum limit. At the same time we want to assess other systematic effects such as the nonperturbative contributions to the lattice QCD determination of Z SMOM m that have not been removed by our extrapolation to zero valence quark mass, and missing higher order perturbative contributions to Z MS/SMOM m . This can be done by comparing results at different µ but, the simplest way to pick out these effects is to run all the results to a common scale, µ ref . We take µ ref = 3 GeV, so that we run up from 2 and 2.5 GeV and down from 4 and 5 GeV. The running is done by integrating the evolution equations numerically in the MS scheme using 4-loop expressions [56,57,61,62]. The result of this is a multiplicative factor R(µ ref , µ) such that Values of R(3 GeV, µ) are given in Table V. Note that, because the uncertainty comes from the uncertainty in α s the uncertainty is 100% (anti-)correlated between the different µ values. The uncertainties in R and Z MS/SMOM m are also 100% (anti-)correlated for the same reason. At this point we also need to include an uncertainty coming from the relative determination of the lattice spacing on coarse, fine and superfine lattices. The uncertainty in w 0 /a, given in Table III, means that the µ values on each set may not match and this gives an additional uncertainty in the running of the mass to the 3 GeV reference point, including in the values obtained at µ = 3 GeV. This gives an additional (correlated) uncertainty of 0.0003 on the coarse lattices, 0.0002 on fine lattices and 0.0008 on superfine lattices. There is an additional correlated 0.1% uncertainty on all points coming from the effect on µ of the uncertainty in the value of w 0 .
We then have results for m s (µ ref ) and m c (µ ref ) that come from lattice calculations with different values of the lattice spacing and different values of µ. We fit these to a function that allows for discretisation effects that depend on a and other systematic effects that depend on µ. It is important to include the correlations between the points: our results at different values of a are correlated through their dependence on the value of w 0 which is used to determine the lattice spacing and our results at different µ for a given ensemble are correlated through the statistical uncertainties in the values of Z MS/SMOM m (see Table IV).
Our results are plotted in Figure 10. Discretisation effects are clearly evident with the slope in a 2 becoming larger with larger µ, not surprisingly. Results at different µ come together on the finer lattices.
A key point, as we have emphasised, is to provide constraints on nonperturbative µ dependence (from condensate terms) that would survive the continuum limit from our lattice QCD calculation but is not part of Z m . To understand how big these terms might be, we turn to the OPE (for more details, see Appendix B). The analysis of the quark propagator is given in [36,63,64]. To lowest order in inverse powers of the momentum (p ≡ µ) and α s this gives (rather than eq. (4)) Here Z q is the perturbative contribution from the leading (unit) operator, A 2 is the square of the Landau gauge gluon field and X denotes vacuum expectation values of dimension 4 operators such as mψψ (which vanishes at zero quark mass), ψAψ and G 2 . The coefficients of the higher dimension operators in the OPE are obtained by matching scattering amplitudes for both sides of the OPE from, for example, low-momentum gluons. Repeating this procedure for the scalar vertex function for the symmetric kinematic point that we use in this calculation (and which allows an OPE treatment), yields 1 12 Tr rather than eq. (6). From eqs (24) and (25) we see that the leading nonperturbative contribution to our determination of Z m is −(2π/3)α s (µ) A 2 /µ 2 . The value of A 2 is not well-known [37,38] and so we simply allow for it to be of size O(1GeV) 2 in our fits to obtain the continuum limit of the MS quark mass. We must also allow for the higher dimension condensates denoted X above, about which even less is known. To do this we include terms at 1/µ 4 and 1/µ 6 in our fits, again allowing the operator vacuum expectation values (summed over all the operators that could appear) to be O(1GeV) 2n .
In fitting our results to obtain physical values for the quark masses we must then allow for: lattice spacing artefacts, non-perturbative effects, sea quark mass effects and missing terms in the perturbative matching to MS. To allow for all of these, we fit our results to the following form: cond,a 2 (Λa/π) 2 .
Here m(µ ref ) is the physical result. The coefficients c µ 2 a 2 allow for discretisation effects set by the scale µ (coming from Z m ) and the c Λ 2 a 2 allow for those set by the scale Λ in the tuning of the bare quark masses, independent of µ. We take Λ to be 500 MeV in the case of the s quark mass, but 1 GeV in the c quark mass case, since it could be set by m c itself. We take the prior on all the c a 2 coefficients to be 0.0 ± 1.0. c α allows for systematic uncertainties in the continuum limit from missing α 3 s terms in the matching of SMOM to MS. We take the prior on c α to be 0.0±0.2, allowing for a size four times larger than c 1 or c 2 . The coefficients c cond allow for nonperturbative condensate effects that have not been removed by extrapolating the valence quark masses to zero. We include three such terms with inverse powers of µ of 2, 4 and 6 since we expect these to be the most significant. The results that we give in Sections IV and IV B show that we need to allow for gauge-noninvariant condensates of size as large as (1 GeV) 2n . We take this to be the generic size of the condensate and give each one a coefficient with prior 0±2 (consistent with the combination of eq. (24) and (25)).
We also allow for a-dependence in each condensate with Λ = 500 MeV and a prior on c cond,a 2 of 0.0 ± 1.0. All h sea and k sea coefficients allow for any small remaining dependence on the sea quark masses, either explicitly in condensate terms or elsewhere with We take the priors on these coefficients to be 0.0 ± 0.2 (consistent with results in Section IV A). The fit is strongly constrained by the number of different correlations included between results from different µ values and different a values. We obtain a χ 2 /dof of 0.8 for both the fits for m c and for m s . We can also do both fits simultaneously, requiring all coefficients to be the same except those for the (aΛ) n terms and then we obtain a χ 2 /dof of 0.7. If we drop the condensate terms from the separate fits the χ 2 /dof increases to 2, indicating that these are important. The χ 2 /dof for the simultaneous fit without condensates increases to 6. We find the total condensate contribution at a = 0 and δ sea = 0 to be relatively small, at -1.0(5)% at µ = 2 GeV and -0.3(1)% at µ = 5 GeV for the separate fits. The simultaneous fit has somewhat more significance, at -1.4(4)% for the condensate contribution at µ = 2 GeV. The condensate contribution is a combination of a negative term at 1/µ 2 (as expected from above), a positive term at 1/µ 4 and a relatively unconstrained term at 1/µ 6 . Figure 11 demonstrates the robustness of our fit, by showing the impact on the final value of numerous modifications. These include missing out sets of results; doubling prior widths on various fit coefficients and changing the numbers of terms used to describe discretisation effects, condensate contributions or missing pieces of the perturbative matching. Effects are relatively minor and generally well within our uncertainties.
The values we obtain for the physical result, m(µ ref ) at the reference scale of 3 GeV (using separate fits to each mass) are: The error budget for the two numbers, evaluated from the fit, is given in Table VI. The uncertainties are dominated by those from the tuned bare quark masses (especially for m s ) but with sizeable contributions from the continuum extrapolation, possible missing α 3 s terms in the SMOM to MS matching and condensate effects.

VI. CONCLUSIONS
Lattice QCD is the method of choice for determining quark masses because it gives direct access to those parameters in the QCD Lagrangian and allows them to be tuned very cleanly against hadron masses measured in experiment. As emphasised in Section I, the key complication in determining quark masses is in providing the matching factor from the quark mass in a particular lattice QCD regularisation scheme to the preferred MS continuum regularisation scheme. Here we compare two accurate methods for providing this matching factor directly for the c quark mass: one is to take the continuum limit of time-moments in the current-current correlator method and the second is to use an intermediate momentum-subtraction scheme whose definition on the lattice translates directly to the continuum. Both methods then use continuum QCD perturbation theory for the final matching step. The two methods are very different in approach; one uses gauge-invariant meson correlators (2-point functions) in position space that are extrapolated to the continuum limit before matching to perturbation theory; the other uses gauge-noninvariant 2-point and 3-point functions in momentum space, obtaining a renormalisation factor at each value of the lattice spacing. Both methods have mechanisms for testing and estimating systematic uncertainties within them and so both are capable of yielding a complete error budget for the final result. A comparison of the two methods is important to make sure that the uncertainties are being fully controlled. The best comparison in this respect is a direct one between the two different methods for the same lattice QCD quark formalism on the same gauge field configurations. This is the comparison that we provide here, for the first time.
to date of each method. We compare earlier results from the improved current-current correlator method (method c) from [12] to those obtained here using the RI-SMOM intermediate scheme [24], which improves on earlier momentum-subtraction schemes in having smaller nonperturbative and perturbative uncertainties 2 .
The lattice QCD results that we give here are for the matching factor, Z SMOM m , determined nonperturbatively on the lattice between the HISQ quark mass and that in the symmetric MOM (RI-SMOM) scheme at multiple  [22], 'HPQCD HISQ JJ + mc/ms' (red) is from [12], 'ETMC RI-MOM' is from [20], 'MP(hotQCD) HISQ JJ + mc/ms' is from [16], 'RBC/UKQCD DW RI-SMOM' is from [67] (note that this uses a different method to fix Zq than is used here), 'BMW clover RI-MOM' is from [68,69], 'HPQCD HISQ JJ + mc/ms' (blue) is from [13,17]. The light grey shaded band shows the evaluation of 96(+8, −4) GeV from the Particle Data Group [66]. The dark grey shaded band gives the weighted average (allowing for correlations) of the n f = 2 + 1 + 1 results as described in the text. different scales between 2 and 5 GeV. We obtained results at different valence (Section IV B) and sea (Section IV A) quark masses and different spatial volumes (Section IV C) to understand in detail what the sources of lattice systematic uncertainty are. We find that variations with sea mass and volume are barely discernible over a large range and valence mass effects are small (in contrast to those seen in the RI-MOM intermediate scheme). We combine results for the matching factor at 3 different values of the lattice spacing with the tuned c and s HISQ masses at that value of the lattice spacing obtained from the HPQCD current-current correlator calculation [12]. We are then able, by including an SMOM to MS matching factor (where we include the impact of having a non-zero c mass in the sea), to extrapolate the resulting masses in the MS scheme to the a = 0 continuum limit (see Section V). By having results at a range of values of the scale, µ, we are able to include a systematic uncertainty for remaining nonperturbative (condensate) contributions that depend on inverse powers of µ. These are expected to be much smaller than in the RI-MOM case, but we demonstrate their existence at several points in our calculation and they cannot be ignored. We find a residual effect of 0.2 % in our error budget from these nonperturbative contributions (see Table VI).
The values we obtain at the reference scale of 3 GeV for the c and s quark masses are: These results are to be compared to those from the current-current correlator method [12] using the same formalism on the same gluon field configurations. The m c value of 0.9851(63) GeV is obtained directly in the current-current correlator method; the value of m s of 0.0845(7) GeV (adjusting the value quoted in [12] from n f = 3 to 4) uses in addition a fully nonperturbative determination of the mass ratio m c /m s [12]. There is good agreement between the two sets of results, within the uncertainties quoted 3 . The uncertainties are small ( 1%) in both approaches and they are not strongly correlated between them. This is because the dominant sources of uncertainty are very different in the two cases: for the RI-SMOM calculation a key source of error is that from the determination of the tuned lattice quark masses, whereas the current-current correlator method is less sensitive to those and a larger source of uncertainty is that from missing higher order terms in the continuum QCD perturbation theory for that quantity. Note that α s can be determined in the same calculation in the current-current correlator approach and the correlation between α s and m c determined [12]. For this RI-SMOM calculation we must take a value of α s from elsewhere when we need one for the SMOM to MS conversion. Both the JJ and RI-SMOM mass determinations have uncertainties from lattice discretisation effects but they appear in different ways: in the RI-SMOM method it is the quark mass itself that is extrapolated to the continuum limit; in the current-current correlator method it is the time-moments of the correlator that are extrapolated. The agreement between the two methods is then a strong additional indication that the separate sources of systematic error are well controlled.
We can run our new RI-SMOM values for m c and m s given in eq. (29) to other scales for comparison with other results. The scale used for the m c is often that of m c itself and for m s , 2 GeV [66]. We give these values below using 4-loop perturbative QCD running in the MS scheme to run down from the higher scales at which the masses were determined: The information that the JJ and RI-SMOM methods agree means that an average of the two results should have a reduced uncertainty. We must take care in the averaging to allow for the correlations between the two methods and we do this by dividing the uncertainty in each case into correlated and uncorrelated pieces and then fitting the two results to a constant. The correlated portion comes from the tuning of the quark mass and determination of the lattice spacing (as given in the error budget) and is taken to be 100% correlated between the two methods. The breakdown in the uncertainty is then: In Figure 12 we compare our results for m c (m c , n f = 4) to previous lattice QCD results. This graphically shows the agreement between the RI-SMOM results here and the directly comparable current-current correlator results (the top and third from top values in the figure). We also include a result (second from the top in the figure) obtained recently by the Fermilab Lattice, MILC and TUMQCD collaborations from a new method they have developed using the minimal renormalon-subtracted (MRS) scheme [22,71]. This uses Heavy Quark Effective Theory to map out the heavy quark mass dependence of pseudoscalar heavy-light meson masses calculated in lattice QCD. A well-defined quark mass for this expansion is obtained by identifying and removing the leading renormalon from the perturbative expansion for the pole mass in terms of the MS mass. This MRS mass then has available a high-order continuum QCD perturbative matching to the MS scheme. Application of this method also yields 1%-accuracy and agrees, well within its uncertainty, with our JJ and RI-SMOM results.
Although the majority of results in Figure 12 have been obtained using the HISQ formalism there are results using other formalisms that demonstrate good agreement, for example the results from JLQCD using domain-wall quarks [15]. The results are divided into those obtained on gluon field configurations that include u, d, s and c quarks in the sea (as here) and those that include u, d and s quarks in the sea. Results obtained on gluon field configurations that include only u and d quarks in the sea are not shown, because it is not clear how to connect them in perturbative QCD (adding an s sea quark) to the values shown here. We see good agreement between almost all the results. The majority of accurate previous results used the current-current correlator method; the RI-MOM intermediate scheme has larger sources of systematic error for the reasons discussed in Section I. The current-current correlator results are tagged with 'a', 'b' or 'c' to denote different implementations of the method; this will be discussed further below.
We can provide a new world average for m c (m c ) allowing for correlations between the two HPQCD results by using the average given in eq. (33) and then combining in a weighted average with the ETMC result [20] and the Fermilab Lattice/MILC/TUMQCD result [22]. The ETMC result (using the RI-MOM approach in the Twisted Mass formalism), 1.348(42) GeV, is nearly 2σ from the HPQCD results and uncorrelated with it. The Fermilab Lattice/MILC/TUMQCD result is 1.273 (10) GeV and is correlated with the HPQCD JJ result because it uses the HPQCD determination of α s obtained concurrently with m c in [12]. The correlation coefficient between α s and m c is given there as 0.16. Since the uncertainty in the Fermilab Lattice/MILC/TUMQCD result is strongly dominated by that from α s we apply a correlation coefficient of 0.16 between that result and the HPQCD average. The HPQCD average will be slightly less correlated with it than the JJ result alone. This allows, however, for some further correlation through the fact that all of these calculations use some of the same sets of gluon field configurations and are done with the same quark formalism, although different lattice QCD quantities are calculated and the lattice spacing was fixed and quark masses tuned in a different way.
The weighted average of lattice QCD results including 4 flavours of sea quarks is then m c (m c , n f = 4) 2+1+1 av. = 1.2753(65) GeV, (34) shown as the dark shaded band on Figure 12. The average has a poor χ 2 /dof of 3 because of the tension between the ETMC result, which has almost no impact on the average, and the other three values. Our dark shaded band is a somewhat narrower band than the evaluation given in the Particle Data Tables [66], shown as the lighter shaded band.
B. World average for ms(2 GeV) Figure 13 provides a comparison of lattice QCD results for the s quark mass in the MS scheme at a scale of 2 GeV. Results are given for calculations with 4 flavours in the sea (as here) and also 3 flavours in the sea. There is very little difference (0.2 MeV, with n f =3 larger) between these from perturbative QCD. Again results for 2 flavours in the sea are not shown since they cannot be connected perturbatively to the more realistic 3 and 4 flavour results. There is reasonably good overall agreement between results using current-current correlator methods or the MRS scheme and m c /m s ratios and those using RI-MOM and RI-SMOM intermediate schemes direcetly for m s , as well as between results using a variety of quark formalisms. Our new RI-SMOM result, however, shows a 2.7σ tension with the earlier RI-SMOM result from RBC/UKQCD [67] using domain-wall quarks. In [67] the RI-SMOM implementation is slightly different, using the vector current vertex to fix Z q . Two values of the lattice spacing are used but only one value of µ (3 GeV) and possible nonperturbative effects are not included in the analysis or allowed for as a systematic uncertainty.
We can determine a new world average for the results from 4 flavour calculations allowing for correlations between the two different HPQCD results and the Fermilab/MILC/TUMQCD MRS result. We combine them at 3 GeV where we have correlated the breakdown of errors in the HPQCD results in eq. (31). There we take the uncertainties associated with fitting/scale-setting from the Fermilab/MILC/TUMQCD result to be 100% correlated with the HPQCD results (an overestimate given the different methods used) but take α s and statistical uncertainties to be uncorrelated, ignoring the relatively small correlations between the uncertainty coming from α s (which does not dominate in this case) and a part of the HPQCD JJ errors coming from m c . We combined this with the uncorrelated ETMC result run to 3 GeV. This average gives 0.08393(43) GeV with a χ 2 /dof of 2.5. Inflating the uncertainty by √ 2.5 to take account of this more general tension, and running down to 2 GeV gives m s (2 GeV, n f = 4) 2+1+1 av. = 0.09291(78) GeV. (35) This is given as the dark shaded band in Figure 13 to be compared to the light shaded band of the evaluation in the Particle Data Tables [66]. The Particle Data Tables result seems unduly pessimistic about our level of knowledge of the s quark mass and has a high central value, given the accuracy of lattice QCD results now available.

C. Future
Accurate though these results are, it is worth asking what the prospects are for reducing uncertainties further in the future. The original current-current correlator method (method a) [4] used lattice QCD results for quarks tuned to the c mass only and so one of the largest sources of uncertainty was from missing higher order terms in the perturbative series for time-moments because the scale of α s was related to m c . Subsequently (in method b) [13] heavier quark masses were included, giving access to m b but also reducing the perturbative uncertainty because the combined fit now included α s terms evaluated at higher scales. The newest variant, method c [12] also includes results for quarks with heavier masses than c. By making use of quark mass ratios, m c is then determined through a perturbative series for the time-moments in which the scale of α s is set by the heavier quark masses. This method then offers the potential to reduce the perturbative uncertainty by working on yet finer lattices where a given value of quark mass in lattice units corresponds to a higher quark mass. This will also reduce the sizeable uncertainty from the a → 0 extrapolation.
For the RI-SMOM intermediate scheme method used here, the largest sources of uncertainty are those from the bare tuned lattice quark masses which in turn depend on the determination of the lattice spacing. Working on finer lattices could cut this uncertainty significantly since the lattice spacing is fixed from f π in the continuum and physical u/d mass limit [55]. The impact of missing higher order terms in the SMOM to MS matching would be reduced by going to higher µ values and this is also possible on finer lattices. At the same time this would also reduce the impact of nonperturbative contributions since they fall rapidly with µ.
In conclusion, lattice QCD now has three very different methods for determination of quark masses at an accuracy of 1% or better. They yield consistent results for m c (m c ) and m s (2GeV). We have given here a particularly strong test of the consistency of the current-current correlator and RI-SMOM methods. Lattice QCD calculations are then well on the way to providing the accuracy needed for stringent future tests of the Standard Model. lana and E. Royo-Amondarain for gauge-fixing the superfine configurations. Computing was done on the Darwin supercomputer at the University of Cambridge High Performance Computing Service as part of the DiRAC facility, jointly funded by the Science and Technology Facilities Council, the Large Facilities Capital Fund of BIS and the Universities of Cambridge and Glasgow. We are grateful to the Darwin support staff for assistance. Funding for this work came from the National Science Foundation, the Royal Society, the Science and Technology Facilities Council and the Wolfson Foundation.
Appendix A: Perturbative matching from SMOM to MS for non-zero sea quark mass We have defined our SMOM scheme to have massive sea quarks so that we can work at physical sea quark masses. We then need a perturbative matching from SMOM to the MS scheme that allows for massive sea quarks. Sea quarks do not appear in the matching calculation until O(α 2 s ) and even then they contribute only a very small part to that coefficient. Given the very small changes that we expect, we need only to consider the case of the most massive sea quark that we have, i.e. the c quark. We work close to the physical mass for c and this in turn (in the MS scheme) is a sizeable fraction of µ for the µ values that we use. We might therefore expect the α 2 s coefficient from the c quark in the sea to be significantly different from that for a massless c quark. Our results show that indeed this is true for µ close to our lower value of 2 GeV but even so this has very little impact on Z MS/SMOM m . The two-loop coefficient c 2 of eq. (8) can be decomposed into two terms: a contribution which is free of as a result of having a massive c quark in the sea. The second column gives the change in Cn f , i.e. ∆Cm = Cn f (m) − Cn f (m = 0) (for n f = 1); the third column gives the change in the α 2 s coefficient c2 (eq. (8) and Table I) for five different values of the scale parameter, µ, which is given in the first column. In the fourth column we give values for αs in the MS scheme at each value of µ. These are obtained by four-loop running in the MS scheme from a value of 0.2128(25) at a scale of 5 GeV with n f = 4 [12]. -0.5 -0.002 0.2128 (25) internal sea quarks, C n f =0 , and a contribution which depends on them, C n f (m). It can be written as where C F = 4/3 and T F = 1/2 are the usual colour factors. The two-loop coefficient C n f =0 of Z MS/SMOM m as well as the piece C n f (m = 0) which is proportional to the number of massless sea quarks has been determined in [25,26]. The result for C n f (m = 0) reads It can be derived with the help of the QCD decoupling functions. In order to access the complete mass dependence, we also have calculated the exact result, which is valid in any mass region. It is plotted in Figure 14. We have checked by computing power corrections in µ 2 /m 2 to eq. (A3) that the expansion coincides as expected with the exact result for large values of m.
If we take m c (3 GeV) from [12], then the value of the ratio µ 2 /m 2 (µ) for µ = 2, 2.5, 3, 4 and 5 GeV reads µ 2 /m 2 (µ) = 3.4, 5.9, 9.3, 18.6 and 31.5. Table VII gives the resulting shifts ∆C m in C n f from the m = 0 result to the result that we need, which includes a massive c quark, at these five values of µ.
We also give the resulting shift ∆c 2 in c 2 of eq. (A1). We can see from Table VII that the only significant effect is for µ = 2 GeV where the shift is about 40% of the α 2 s coefficient. This coefficient is very small, however, and so the impact of this shift is also very small, less than about -0.2% on Z MS/SMOM m . We use the results in Table VII to shift the values of Z MS/SMOM m used in order to convert our SMOM results to the MS scheme in Section V.