Up, down, strange and charm quark masses with Nf = 2+1+1 twisted mass lattice QCD

We present a lattice QCD calculation of the up, down, strange and charm quark masses performed using the gauge configurations produced by the European Twisted Mass Collaboration with Nf = 2 + 1 + 1 dynamical quarks, which include in the sea, besides two light mass degenerate quarks, also the strange and charm quarks with masses close to their physical values. The simulations are based on a unitary setup for the two light quarks and on a mixed action approach for the strange and charm quarks. The analysis uses data at three values of the lattice spacing and pion masses in the range 210 - 450 MeV, allowing for accurate continuum limit and controlled chiral extrapolation. The quark mass renormalization is carried out non-perturbatively using the RI-MOM method. The results for the quark masses converted to the bar{MS} scheme are: mud(2 GeV) = 3.70(17) MeV, ms(2 GeV) = 99.6(4.3) MeV and mc(mc) = 1.348(46) GeV. We obtain also the quark mass ratios ms/mud = 26.66(32) and mc/ms = 11.62(16). By studying the mass splitting between the neutral and charged kaons and using available lattice results for the electromagnetic contributions, we evaluate mu/md = 0.470(56), leading to mu = 2.36(24) MeV and md = 5.03(26) MeV.


Introduction
The precise knowledge of the quark masses and of the hadronic parameters in general plays a fundamental role both in testing the Standard Model (SM) and in the search for new physics (NP). Despite its unquestionable successes in describing experimental data the SM does not provide any explanation for the quark masses. On the theoretical side, the understanding of the hierarchical pattern of the quark masses remains an open and fascinating challenge. On the phenomenological side, since several important observables depend on the quark masses, a precise determination of these values is crucial to constrain the SM and through comparisons between theory and experiments to search for NP.
In the determination of the quark masses lattice QCD (LQCD) plays a primary role as it is a non-perturbative approach based on first principles. It consists in simulating QCD by formulating the Lagrangian on a discrete and finite Euclidean space-time which allows for a numerical computation of the path integral via MonteCarlo methods. The finite volume, the lattice spacing and generally the lower bound on the simulated light quark masses, which are limited by the currently available computing power, introduce errors which have to be well under control and accounted for.
Thanks to the increased computational power as well as to the algorithm and action improvements of the last decade, LQCD simulations have made significant progresses reaching a remarkable level of precision. In particular, this is due to the so-called unquenched calculations, where the contribution of loops of dynamical sea quarks is taken into account. As a matter of fact, most of the recent lattice determinations of quark masses have been performed with either two (up and down) [1,2] or three (up, down and strange) [3]- [10] dynamical sea quarks.
In this paper we present an accurate determination of the up, down, strange and charm quark masses using the gauge configurations produced by the European Twisted Mass (ETM) Collaboration with four flavors of dynamical quarks (N f = 2 + 1 + 1), which include in the sea, besides two light mass degenerate quarks, also the strange and charm quarks with masses close to their physical values. Such a setup is the closest one to the real world, adopted till now only by the ETM [11,12,13,14] and the MILC [15] Collaborations.
The simulations have been carried out at three different values of the inverse bare lattice coupling β, namely β = 1.90, 1.95 and 2.10, to allow for a controlled extrapolation to the continuum limit. For β = 1.90 and β = 1.95 two different lattice volumes have been considered. We also used non-perturbative renormalization constants evaluated in the RI-MOM scheme, whose calculation is discussed in Appendix A. The fermions were simulated using the Wilson Twisted Mass Action [16,17] which, at maximal twist, allows for automatic O(a)-improvement [18,19]. In order to avoid the mixing in the strange and charm sectors we adopted the non-unitary set up described in Ref. [19], in which the strange and charm valence quarks are regularized as Osterwalder-Seiler (OS) fermions [20]. For the links the Iwasaki action [21] was adopted, because it proved to relieve simulations with light quark masses allowing to bring the simulated pion mass down to approximately 210 MeV.
Since simulations were not performed at the physical point for the up and down quark masses, a chiral extrapolation is needed. In order to estimate the associated systematic error we studied the dependence on the light quark mass by using different fit formulae based on the predictions of Chiral Perturbation Theory (ChPT) as well as on polynomial expressions.
To account for finite size effects (FSE) we used the resummed asymptotic formulae developed in Ref. [22] for the pion sector, which include the effects due to the neutral and charged pion mass splitting (present in the twisted mass formulation), and the formulae of Ref. [23] for the kaon sector. We checked the accuracy of these predictions for FSE on the lattice data obtained at fixed quark masses and lattice spacings, but different lattice volumes.
As for the continuum limit, in order to lower the impact of discretization effects as much as possible and to keep the continuum extrapolation under control we tried two different procedures, which both use f π to set the scale. The first one involves the Sommer parameter r 0 [24] in units of the lattice spacing a, i.e. r 0 /a, as the intermediate scaling variable, while in the second one we used the mass of a fictitious pseudoscalar (PS) meson made of two strange-like quarks (or a strange-like and a charm-like quark), aM s s (or aM c s ), trying to exploit cancellation of discretization effects in ratios like M K /M s s (or M Ds /M c s ). In particular for the kaon and D s (D) meson masses these ratios lead to a significant reduction of discretization effects. Of course, in order to determine the lattice scale, the continuum limit of M s s (or M c s ) has to be performed eventually. The fact that we obtain compatible predictions from the two procedures strengthens the validity of our results and shows that the impact of the discretization effects is safely kept under control.
As described in Appendix A, by using dedicated ensembles of gauge configurations produced with N f = 4 degenerate flavors of sea quarks [25], we computed the quark mass renormalization constants (RCs) Z µ = 1/Z P in the RI-MOM scheme using two different methods, labelled as M1 and M2. The first method (M1) aims at removing O(a 2 p 2 ) effects, while in the second method (M2) the renormalization constants are taken at a fixed reference value of p 2 . The use of the two sets of renormalization constants is expected to lead to the same final results once the continuum limit for the physical quantity of interest is performed.
Summarizing, our analysis has followed eight branches differing in the choice of the scaling variable (either r 0 /a or aM s s ), the fitting procedures (either ChPT or polynomial expansion) and the method (either M1 or M2) used to determine the values of the RCs Z P .
First we calculated the up/down average quark mass from the analysis of the pion mass and decay constant. Then, using either r 0 or M s s (M c s ) as well as the lattice spacing and the light quark mass determined from the pion sector, we extracted the strange and charm quark masses from the analysis of K-and D-meson correlators, respectively. The differences among the results obtained within the various branches of the analysis have been used to estimate the systematic uncertainties.
where the errors are the sum in quadrature of the statistical and systematic uncertainties. By studying the light-quark mass dependence of the squared kaon mass we calculated also the leading strong isospin breaking (IB) effect on the charged and neutral kaon masses,M K 0 andM K + , which occurs in the pure QCD sector of the SM due to the quark mass difference (m d − m u ). Adopting the recent FLAG estimatê M K + −M K 0 = −6.1(4) MeV [26], we find which is independent of both the renormalization scheme and scale. Combining Eqs.
(1-2) we obtain the following predictions for the up and down quark masses: which are independent of both the renormalization scheme and scale. We also quote our results for the ratios which provide information on the relative size of SU(3) and SU(2) symmetry breaking effects.

