Light quark masses in N f = 2 + 1 lattice QCD with Wilson fermions

We present a lattice QCD determination of light quark masses with three sea-quark flavours (Nf = 2+1). Bare quark masses are known from PCAC relations in the framework of CLS lattice computations with a non-perturbatively improved Wilson-Clover action and a tree-level Symanzik improved gauge action. They are fully non-perturbatively improved, including the recently computed Symanzik counter-term bA − bP. The mass renormalisation at hadronic scales and the renormalisation group running over a wide range of scales are known non-perturbatively in the Schrödinger functional scheme. In the present paper we perform detailed extrapolations to the physical point, obtaining (for the four-flavour theory) mu/d(2 GeV) = 3.54(12)(9) MeV and ms(2 GeV) = 95.7(2.5)(2.4) MeV in the MS scheme. For the mass ratio we have ms/mu/d = 27.0(1.0)(0.4). The RGI values in the three-flavour theory are Mu/d = 4.70(15)(12) MeV and Ms = 127.0(3.1)(3.2) MeV.


Introduction
The lattice regularisation of QCD provides a well-defined procedure for the determination of the fundamental parameters of the theory (i.e. the gauge coupling and the quark masses) from first principles. The aim of the present work is the determination of the three lightest quark masses (i.e. those of the up, down and strange flavours), in a framework in which the up and down quarks are degenerate and all heavier flavours (i.e. charm and above), if present in the theory, would be quenched (valence) degrees of freedom. This is known as N f = 2 + 1 lattice QCD. Moreover, QED effects are ignored.
The three-flavour theory adopted in this paper is presumably sufficient for determining light quark masses due to the decoupling of heavier quarks [1][2][3][4]. Indeed, lattice world averages of light quark masses m u/d , m s do not show a significant dependence on the number of flavours at low energies for N f ≥ 2 within present-day errors [5]. This also holds for the more accurately known renormalisation group independent ratio m u/d /m s . Recently, heavy-flavour decoupling has been substantiated also non-perturbatively [6].
This paper is based on large-scale N f = 2 + 1 flavour ensembles produced by the Coordinated Lattice Simulation (CLS) effort [7,8]. The simulations employ a tree-level Symanzik-improved gauge action and a non-perturbatively improved Wilson fermion action; see references [9][10][11][12]. The sea quark content is made of a doublet of light degenerate quarks m q,1 = m q,2 , plus a heavier one m q, 3 . At the physical point m q,1 , m q,2 = m u/d ≡ 1 2 (m u + m d ) and m q,3 = m s .
The bare quark masses produced by CLS [7,8] need to be combined with renormalisation and improvement coefficients in order to obtain renormalised quantities with O(a 2 ) discretisation effects. We use ALPHA collaboration results for the quark mass renormalisation and Renormalisation Group (RG) running [13] in the Schrödinger functional (SF) scheme. Symanzik improvement is implemented for the removal of discretisation effects from correlation functions, leaving us with O(a 2 ) uncertainties in the bulk and O(g 4 0 a) ones at the time boundaries. We find that correlation functions extrapolations to the continuum limit are compatible with an O(a 2 ) overall behaviour. The counter-terms required for the improvement of the axial current are known from refs. [14][15][16]. The present work has combined all these elements, obtaining estimates of the up/down and strange quark masses, as well as their ratio. These are expressed as renormalisation scale-independent and schemeindependent quantities, known as Renormalisation Group Invariant (RGI) quark masses. Of course we also give the same results in the MS scheme at scale µ = 2 GeV.
The bare dimensionless parameters of the lattice theory are the strong coupling g 2 0 ≡ 6/β and the quark masses expressed in lattice units am q,1 = am q,2 and am q,3 , with a the lattice spacing. They can be varied freely in simulations. Having chosen a specific discretisation of the QCD action, its parameters must be calibrated so that three "input" hadronic quantities (one for each bare mass parameter and one for the lattice spacing) attain their physical values. Other physical quantities can subsequently be predicted. Such input quantities are typically very well-known from experiment, but they also need to be precisely computed on the lattice; examples are ground state hadron masses and decay constants (m π , f π , m K , f K , . . .). Since the majority of numerical large-scale simulations do not yet include the small strong-isospin breaking and electromagnetic effects, the physical input quantities have to be corrected accordingly. Following ref. [8], When simulations approach the point of physical mass parameters while the lattice spacing is lowered, computational demands rapidly increase. In the present work results are obtained at non-zero lattice spacings and at quark masses which correspond to unphysical meson and decay constant values. Thus our data need to be extrapolated to the continuum limit and extra/interpolated to the physical quark mass values. This is achieved with a joint chiral and continuum extrapolation. The present work pays particular attention to these extrapolations and interpolations and the ensuing sources of systematic error.
So far we did not specify the reference scale f ref . In ref. [8] the three-flavor symmetric combination f phys πK = 2 3 (f phys K + 1 2 f phys π ) = 147.6(5) MeV, (obtained from the physical input of eqs. (1.1)) was used for calibration, and for the determination of the hadronic gradient flow scale t 0 [17]. 1 An artificial (theoretical) hadronic scale with mass dimension -2, t 0 is precisely computable with small systematic effects [18,19], and thus well-suited as intermediate scale on the lattice. Its physical value determined from CLS ensembles [8]  where the first error of 8t phys 0 is statistical and the second systematic. The theoretical framework of our work is explained in Section 2. The definitions of bare current quark masses, their renormalisation parameters and the O(a)-improvement counter-terms are provided in standard ALPHA-collaboration fashion. There is also a fairly detailed exposition of how the so-called "chiral trajectory" (a Line of Constant Physics -LCP) is traced by N f = 2 + 1 CLS simulations. In Section 3 we outline the computations leading to renormalised current quark masses as functions of the pion squared mass. These are computed in the SF renormalisation scheme at a hadronic (low energy) scale. In Section 4 we perform the combined chiral and continuum limit extrapolations in order to obtain estimates of the physical up/down and strange quark masses. Details of the ansätze we have used are provided in Appendices A and B. Our final results are gathered in Section 5. Preliminary results have been presented in [20].