Simulation details
The present work is based on the N f = 2 + 1 + 1 gauge field configurations generated by the ETMC [11,13] using the following action where the gluon action S g is the Iwasaki one [21]. For the fermions we have adopted the Wilson twisted-mass action, given explicitly for the mass-degenerate up/down quark doublet by [16] S tm = a 4 x ψ(x) and for the strange and charm doublet by [17] S h tm = a 4 x ψ(x) where ∇ µ and ∇ * µ are nearest-neighbor forward and backward covariant derivatives, µ is the light quark mass and m 0 is the "untwisted" mass. The latter is tuned to its critical value m cr as discussed in Ref. [11] in order to guarantee the automatic O(a)-improvement at maximal twist [18,19]. Finally in Eq. (8) the twisted masses µ σ and µ δ are related to the renormalized strange and charm sea quark masses via the relation [19] m sea c,s = 1 Z P µ σ ± Z P Z S µ δ (9) with Z P and Z S being the pseudoscalar and scalar renormalization constants, respectively.
The twisted-mass action (6) leads to a mixing in the strange and charm sectors [17,12]. In order to avoid the mixing of K-and D-meson states in the correlation functions, we adopted a non-unitary set up [19] in which the strange and charm valence quarks are regularized as Osterwalder-Seiler (OS) fermions [20]. Thus, while we keep the light sector unitary, the action in the strange and charm sectors (f = s, c) reads as where r f = ±1. When constructing meson correlation functions (including the pion) the Wilson parameters of the two valence quarks are always chosen to have opposite values. This choice guarantees that the squared PS meson mass, M 2 P S , differs from its continuum counterpart only by terms of O(a 2 µ) [18,27].
The details of our lattice set up are collected in Table 1, where the number of gauge configurations analyzed (N cf g ) corresponds to a separation of 20 trajectories. At each lattice spacing, different values of the light sea quark masses have been considered. The light valence and sea quark masses are always taken to be degenerate. The masses of both the strange and the charm sea quarks are fixed, at each β, to values close to the physical ones [11]. We have simulated three values of the valence strange quark mass and six values of the valence heavy quark mass, which are needed for the interpolation in the physical charm region as well as to extrapolate to the b-quark sector for future studies. In particular, for the light sector the quark masses were simulated in the range 3 m phys µ 12 m phys , for the strange sector in 0.7 m phys s µ s 1.2 m phys s , while for the charm sector in 0.7 m phys c µ c 2.5 m phys c . Quark propagators with different valence masses are obtained using the so-called multiple mass solver method [29,30], which allows to invert the Dirac operator for several quark masses at a relatively low computational cost.  Table 1: Values of the simulated sea and valence quark bare masses for each ensemble used in this work.
We studied the dependence of the PS meson masses and of the pion decay constant on the renormalized light quark mass fitting simultaneously the data at different lattice spacings and volumes. In particular, we anticipate that the values of the lattice spacing found in our pion analysis are a = 0.0885(36), 0.0815(30), 0.0619(18) fm at β = 1.90, 1.95 and 2.10, respectively, so that the lattice volume goes from 2 to 3 fm. In Table 2 we provide for each ensemble the central values of the pion mass (covering the range 210 ÷ 450 MeV), of the lattice size L and of the product M π L.
The statistical accuracy of the meson correlators is significantly improved by using the so-called "one-end" stochastic method [31], which includes spatial stochastic sources at a single time slice chosen randomly. Statistical errors on the meson masses ensemble β L(fm) M π (MeV) M π L A30. 32 1  Table 2: Central values of the pion mass M π , of the lattice size L and of the product M π L for the various ensembles used in this work. The values of M π are extrapolated to the continuum and infinite volume limits, according to the ChPT fit (16), described in Section 3.1.
are evaluated using the jackknife procedure, while statistical errors based on data obtained from independent ensembles of gauge configurations, like the errors of the fitting procedures, are evaluated using a bootstrap sampling with O(100) events to take properly into account cross-correlations. In Table 3 we present the values of the RCs Z P corresponding to the two methods M1 and M2, described in Section 1 (see also Appendix A), and the values of r 0 /a used to convert the data at different values of lattice spacing to the common scale given by the Sommer parameter r 0 . For each β the values of r 0 /a have been calculated at the various values of the light quark mass [11,13] and then extrapolated to the chiral limit, assuming either a linear or a quadratic dependence in aµ sea . Our results for r 0 /a are consistent within the errors with the findings of Refs. [14,32], where the extrapolation to the chiral limit was performed using only a linear dependence on aµ sea . The errors reported in Table 3 represent the sum in quadrature of the statistical uncertainty and of the systematic error associated to the two different chiral extrapolations.
Since the renormalization constants Z P and the values of r 0 /a have been evaluated using different ensembles of gauge configurations, their uncertainties have been taken into account in the fitting procedures as follows. First we generated randomly a set of  values of (r 0 /a) i and (Z P ) i for the bootstrap event i assuming gaussian distributions corresponding to the central values and the standard deviations given in Table 3. Then we added in the definition of the χ 2 the following contribution where (r 0 /a) f it i and (Z P ) f it i are free parameters of the fitting procedure for the bootstrap event i. The use of Eq. (11) allows the quantities r 0 /a and Z P to slightly change from their central values (in the given bootstrap event) with a weight in the χ 2 given by their uncertainties. This procedure corresponds to impose a gaussian prior for Z P and r 0 /a.
Before closing this section we have collected in Table 4 the time intervals (conservatively) adopted for the extraction of the PS meson masses (and of the pion decay constant) from the 2-point correlators at each β and lattice volume in the light, strange and charm sectors.

Average up and down quark mass
For each ensemble we computed the 2-point PS correlators defined as where P 5 (x) = u(x)γ 5 d(x) 1 . As it is well known at large time distances one has so that the pion mass and the matrix element Z π = | π|uγ 5 d|0 | 2 can be extracted from the exponential fit given in the r.h.s. of Eq. (13). The time intervals used for the pion case can be read off from Table 4. For maximally twisted fermions the value of Z π determines the pion decay constant without the need of the knowledge of any renormalization constant [16,18], namely Then we have studied the dependence of the pion mass and decay constant on the renormalized light quark mass through simultaneous fits based either on ChPT at next-to-leading order (NLO) or on a polynomial expansion in m . This was done following two procedures that differ for the choice of the scaling variable. In the first one we used r 0 /a, while in the second one the fictitious meson mass aM s s is adopted in order to reduce the impact of discretization effects of the PS meson masses.

Analyses in units of r 0 (analyses A and B)
Since the chiral extrapolation is an important source of uncertainty in our analysis, we have fitted the dependence of both M 2 π and f π on the renormalized light quark mass m using two different fitting functions: the one predicted by ChPT at NLO and a polynomial expansion. These two choices correspond to expanding the squared pion mass and decay constant either around the chiral point m = 0 up to higher masses including the effects of chiral logarithms, or around a non-vanishing mass m = m * down to the physical pion point without reaching the chiral limit, where non-analytic terms arise in the expansion. The ChPT approach at NLO is expected to be more accurate in the region of low m , but to suffer from possible higher order corrections at large values of m , where the polynomial expansion is expected to be more accurate.
Both solutions are in principle legitimate to perform the chiral extrapolation. Since both fits turn out to describe our lattice data nicely, the spread between the results obtained using NLO ChPT and those corresponding to the polynomial expansion represents our uncertainty on the chiral extrapolation and it will be used to estimate the corresponding systematics. This is reasonable also because the polynomial ansatz might underestimate the curvatures of f π and M 2 π /m at small values of m (as it does not contain any chiral logarithm), while the NLO ChPT fit applied to the range of our pion data (see Table 2) might overestimate the curvatures in the small m region, as suggested by the results of NNLO fits (see later Section 3.3) and indicated also by the findings of Refs. [33,34] at N f = 2 and of Refs. [35,36] at N f = 2 + 1.
Let us consider the SU(2) ChPT approach in units of r 0 which hereafter will be referred to as analysis A. The ChPT predictions at NLO can be written in the following way where P 1 -P 4 are free parameters and with B and f being the SU(2) low-energy constants (LECs) entering the LO chiral Lagrangian, which have been left free to vary in our fits. In Eqs. (16-17) the parameters P 1 and P 3 are related to the NLO LECs 3 and 4 by with M phys π being the value of the pion mass at the physical point, while the quantities K F SE M 2 and K F SE f represent the finite size effects (FSE) for the squared pion mass and the pion decay constant, respectively. They will be discussed in a while.
For the moment notice the presence of the terms proportional to a 2 log ξ in Eqs. (16)(17). These terms originate from the mass splitting between the charged and the neutral pions, which is a discretization effect appearing within the twisted mass formulation. Its impact on the ChPT expansion of M 2 π and f π (see Ref. [37] and references therein) has been worked out in Ref. [28], where a power counting scheme was adopted in which a 2 Λ 4 QCD ≈ 2Bm . We have expanded the resulting formulae up to O(a 2 ), leading to Eqs. (16)(17) with the presence of the parameter c 2 which is directly related to the neutral and charged pion mass splitting at LO by In the χ 2 -minimization procedure we have given to c 2 a prior based on the values found in Ref. [32] by analyzing charged and neutral pion data for a set of ETMC ensembles consistent with the one considered in this work 2 . In Ref. [32] two different determinations of c 2 are reported, one in which the chiral limit is performed through a constant fit in M 2 π and the other one in which the fit was assumed to be linear. In the present work we have used an average of the two determinations including the spread in the error, which in units of r 0 reads as c 2 r 4 0 = −1.7 ± 0.6. On the theoretical side the impact of FSE on M π and f π has been studied within ChPT at NLO in Ref. [38] and using a resummed asymptotic formula in Ref. [23], where both leading and sub-leading exponential terms are taken into account and the chiral expansion is applied to the π−π forward scattering amplitude. When the leading chiral representation of the latter is considered, the resummed approach coincides with the NLO result of Ref. [38]. Viceversa at NNLO the resummation technique includes only part of the two-loop effects as well as of higher-loop effects. The resummed approach was positively checked against a full NNLO calculation of the pion mass in Ref. [39], showing that the missing two-loop contributions are actually negligible for M π L ∼ > 2 and L ∼ > 2 fm. Finally, we considered that the cutoff effects, giving rise to the splitting between charged and neutral pions, enter also the determination of FSE, as explicitly worked out within the resummed approach in Ref. [22].
Thus, as far as FSE are concerned, we have investigated three different approaches: the NLO ChPT predictions of Ref. [38] (which will be labelled hereafter as GL), the resummed formulae of Ref. [23] including higher order corrections (labelled as CDH) and the formulae developed in Ref. [22] which accounts for the π 0 − π + mass splitting (labelled as CWW).
The predictions of both CDH and CWW approaches require the knowledge of the LECs 1 − 4 and eventually of the splitting parameter c 2 . The LECs 3 and 4 , which are related to the ξ -dependent NLO terms in M 2 π and f π [see Eqs. (16)(17)(18)(19)], have been treated as free parameters in our fitting procedures, while for 1 and 2 we used the values given in Ref. [22]. The CWW corrections depend also on the neutral pion mass M π 0 , which was estimated at LO through Eq. (20) using (M π + ) LO = 2Bm . We have checked that such values of M π 0 are consistent with those extracted directly from the neutral PS correlator in Refs. [14,32].
In order to check how well the finite volume corrections predicted by the three chosen approaches are working, we have used the two ensembles A40.32 and A40.24 (see Table 1), which correspond to the same quark mass and lattice spacing, but different lattice volumes. Notice that the ensemble A40.24 has both the lowest value of the quantity M π L (see Table 2) and the largest pion mass splitting, being M π 0 /M + π ≈ 0.5 [14,32]. Therefore FSE are expected to be maximal for this ensemble. The and in analogous way for K F SE f, [32] and K F SE f, [24] in the case of the decay constant f π . Taking the ratio of the above relations we see that for an ideal correction the ratio of the multiplicative factors K F SE should match the ratio of the uncorrected lattice data independently of the infinite volume values. The more accurate the correction is, the more the prediction for (K F SE M 2 , [32] /K F SE M 2 , [24] ) matches the lattice data (M 2 [32] /M 2 [24] ). The corresponding numerical results are reported in Tables 5 and 6   for the ensembles A40.32 and A40.24, obtained within the approaches GL, CDH and CW W (see text), compared with the corresponding ratio of lattice data.
From these tables one can see that the corrections calculated using the CWW approach are well compatible with the lattice data for both the pion mass and the decay constant. It is also possible to see how large the relative contribution of the various FSE corrections is.
In table 7 we collected the values of the coefficients (K F SE M 2 , [24] − 1) and (K F SE f, [24] − 1), representing the FSE correction for the ensemble A40.24, which, as already noted, is affected by the largest FSE correction in the whole set of ensembles. By comparing CDH and CWW predictions it can also be seen that the O(a 2 ) term related to the pion mass splitting, though not negligible, is not the dominant one and appears to be at the percent level. In what follows, the pion data will be corrected for FSE using the CWW formulae unless explicitly stated otherwise.
The dependence of our lattice data for r 0 M 2 π /m and r 0 f π on the renormalized quark mass r 0 m is shown in Figs. 1 and 2  From Figs. 1 and 2 it can be seen that the impact of discretization effects using the values of r 0 /a is almost completely negligible in the case of r 0 f π , while it is at the level of 10% in the case of r 0 M 2 π /m (using the difference between the continuum results and the ones at the finest lattice spacing).
The value of the physical average up/down quark mass, m ud , can be extracted from the ratio M 2 π /f 2 π using as input its experimental value, obtained from the central values of Ref. [40] (see also Ref. [26]) The numerical results for m ud as well as those for the lattice spacing and the relevant LECs will be collected and discussed in Section 3.3.
As anticipated in Section 1, we studied the chiral extrapolation also by replacing the NLO ChPT ansatz with a simple polynomial expansion in the renormalized light quark mass, namely where B, f and P 1 -P 6 are free parameters. This analysis will be referred to as analysis B. Since the calculation of K F SE M 2 and K F SE f is based on ChPT, the FSE corrections have been taken from the analysis A and applied directly to the lattice data.
The chiral extrapolations of our lattice data for r 0 M 2 π /m and r 0 f π , obtained using the polynomial fits (23)(24), are shown for each lattice spacing and in the continuum limit in Figs. 3 and 4, respectively.
Notice that the impact of discretization effects on r 0 M 2 π /m obtained using the polynomial fit (see Fig. 3) is very similar to the one found in the case of the NLO ChPT prediction (see Fig. 1), while in the case of r 0 f π , at variance with the NLO ChPT fit (see Fig. 2), the polynomial expansion exhibits visible cutoff effects (see Fig. 4) though limited at the level of few percent only. Nevertheless, both the NLO ChPT and the polynomial fits describe quite well the lattice data for the pion mass and the decay constant, yielding only slightly different results, at the percent level, at the physical pion point.
Continuum Limit Figure 1: Chiral and continuum extrapolation of r 0 M 2 π /m based on the NLO ChPT fit given by Eq. (16). Lattice data have been corrected for FSE using the CWW approach [22] and correspond to the RCs Z P calculated with the method M1 (see text).
Continuum Limit Continuum Limit Continuum Limit Figure 4: The same as in Fig. 3, but for the decay constant r 0 f π .

Analyses in units of M s s (analyses C and D)
The results presented in Figs. 1 and 3 show that the impact of discretization effects using r 0 as the scaling variable is at the level of 10% for the squared pion mass.
In order to keep the extrapolation to the continuum limit under better control we repeated the analyses A and B adopting a different choice for the scaling variable, namely instead of r 0 we introduced the mass M s s of a fictitious PS meson made of two strange-like valence quarks 3 . The PS mass M s s has a very mild dependence on the light-quark mass and is affected by cutoff effects similar to the ones of a K meson. Thus, we tried to improve the continuum extrapolation by considering the ratio M 2 π /M 2 s s which may exploit a partial cancellation of discretization effects. To construct the meson mass ratio we first performed a slight interpolation in the strange valence quark mass to get the quantity aM s s at a common (but arbitrary) value r 0 m s = 0.22 for each β and light quark mass. Since, as expected, we found no significant dependence of aM s s on the light quark mass, we performed a constant fit in aµ to obtain the values of aM s s at each β. In this way we find The values of aM s s have been used to bring to a common scale all lattice quantities, covering the role that in analysis A and B was played by r 0 /a. The (quite small) errors on aM s s are propagated via the bootstrap sampling.
The new analyses, which will be referred to as analyses C and D, proceed in the same way as in the previous Section, namely in the case of the NLO ChPT fit (analysis C) one employs the ansatz where again the parameters P 1 and P 3 are related to the NLO LECs 3 and 4 through Eq. (19). In the case of the polynomial fit (analysis D) one fits the data with the analogue of Eqs. (23) and (24)  The comparison of Figs. 1 and 5 clearly shows that, when M s s is chosen as the scaling variable, the discretization effects on the squared pion mass are significantly reduced from 10% down to 4.5%. At the same time the discretization effects on the pion decay constant, which are almost negligible in units of r 0 (see Fig. 2), are kept to be within 4% when M s s is used as the scaling variable (see Fig. 6).
Continuum Limit Continuum Limit Figure 6: The same as in Fig. 5, but for the pion decay constant f π in units of M s s .