Theoretical framework
We review our strategy for computing light quark masses with improved Wilson fermions.
In what follows equations are often written for a general number of flavours N f . In practice N f = 2 + 1. Flavours 1 and 2 indicate the lighter fermion fields, which are degenerate; at the physical point their mass is the average up/down quark mass. Flavour 3 stands for the heavier fermion, corresponding to the strange quark at the physical point.

Quark masses, renormalisation, and improvement
The starting point is the definition of bare correlation functions on a lattice with spacing is a and physical extension L 3 × T : where the pseudoscalar density and axial current are The indices i, j = 1, 2, 3 label quark flavours, which are always distinct (i = j). The bare current (or PCAC) quark mass is defined via the axial Ward identity at zero momentum and a plateau average between suitable initial and final time- with the source P ji positioned either at y 0 = a or y 0 = T − a. 2 The mass-independent improvement coefficient c A is determined non-perturbatively [14]. The average of two renormalised quark masses is then expressed in terms of the PCAC mass m ij as follows: (2.5) where M q ≡ diag(m q,1 , m q,2 , · · · , m q,N f ) is the matrix of the sea quark subtracted masses, characteristic of Wilson fermions. Given the bare mass parameter m 0,i ≡ (1/κ i − 8)/(2a), with κ i the Wilson hopping parameter, these are defined as where m cr ∼ 1/a is an additive mass renormalisation arising from the loss of chiral symmetry by the regularisation and κ cr is the critical (chiral) point. The average of two subtracted masses is then denoted by m q,ij ≡ 1 2 (m q,i + m q,j ) in eq. (2.5). The axial current normalisation Z A (g 2 0 ) is scale-independent, whereas the current quark mass renormalisation parameter 1/Z P (g 2 0 , aµ) depends on the renormalisation scale µ. The renormalisation condition imposed on the pseudoscalar density operator P ji defines the renormalisation scheme for the quark masses. The schemes used in the present work (SF and MS) are mass-independent. Pertinent details will be discussed in latter sections.
The improvement coefficients b A − b P andb A −b P of eq. (2.5) cancel O(a) massdependent cutoff effects; they are functions of the bare gauge coupling g 2 0 . The corresponding counter-terms of eq. (2.5) contain the subtracted masses am q,ij and Tr[aM q ], which require knowledge on the critical mass m cr . This can be avoided by substituting these masses with current quark masses and their sum. Their relationship is [21], where Z(g 2 0 ) ≡ Z P /(Z S Z A ) and r m (g 2 0 ) are finite normalisations. Z S is the renormalisation parameter of the non-singlet scalar density S ij ≡ψ i ψ j and r m /Z S is the renormalisation parameter of the singlet scalar density, which indirectly defines r m ; cf. ref. [21]. In the above we neglect O(a) terms, as they only contribute to O(a 2 ) in the b-counter-terms of eq. (2.5). Substituting am q,ij → am ij in the latter expression, we obtain To leading order in perturbation theory the difference b A − b P is O(g 2 0 ) and equalsb A −b P . However, non-perturbative estimates are likely to differ significantly, especially in the range of couplings considered here (1.56 g 2 0 1.76). We will employ non-perturbative estimates of b A − b P and Z; cf. ref. [22]. The term multiplying M sum contains (1 − r m )/r m and (b A −b P ). In perturbation theory r m = 1 + 0.001158 [21]. A first non-perturbative study of the coefficientsb A andb P produced noisy results with 100% errors [25]. Given the lack of robust non-perturbative results and the fact that the term in curly brackets is O(g 4 0 ) in perturbation theory, it will be dropped in what follows.
Once the quark mass averages m 12R and m 13R are computed say, in the SF scheme at a scale µ had , the three renormalised quark masses can be determined. Since we are working in the isospin limit (m q,1 = m q,2 ), the lighter quark mass is given by m 12R . Then one can isolate m 13R from the ratio m 13R /m 12R in which, as seen from eq. (2.8), the M sum counter-term cancels out.
We will be quoting results also for the scheme-and scale-independent renormalisation group invariant (RGI) quark masses M 12 and M 13 (corresponding to the current masses m 12 and m 13 ) as well as the physical RGI quark masses M u/d and M s derived from them. They are conventionally defined in massless schemes [38] by for each quark flavour i. In our opinion, M i is better suited for comparisons either to experimental results or other theoretical determinations. Equation (2.10) is formally exact and independent of perturbation theory as long as the renormalised parameters (g R , m iR ) and the continuum renormalisation group functions (i.e. the Callan-Symanzik β-function and the mass anomalous dimension τ ) are known non-perturbatively with satisfacory accuracy [13,[32][33][34][35]. Their computation in the SF scheme with N f = 3 massless quarks has been carried out in ref. [13]. Our determination of the renormalised quark masses is based on the bare current mass averages m ijR ; cf. eqs. (2.5) and (2.8). The analogue of these expressions for the RGI mass averages is given by Note that the ratio M/m R (µ had ) is flavour-independent; cf. eq. (2.10). In ref. [13] it has been computed in the SF scheme for the N f = 3 massless flavours at µ had = 233(8) MeV.

The chiral trajectory and scale setting
Our aim is to stay on a line of constant Physics within systematic uncertainties of O(a 2 ), as we vary the bare parameters of the theory (i.e. the gauge coupling g 0 and the N f = 2+1 quark masses). In particular, if the improved bare gauge coupling is kept fixed in the simulations, so does the lattice spacing, with any fluctuations being attributed to O(a 2 )-effects. The problem is that b g (g 2 0 ) is only known to one-loop order in perturbation theory [39,40] Thus, following refs. [41,42], we vary the quark masses at fixed g 2 0 , ensuring that the trace of the quark mass matrix remains constant: Tr[M q ] = 2m q,1 + m q,3 = const . (2.13) In this way the improved bare gauge couplingg 2 0 is kept constant at fixed β for any b g . 3 This requirement leads to an unusual but unambiguous approach to the physical point, shown in the (M 12 , M 13 )-plane in the left panel of figure 1. Initially, one starts at the symmetric point (am q,1 = am q,2 = am q,3 = am sym q ) for some fixed coupling β = 6/g 2 0 , and tunes the mass parameter of the simulation in such a way that Tr  12 . Coloured (gray) points correspond to mass-shifted (-unshifted) points in parameter space, cf. the discussion in the text. a good approximation 4 . This is achieved by varying am sym q until (m 2 K + 1 2 m 2 π )/f ref takes its physical value. Since it is proportional to Tr[M q ] at leading order in chiral perturbation theory (χPT), it suffices as tuning observable. In subsequent simulations, one successively lifts the mass-degeneracy towards the physical point by decreasing the light quark masses while maintaining the constant-trace condition. By doing so the physical strange quark mass is approached from below as in figure 1 (left panel). We call this procedure "the determination of the chiral trajectory".
Note that the improved renormalised quark mass matrix M R is given by [21] Tr Since the d m -counter-term is proportional to squared bare masses, a constant Tr[M q ] does not correspond to a constant Tr[M R ]; the latter requirement is violated by O(a) effects. This is an undesirable feature, as it implies that the chiral trajectory is not a line of constant-physics. In practice these violations have been monitored in ref.
and t 0 is the gluonic quantity of the Wilson flow [17]; it has mass dimension -2. Here m π and m K are the lightest and strange pseudoscalar mesons respectively; at the physical point these are the pion m phys π and kaon m phys K . Keeping φ 4 constant is a Symanzik-improved constant physics condition. But φ 4 is proportional to the sum of the three quark masses only in leading-order (LO) chiral perturbation theory (χPT). Thus, the improved bare couplingg 2 0 now suffers from O(am q Tr[M q ]) discretisation effects due to higher-order χPT contributions. In practice, these turn out to be small, as can be seen from ref. [8] (Fig. 4, lowest rhs panel), where Tr[M R ] has been computed, at constant φ 4 . The violations appear to be at most 1% and thus the variation of the O(a) b g -term ing 2 0 can be ignored. Obviously, one must also ensure, through careful tuning, that the chosen φ 4 = const. trajectory passes through the point corresponding to physical up/down and strange renormalised masses (i.e. quark masses that correspond to the physical pseudoscalar mesons m phys π and m phys K ). This is done by driving φ 4 to its physical value φ phys 4 = 8t 0 [(m phys K ) 2 + (m phys π ) 2 /2] through mass shifts [8]. The aim is to express the computed quantities of interest (in our case the quark masses) as functions of with φ 4 held fixed at φ phys 4 , and eventually extrapolate them to φ phys 2 = 8t 0 (m phys π ) 2 . The determination of the redefined chiral trajectory is not straightforward. One needs to know the value of φ phys 4 . The latter is obtained from t 0 and the pseudoscalar masses (corrected for isospin-breaking effects) quoted in eq. (1.1). But since the value of t 0 is only approximately known, one starts with an initial guesst 0 , which provides an initial guess φ 4 . At each β, the symmetric point with degenerate masses (κ 1 = κ 3 ) is tuned so that the computed t 0 /a 2 , am π and am K combine as in eq. (2.15) to give a value close toφ 4 . The other ensembles at the same β have been obtained by decreasing the degenerate (lightest) quark mass m q,1 = m q,2 , while increasing the heavier mass m q,3 so as to keep Tr[M q ] constant. Thus they do not correspond exactly to the sameφ 4 . Small corrections of the subtracted bare quark masses (or hopping parameters) are introduced, using a Taylor expansion discussed in sect. IV of ref. [8], in order to shift φ 4 to the reference valueφ 4 and correct analogously the measured PCAC quark masses and other quantities of interest such as the decay constants. The procedure is repeated for each β and the same valueφ 4 at the starting symmetric point.
All shifted quantities are now known atφ 4 as functions of φ 2 . Defining the combination of decay constants the dimensionless 8t 0 f πK is computed for all φ 2 and extrapolated toφ 2 = 8t 0 (m phys π ) 2 . The extrapolated t 0 f πK , combined with the experimentally known f phys πK , gives a better estimate oft 0 , and thus ofφ 4 . As described in sect. V of ref. [8], this procedure can be recursively repeated and eventually the physical value of t 0 is fixed through f phys πK ; the value in eq. (1.2) from ref. [8] Table 1: Details of CLS configuration ensembles, generated as described in ref. [7]. In the last column, ensembles are labelled by a letter, denoting the lattice geometry, a first digit for the coupling and a further two digits for the quark mass combination.
In analogy to the definitions (2.15) and (2.16), we also define rescaled dimensionless bare current quark masses and their renormalised counterparts at scale µ had The redefined chiral trajectory is shown in figure 1 (right panel), where the light-light and heavy-light dimensionless mass averages (φ 12R and φ 13R respectively) are plotted as functions of φ 2 . Extrapolating in φ 2 to φ phys 2 amounts to the simultaneous approach of the light and heavy quark masses to the corresponding physical up/down and strange values. All other physical quantities are then also at the physical point. Sect. 4 is dedicated to these extrapolations.

Quark mass computations
We base our determination of quark masses on the CLS ensembles for N f = 2 + 1 QCD, listed in Table 1. The bare gauge action is the Lüscher-Weisz one, with tree-level coefficients [11]. The bare quark action is the Wilson, Symanzik-improved [10] one. The Clover term coefficient c sw has been tuned non-perturbatively in ref. [12]. Boundary conditions are periodic in space and open in time, as detailed in ref. [44].
For details on the generation of these ensembles see ref. [7]. As seen in Table 1, results have been obtained at four lattice spacings in the range 0.05 a/fm 0.086. For each lattice coupling β = 6/g 2 0 , gauge field ensembles have been generated for a few 5 values of the Wilson hopping parameters κ 1 = κ 2 and κ 3 . The light pseudoscalar meson (pion) varies between 200 MeV and 420 MeV. The heaviest value corresponds to the symmetric point where the three quark masses and the pseudoscalar mesons are degenerate. The strange meson (kaon) varies between 420 MeV and 470 MeV. Given that our lightest pseudoscalars are relatively heavy (200 MeV), the chiral limit ought to be taken with care.
The bare correlation functions f ij P , f ij A of eqs. (2.1) are estimated with stochastic sources located on time slice y 0 , with either y 0 = a or y 0 = T − a. From them the current quark masses m 12 , m 13 are computed as in eq. (2.4), with the O(a)-improvement coefficient c A determined non-perturbatively in ref. [14]. The exact procedure to select the plateaux range in the presence of open boundary conditions has been explained in refs. [7,43,45].
Having obtained the bare current quark masses m 12 , m 13 at four values of the coupling g 2 0 , we construct the renormalised dimensionless quantities m 12R (µ had ) and m 13R (µ had ); cf. eq. (2.8). For this we need the ratio Z A (g 2 0 )/Z P (g 2 0 , µ had ) and the Symanzik b-counterterms. Results for the axial current normalisation Z A (g 2 0 ) are available in ref. [46], from a separate computation based on the chirally rotated Schrödinger Functional setup of refs. [47][48][49]. The computation of Z P (g 2 0 , µ had ) in the SF scheme, for µ had = 233 (8) MeV, was carried out in ref. [13] for a theory with N f = 3 massless quarks and the lattice action of the present work. The Z P results, shown in eqs. (5.2) and (5.3) of ref. [13], are in a range of inverse gauge couplings which covers the β ∈ [3.40, 3.85] interval of the large volume simulations of ref. [8], from which our bare dimensionless PCAC masses are extracted.
Besides the ratio Z A /Z P , we also need the improvement coefficient (b A −b P ), which multiplies the O(a) counter-term proportional to am ij in eq. (2.8). To leading order in perturbation theoryb A −b P = −0.0012g 2 0 . Non-perturbative estimates based on a coordinate-space renormalisation scheme have been provided for N f = 2 + 1 lattice QCD in ref. [25]. More accurate non-perturbative results have been subsequently obtained by the ALPHA Collaboration, using suitable combinations of valence current quark masses, measured on ensembles with N f = 3 nearly-chiral sea quark masses in small physical volumes [16,22]. Also these simulations have been carried out in an inverse coupling range that spans the interval β ∈ [3.40, 3.85] of the large volume CLS results of ref. [8]. They are expressed in the form of ratios R AP and R Z , from which (b A − b P ) and Z are estimated; thus (b A −b P ) = R AP /R Z . In ref. [22], results are quoted for two values of constant Physics, dubbed LCP-0 and LPC-1. In LCP-0, R AP and R Z are obtained with all masses in the chiral limit. In LCP-1, one valence flavour is in the chiral limit (so it is equal to the sea quark mass), while a second one is held fixed to a non-zero value. The physical volumes are always kept fixed. In ref. [22], eqs. (5.1), (5.2) and (5.3) refer to LCP-0 results, while those in eqs. (5.1), (5.4) and (5.5) refer to LCP-1; differences are due to O(a) discretisation effects.
We have opted to use the LCP-0 values ofb A −b P in the present work. The covariance matrices of the fit parameters of R AP as well as those of R Z are provided in ref. [22]. We assume that the covariance matrix between fit parameters of R AP and R Z is nil. This is justified a posteriori, by repeating the analysis with LCP-1 values, as a means to estimate the magnitude of systematic errors arising from our choice. Moreover, we have also compared our LCP-0 results to those obtained from different fit functions, used in the preliminary analysis of ref. [16], as well as from the perturbative estimateb A −b P . We find that the contribution arising from such variations is below ∼ 1% of the total error on renormalised quark masses at the physical point.
As discussed in Section 2, the complicated Symanzik counter-term in curly brackets, multiplying aM sum in eq. (2.8), is O(g 4 0 a) in perturbation theory. As there are no robust non-perturbative estimates of its magnitude at present, we will drop this term, assuming  Rescaled dimensionless current quark masses φ 12 and φ 13 , renormalised in the SF scheme at µ had , for each CLS ensemble used in our analysis. Note that for simulation points H102, H105, C101 more than one independent ensembles exist, which have been run with different algorithmic setups; we keep those separate before fits. All points have been shifted to the target chiral trajectory as described in the text, and the quoted errors contain both statistical uncertainties and the contribution from renormalisation and the mass shift.
that the O(g 4 0 a) effects it would remove are subdominant compared to O(a 2 ) uncertainties. As already explained in subsect. 2.2, our analysis is based on the rescaled dimensionless quantities defined in eqs. (2.15), (2.16), and (2.20). At each β value and for each gauge field configuration, we have results for t 0 /a 2 , am 12 and am 13 from refs. [7,8], from which φ 12 and φ 13 are obtained. The error analysis is carried out using the Gamma method approach [50][51][52][53] and automatic differentiation for error propagation, using the library described in ref. [54]. This takes into account all the existing errors and correlations in the data and ancillary quantities (renormalisation constants, improvement coefficients, etc.), and estimates autocorrelation functions (including exponential tails) to rescale the uncertainties correspondingly. Following [8], the estimate of the exponential autocorrelation times τ exp used in the analysis is the one quoted in [7], viz., We have checked that without attaching exponential tails statistical errors are 40 to 70% smaller in our final results. The full analysis has been crosschecked by an independent code based on (appropriately) binned jackknife error estimation.
The starting values for φ 12R and φ 13R on which the analysis is based are shown in Table 2, where renormalised quark masses are in the SF scheme at a scale µ had = 233(8) MeV. By suitably fitting these quantities as functions of φ 2 , and extrapolating to φ phys 2 , we obtain the results for physical up/down and strange quarks at scale µ had , as detailed in sect. 4. Only then do we convert them to the RGI masses, by multiplying them with the RG-running factor [13] M m R (µ had ) = 0.9148(88) , with the error added in quadrature; cf. eq. (2.11). Before presenting our chiral fits in section 4, we conclude this section with a comment on finite-volume effects. Current quark masses are not expected to be affected by finitevolume corrections, since their values are fixed by Ward identities. On the other hand, meson masses, decay constants, and the ratio t 0 /a 2 are expected to suffer from such effects. Standard SU(3) χPT NLO formulae are available for masses and decay constants [55]; t 0 /a 2 does not suffer from finite-volume effects up to NNLO corrections [19]. In particular, the χPT-predicted effects for meson masses are below the percent level, since the lattice spatial size in units of the inverse lightest pseudoscalar meson mass is in the range [3.9, 5.8].
On the other hand, by directly comparing the values in Table 2 obtained at the same lattice spacing and sea quark masses but different volumes (cf. ref. [43]), it is seen that the finitevolume effects on t 0 /a 2 and m 2 π are comparable and come with opposite signs. As a result, they largely cancel in φ 2 , the variable in which chiral fits are performed. Decay constants, which generally suffer from larger finite-volume effects than meson masses, enter our computation indirectly only -firstly through NLO terms in chiral fits, where the finite-volume correction is sub-leading, and secondly through the physical value of √ 8t 0 determined in [8], where these corrections have already been taken into account. We therefore expect that the quantities most affected by finite-volume effects are the rescaled current quark masses φ 12 , φ 13 , due to the presence of √ 8t 0 /a in their definition. As mentioned above, these are much smaller than our statistical uncertainty, cf. Table 2.
In the rest of our analysis we will therefore neglect this source of uncertainty.

Extrapolations to physical quark masses
Having obtained the dimensionless renormalised current mass combinations φ 12R and φ 13R at each β as functions of φ 2 , we now proceed with the determination of the physical values φ ud and φ s . This is done in fairly standard fashion through fits and extrapolations. To begin with, we note that the two lighter degenerate quark masses are simply given by φ 12 , whereas the heavier strange one is obtained from the difference 6 φ h = 2 φ 13 − φ 12 . (4.1) It is then possible to perform simultaneous fits of φ 12 and φ h as functions of φ 2 and the lattice spacing, subsequently extrapolating the results to φ phys 2 of eq. (2.19) and the continuum limit, so as to obtain φ ud and φ s . Variants of this method consist in simultaneous fits and extrapolations of either φ 13 or φ 12 on one hand and their ratio φ 12 /φ 13 on the other. These turn out to be advantageous, as does a certain combination of ratios involving φ 12 , φ 13 , φ 2 , and φ 4 , for reasons discussed below. We recall in passing that in the ratio φ 12 /φ 13 all renormalisation factors cancel.
We use fits based in chiral perturbation theory (χPT fits) which are expected to model the data well close to the chiral limit φ 2 = 0. Recall that we have performed N f = 2 + 1 simulations on a chiral trajectory; starting from a symmetric point where all quark masses are degenerate, we increase the mass of the heavy quark while decreasing that of the light one, until the physical point is reached. Since both masses are varying, it is natural to use SU(3) L ⊗ SU(3) R chiral perturbation theory, which bears explicit dependence on both masses. This works when all three quark masses in the simulations are light enough for say, NLO χPT with three flavours to provide reliable fits. In ref. [56] it is stated that this is the case for their data, obtained with domain wall fermions, as long as the average quark mass satisfies am avg < 0.01. As seen in Table 2 of ref. [8], our PCAC dimensionless quark masses am 12 and am 13 also satisfy this empirical constraint. The real test comes about a posteriori, when the SU(3) L ⊗ SU(3) R NLO ansätze are seen to fit our results well.
In Appendices A and B, ansätze for NLO χPT and discretisation effects are adapted to our specific parametrisation in terms of φ 2 and φ 4 . For the current quark masses these are where p 1 , p 2 , p 3 (4.4) As discussed in Appendix B, the form of the cutoff effects respects the constraint φ 12 /φ 13 = 1 at the symmetric point m q,1 = m q,3 , which is exact at all lattice spacings by construction. For the combination defined in eq. (A.14), we have An alternative to NLO χPT fits is the use of power series, based simply on Taylor expansions around the symmetric point m q1 = m q2 = m q3 , for which φ sym 2 = 2φ phys 4 /3: Note that imposing the constraint φ 12 = φ 13 at the symmetric point implies that s 0 and S 0 are common fit parameters. These expansions are expected to give reliable results in the higher end of the φ 2 range, underperforming close to the chiral limit. They are thus complementary to the chiral fits, which are better suited for the small-mass regime. In this sense the two approaches may provide a handle to estimate the systematic uncertainties due to these fits and extrapolations. We explore various fit variants, in order to unravel the presence of potentially significant systematic effects. They are encoded as follows: • Fitted quantities and ansätze: [chi12] Fit of φ 12 data only, using the χPT ansatz.
[tay12] Fit of φ 12 data only, using the Taylor expansion ansatz.
[tay13] Fit of φ 13 data only, using the Taylor expansion ansatz.
• Cuts on pseudoscalar meson masses: [420] Fit all available data, including the symmetric point; i.e. data satisfies m π 420 MeV.
Any given fit will thus be labelled as [a2] we find that they are sensitive to the presence of a discretisation term ∝ a 2 /t 0 , albeit within ∼ 1 − 2σ. This difference is attenuated when the more stringent mass cutoff [chi12][300] is enforced, mainly because the error increases as less points are fitted. The same qualitative conclusions are true for the Taylor expansion fits [tay12] of the light quark mass. On the other hand, the lower panel of Fig. 2 shows that the average quark mass M phys 13 is not sensitive to the details of the fit ansätze. This is not surprising, given that our simulations have been performed in a region of rather heavy pions 220 MeV ≤ m π ≤ 420 MeV, with data covering the physical point M phys 13 , while M u/d requires long extrapolations. The conclusion is that independent fits are reliable for φ 13 but less so for φ 12 , and so we discard their results.
Combined fits to φ 12 and φ 13 : Fig. 3  give results which are sensitive to the ansatz employed for the cutoff effects. This is more    A few general points concerning the fit analysis deserve to be highlighted: • Unsurprisingly, the inclusion of a second discretisation term ∝ φ 2 (a 2 /t 0 ) in the fits contributes to an increase of the error. This term is often compatible with zero, • Within large uncertainties, the coefficients of the leading cutoff effects (i.e. those ∝ a 2 /t 0 ) depend on the fitted observable, and are larger for φ 13 than for φ 12 .  are generally in good agreement, signalling the consistency of the approach. It is also interesting to note that the coefficient of the quadratic term is very small and always compatible with zero within 1σ (save for two cases where it vanishes within 2σ).
• NLO χPT appears to be suffering around and above 400 MeV.

Final results and discussion
Following the analysis of sect. 4, we quote as final results those obtained from the following procedure: • The central values are those of a combined fit to the ratio φ 12 /φ 13 and the quantity 2φ 13 /φ K + φ 12 /φ 2 , using NLO χPT, with pseudoscalar meson masses less than 360 MeV and a discretisation term proportional to a 2 /t 0 (i.e., fit [chirr][a1][360]). The error from this fit will appear as the first uncertainty in the results below. The fit is illustrated in Fig. 4. fits. This is the second error of the results below. Recall that [chirc] are combined fits to φ 13 and the ratio φ 12 /φ 13 , using NLO χPT.
• All results have been obtained using the Symanzikb-parameters computed in the LCP-0 case (see discussion in sect. 3). Using LPC-1 results instead, has very marginal effects on the error. The quark mass ratio is obtained from Dependence on renormalisation is only implicit, from the joint fit with φ 13 . The same procedures as above yield The above results for RGI masses refer to the N f = 2 + 1 theory. It is customary in phenomenological studies to report light quark masses measured in the N f = 2+1 lattice theory in the MS scheme at 2 GeV, referred to the more physical QCD with four flavours. This entails using N f = 3 perturbative RG-running from 2 GeV down to the charm threshold, followed by N f = 4 perturbative RG-running back to 2 GeV; see for example ref. [5]. We use 4-loop perturbative RG-running and the value for the Λ MS QCD parameter computed by the ALPHA Collaboration in ref. [33]  The mass ratio is obviously the same as in eq. (5.3). We note in passing that switching to the four-flavour theory has a very small effect on MS results, since at 2 GeV the matching factor is m R (N f = 4)/m R (N f = 3) = 1.002. The error budget for our computation is summarised in Table 3 and Fig 5. Uncertainties are completely dominated by our chiral fits. We have separated these errors into two contributions; see first two lines of Table 3. The first error is that of our best fit [chirr][a1] [360], and includes the statistical errors as well as the error from combined fits in φ 2 and a. The second uncertainty is the one arising upon varying the fit ansätze and their φ 2 range. All other errors are clearly seen to be subdominant. It is worth noting 7 In converting our results to MS we have taken into account the uncertainty in the matching factor coming from the error on Λ MS QCD , as well as the covariance of the latter with our determination of Mi.

A Chiral perturbation theory expansions
We adapt standard χPT expressions to our specific parametrisation of the data, stemming from our choice of chiral trajectory. We start, for example, from eqs. (B5) and (B6) of ref. [56], which are NLO chiral expansions of the light and strange pseudoscalar mesons m π and m K in terms of light and strange quark masses m 1 = m 2 and m 3 . These series are inverted, so that quark masses are functions of meson masses. The PCAC quark mass combinations m 12 = (m 1 + m 2 )/2 and m 13 = (m 1 + m 3 )/2 are then formed and everything is re-expressed in terms of the dimensionless quark masses φ 12 , φ 13 and the dimensionless quantities φ 2 and φ 4 , so that we arrive at where B 0 , f 0 , L k (k = 4, 5, 6, 8) are standard χPT parameters and The NLO LECs L i and B 0 are implicitly renormalised at scale Λ χ . It is also useful to consider the ratio which does not require renormalisation. Note that at the symmetric point (φ 2 = 2φ 4 /3) the current quark masses of eqs. (A.1) and (A.2) respect the constraint φ 12 = φ 13 , while the ratio (A.4) is exactly 1. Note that the sum of ratios has the remarkable advantages of depending on just one combination of NLO LECs, and of being free of polynomial dependence on φ 2 .
We next rewrite eqs. (A.1) and (A.2) in forms which are suitable for combined fits, with common coefficients for φ 12 and φ 13 , obtaining where the coefficients p 1 , p 2 , and p 3 relate to LECs as follows: with f πK given by eq. (2.17). The chiral logarithms are The following points should be kept in mind: • We are using only configurations along the φ 4 = constant chiral trajectory. Terms proportional to φ 4 are thus reabsorbed into constant fit terms.
• Our expressions are linear in fit parameters, rather than non-linear factors in which LECs appear explicitly. Determination of LECs is beyond the scope of the present work.
• By replacing f 2 0 by f 2 πK in the above definition of K, the coefficients of chiral logarithms are completely fixed relative to the LO value; cf. eqs. (A.1) and (A.2). This eliminates one fit parameter, pushing its effect to NNLO LECs. In practice, the fact that terms with φ 4 are reabsorbed into the LO terms nullifies the effect in some fits, e.g., those for φ 12 and φ 13 . A second advantage of this choice is that the resulting ansätze are fully linear in the fit parameters. See also eq. (2.5) in [8] and comments therein on the reasons that f 0 ≈ f πK and for preferring f πK to f 0 .
• We conveniently set the renormalisation scale to Λ χ = 1/ √ 8t 0 476 MeV, simplifying the chiral logs. There is no need to reabsorb ln(8t 0 Λ 2 χ ) terms in fit parameters. This is an unconventional choice, as common practice consists in providing results for LECs at Λ χ = m ρ or Λ χ = 4πf 0 . Consequently, NLO LECs eventually obtained with our methodology may only be compared to results in the literature after some extra work.
Using the above expressions and consistently neglecting higher mass orders, we obtain for the ratio (A.4) of PCAC masses For the combination (A.5) we have (A.14) With φ 4 held constant, the quantities of eqs. (A.6), (A.7), (A.13), and (A.14) are functions of φ 2 only. We use these expressions to fit our data, after adding O(a 2 ) terms which model leading discretisation effects that have been neglected throughout this Appendix.
expressing the current quark masses on the rhs of the above equation in terms of φ 12 and φ 13 , followed by using their LO χPT relations to φ 2 and φ 4 . In particular, with β 0 ≡ 1/(2B 0 √ 8t 0 ), we see from eqs. (A.6), (A.7) that to LO: Inserting the above LO expressions in eq. (B.2) we obtain, after some straightforward algebra, that for two light quarks the discretisation function has the form 2 ) contributions of f 12 and f 13 . Moreover, standard power-counting schemes in Wilson χPT [64,65] suggest that terms of O(a 2 ) enter at the same order as O(m 2 π ), which would imply that the terms of O(a 2 φ 2 ) should also be dropped. We will nevertheless keep this term and explore its impact.
For the ratio of φ 12 and φ 13 we have that The coefficients D 0 , D 1 , D 2 , . . . depend on the c i 's. The factor 2φ 4 −φ 2 in the discretisation term vanishes at the symmetric point φ 2 = 2φ 4 /3. This confirms that at the symmetric point the ratio φ 12 /φ 13 is 1 by construction, for any lattice spacing. In analogy to the arguments exposed above for f 12 and f 13 , we drop the D 2 φ 2 2 term in our fits. Moreover, the variation of the denominator (2φ 4 − φ 2 ) 2 is relatively mild, ranging between ∼ 2 and ∼ 4.6 as φ 2 varies between ∼ 0.1 and ∼ 0.8 in our simulations. To simplify matters, we reabsorb this O(1) term in re-definitions of D 0 and D 1 . Finally, for the combination of eq. (A.14), we straightforwardly parametrise the discretisation errors in a way analogous to f 12 and f 13 ; see eq. (4.5).