Results for the pion sector
In this section we present the results of the four analyses (A, B, C, D) carried out in the pion sector. We have adopted the values of the RCs Z P corresponding to the methods M1 and M2, so that we end up with eight analyses, which will be referred to as analyses A1, B1, C1, D1 and A2, B2, C2, D2, respectively. Using the experimental value of the ratio M 2 π /f 2 π [see Eq. (22)], the average up/down quark mass m ud is determined, so that the quantity (r 0 f π ) is calculated at the physical point within the analyses A1 (A2) and B1 (B2). Then, using the experimental value of f π as input, the Sommer parameter r 0 is extracted and this in turn allows to get the values of the lattice spacing at each β using the determinations of r 0 /a collected in Table 3.
The analyses C1 (C2) and D1 (D2) proceed in the same way: the average up/down quark mass m ud is determined through the experimental value of the ratio M 2 π /f 2 π , while the mass M s s is obtained by combining the value of f π /M s s , calculated at the physical point, and the experimental value of f π . However, in order to determine the lattice spacing at each β we did not use the quantities aM s s given in Eq. (25), since they are affected by discretization effects larger than those occurring in the values of r 0 /a. Thus we proceeded as follows. First we converted the results (25) for aM s s to r 0 M s s using the values of r 0 /a from Table 3, and then we performed a simple fit of the form r 0 M s s = P 1 + P 2 a 2 /r 2 0 . Finally, we determined the values of the lattice spacing at each β by combining the values of a/r 0 with the continuum extrapolation of r 0 M s s and the value of M s s obtained from the experimental value of f π .
For convenience the results obtained for the quark mass m ud , the scaling variables r 0 and M s s , the values of the lattice spacing and the LECs B, f , 3 and 4 , are collected in Tables 8 and 9.
It is quite reassuring to find that different ways of handling both the chiral extrapolation and the discretization effects produce consistent results.
Combining the results reported in Tables 8 and 9 provides us with the final determinations and the estimates of the various sources of systematic uncertainties. For each quantity we have a set of N results (where N = 4 or N = 8 depending on the specific quantity) coming from the various analyses A1 -D2. We assign to all analyses the same weight and therefore we assume that the observable x has a distribution f (x) given by f ( is the distribution provided by the bootstrap sample of the i-th analysis and characterized by central value x i and standard deviation σ i . Thus we estimate the central value and the error for the observable x through the mean value and the standard deviation of the distribution f (x), which are given by The second term in the r.h.s. of Eq. (28), coming from the spread among the results of the different analyses, corresponds to a systematic error which accounts for the uncertainties due to the chiral extrapolation, the cutoff effects and the RCs Z P . Finally we (13) 3.87 (17) 3.66 (10) 3.75 (13) r 0 (GeV −1 ) 2.39 (6) 2.42 (7) -- (12) 0.477 (14) -- (10) a(β = 1.90)(fm) 0.0886 (27) (12) 3.78 (16) 3.55 (9) 3.63 (12) r 0 (GeV −1 ) 2.40 (6) 2.42 (7) -- The first error includes the statistical one as well as the error associated with the fitting procedure. This error is larger than the typical statistical error of the lattice data, being amplified by the chiral and continuum extrapolations. For m ud we get a (stat+fit) error equal to 3.5%.
In order to separate in Eq. (29) the uncertainties related to the chiral extrapolation, the discretization effects and the choice of the RCs Z P we split the contribution coming from the second term in the r.h.s. of Eq. (28) into those related to the differences of the results obtained using r 0 or M s s (labelled as Disc), chiral or polynomial fits (labelled as Chiral) and the two methods M1 and M2 for the RCs Z P (labelled as Z P ). In this way we found them to be at the level of 1.6%, 1.6% and 1.4%, respectively.
For the FSE we considered the difference between the result obtained using the most accurate correction, i.e. the CWW one, and the one corresponding to no FSE correction at all. This gave rise to an error on m ud equal to 1.1%.
Our determination (29) for m ud is the first one obtained at N f = 2 + 1 + 1. The recent lattice averages, provided by FLAG [26] and based on the findings of Refs. [1,4,5,41,42], are: m ud = 3.6(2) MeV at N f = 2 and m ud = 3.42 (9) MeV at N f = 2 + 1. The comparison of these results with our finding (29) shows that the partial quenching of the strange and/or charm sea quarks is not yet visible at the (few percent) level of the present total systematic uncertainty.
As it is known (see the findings of Refs. [33,34] at N f = 2 and of Refs. [35,36] at N f = 2 + 1), a precise determination of the NLO LECs 3 and 4 requires refined analyses addressing the impact of the choice of pion mass range used for the chiral extrapolation as well as the effects of NNLO corrections. Such analyses are beyond the scope of the present work. Here we mention only that we have performed NNLO fits in the whole mass range covered by our data (M π < 450 MeV) as well as NLO fits restricted to pion masses smaller than 300, 350 or 400 MeV. The results of these fits (see Table 10) indicate that the curvatures of M 2 π /m and f π are within the range already selected by the polynomial and the NLO ChPT fits performed in the full range of simulated pion masses. In particular, for the average up/down quark mass 3.72 (13) 3.87 (17) 3.77 (21)   m ud , whose determination is one of the main goals of the present work, and for the LECs B and f , we have found results always in between those obtained with the polynomial and the NLO ChPT fits.
It is interesting to show in detail the impact of the various approaches used to calculate the FSE for the various quantities extracted from the pion analysis. The results obtained within the eight analyses A1 -D2 are quite similar to each other. In Table 11 we have reported the findings corresponding to the analysis A1.
3.68 (14) 3.76 (14) 3.73 (13)   From Table 11 it can be seen that, though the FSE corrections in some particular ensemble can be as large as 4.9% and 6.3% for the pion mass and decay constant, respectively (see Table 7 for the ensemble A40.24), the overall final impact on m ud , r 0 and the LECs B, f , 3 and 4 is limited to be well below the (stat+fit) error.
Before closing this Section, we notice that a set of ETMC data consistent with the ones considered in this work have been analyzed in Ref. [13] adopting ChPT at NLO for the chiral extrapolation, but without accounting for the effect of the charged/neutral pion mass splitting and without involving the determinations of the RCs Z P . The findings of Ref. [13] concerning both the lattice spacings and the LECs B, f , 3 and 4 nicely agree with our results of Tables 8-9 within one standard deviation. This indicates that the role played in our analyses by the pion mass splitting and by the RCs Z P is well under control.

Strange quark mass
In this Section we present our determination of the strange quark mass m s . The analysis follows a strategy similar to the one presented for the pion sector. As a preliminary step, however, we performed an interpolation of the lattice kaon data to a fixed value of the strange quark mass in order to arrive iteratively at the physical one (see next Section).
As in the pion sector, we handled discretization effects by performing a first analysis which uses r 0 /a as scaling variable, and a second one in which the fictitious PS meson mass aM s s is used to build the ratios M K /M s s , which are expected to have milder lattice artifacts. For both approaches we considered two different chiral extrapolations in the light quark mass m , namely either the predictions of SU(2) ChPT or the polynomial expansion. All these analyses are then repeated with the two sets of values of the RCs Z P obtained within the methods M1 or M2. In this way, as in the pion sector, there are eight different branches of the analysis. In all cases the quark masses are converted directly to physical units using the values of the lattice spacing found in the pion sector.
To determine the strange quark mass we made use of several quantities extracted from the pion sector, like the lattice spacing, the LECs B and f , the Sommer parameter r 0 and the results for the average up/down quark mass m ud . In order to preserve the physical correlations, in each of the eight kaon analyses we adopted the inputs coming from the corresponding pion fit. For instance, if SU(2) ChPT is used for the pion, then the same approach is applied to the kaon as well. The uncertainties on the input quantities are propagated through the bootstrap sampling for each of the branches of the kaon analysis. Combining the results from all the eight analyses we obtained our final result for m s and the estimates of the various sources of systematic uncertainty.

Chiral extrapolation in units of r 0 (analyses A and B)
The analysis is performed iteratively. We start from an initial guess for the physical strange quark mass m s . Then, adopting a quadratic spline, the lattice data for the kaon masses are interpolated in the strange quark mass to the (guessed) physical value m s and brought to a common scale using r 0 /a. A combined fit is performed to extrapolate M 2 K in the light quark mass and in the (squared) lattice spacing to the physical point and to the continuum limit. Afterwards the value obtained for the kaon mass, converted in physical units using the value of r 0 obtained from the pion analyses, is compared with the experimental one. If the latter is not reproduced, a new guess for m s is done and the whole process is repeated again.
The experimental value of the kaon mass to be matched is the one in pure QCD corrected for leading strong and electromagnetic isospin breaking effects according to where ε = 0.7(3), ε K 0 = 0.3(3) and ε m = 0.04(2) [26].
For the analysis A we used the SU(2) ChPT predictions at NLO, which assume the chiral symmetry to be satisfied only by the up and down quarks and read as Alternatively we considered a polynomial fit (analysis B) according to the following expression Notice that for the squared kaon mass SU (2) ChPT predicts the absence of chiral logarithms at NLO, so that the expressions (33) and (34) actually correspond to a linear and a quadratic fit in m , respectively. The data for the kaon mass have been corrected for FSE using ChPT formulae. The absence of the chiral log at NLO makes the corresponding FSE correction (GL) vanishing identically, i.e. K F SE M 2 K = 1. The first non-vanishing correction appears at NNLO and it was calculated in Ref. [23]. The pion mass splitting is expected to give a contribution to the FSE also for the kaon mass. However explicit calculations are not yet available 4 . In Table 12 the relative size of the FSE correction for the kaon mass is presented, together with a comparison to the lattice data. It can clearly be seen that: i) FSE on the kaon mass are definitely smaller compared to the pion case (see Table  5), and ii) even if the contribution from the pion mass splitting is neglected, the CDH predictions appear to work quite well, reproducing the observed ratio of lattice data.     Continuum Limit Figure 8: The same as in Fig. 7, but in the case of the polynomial fit (34).
In both cases the lattice data are reproduced quite well by the fitting formulae.
Notice the size of discretization effects, which can be quantified at the level of 10% taking the difference between the results at the finest lattice spacing and the ones in the continuum limit.

Chiral extrapolations in units of M s s (analyses C and D)
Following the same strategy adopted in the pion analyses, the kaon masses simulated at different β values can be brought to a common scale by constructing the ratios M 2 K /M 2 s s , which are expected to suffer only marginally by discretization effects. The values of aM s s for each β are given in Eq. (25). The light quark mass m is expressed directly in physical units by using the values of the lattice spacing found in the corresponding pion analysis.
As for the analyses in units of r 0 , we considered two different chiral extrapolations, adopting formulae similar to Eqs. (33) and (34), but expressed in units of M s s . After the chiral extrapolation and the continuum limit are carried out, the result for M K /M s s can be combined with the value of M s s obtained in the corresponding pion analysis in order to compare with the experimental kaon mass (32).
The dependencies of M 2 K /M 2 s s on the renormalized light quark mass at the three values of β as well as in the continuum limit are shown in Fig. 9 using the SU(2) ChPT prediction (analysis C). Results of the same quality are obtained within the analysis D, which makes use of the polynomial fit for the chiral extrapolation.
In the case of the kaon mass the use of the hadron scale M s s turns out to be an extremely efficient choice for an almost total cancellation of the discretization effects, namely from 10% (see Figs. 7 and 8) to about 0.4% (see Fig. 9). This allows us to keep the extrapolation to the continuum limit under a very good control in the whole range of values of the renormalized light quark mass.

Results for the kaon sector
Our results for the strange quark mass m s are those reproducing after the chiral and continuum extrapolations the experimental value of the K-meson mass given in Eq. (32). The results of the eight analyses for the strange quark mass, given in the MS scheme at a renormalization scale of 2 GeV, are shown in Table 13.
The chiral extrapolation error has been evaluated from the spread among the results obtained using the chiral and the polynomial fits in units of either r 0 or M s s . This corresponds in the error budget to a 0.6% systematic uncertainty. The discretization error has been calculated from the spread among the results obtained in units of r 0 or M s s and represents a 1.1% uncertainty on m s .
The two different sets of values of Z P , calculated using the methods M1 and M2, introduce an additional uncertainty of 1.4%.
The difference of the results for the strange quark mass obtained without correcting for the FSE and the one obtained using the CDH approach [23] has been conservatively taken as the estimate of the corresponding systematic uncertainty, which turns out to be equal to 0.5%.
The largest uncertainty, equal to 3.6%, comes from the statistical error plus the uncertainties due to the fitting procedure. The latter is the dominating one and it mainly depends on the distance between the lowest simulated quark mass and the physical point m ud in the chiral extrapolation.
Our determination (35) for m s is the first one obtained at N f = 2 + 1 + 1. The recent lattice averages, provided by FLAG [26] and based on the findings of Refs. [1,2,5,41,42,44], are: m s = 101 (3) MeV at N f = 2 and m s = 93.8 (2.4) MeV at N f = 2+1. The comparison of these results with our finding (35) shows that the partial quenching of the strange and/or charm sea quarks is not yet visible at the (few percent) level of the present total systematic uncertainty.

The ratio m u /m d
The light quark mass dependence of the squared kaon mass can be used to calculate the mass difference between the u and d quark masses, leading to an estimate of the ratio m u /m d . In the limit of vanishing electromagnetic interactions the difference between the neutral and charged squared kaon masses can be expanded in terms of the (small) quark mass difference (m d − m u ) as (see Ref. [45] and references therein) The slope (∂M 2 K /∂m ) m =m ud is defined in the isospin symmetric limit and therefore it can be computed directly using our ensembles by taking the derivative of the continuum and infinite volume limits of our fitting formulae, like Eqs. (33)(34), with respect to m ,obtaining We observe that in Ref. [45], using a different method based on the insertion of the isovector scalar density, the slope was found to be equal to (∂M 2 K /∂m ) m =m ud = 2.57 (8) GeV at N f = 2.
The charged and neutral kaon masses,M K 0 andM K + , are those defined in pure QCD. For them we adopt the FLAG estimateM K + −M K 0 = −6.1(4) MeV [26]

Determinations of the strange and charm sea quark masses
As discussed in Section 2, within the twisted mass formulation adopted in the present work the (renormalized) strange and charm sea quark mass are related to the bare twisted parameters µ σ and µ δ by Using the results found for the RCs Z P and Z P /Z S (see Appendix A for the latter), it turns out that the values of m sea s obtained from Eq. (41) are plagued by large uncertainties that can reach the 20% level, mainly because of a large cancellation between the two terms in the r.h.s. of Eq. (41). Moreover, the definition (41) is affected by the lattice artifacts that unavoidably enter the determination of the RCs.
A more accurate determination of m sea s can be obtained using the results of Refs. [11,12,13,14,46], where for all the ensembles used in the present work the kaon mass has been determined in the twisted-mass unitary setup, in which the valence quarks are described by the same action (8) adopted for the sea quarks.
In terms of the valence (m and m s ) and strange sea (m sea s ) quark masses the OS kaon masses, computed in the present study, can be represented as M OS .
This ratio is equal to the ratio m sea s /m s in ChPT at LO and it is equal to unity when m sea s = m s up to lattice artifacts corresponding to the difference of the discretization effects in the unitary and OS setups. Therefore, for each ensemble a smooth local interpolation (carried out with quadratic splines) allows us to find the value of the valence strange quark mass m s that makes the ratio R sea (m , m s , m sea s ) equal to unity. The results of the above procedure are shown in Fig. 10, where it can be seen that the matching mass can be determined with good precision and it is almost independent on the values of the light quark mass for fixed β.
where each error includes also the spread of the matching mass with respect to the light quark mass (see Fig. 10). The results (44) differ from the determination (35) of the physical strange quark mass by ≈ 10% at most, with the largest difference at β = 1.95. We tried to estimate the effect of the mistuning of the strange sea quark mass using the SU(3) ChPT predictions developed in Refs. [47]- [49] for arbitrary values of sea and valence quark masses. For the squared pion and kaon masses one gets at NLO where and B 0 and f 0 are the LO SU(3) LECs, while L r 4 (µ) and L r 6 (µ) are the NLO LECs evaluated at the renormalization scale µ. For the pion decay constant one gets Using from the results quoted in Ref. [26] the values B 0 /f 0 = 19 (2) and 2L r 6 (µ) − L r 4 (µ) = 0.14 (12) · 10 −3 , L r 4 (µ) = 0.09 (34) · 10 −3 at µ = M ρ = 0.770 GeV, the corrections (45), (46) and (48) are below the 1% level at our simulated quark masses and at the physical point.
We have also verified that by including the corrections (45), (46) and (48) in the lattice data the changes observed in the predictions of our analyses for m ud and m s are smaller than the other systematic uncertainties.
We close this subsection by presenting the estimate of the charm sea quark mass  Table 14). It follows that, while there is a good agreement within the errors at β = 1.90 and 1.95, a ≈ 18% mistuning is present at β = 2.10. Since scaling distortions are not visible in our data, we expect that the mistuning of the charm sea quark mass has a negligible effect with respect to the one of the strange sea quark and, therefore, it does not affect our determination of the quark masses in a significant way.

Determination of the ratio m s /m ud
The results for the strange quark mass m s and for the average up/down quark mass m ud (see Tables 8, 9 with a total uncertainty of 5.3%. In order to reduce the uncertainty we have investigated an alternative approach, which leads to a more precise determination of the ratio m s /m ud . Using the kaon and pion lattice data we define the quantity R(m s , m , a 2 ) as which, by construction, is independent on the values of Z P as well as of the lattice spacing up to cutoff effects. In ChPT at LO the ratio R(m s , m , a 2 ) is equal to unity. At the physical point one gets [(2M 2 K − M 2 π )/M 2 π ] phys 25.8 and, adopting the estimate (50) for (m s /m ud ), one has R phys ≡ R(m s , m ud , 0) 0.96. Therefore, the dependence of R(m s , m , a 2 ) on the strange and light quark masses is expected to give rise to small corrections only. This is a very useful feature, since the mild dependence on the strange quark mass allow us to interpolate the ratio R(m s , m , a 2 ) at the physical value (35) with a small sensitivity to the error on m s , while the mild dependence on the light quark mass m represents a way to reduce the uncertainty introduced by the chiral extrapolation. In this way a precise determination of the mass ratio m s /m ud can be obtained as where R phys is computed on the lattice. In Fig. 11 the lattice data for R(m s , m , a 2 ), interpolated at the physical strange mass (35) and corrected for FSE (using the CWW predictions [22] for M π and the CDH ones [23] for M K ), are shown versus the light quark mass m for our ensembles. As expected, the dependence on the light quark mass is found to be quite mild and the ratio R(m s , m , a 2 ) is close to unity at all the simulated quark masses. 10   We performed the chiral and continuum extrapolations through a simple fit of the form R(m s , m , a 2 ) = R 0 + R 1 m + R 3 a 2 .
The results are presented in Fig. 11 for each β value and in the continuum limit. It can be seen that discretization effects are quite small, being the difference between the result at the finest lattice spacing and the one in the continuum less than 1%.
which has an accuracy at the level of 1.2%. For comparison, the updated FLAG averages [26] are m s /m ud = 28.1 (1.2) at N f = 2 and m s /m ud = 27.5(4) at N f = 2 + 1.

Charm quark mass
In this Section we present our determination of the mass of the charm quark obtained by analyzing both the D-and D s -meson masses, following a strategy similar to the one presented for the K-meson.
The lattice data for the D-and D s -meson masses are interpolated to the physical strange and charm quark masses using a quadratic spline. The physical strange quark mass is the one determined in the previous Section, while the physical charm quark mass is defined such that the experimental value of the D-or D s -meson mass is reproduced. Then the dependence of M D and M Ds on the light quark mass and on the lattice spacing is studied at fixed strange and charm quark masses, and the continuum limit and the chiral extrapolation to the physical point m ud of the light quark mass are performed. The analysis based on the D s -meson masses is expected to have smaller systematic uncertainty associated to the chiral extrapolation because of the milder light quark dependence, which occurs only through the sea effects. Therefore, our final result for the charm quark mass is derived from the D s -meson analysis and the value obtained from fitting the D-meson mass is used as a consistency check.
As For the chiral extrapolation in the light quark mass, the Heavy Meson ChPT (HMChPT) predicts no chiral logarithms at NLO for both D-and D s -meson masses. Therefore, we have adopted either a linear or a quadratic expansion in m and the latter is used only to estimate the uncertainty related to the chiral extrapolation.

Fit in units of r 0
Our analyses follow closely the strategy already applied to the kaon case. We start from an initial guess for the physical charm quark mass m c and consider the value of the physical strange quark mass m s given in (35). After a smooth interpolation in the strange and charm quark masses, the D-and D s -meson masses, extracted from the corresponding 2-point correlators, are brought to a common scale using r 0 /a. The light quark mass is directly converted in physical units using the values of the lattice spacing obtained in the pion sector.
As discussed in the previous Section, the dependence of both r 0 M D and r 0 M Ds on the light quark mass m is well described by a simple polynomial dependence, namely where P 0 -P 3 and P 0 -P 3 are free parameters. For both D and D s mesons we have investigated either a linear (i.e. with P 2 = P 2 = 0 in Eqs. (55)(56)) or a quadratic fit. As in the previous analyses, the prior information on Z P and r 0 /a is introduced through the contribution to the χ 2 given in Eq. (11). Moreover, since the results for the D-and D s -meson masses corresponding to the ensembles A40.24 and A40.32 (which differ only for the lattice volume) almost coincide, we did not apply any FSE correction.
The dependence of M Ds on the light quark mass m for each β value and in the continuum limit is illustrated in Figs. 12-13, adopting a linear or a quadratic fit, respectively. It can be seen that the discretization effects, which can be quantified by the difference between the results at the finest lattice spacing and those in the continuum limit, are found to be of the order of 3%.

Fit in units of M c s
The impact of discretization effects can be reduced using the reference meson mass      Continuum Limit Figure 13: The same as in Fig. 12, but in case of the quadratic fit of Eq. (56).      The results for the charm quark mass in the MS(2 GeV) scheme, obtained form the D s -meson analysis, are shown in Table 14. Each entry in the Table is already the average value evaluated according to Eq. (28) of the results of analyses which differ only for the choice of the set of the input parameters: those coming from pion and kaon analysis A (chiral extrapolation) and B (polynomial extrapolation) when r 0 is used as scaling variable, and those coming from the analyses C and D when M c s is considered.

Results for the charm mass
It is interesting to compare the results for m c obtained analyzing the D s -meson mass with those obtained using the D-meson mass. The latter are presented in Table  15. It can be clearly seen that there is indeed a full compatibility, the differences being much smaller than the quoted uncertainties.    Table 14, but using the data for the D-meson mass.
The quality of the chiral and continuum extrapolation performed on the D-meson mass is illustrated in Fig. 16 in the case of the quadratic fit in m .
The results from Table 14 corresponding to the linear fit in the light quark mass have been averaged to get our final result for m c and its systematic uncertainties. The results of the quadratic fit in m have not been included in the average, but they have been considered to estimate the uncertainty related to the chiral extrapolation with a total uncertainty equal to 3.1% of the central value.
The strategy followed to separate the various sources of the systematic error is the same as the one used in the pion and kaon cases.
The first error in Eq. (60) includes the statistical uncertainties combined with the systematic error associated with the fitting procedure, the physical strange quark mass and the scale setting. This error is the dominant one and corresponds to a 2.7% of the central value.
The systematic uncertainty due the chiral extrapolation, estimated from the difference between the results of the linear and quadratic fits, is equal to 0.15%.
The difference among the results obtained using r 0 or M c s is used to estimate the uncertainty coming from the discretization effects, which results to be of the order of 0.45%.
The effect of choosing the values of Z P obtained either from the method M1 or M2 gives rise to a systematic uncertainty of 1.5%.
Our determination (60) for m c (m c ) is the first one obtained at N f = 2 + 1 + 1 and it is consistent with the result m c (m c ) = 1.28(4) GeV obtained in Ref. [1] at N f = 2, with the finding m c (m c ) = 1.273(6) GeV of Ref. [9] at N f = 2 + 1 as well as with the PDG value m c (m c ) = 1.275(25) GeV [40].

Determination of the ratio m c /m s
The results for the strange and charm quark masses given in Tables 13 and 14 with a total uncertainty of 4.7%.
In order to improve the precision of this determination we followed an approach similar to the one used in the case of the mass ratio m s /m ud discussed in Section 4.6. Using the lattice data for the masses of the η c and D s mesons, we define the quantity R(m c , m s , m , a 2 ) as which by construction is independent of the values of Z P and of the lattice spacing up to cutoff effects. In Eq. (62) the mass M ηc of the η c meson corresponds to the (fermionic) connected diagram only, or in other words it is the mass of a fictitious cc PS meson with mc = m c and rc = −r c . Let us explain the choice of the ratio (62). For a PS meson made of two valence quarks with (renormalized) masses m 1 and m 2 , in which one of the two quarks is around the charm mass, the meson mass M 12 can be written up to cutoff effects as where the function r(m 1 , m 2 ) includes higher order contributions in the quark masses. Therefore, up to cutoff effects, the ratio R(m c , m s , m , a 2 ) has a leading term, which is a constant and receives corrections only from terms in m c appearing in Eq. (63) at orders higher than the linear one 5 .
As in the case of the ratio R(m s , m , a 2 ) defined in Section 4.6, the useful features of R(m c , m s , m , a 2 ) are that: i) its interpolation at the physical charm and strange quark masses is only slightly sensitive to the uncertainties on m c and m s , and ii) its dependence on the light quark mass m is expected to be mild, so that the uncertainty introduced by the chiral extrapolation is largely reduced. A determination of the mass ratio m c /m s is then obtained from where R phys ≡ R(m c , m s , m ud , 0) is computed from lattice data.
In Fig. 17 and the results are shown in Fig. 17 as the solid lines at each β value and in the continuum. It can be seen that the dependence on the light quark mass is very mild, allowing to get a precise chiral extrapolation to the physical point, namely R phys = 0.1772 (24) stat+f it (2) Z P .
From the PDG [40] one gets: M ηc = 2.9837(7) GeV and M D ± s = 1.9690 (14) GeV. The disconnected contribution to the physical η c meson, which is neglected in the present calculation, can be estimated from the annihilation rate into gluons, leading to an estimate of 2.5 MeV (see Ref. [8] and references therein). Assuming a 50% error on the latter, the "connected" η c mass to be used in Eq.  (40) [50] at N f = 2, and m c /m s = 11.85 (16) [8] at N f = 2 + 1.

Conclusions
We have presented results for the up, down, strange and charm quark masses, obtained with N f = 2 + 1 + 1 twisted-mass Wilson fermions. We have used the gauge configurations produced by the ETMC, which include in the sea, besides two light mass degenerate quarks, also the strange and the charm quarks with masses close to their physical values. Such a setup is the closest one to the real world, adopted till now only by the ETM [11,12,13,14] and the MILC [15] Collaborations.
The analysis includes data at three values of the lattice spacing and pion masses in the range 210 ÷ 450 MeV, allowing for accurate continuum limit and controlled chiral extrapolation. In order to estimate the systematic error associated with the chiral extrapolation we studied the dependence on the light quark mass by using different fitting formulae based either on the predictions of ChPT or on polynomial expressions.
As for the continuum limit, in order to lower as much as possible the impact of discretization effects and to keep the continuum extrapolation under control we investigated two different procedures, which both use f π to set the scale. The first one involves the Sommer parameter r 0 as the intermediate scaling variable, while in the second one we used the mass of a fictitious pseudoscalar meson made of two strangelike quarks (or a strange-like and a charm-like quark), M s s (or M c s ), trying to exploit cancellation of discretization effects in ratios like M K /M s s (or M Ds /M c s ). For the kaon and D s (D) meson masses these ratios really lead to a significant reduction of discretization effects.
To account for FSE we used the resummed asymptotic formulae developed in Ref. [22] for the pion sector, which include the effects due to the neutral and charged pion mass splitting (present in the twisted mass formulation), and the formulae of Ref. [23] for the kaon sector. We checked the accuracy of these predictions for FSE on the lattice data obtained at fixed quark masses and lattice spacings, but different lattice sizes.
The quark mass renormalization has been carried out non-perturbatively using the RI-MOM method, adopting dedicated ensembles of gauge configurations produced by ETMC with N f = 4 degenerate flavors of sea quarks.
The main results obtained in this paper for the up, down, strange and charm quark masses and for some important quark mass ratios have been collected in Section 1, see Eqs.

A Renormalization constants
In order to obtain results for the quark mass m f (f = u, d, s, c) in the MS scheme at a given renormalization scale, chosen to be 2 GeV in the present study, a necessary step is the evaluation of the quark mass renormalization constant (RC) in a suitable intermediate lattice renormalization scheme, which here we take to be the RI-MOM scheme [51].
In the lattice framework employed in the present paper, which technically is a mixed action setup based on twisted mass Wilson fermions (see Section 2), what we really need is the renormalization constant of the valence quark mass µ f appearing in the valence fermion action (10). As discussed in Refs. [19,18,16], such a renormalization constant, Z µ , is independent of the flavor f and the sign of r f in Eq. (10), as well as of all the sea Wilson parameters. The RC Z µ can be conveniently chosen and evaluated as i.e. as the inverse of the renormalization constant Z P of the pseudoscalar, flavor nonsinglet density P f f =q f γ 5 q f , where q f and q f are two distinct valence flavors of maximally twisted Wilson fermions with action as in Eq. (10) and r f = −r f . Since RI-MOM is a mass-independent scheme [52], the RCs of operators with non-vanishing anomalous dimension must be defined in the massless limit of the UVregulated theory, i.e. QCD with N f = 4 massless quark flavors. For this purpose the ensembles with fixed (non small) strange and charm sea quark masses, summarized in Table 1 and employed to compute physical observables (the so-called "production ensembles"), are not well suited. Rather one needs to produce dedicated ensembles with N f = 4 "moderately light" and, for simplicity, degenerate dynamical quarks in a lattice setup whose chiral limit coincides with the one of the lattice formulation chosen for the "production ensembles". Doing so for a sequence of progressively smaller dynamical quark mass values allows for a controlled extrapolation of massive RCestimators to the desired chiral limit.
With an eye to Section 2, a moment of thought reveals that in the chiral limit the relevant lattice regulated theory is unique (up to a choice of sign for the Wilson parameters r u , r d , r s , r c ) and corresponds to the Iwasaki action in the pure gauge sector and the standard Wilson action for N f = 4 massless fermions in the quark sector 6 .
Since RC-dedicated simulations have eventually to be performed outside the chiral limit, different numerical strategies are conceivable. The simplest and most attractive one for a RI-MOM scheme computation of RCs probably amounts to working with two degenerate maximally twisted doublets with twisted masses µ u,d,s,c = µ, which is obtained by setting m 0 = m cr , r d = −r u and r c = −r s . In such a setup RC-estimators are free of O(a) lattice artifacts at arbitrary values of the twisted mass µ and momenta p [53].
However, for at least two (β = 1.90 and β = 1.95) of the three β-values considered in this paper, the strategy outlined above could not be carried out due to numerical difficulties in implementing the maximal twist setup. In fact, in the region of small PCAC quark masses, which must be accessed when m 0 approaches m cr , Monte Carlo simulation instabilities were observed leading to very large autocorrelation times [54]. Hence we opted for an alternative strategy to achieve O(a) improvement, already proposed in Ref. [18], that does not require to work at maximal twist.
The method is based on averaging results obtained at opposite values of the PCAC quark mass and thus requires a doubling of the (in any case reasonably low) CPU time cost for producing ensembles at non-zero standard and twisted quark mass. Naturally, as it is customary in RI-MOM scheme studies of RCs, we have to consider several values of the valence mass parameters for each given choice of sea mass parameters [54,25] in order to have stable and reliable valence quark mass chiral extrapolations.
The corresponding RC computational setup, which can be viewed as a partially (un)quenched setting for Wilson tmLQCD with N f = 4 mass degenerate quark flavors at generic twist angle(s), is outlined in Section A.1, where also the choice of the relevant quark mass parameters is discussed.
In Section A.2 we recall why O(a) artifacts get canceled in correlation functions of parity-even (multi-)local operators upon averaging results obtained at opposite values of the PCAC quark mass. This is sufficient to prove the O(a) improvement of Z q and, with little more effort, of Z P and the other RCs of bilinear quark operators. In Section A.3 we report on the numerical parameters of our RC-dedicated simulations and the analysis procedure we followed. The latter is illustrated in its key aspects for a few typical examples. Our final results for the RCs in the RI-MOM scheme for the three β values considered here are given in Tables 17 and 18 together with few remarks on the conversion to the MS scheme at the 2 GeV scale. the parameters r val,sea f are restricted to ±1, while the twisted mass parameters aµ val,sea f are assumed to be non-negative.
In the partially quenched situation of interest, with all flavors mass-degenerate, a convenient and chiral covariant choice of renormalized quark mass parameters is given by [19] M sea,val = Z −1 P M sea,val Here Z A stands for the RC of the (flavor non-singlet) axial current for untwisted Wilson fermions and m sea P CAC denotes the standard PCAC quark mass computed in the unitary setup, while m val P CAC is the analogous quantity that is obtained from correlators defined in terms of valence quark fields, with valence mass parameters possibly different from their sea counterparts. More precisely, the angles θ sea which are trivially defined in terms of bare parameters of the Lagrangians (69) and (71). This is due to loop effects induced by the chiral-breaking Wilson terms. For the same reason the twist angles ω sea,val f differ from their tree-level counterparts (76)

A.2 O(a) improvement via θ-average
In the physical quark basis the parity P entails the standard fermion field transforma- where L 4 is the formal Lagrangian of the partially quenched (Euclidean) continuum QCD, 8 We define m val cr,f as the value of m val 0,f where m val P CAC = 0 at given sea quark masses. 9 In this Section to lighten notation we shall omit the flavor labels. In the physically interesting case where O is P -even, one can check by inserting the expression (80) of L 5 into Eq. (78) that terms linear in a appear either as vev's of P -odd operators (in the target continuum L 4 -theory) times a factor sin(ω sea ) or sin(ω val ), or as vev's of P -even operators multiplied by a factor cos(ω sea ) or cos(ω val ). The former terms vanish by parity, which is a symmetry of the target L 4 -theory, while the latter are in general non-zero. They however get canceled if the lattice correlator O latt M,ω of Eq. (78) is averaged with its counterpart O latt M,π−ω , as cos(π−ω sea,val ) = − cos(ω sea,val ). Notice that, in view of the twist angle definition (73), the average over ω sea,val and π − ω sea,val corresponds to averaging over θ sea,val and −θ sea,val .
This O(a) improvement property can be viewed [54,25] as a consequence of the formal invariance of the lattice mixed action (68) rewritten in the physical quark basis (74) under the spurionic transformation D d × (M sea,val 0,class → −M sea,val 0,class ) × P × (θ sea,val 0,class → −θ sea,val 0,class ). Since D d × (M sea,val 0,class → −M sea,val 0,class ) just counts the parity of the dimension of Lagrangian terms, the above spurionic invariance implies that in the SLEL of the vev's of multiplicatively renormalizable P -even (multi-)local operators all the lattice artifact contributions with odd powers of a appear with a coefficient also odd in θ sea and θ val . Hence they get canceled upon averaging vev's taken at opposite values of θ sea and θ val . We remark that by definition (see Eqs. (72) and (76)) a sign change in θ sea,val is equivalent to a sign change in θ sea,val class . In the following this way of removing the cutoff effects of first order (as a matter of fact of all odd integer orders) in a will be referred to as θ-average.
In the RI-MOM scheme the quark wave function renormalization constant, Z q , at the scale p 2 is defined by the condition Tr 10 In L 5 the occurrence of terms likeq sea,val (γ 5 γ · D)q sea,val is forbidden by charge conjugation invariance, while F ·F terms are ruled out by the P × (u ↔ d) × (c ↔ s) symmetry (thanks to r sea u,c = −r sea d,s ). where are the (momentum space) lattice propagators of the valence quark field of flavor f in the chiral limit expressed in the physical (q) and twisted (χ) basis, respectively. In practice one imposes the condition (81) at non-zero quark mass obtaining Z qestimators that must be subsequently extrapolated to the chiral limit. Applying to the massive quark propagator S q (p) the arguments on leading cutoff effects developed in the introductory part of Section A.2 and noting the P -invariance of S −1 q (p), it follows that in the lattice expression (81) the cutoff effects linear in a get canceled if S −1 q (p) is replaced by its θ-average, i.e. by the average of S −1 q (p) evaluated at (M, ω) and its analog evaluated at (M, π − ω). The θ-average procedure guarantees O(a) improvement already at the level of the RC-estimators in the massive theory.
A.2.2 O(a) improvement of Z P , Z S and Z T With the usual notation, according to which RCs are denoted by the name they would have in the twisted basis, where the fermionic sector of the Lagrangian is given by Eqs. (69) and (71), the formulae that define in the chiral limit the RCs of quark bilinear operators in the RI-MOM scheme read where Γ = 1 1, γ 5 , γ µ , γ µ γ 5 , σ µν , P Γ is a Dirac projector satisfying Tr(ΓP Γ ) = 1, while χ val 1 and χ val 2 are valence quark fields with flavor indices f = 1 and f = 2, and parameters r val 1 and r val 2 , respectively. For the case of r val 2 = −r val 1 in the valence fermion Lagrangian (71), passing to the physical quark basis we have ω val 2 = −ω val 1 . If Γ = γ 5 in Eq. (84) we thus find the identity From this identity, applying the arguments developed in the introductory part of Section A.2 to a 8 x,y e −ip(x−y) q val 1 (x)(q val 1 γ 5 q val 2 )(0)q val 2 (y) latt M,ω , as well as to S −1 q 1 (p) and S −1 q 2 (p), and noting that the Dirac trace of their combination in the r.h.s. is a parity invariant form factor, we conclude that taking the θ-average of the lattice expression in Eq. (85), improved estimators of Z q /Z P for all values of M val and M sea are obtained. Once an O(a) improved determination of Z q is available, Z P can be extracted with only O(a 2 ) artifacts by appropriate chiral extrapolations. The argument for the O(a) improvement via θ-average of the lattice estimators of Z S and Z T is identical to the one given above for Z P , because for Γ = 1 or Γ = σ µν and r val 2 = −r val 1 we find identities completely analogous to Eq. (85) -of course with γ 5 and P γ 5 replaced by the relevant Dirac matrix Γ and associated projector P Γ .

A.2.3 O(a) improvement of Z V and Z A
As the massive lattice estimators for Z q /Z V,A in the RIMOM approach are provided by Eq. (84) with Γ = γ µ or Γ = γ µ γ 5 , passing from the twisted to the physical quark basis, in the case of r val 2 = −r val 1 , identities analogous to Eq. (85) are obtained, but (owing to anti-commutation of Γ with the γ 5 occurring in the equation relating χ val f and q val f ) with a more complicated r.h.s. If, for instance, Γ = γ µ , upon setting −ω val 2 = ω val 1 ≡ ω val , we find x,y e −ip(x−y) q val 1 (x)V 12,µ (0)q val 2 (y) latt M,ω S −1 q 2 (p)P γµγ 5 + +iCS Tr S −1 q 1 (p) x,y e −ip(x−y) q val 1 (x)A 12,µ (0)q val 2 (y) latt M,ω S −1 q 2 (p)P γµ (86) with C = cos ω, S = sin ω, V 12,µ =q val 1 γ µ q val 2 and A 12,µ =q val 1 γ µ γ 5 q val 2 . Looking at the r.h.s. of this identity, we note that the expressions with pre-factors C 2 and S 2 (±iCS) are parity-even (parity-odd) form factors. Then applying the Symanzik analysis arguments developed in the introductory part of Section A.2 to x,y e −ip(x−y) q val 1 (x)V 12,µ (0)[A 12,µ (0)]q val 2 (y) latt M,ω , as well as to S −1 q 1 (p) and S −1 q 2 (p), we see that • to order zero in a, the terms with pre-factors ±iCS vanish by parity, while those with pre-factors C 2 and S 2 are non-zero (and coinciding in the limit of unbroken chiral symmetry); • to first order in a the contributions that do not vanish by parity are those obtained either by inserting the P -even piece of L 5 in the terms on the the r.h.s. of Eq. (86) with pre-factors C 2 and S 2 or by inserting the P -odd piece of L 5 in the terms with pre-factors ±iCS.
Taking also into account the ω-and ω sea -dependence of L 5 (see Eq. (80)), one checks that, while the contributions of zero (actually even integer) order in a are even under ω → π −ω, all the contributions of first (actually odd integer) order in a are odd under ω → π − ω. Hence by taking the θ average of the lattice expression (86), the lattice artifacts of odd order in a get canceled, leaving out the contributions of order a 2n . This proves the O(a) improvement by θ average of the massive lattice estimators of Z V , from which the RC is extracted after chiral extrapolations. Identical arguments clearly hold as well for Z A , because for Γ = γ µ γ 5 we find an identity completely analogous to Eq. (86) with the axial and vector operators and the associated projectors properly interchanged.
We have focused here on the choice r val 2 = −r val 1 for the valence parameters of quark bilinear operators, because this is the case with smallest statistical fluctuations in the numerical evaluation of RCs and which the results quoted in the following refer to. The discussion of the alternative (and computationally more noisy) choice r val 2 = r val 1 could be carried out along similar lines 11 , finding again that upon θ average the estimators of the RCs of all quark bilinear operators are O(a) improved.

A.3 Numerical details and results
In Table 16 we report the information on the relevant simulation parameters for the three ensembles we have considered in this paper. Except for the θ-average, which is implemented in order to achieve the O(a) improvement out of the maximal twist, the other parts of the analysis follow closely the procedure described in Ref. [53].
For each ensemble in the table we compute the RC-estimators at values of momenta, p µ = (2π/L µ ) n µ , with components lying in the following intervals ( [2,5] , [2,5] , [2,5] , [4,9]) , for β = 1.90 and 2.10 (87) and L µ denoting the lattice size in the direction µ. Anti-periodic boundary conditions on the quark fields in the time direction are implemented by a shift of the time component of the four-momentum by the constant ∆p 4 = π/L 4 . The final analysis of the RC estimators has been performed at four-momenta that pass the "democratic" momentum cut defined by 11 If r val 2 = r val 1 , however, the identities obtained when passing from the twisted to the physical quark basis have a simple r.h.s. for the case of Z V,A and a more complicated one (like the one of Eq. (86)) for Z P,S,T .
As a typical example, in Fig. 18 we show the effect of the subtraction of the Goldstone pole in the amputated two-point correlators for the ensembles B4m and B4p (the most critical ones at β = 1.95). The quantities V P and V sub P are defined according to Eqs. (3.4) and (3.12)-(3.13) of Ref. [53].  Figure 18: Amputated pseudoscalar density two-point correlators before (V P , red squares) and after (V sub P , blue dots) the Goldstone pole subtraction, at β = 1.95 and (ap) 2 ≈ 1.5. Panels (a) and (b) correspond to data from ensembles B4m and B4p, respectively.  Figure 19: Chiral extrapolations of θ-averaged B4m and B4p data at β = 1.95 and (ap) 2 ≈ 1.5. Panel (a): valence quark mass extrapolation of Z q data. Panel (b): sea quark mass extrapolation of Z q and Z P data. The blue dots have been displaced by 0.2 from the red squares for better visibility.
In Fig. 19(a) we plot the valence mass extrapolation of Z q , while Fig. 19(b) shows the sea quark mass extrapolation in the chiral limit of θ-averaged B4m and B4p data for the cases of Z q and Z P at β = 1.95 and (ap) 2 ≈ 1.5.
Our final estimates of the RCs of all the quark bilinear operators in the RI-MOM(3 GeV) and MS(2 GeV) schemes are collected in Tables 17-20. They have been obtained after subtracting the perturbative cutoff effects up to O(a 2 g 2 boost ), where as usual g 2 boost = 6/(β P ), with P the average plaquette value. Perturbative estimates of the discretization effects on the Green function of the operators of interest can be found in Ref. [55]. Finally the RCs Z q , Z P and Z S , obtained in the RI-MOM(3 GeV) scheme, have been converted to the MS(2 GeV) scheme using the appropriate N 3 LO formulae from Ref. [56], while for the RC Z T the N 2 LO formula of Ref. [57] has been employed.
Recently, the RCs Z V , Z A , Z P and Z S have been computed perturbatively up to three loops in Ref. [58] at β = 1.95 and 2.10. The comparison with our nonperturbative results of Tables 17 and 20 shows a remarkable, fair agreement within the quoted errors.
In Tables 17-19 we have given our results for the RCs derived from two different methods, M1 and M2, that differ for the way one deals with the residual (ap) 2 discretization effects [53]. The method M1 consists in extrapolating the RCs linearly to (ap) 2 → 0, after fitting the functions Z q,Γ (µ = 3 GeV; (ap) 2 ) in the wide momentum interval (ap) 2 ∈ [1.5, 2.2]. The slopes of the fits, λ q,Γ = dZ q,Γ (µ = 3 GeV; (ap) 2 )/d(ap) 2 at each value of β exhibit only a very mild dependence on the coupling constant. Following the discussion of Ref. [53] (see Section 3.2.2 and in particular the arguments leading to Eq. (3.24) of that reference), we assume a simple linear dependence of λ q,Γ on β, and perform a simultaneous extrapolation of Z q,Γ (µ = 3 GeV; (ap) 2 ) towards (ap) 2 → 0 for the three values of β (see for instance    Table 19: RCs Z q , Z P , Z S and Z T obtained in the MS(2 GeV) scheme using the methods M1 and M2.
artifacts occurring in the RCs of the method M2 will be removed once the continuum limit of the physical quantities of interest is taken, as shown in Fig. 21 Table 20: RCs Z V , Z A and Z P /Z S obtained with the methods M1 and M2. We also present the accurate results for Z V obtained using the Ward-Takahashi identity (WTI) (for more details see Section 2.3 of Ref. [53]).
the squared pion mass.