Quark masses using twisted mass fermion gauge ensembles

We present a calculation of the up, down, strange and charm quark masses performed within the lattice QCD framework. We use the twisted mass fermion action and carry out simulations that include in the sea two light mass-degenerate quarks, as well as the strange and charm quarks. In the analysis we use gauge ensembles simulated at three values of the lattice spacing and with light quarks that correspond to pion masses in the range from 350 MeV to the physical value, while the strange and charm quark masses are tuned approximately to their physical values. We use several quantities to set the scale in order to check for finite lattice spacing effects and in the continuum limit we get compatible results. The quark mass renormalization is carried out non-perturbatively using the RI'-MOM method converted into the $\overline{\rm MS}$ scheme. For the determination of the quark masses we use physical observables from both the meson and the baryon sectors, obtaining $m_{ud} = 3.636(66)(^{+60}_{-57})$~MeV and $m_s = 98.7(2.4)(^{+4.0}_{-3.2})$~MeV in the $\overline{\rm MS}(2\,{\rm GeV})$ scheme and $m_c = 1036(17)(^{+15}_{-8})$~MeV in the $\overline{\rm MS}(3\,{\rm GeV})$ scheme, where the first errors are statistical and the second ones are combinations of systematic errors. For the quark mass ratios we get $m_s / m_{ud} = 27.17(32)(^{+56}_{-38})$ and $m_c / m_s = 11.48(12)(^{+25}_{-19})$.

In the analysis we use gauge ensembles simulated at three values of the lattice spacing and with light quarks that correspond to pion masses in the range from 350 MeV to the physical value, while the strange and charm quark masses are tuned approximately to their physical values. We use several quantities to set the scale in order to check for finite lattice spacing effects and in the continuum limit we get compatible results. The quark mass renormalization is carried out nonperturbatively using the RI'-MOM method converted into the MS scheme. For the determination of the quark masses we use physical observables from both the meson and the baryon sectors, obtaining m ud

I. INTRODUCTION
Quark masses are essential inputs of the Standard Model (SM) and play a primary role for the description of a large number of physical processes that can provide insights into the dynamics of the SM as well as in the search of beyond the Standard Model physics. The quark masses together with the strong coupling constant can be regarded as the fundamental parameters of Quantum chromodynamics (QCD), the renormalizable theory of the strong interactions. Therefore, their determination plays a crucial role in the phenomenological description of the plethora of complex phenomena governed by strong nuclear forces taking place in the universe as well as at particle colliders. Lattice QCD provides a non-perturbative approach based on first principles and systematically improvable for determining the quark masses and the strong coupling. In this approach the QCD Lagrangian is defined on a discrete Euclidean 4-dimensional space-time lattice of large but finite extent, which allows for numerical simulation of the theory via Monte Carlo methods. The finite volume and the non-vanishing lattice spacing introduce systematic artefacts, which can be theoretically understood, kept under numerical control and extrapolated away in order to extract the physical quantities of interest.
Theoretical progress in lattice field theory and improvement in numerical algorithms, accompanied with a continuously increasing computational power, are allowing us to perform simulations using physical values of the light-quark masses. However, most of these simulations are still carried out using a single lattice spacing and volume, although this is rapidly changing as more lattice QCD collaborations gain access to larger computational resources and can produce multiple ensembles of gauge configurations generated with physical values of the light quark masses. Such ensembles will be referred to as physical point ensembles. In this work we include two physical point ensembles at two different lattice spacings. In order to take the continuum limit we also employ additional ensembles at a coarser lattice resolution with larger than physical pion masses. Globally we thus use ensembles with three values of the lattice spacing and spanning pion masses in the range from about 350 MeV to 135 MeV, which enable us to perform a combined chiral and continuum extrapolations. In order to study systematic effects in the determination of the quark masses, we use two sets of observables to set the scale and to evaluate the quark masses. One set of observables is based on quantities from the meson sector of QCD while the other set relies on baryonic observables. In the former case we use the pion mass and decay constant to set the scale and to determine the average up/down quark mass. The mass of the strange and charm quarks are extracted using the kaon and D-meson masses, respectively. In the latter case, instead, the masses of the pion and nucleon are employed to set the scale and fix the average up/down quark mass, while the Ω − and the Λ c masses determine the strange and charm quark masses, respectively. In this way we obtain a valuable consistency check with respect to the results coming from the mesonic sector.
For the renormalization of the quark mass we employ a dedicated set of gauge ensembles with four mass-degenerate sea quarks (having mass around half of the strange quark mass). Such a set ensures a good control of the extrapolation to the massless limit. We perform the computation in an intermediate mass-independent scheme, which is finally converted to the standard MS scheme.
The paper is organized as follows: In Section II we describe the gauge ensembles used in this study and explain our methodology. In Sections III and IV we present the methods used to set the lattice spacing a and to carry out a non-perturbative computation of the renormalization constant Z P including a detailed discussion on the control of hadronic contaminations and other systematic errors. In Sections V and VI we describe the extraction of the quark masses and their ratios using inputs from the mesonic and baryonic sectors, respectively. In Section VII, we discuss our final results and give our conclusions and outlooks.

II. METHODOLOGY
In twisted-mass lattice QCD [1] the discretized Dirac operator in the physical quark basis is written as where r f = ±1, ∇ µ and ∇ * µ are nearest-neighbor forward and backward covariant derivatives, µ f sets the mass of the quark field q f of flavor f , c SW is the coefficient of the clover-term Q µν [2] and m cr is the critical value of the "untwisted" mass, m 0 , obtained by requiring the vanishing of the partially conserved axial current (PCAC) mass, as discussed in Ref. [3]. This condition, referred to as maximal twist, guarantees automatic O(a)-improvement of physical observables [4,5]. In the twisted-mass fermion formulation at maximal twist, the renormalized quark masses are thus given by where a is the lattice spacing and Z P is the pseudo-scalar renormalization constant. Therefore, the determination of both Z P and the lattice spacing, combined with inputs from known physical quantities depending on the quark masses, enables us to extract m f . Having gauge ensembles with at least three different lattice spacings at several pion masses allows us to take the continuum limit and by performing a chiral extrapolation we can determine the quark masses at the physical point.

A. Gauge ensembles
We use the twisted-mass fermion discretization scheme [1,4] with the inclusion of a clover-term [2]. As already explained, twisted-mass fermions (TMF) provide an attractive formulation for lattice QCD simulations allowing for automatic O(a) improvement of physical observables as well as renormalization constants [4,6]. This is an important property since quantities of interest have lattice artifacts of O(a 2 ) and are thus closer to the continuum limit. A clover-term is added to the TMF action to suppress O(a 2 ) breaking effects between the neutral and charged pions, which eventually leads to the stabilization of simulations with light quark masses close to the physical pion mass. For more details on the TMF formulation see Refs. [7,8] and on the simulation and tuning strategies see Refs. [9,10].
In this study we analyze ten gauge ensembles simulated at three values of the lattice spacing and at several values of the pion mass, spanning a range from the physical pion mass up to 350 MeV. Some parameters of these ensembles and the values of few key physical quantities are listed in Table I. More details are given in Ref. [11]. With respect to Ref. [11] the ensemble cC211. 20.48 has been added in order to investigate the light quark mass dependence at the finest lattice spacing.
The ensembles are generated with two mass-degenerate light quarks and the strange and charm quarks in the sea (N f = 2 + 1 + 1 ensembles). The strange and charm sea quark mass parameters, aµ σ and aµ δ (see, e.g., Eq. (8) of Ref. [12]) have been adjusted so as to reproduce the phenomenological conditions m c /m s 11.8 and m Ds /f Ds 7.9 [13], which are easy to implement with few percent level precision even using simulations with larger than physical pion masses and on lattices of linear size L 2.5 fm, as detailed in Ref. [10]. The condition on m Ds /f Ds is sensitive to the charm quark mass while the one on m c /m s fixes the strange quark mass. In this way the charm and strange sea quark mass parameters have been tuned, separately for each lattice resolution (or β), to bare values that a posteriori turn out to yield values for the renormalized sea quark masses that are consistent within statistical errors of few percents with those we determine in this paper, as discussed in the following, at the physical pion mass point, on large volumes and in the continuum limit.

Ensemble
L 3 × T MDUs aµ amπ afπ mπL mN /mπ mπ [MeV] β = 1.726, cSW = 1.74, aµσ = 0.1408, aµ δ = 0.1521, w0/a = 1.8352 (35) cA211.53. 24  Parameters of the N f = 2 + 1 + 1 ensembles analyzed in this study. In the first column we give the name of the ensemble, in the second the lattice volume, in the third the number of molecular dynamics units simulated per ensemble, in the fourth the twisted-mass parameter, aµ , for the average up/down (light) quark, in the fifth and in the sixth the pion mass amπ and decay constant afπ in lattice units from Ref. [11], in the seventh the pion mass times the lattice spatial length, mπL, in the eighth the ratio mN /mπ as determined in Section VI and, finally, in the last column the pion mass in physical units, using our determination of the gradient-flow scale w0 obtained in Ref. [11] (see later Eq. (32)). We also include for each set of ensembles with the same lattice spacing the coupling constant β, the clover-term parameter cSW , the parameters of the non-degenerate operator aµσ and aµ δ , related to the renormalized strange and charm sea quark masses [5], and the value of the gradient-flow scale w0/a determined at the physical pion mass in Ref. [11].

B. Osterwalder-Seiler fermions
A naive use of the twisted-mass action for non-degenerate strange and charm quarks would lead to an undesired O(a 2 ) mixing of the strange and charm flavours in the correlation functions of interest to determine physical quantities [14,15]. In order to avoid such a mixing in the correlation functions, we adopt a non-unitary lattice setup [5] where the twisted-mass action for non-degenerate strange and charm quarks is employed only in the sea sector, while the valence strange and charm quarks that enter the correlation functions are regularized as exactly flavour-diagonal Osterwalder-Seiler fermions [16]. Thus, the valence action in the strange and charm sectors (f = s, c) reads where D(µ f ) is the twisted mass Dirac operator in Eq. (1) with the same values of m 0 and c SW as in the sea sector action used for ensemble generations. As the renormalized strange and charm sea quark masses are matched with few percent relative accuracy to their valence counterparts, no significant unitarity violation is expected in our continuum limit results. When constructing meson correlation functions, the Wilson parameters of the two valence quarks are always chosen to have opposite values. This choice guarantees that squared pseudoscalar meson masses, generically indicated by m 2 P S , differ from their continuum counterparts only by terms of order O(a 2 µ) [4,7]. As we said above, in our lattice setup the (valence) flavour conservation is guaranteed in all correlation functions and automatic O(a)-improvement is maintained. Of course we need to fix the valence strange and charm quark masses, µ s and µ c , by imposing suitable renormalization conditions. For this purpose in the present work we use two different sets of observables. Namely, in one case we use the mass of the physical masses of kaon and D (or D s ) mesons and in the other case the masses of the Ω and Λ c baryons. These two different choices will lead to two different determinations of the strange and charm quark masses, which will enable us to check the consistency of our results when using physical inputs from the mesonic and baryonic sectors.

III. SCALE SETTING
As already mentioned, we have ensembles at three different lattice spacings. We will refer to the ensembles in Table I that start with cA in their names as A ensembles, those starting with cB as B ensembles and those with cC as C ensembles (the label c stands for clover). Each of these groups have the same lattice spacing, with the A ensembles having the largest lattice spacing and C ensembles the smallest one. In what follows we will use different quantities to determine the three lattice spacings. This will allow us to check consistency while taking the continuum limit when different inputs are used. In the pion sector, the pion mass and decay constant are used as input. Within this approach one also determines the value of the gradient-flow scale w 0 . We use the iso-symmetric values of the pion mass and decay constant, given respectively by [17], m isoQCD π = 135.0(2) MeV and f isoQCD π = 130.4(2) MeV .
We also compute the value of w 0 /a for each ensemble (see Table I) and extrapolate to the physical pion mass and continuum limit. We find w 0 = 0.17383(63) fm [11] and using this value one determines the lattice spacings shown in Table II. We refer to this determination of the lattice spacings as coming from the "pion" sector. Details are given in Ref. [11]. Another quantity used for the determination of the lattice spacings is the mass of the nucleon [18,19]. Details on the extraction of the nucleon mass are given in Section VI. In order to fit the pion mass dependence of the nucleon mass, we use the well established SU(2) chiral perturbation theory result to one-loop [20,21] The three values of the lattice spacing, which will be denoted by a A , a B and a C , can be determined from the lattice data for the nucleon and pion masses by rewriting Eq. (5) as where (a i m N ) and (a i m π ) are our lattice QCD results and i = A, B, C. The three quantities a i as well as the nucleon mass in the chiral limit, m 0 N , are treated as fitting parameters, while the value of c 1 is fixed by requiring the reproduction of the physical value of the nucleon mass, m isoQCD N , at the physical pion point (4), namely We restrict ourselves to using ensembles for which the pion mass is below 260 MeV since chiral perturbation theory to higher orders has larger ambiguities. The simulations of the gauge ensembles use mass-degenerate up and down quarks and include no electromagnetic effects. Thus, we use the average value of the proton and neutron mass as our input for fixing the lattice spacings, namely we assume m isoQCD N = 0.9389 GeV in Eq. (7). We also use the physical value for the axial charge, g A = 1.27641(56) [22] and for consistency the physical value of f π from Eq. (4). The ratio g A /f π appears in the m 3 π term and any residual correction due to strong isospin breaking and electromagnetism is neglected.
The result of the fit to the mass of the nucleon m N is depicted in Fig. 1 and describes very well the data, yielding  Table II. In Fig. 1 we also show the ratio m N /m π and the resulting fit using the parameters extracted from the fit to the nucleon mass. As it can be seen, the data for m N /m π are well described. In the left panel we show the nucleon mass as a function of the pion mass mπ squared. In the right panel we show the dimensionless quantity mπ/mN as a function of mπ, using the lattice spacing extracted from the nucleon mass. The values of mπ/mN are listed in Table I, while amN and amπ are determined in Section VI. The solid line shows the fit to the lattice QCD data using Eq. (5). The value of χ 2 /d.o.f. is 0.19, where the number of degrees of freedom is five.
The values of the lattice spacing extracted from the pion sector and from the nucleon mass are shown in Table II and they differ by O(a 2 ) effects. Fitting their difference as a function of a 2 is shown in Fig. 2. We observe that in the continuum limit the difference vanishes, as expected for our O(a)-improved formalism. In what follows we will use both determinations to extract the quark masses. This provides a cross-check for our procedure and for the magnitude of any residual lattice spacing effect. MeV. The solid line shows the linear fit in a 2 to the results extracted by using ensembles with mπ < 260 MeV (full symbols), which is largely consistent with zero in the continuum limit.

IV. COMPUTATION OF ZP
In order to determine the renormalized quark masses it is crucial to perform an accurate evaluation of the mass renormalization factor Z m , that in the maximally twisted-mass formulation is given by Z m = 1/Z P (see Eq. (2)). For this reason the details of the procedure we have followed to compute Z P will be given in this Section.
For the calculation of Z P we employ the non-perturbative RI'-MOM renormalization scheme [23], which is a mass-independent scheme since the renormalization constants are defined in the massless limit. The choice of this intermediate scheme is convenient in that the scale evolution for the renormalization constants of the operators with non-trivial anomalous dimension is controlled by the renormalized gauge coupling alone. This requires however simulations close enough to the chiral limit, which is not the case of the N f = 2 + 1 + 1 ensembles of Table I, mainly due to the presence of the heavy sea charm quark. In order to safely take the chiral limit in the computation of the renormalization constants, we have thus separately produced gauge field configurations with four mass-degenerate quarks (N f = 4) at the same value of the coupling β and with the same clover term included in the fermionic action as for our A, B and C ensembles of Table I. This ensures that in the chiral limit the same massless N f = 4 QCD theory underlies both the ensembles used for computing hadronic observables and setting the energy scale (see Table I) and the ensembles dedicated to the evaluation of the renormalization constants, about which details are given in Table III. The four degenerate quarks are taken with masses from ∼ 8 to ∼ 16 times larger than the average up-down quark mass, which simplifies both the simulations and the tuning to maximal twist. The values of the critical mass m cr have been chosen in order to satisfy the maximal twist condition (which is convenient to reduce lattice artifacts) to a very good accuracy level, as it can be deduced from the smallness of the PCAC masses in Table III (actually the tuning is even slightly better than the one corresponding to the N f = 2 + 1 + 1 ensembles of Table I   In the RI'-MOM scheme we obtain the renormalization constant of the (flavour non-singlet) pseudoscalar density operator O P = iqγ 5 q with r q = −r q = −1 and µ q = µ q via the following condition 1 where V P is the pseudoscalar vertex function between quark and antiquark states with momentum p and Z q is the renormalization constant of the quark field, defined as Here S q is the quark propagator at momentum p, which is identified with the renormalization scale µ 0 . In this work we adopt the alternative definition of Z q first proposed in Ref. [6], where the sum is over the N p non-vanishing components of the lattice momentum ap µ = sin(ap µ ). The prescription of Eq. (10), unlike the naive RI'-MOM definition (9), has no lattice artifacts at tree-level and beyond tree-level it exhibits quite small O(a 2 ) cutoff effects.
The subtraction of the Goldstone pole in the vertex function V P requires a good control of the vertex mass dependence. Therefore, we find it more suitable to adopt a partially quenched (PQ) setup, in which propagators and vertices are computed for multiple values of valence quark masses µ val at fixed sea mass µ sea , namely we use In order to reduce the effect of Lorentz non-invariant cut-off effects, we filter the momenta selecting the ones that are isotropic ("democratic") in the spatial directions, thus satisfying: The above constraint ensures that unwanted hypercubic lattice artifacts are suppressed [6]. We further improve the Z P estimator by using results from lattice perturbation theory (for more details see, e.g., Refs. [24,25]). In summary, we calculate the lattice artifacts at one-loop level and to all orders in the lattice spacing, O(g 2 a ∞ ). The perturbative corrections to the Green function of the pseudoscalar operator, as well as to the quark propagator, are evaluated for each momentum ap at which the renormalization constants are computed. It should be noted that each value of ap requires a separate calculation of the O(g 2 a ∞ ) correction, as the perturbative contributions are not analytical and require the numerical evaluation of one-loop integrals. Such contributions also include the leading order terms, O(g 2 a 0 ), that have to be separated from the pure O(a n ) terms (n = 0). The computation of perturbative lattice artifacts to O(g 2 a ∞ ) done in Ref. [25] are adapted for the case of the specific definition of Z q in Eq. (10). Thus, the one-loop perturbative corrections that we use are defined as follows

∆V
(1) where Z (1) q and V (1) P are the one-loop contributions to the quark field renormalization constant and to the pseudoscalar vertex, respectively. Using the above quantities, we extract improved non-perturbative estimates for Z q and Z P by modifying the renormalization conditions as: A. Analysis method and safety checks against hadronic contaminations The improved Z P estimators, namely Z impr P obtained from Eq. (16), are evaluated at different values of p 2 = µ 2 0 and extrapolated to the chiral limit, a step which is discussed in detail in the following Section IV B. Then, the chirally extrapolated lattice estimators of Z P are evolved to a common reference scale p 2 = µ 2 ref in the RI'-MOM scheme using the anomalous dimension known up to N 3 LO according to Ref. [26] and adopting Λ QCD (N f = 4) = 294 (12) MeV from Ref. [13]. The curves obtained at the three β values are reported in Fig. 3   A dependence on p 2 of the Z P (µ 2 ref ) estimators in Fig. 3 is expected due to lattice artifacts, i.e. O(a 2 p 2 ) terms, and possibly also to residual O(a 0 ) hadronic contaminations 2 , which, however, must vanish as 1/p 2 at large p 2 . The plots of Z P (µ 2 ref ) in Fig. 3 show a very good linearity in p 2 within the range p 2 ∈ [15, 24] GeV 2 , which is the one relevant for our determination of Z P and hence of the renormalized quark masses. This fact indicates that lattice artifacts other than a 2 p 2 -terms and possible hadronic contamination effects are negligible within our small statistical errors. Such a property is explicitly checked at each β value by performing a fit of Z P (µ 2 ref = 17 GeV 2 ) (left panel of Fig. 3) with the two Ansätze given by The resulting best fit values for z i are given in Table IV. The coefficients z 2 and z −1 are compatible with zero within statistical errors at all β values, while the coefficients z 1 scale nicely with a 2 (see Tab. II). From this check we see that the systematic uncertainties on Z P are negligible within our small statistical errors. We will comment on the value of z 0 when we present our results in Sec. IV C.
Within the present study of renormalized quark masses, we follow two different methods for determining Z P in the RI' scheme and use data in two different p 2 -ranges. The first method (M1) consists in fitting the Z P (µ 2 ref ) data linearly in p 2 in a given p 2 -range with the aim of removing O(a 2 p 2 ) discretization effects, while in the second method (M2) the same data are fitted to a constant [12]. Method M2 is by construction much less sensitive than M1 to possible small residual hadronic contaminations but at the expense of leaving some O(a 2 ) artifacts in the determination of Z P . The ranges of p 2 used in the present analysis are p 2 ∈ [15, 19] GeV 2 and p 2 ∈ [18, 24] GeV 2 for determining Z P (17 GeV 2 ) and Z P (21 GeV 2 ), respectively.
As an additional check of our determination of Z P (µ 2 ) in the RI'-MOM scheme, we show in Fig. 4 the results for the non-perturbative step scaling function Σ P (µ 2 We see that the lattice QCD data exhibit small discretization errors and agree in the continuum limit with the perturbative counterpart Σ pt Moreover, as will be shown in Section V, using the four Z P determinations corresponding to the methods M1 and M2 and at the two reference scales µ 2 ref = 17 GeV 2 and µ ref = 21 GeV 2 we obtain in the continuum limit consistent final results for the renormalized quark masses.

B. Chiral extrapolation and Goldstone boson pole subtraction
A crucial step in determining the renormalization constant Z P is the extrapolation of its lattice estimators to the chiral limit, where the mass-independent RI'-MOM scheme is defined.

Hadronic contaminations in the pseudoscalar vertex
It is well-known that the pseudoscalar vertex V P receives contributions at the non-perturbative level by hadronic contaminations whose leading term scales as ∼ (p 2 m 2 π ) −1 [23]. Such Goldstone boson pole has to be identified and subtracted from the data. In a unitary lattice setup for QCD with N f = 4 degenerate flavours of mass m q , the lattice estimator of the vertex, v P (p 2 , m q ), is expected to have the form where the dimensionless quantities V P (our target vertex), h, h and h depend in general on a 2 p 2 , a 2 Λ 2 QCD , a 2 m q Λ QCD and a 2 m 2 q , while the ellipses stand for terms suppressed by higher powers of 1/p 2 as p 2 → ∞. We note that terms linear in m q are either hadronic contaminations suppressed as ∼ 1/p 2 at large p 2 or lattice artifacts of the form ∼ a 2 m q Λ QCD , which are numerically tiny for the am q values of interest here. Since close to the chiral limit m 2 π ∼ m q , an equivalent Ansatz for v P (p 2 , m q ) can be written in the form 3 where we separate the hadronic contaminations decreasing, for large p 2 , like 1/p 2 from the vertex of interest V P (p 2 ).

Choice of a partially quenched setup
In a PQ setup, such as the one adopted in the present analysis (see Sec. IV A), the lattice action is power counting renormalizable and the operator vertices evaluated at several values of valence (µ val ) and sea (µ sea ) quark masses approach, as (µ val , µ sea ) → (0, 0), the corresponding unitary vertices from which the RI'-MOM renormalization constants can be computed. As detailed in Sec. IV A above, at all β values we use nine values of µ val for each of the four µ sea values. This allows us to have a good control on the mass dependence of the pseudoscalar vertex and to adopt, at fixed β, p 2 and µ sea values, the following fit Ansatz for the chiral fit in where the quantities V P (p 2 , µ sea ), H, H and H depend, besides on µ sea (to a numerically negligible level, as we shall see below), also on a 2 p 2 , a 2 Λ 2 QCD , a 2 m q Λ QCD and a 2 m 2 q , while the ellipses have the same meaning as in Eq. (18). Noting that the hadronic contaminations in the three-point correlation function of quark, pseudoscalar bilinear and antiquark fields at fixed four-momenta (and in the derived quantity v P ) arise from the time orderings where the quark and antiquark fields are located at time distances both before or both after the pseudoscalar density [23], it follows that the Goldstone boson pole contamination is controlled by the mass [m 2 π ] val of the valence pion that appears as an intermediate state in the aforementioned time orderings. Recalling also that, to leading order in PQ chiral perturbation theory, [m 2 π ] val ∼ µ val [27], we choose to use the equivalent Ansatz where again we separate the hadronic contaminations (suppressed like 1/p 2 as p 2 → ∞) from the vertex of interest and the dimensionful coefficients K, K and K in general depend on µ sea and may be affected by lattice artifacts. From the fit of v P in µ val at fixed p 2 and µ sea , one can determine the coefficients K and K , but it is not possible to disentangle the vertex V P (p 2 , µ sea ) from the hadronic contamination K /p 2 . At this point one can safely take the limit µ sea → 0 and determine V P (p 2 , 0) + K | µsea=0 /p 2 . After taking the full chiral limit, one can check that the residual hadronic contamination K | µsea=0 /p 2 is completely negligible in the range of p 2 used for the extraction of Z P , as already detailed in Section IV A (see the discussion around Tab. IV).

Numerical data and intermediate analysis results
Within the p 2 ranges used in the present analysis, our data for v P (p 2 , µ val , µ sea ) exhibit a tiny linear dependence on µ val , which is compatible with K = O(a 2 ) up to statistical errors. This is checked by fitting the data to the Ansatz of Eq. (21) for each fixed aµ sea and p 2 , which determines the quantities [V P (p 2 , µ sea ) + K /p 2 ], K/p 2 and K /p 2 , and then studying the dimensionless ratio d ≡ K /(w 0 p 2 ) as a function of (a/w 0 ) 2 . At fixed β value, we observe that d changes with the sea quark mass of a given gauge ensemble non monotonically in µ sea by the same amount as the statistical errors. Therefore, we average over the values of d at each µ sea for a fixed β value, obtaining the quantity d sea , which is shown in Fig. 5 for three representative p 2 values, namely 13, 20 and 26 GeV 2 . As one can see, the continuum limit of d is consistent with zero up to statistical errors and/or a small residual term, which even if present (given the numerical values of aµ val 0.02) would alter Z P only to O(10 −4 ). In view of the evidence that K is O(a 2 ) or numerically negligible in the p 2 range of interest here, we can perform the fit in µ val on the data for v P excluding the term linear in µ val , namely we use This procedure has the advantage of yielding small statistical errors at the price of including well-controlled O(a 2 ) artifacts in the numerical estimate of [V P (p 2 , µ sea ) + K /p 2 ] and hence of Z P . The results of the fit on v P for few values of p 2 and the two extreme values of µ sea are shown in Fig. 6 for the cases β = 1.726 and β = 1.836. Besides the very good quality of the fits we remark that the resulting estimates of [V P (p 2 , µ sea )+K /p 2 ] at µ val = 0 indeed shows a very tiny dependence on µ sea as mentioned above. Such a dependence on µ sea turns out to be of the same size as the statistical errors (about ∼ 0.5%), non monotonic in µ sea at fixed β and with different trends at different β's.
This feature is illustrated in Fig. 7, where the resulting estimates of Z P (µ sea , p 2 ), obtained using the RI condition of Eq. (8), are shown at β = 1.726 and β = 1.836. Therefore, as the dependence on µ sea of Z P (µ sea , p 2 ) is not statistically significant, we average them in order to estimate Z P (p 2 ) in the unitary chiral limit.  (24)(35) 0.5159(24)(37) 0.4851(26)(32) 0.5229(24)(34) 1.778 0.4888 (33)(35) 0.5121(26)(37) 0.4877(27)(32) 0.5184(23)(34) 1.836 0.4976 (26)(36) 0.5133(23)(37) 0.4978(27)(33) 0.5169(24)(34) TABLE VI. Results for ZP in the RI'-MOM scheme but evolved to the common reference scale µ 2 ref = 19 GeV 2 . The notation is the same as that in Table V and   In Table V we give the results of Z P determined in the RI'-MOM scheme using methods M1 and M2 (see Sec. IV A) for the two reference scales µ 2 ref = 17 GeV 2 (labelled by "a") and µ 2 ref = 21 GeV 2 (labelled by "b"). We find that at each β value, the parameter z 0 appearing in Eq. (17) is found to be compatible with the results for Z P (µ 2 ref = 17 GeV 2 ) extracted using the method M1, as expected since this method corresponds to a linear fit Ansatz in p 2 . The results are obtained by fitting the data in the momentum ranges p 2 ∈ (15, 19) GeV 2 and p 2 ∈ (18, 24) GeV 2 , respectively, for the reference scales 17 GeV 2 and 21 GeV 2 . The values of Z P given in ) cut-off effects. This implies that using whichever of them leads to equivalent results for the renormalized quark masses and renormalized matrix elements of the pseudoscalar density in the continuum limit. As a check of the good accuracy to which this property is expected to be satisfied, we show in Fig. 8 Table, which is due to the N 3 LO evolution to 19 GeV 2 and is independent from the lattice spacing. Therefore the continuum limit value and its uncertainty are obtained by taking into account only the statistical errors at finite a.
Since quark masses are generally given in the MS scheme at 2 or 3 GeV, we obtain the corresponding renormalization constants at these scales by using the following evolution factors Our evolution function is accurate at N 3 LO [26], i.e. O(α 3 s ), and therefore, we estimate the uncertainty due to higher orders as the last known term raised to the power 4/3 (see the second error in the results of Table VI). When computing conversion factors, which are ratios of evolution functions, we add in quadrature the error coming from the numerator and the denominator. We verified that this procedure provides a good estimate of the uncertainty due to higher orders when applied to the N 2 LO conversion factors in order to estimate the (known) N 3 LO results.

V. MESON SECTOR ANALYSIS
In this Section we describe the determination of the quark masses taking as input the iso-symmetric values of the pion, kaon and D (s) -meson masses.

A. Methodology
For each ensemble, we compute the two point function where is the meson interpolating field with q f being the valence quark field of flavor f ∈ { , s, c}. By we denote the average up/down (light) quark. The correlators for the pion and kaon are the same as those used in Ref. [11]. For all mesons the two valence quarks q f and q f are always taken with opposite Wilson parameters, i.e. r f = −r f = 1, as this choice is known to suppress O(a 2 ) lattice artefacts [4,7]). For the valence mass parameters, we evaluate correlators at µ values equal to its sea counterpart, as well as at three values of the quark mass parameter µ s in the range of the strange quark masses and four values of the quark mass parameter µ c in the range of the charm quark masses. The chosen values of valence quark masses are collected in Table VIII and allow for a precise interpolation to the physical strange and charm quark masses as determined by the kaon and D-meson masses in the isosymmetric QCD. The latter ones, following the FLAG report [17], are given by From the correlator given in Eq. (26), the overlap S = | P S|J f f |0 | 2 can be extracted using an exponential fit at large time distances where m P S is the ground-state mass of a pseudoscalar (PS) meson made of the two valence quarks with flavor f and f . For maximally twisted quarks the value of the matrix element S determines the PS-meson decay constant with no need of any renormalization constant [1], from the formula The slight deviation from maximal twist of the ensemble cA211.12.48 is corrected according to Appendix C of Ref. [11]. The global energy scale is set using the isosymmetric QCD inputs (4) and data at different lattice spacings are connected by exploiting the gradient-flow (GF) quantities w 0 [28], √ t 0 [29] and t 0 /w 0 measured in lattice units. Their values have been already determined quite precisely in Ref. [11], namely 4 Nevertheless, in order to take properly into account all the correlations with the meson data the GF scales are determined again in the present analysis (see the next subsection), obtaining results well compatible with Eqs. (32)(33)(34).

B. Light quark mass
The lattice QCD data on the pion mass and decay constant are computed in an unitary setup, i.e. with µ sea = µ valence = µ , the values used in this Section are reported in the Table VIII. The lattice QCD data on the pion mass and decay constant are analyzed relying on SU(2) chiral perturbation theory (ChPT) using the formulae (m π w 0 ) 2 = 2(Bw 0 )(m w 0 ) 1 + ξ log ξ + P 1 ξ + P 2 a 2 /w 2 where the variable ξ = 2Bm /(16π 2 f 2 ) is related to the quark renormalized mass m = µ /Z P . The parameters P 1 and P 3 are related to the low-energy constants¯ 3 and¯ 4 by The quantities K F SE M 2 and K F SE f represent the finite size effects (FSE) on the squared pion mass and the pion decay constant, respectively. In Ref. [11] it was shown that SU(2) ChPT at NLO [31] adequately describes our lattice data, once discretization effects proportional both to a 2 and to a 2 m are included in f π (see Eq. (36)), while in m π the leading lattice artefact is already directly proportional to a 2 m (see Eq. (35)). The fit parameters are Bw 0 ,¯ 3 , P 2 , f w 0 ,¯ 4 , P 4 and P 5 . We repeat the fit procedure adopting the values of the renormalization constant Z P determined using the methods M1a, M1b, M2a, M2b given in Table VII.
To estimate possible systematics due to the scale setting and the chiral extrapolation, we repeat the analysis using: • the ratio t 0 /w 0 to set the scale; • the GF scale √ t 0 to set the scale; • only a combination of two lattice spacing 5 , namely  ) and (36) and ZP for the M2b method. Different colored bands correspond to different lattice spacings (red for the A ensembles, blue for the B and green for the C). The grey band is the extrapolation to the continuum limit. Note that for w0fπ discretization effects proportional both to a 2 and to a 2 m are visible (see text).  We need now to average the results coming from the different analyses collected in Table IX. To this end we adopt a simple generalization of Eq. (28) of Ref. [12]. For a given observable x we assume that its probability distribution f (x) is given by where w i are weights to be specified and f i (x) are the probability distributions corresponding to the individual analyses (labelled with i = 1, 2, ..., N ). It is not necessary to specify the form of the individual distributions. It suffices to know that x i and σ i are the mean value and standard deviation of the distribution f i (x).
Thus, using Eq. (38) we can represent the combination of the N results of the various analyses in the form Eq. (41) represents the square of a "statistical" error given by the weighted average of the individual variances, while Eq. (42) corresponds to the square of a "systematic" error related to the spread among the results of the different analyses. The total error σ is given by the sum in quadrature of σ stat and σ syst . Given the limited number of data points, we refrain in using the values of χ 2 , shown in Table IX, as a quantitative estimate of the quality of the various fits. Instead, since the results of Table IX suggest the dominance of the "statistical" uncertainties over the "systematic" ones, a reasonable choice for the weights w i is w i ∝ 1/σ 2 i , namely Thus to obtain the value of m ud we combine the values of Table IX  In this Section we present our determination of the strange quark mass m s . For the valence mass parameters, we evaluate correlators for µ values equal to its sea counterpart, as well as at three values of the quark mass parameter µ s in the range of the strange quark masses shown in Table X.  Table VIII. For each ensemble we perform a linear interpolation of the kaon mass to three reference values of (m s w 0 ) ref = 0.064, 0.080, 0.095 using the Ansatz A similar interpolation is also performed for the other GF scales t 0 /w 0 and √ t 0 . Then for each value of (m s w 0 ) ref we extrapolate to the continuum limit and to the isosymmentric QCD point using the value of m = m ud determined in the previous Section and our best fit to the data for m K according to the Ansatz At NLO order of SU (2) ChPT there are no finite volume effects on the kaon mass, and in Ref. [11] it has been shown that the lattice QCD data on the kaon masses agree with this prediction. The fit parameters in Eqs. (49) are P 0 , P 1 , P 2 , P 3 , while the LO low-energy constants f and B are taken from our pion sector fit. The quality of the resulting fit to Eq. (49) is shown in Fig. 10 as an example for the specific determination of Z P . Other determinations yield similar results. The last step of the analysis is an interpolation using Eq. (48) to find the value of m s that reproduces m isoQCD K = 494.2(3) MeV given in Eq. (27). As in the case of the pion, to estimate the systematic errors related to the scale setting, in the chiral extrapolation and the continuum limit we repeat the analyses using two different GF scales, excluding the ensembles with pion mass larger than 190 MeV and the term proportional to P 2 in Eq. (49), and with only pairs of values of the lattice spacing. The results are shown in Table XI.
We use the same procedure as for m ud to obtain the mean value, statistical and systematic errors for m s using Eqs. (40)(41)(42)

D. Charm quark mass
In this Section we present our determination of the mass of the charm quark obtained by analyzing D-and D smeson masses, following a strategy similar to the one presented for the determination of m s . For the valence mass parameters, we evaluate correlators for µ values equal to its sea counterpart, as well as at four values of the quark mass parameter µ c in the range of the charm mass. In the case of the D s meson we also use three values of the quark mass parameter µ s equal to the values used in the kaon analysis (see Table X). The values for the D-meson masses are given in Table XII, while the ones for the D s -meson in Table XIII. The D-and D s -meson correlators are computed using both smeared and local interpolating fields. Using the four combinations of smeared-smeared, smeared-local  Table X are given in the top most panel, using t0/w0 in the second panel, using √ t0 in the third panel, using w0 and limiting mπ < 190 MeV in the fourth panel, using w0 and only the two coarser lattice spacings in the fifth panel, using w0 and only the coarser and finest lattice spacings in the sixth panel and using w0 and only the two finest lattice spacings in the last panel.
To determine the ratio ms/m ud we use the values of m ud from Table IX. and local-local correlators we construct a 2 × 2 matrix and perform a GEVP analysis [32] to extract the mass of the D-and D s -mesons 6 . We employ Jacobi smearing for the quark fields [33], combined with APE smearing of the gauge links [34] used in the Jacobi smearing function. The values of m D and m Ds used in this analysis are reported in Table XII and Table XIII.
Analogously to the case of the analysis for the strange quark mass determination, we interpolate the D and D s For the D s meson we also perform an interpolation to the the mass m s given in Table XI . At each of the reference charm quark masses, we extrapolate to the continuum and to the isospin-symmetric QCD (isoQCD) light quark mass m = m ud using the following polynomials in m m D = P 0 + P 1 m w 0 + P 2 a 2 /w 2 0 , (53) where P j and P s j , j = 0, 1, 2 are fit parameters. For each reference mass (m c w 0 ) ref we compute the masses m D and m Ds in the continuum limit at the isoQCD value of m = m ud given in Table IX. We then perform an interpolation in m c with the Ansatz given in Eq. (52) to compute the value of m c that reproduces the isoQCD masses of the D of D s mesons, given in Eqs. (28) and (29). We note that the analysis is done separately using either the D or the D s meson.
The resulting fits to Eqs. (53) and (54) for the D and D s mesons are shown in Fig. 11 for the case where Z P is determined from the M2b method.
The values in physical units that we obtain for the m c mass are shown in Table XIV ) using the value of ZP extracted from the M2b method and using w0 to set the scale. The notation is the same as that of Fig. 10.
We combine all the values given in Table XIV m c m s = 11.43(9) stat (10) syst = 11.43 (13) , where the charm quark mass is given in the MS at 3 GeV.

VI. BARYON SECTOR ANALYSIS
In the baryon sector, we use the nucleon and pion masses to set the scale and determine the light quark mass. We use the Ω − (sss) and the Λ c (udc) masses to determine, respectively, the strange and charm quark masses. The range of validity of ChPT in the baryon sector is more limited as compared to that in the pion sector and thus we restrict ourselves to using pion masses up to 260 MeV.

A. Methodology
In order to compute the baryon masses we construct the following two-point correlation functions at zero momentum, defined as  Table XII are given in the top most panel, using t0/w0 in the second panel, using √ t0 in the third panel, using w0 and limiting mπ < 190 MeV in the fourth panel, using w0 and only the two coarser lattice spacings in the fifth panel, using w0 and only the coarser and finest lattice spacings in the sixth panel and using w0 and only the two finest lattice spacings in the last panel.
To determine the ratio mc/ms we use the values of ms from Table XI. where J B(qq q ) is the interpolating operator for the baryon B(q, q , q ) with q, q and q ∈ {l, s, c}. In this work, we increase statistics by considering both 1 2 (1 ± γ 0 ) projectors. For the interpolating fields of the nucleon, the Ω and the Λ c we take, respectively where latin indices refer to colour, abc is the antisymmetric tensor and C is the charge conjugation matrix. In order to suppress contributions from excited states we apply Gaussian smearing to each quark field q( x, t). The smeared quark field is given by q smear ( x, t) = y F ( x, y; U (t))q( y, t), where F is the gauge invariant smearing function F ( x, y; U (t)) = (1 + αH) n ( x, y; U (t)), constructed from the hopping matrix understood as a matrix in coordinate, color and spin space, In addition, we apply APE smearing to the spatial links that enter the hopping matrix H. Different Gaussian smearing is applied to the light and strange quarks. The parameters of the Gaussian and APE smearing for each ensemble for the light and strange quarks are given in Table XV. The charm quark interpolating fields are not smeared.
can be written using the spectral decomposition of the two-point correlators as where ∆ j is the mass difference of the j-th excited state with respect to the ground state mass m B . We consider one-, two-and three-state fits by taking K = 0, 1, 2 in Eq. (62). This allows us to check the consistency in our determination of the ground state mass m B . Since statistical errors are larger for baryons as compared to those of mesons and grow rapidly with increasing time separation t, it is important to identify the ground state for as small time separation as possible, so that we can be confident that excited are sufficiently suppressed. Our procedure for identifying m B is as follows: 1. We keep the upper time used in the fit constant. The upper time is chosen so that statistical errors are reliably evaluated.
2. We fit the effective mass keeping two excited states, i.e., we take K = 2 in Eq. (62) and vary the lower time used in the fit t 3st low /a from one to three. We choose the parameters of the fit that has the smallest t 3st low for which χ 2 /d.o.f. 1. This determines m 3st B . 3. We then fit the effective mass including one excited state, i.e., we set K = 1 in Eq. (62) and vary t 2st low for t 2st low > t 3st low until the extracted mass m 2st B satisfies the criterion |m 2st 4. Having determined m 2st B we make a single state fit, i.e., we set K = 0 in Eq. (62) and vary the lower value of t. We choose t 1st low > t 2st low and take the smallest value that satisfies |m 1st We illustrate our analysis for the extraction of the masses by giving representative examples for the nucleon, Ω − and Λ c . In all cases we use correlation functions with smeared sources and for the ensembles listed in Table XVI. In Fig. 12 we show an example of the results obtained using the nucleon correlators for the cC211.06.80 ensemble and in Table XVI we give the number of configurations and source positions used, the fit ranges for the one-, two-and three-state fit, as well as the extracted nucleon mass and the χ 2 /d.o.f.. As can be seen, the mass of the first excited state converges to a value compatible with the mass of the Roper for the physical point ensembles. The nucleon-pion state, although it has lower energy, it is volume suppressed. In our two-state fits to extract the energy of the first excited state we find that the coefficient of the second exponential compared to that of the ground state is of order 1. This is to be contrasted with the chiral perturbation theory analysis of Ref. [35] which predicts a few percent for two-particle states. This indicates that the contribution of two-particles is suppressed. TABLE XVI. We give for each ensemble the resulting values of amN , when using one-state (fourth main column), two-state (fifth main column) and three-state (sixth main column) fits,χ 2 ≡ χ 2 /d.o.f. is the reduced χ 2 , n conf is the number of configurations analyzed, nsrcs is the number of two-point functions generated per configuration at different source positions and [t low , tmax] is the fitting range. We also show the values for the pion mass computed on the same statistics, amπ, noticing that they are compatible with those given in Table I. We analyze in a similar way the effective mass defined by the Ω correlator given in Eq. (57). In Fig. 13 we show an example of the effective mass m eff Ω for the cB211.072.64 ensemble at µ s = 0.017, 0.0195 and 0.022. As can be seen, we obtain accurate results that allow us to perform a fit including up to the second excited state. We fix the maximum time for these fits to be t max /a = 34. The convergence of the effective mass for Ω − as we vary t low is demonstrated when using one-, two-and three-state fits. In a similar manner, the convergence of the first excited energy E 1 Ω is demonstrated by varying t low . We employ the criterion described above to choose the value of m Ω from the one-state fit at each µ s . We note that for all the three values of µ s we find the same t low . The masses extracted are given in Table XVII, where we also quote the reduced-χ 2 ,χ 2 ≡ χ 2 /d.o.f., of the various fits.
The analysis of the two-point correlator for the Λ c proceeds in an analogous manner. We illustrate the results for the cB211.072.64 ensemble in Fig. 14 for two different values of the charm mass parameter µ c . From the study of the Ω − mass we find that there is strong correlation among the data for the three values of µ s as demonstrated in Fig. 16 and, thus, for Λ c we opt to use two different values of µ c in the interpolation. Since Λ c is heavier and decays faster, a three-state fit is not possible and we limit ourselves to comparing one-and two-state fits. The masses extracted, the  (5)  statistics used and the value ofχ 2 are given in Table XVIII In what follows we will use the values of m N , m Ω and m Λc extracted from the one-state fit given in Tables XVI,  XVII and XVIII, respectively, to determine the light, strange and charm quark masses. In order to estimate the systematic error due to the fit range we will also use the values for the masses extracted from the one-state fit at t low /a + 1.

B. Light quark mass
We use the ChPT expression of Eq. (5) to extrapolate to the physical point. To one-loop order in ChPT (up to which the nucleon mass is expanded in Eq. (5)) we can substitute the pion mass by m 2 π = 2Bm ud (1 + c 2 a 2 ) to obtain the expansion consistent to the order we are working and including O(a 2 ) effects both in the pion expansion in Eq. (63) with the coefficient c 2 and in the nucleon expansion in Eq. (5). We thus have two fit parameters, B and c 2 , while the lattice spacings, m 0 N and c 1 are determined from Eqs. (5)- (7). The fit procedure is performed for the four values of the renormalization constant from Table VII, checking for consistency and estimating systematic effects in the determination of Z P . Since the values we obtain using the different methods of extracting Z P are in very good agreement we average over them. The statistical error in the lattice spacing is taken into account in the jackknife analysis. The fit results are reported in Table XIX and depicted in Fig. 15 FIG. 14. The same as in Fig. 12 but for the case of Λc using the cB211.072.64 ensemble for the two values of µc given in the figure legend.
using the nucleon and pion mass, given in the MS scheme at 2 GeV, is obtained by averaging the values in Table XIX. The systematic error is computed as in Eq. (42) but in the sum we only take into account the mean values not included in the computation of the average. Namely, the systematic error reflects the choice of the fitting range estimated by increasing t low by one unit and the sensitivity due to the chiral extrapolation estimated by using ensembles with pion mass smaller than 190 MeV. We will follow this procedure also for the computation of the systematic errors also for the strange and charm quark masses.

C. Strange and charm quark masses
We determine the strange and charm quark masses using the experimental value of the Ω (sss) and Λ c (udc) masses and the lattice spacings determined from the nucleon mass. Namely, we use m  [36]. We use the renormalization constants Z P given in Table VII.
We parametrize the Ω − and Λ c mass dependence on the strange and charm quark mass by expanding aroundm s andm c , that we chose to be in the same ballpark of the physical quark masses. In particular we usem s = 95 MeV andm c = 1.2 GeV, and we interpolate around these reference points using We first discuss our procedure for the determination of the strange quark mass from the Ω mass. We then apply the same procedure for determining m c using the Λ c mass.

Strange quark mass
The knowledge of the Ω − mass at three values of the valence strange quark mass parameter µ s allows us to determine the Ω mass as a function of µ s using the linear Ansatz of Eq. (65). We show a representative example of the resulting fit in Fig. 16 for the ensemble cB211.072.64. The same analysis is carried out for all the ensembles listed in Table XVII, where we give the values of A Ω and B Ω , defined in Eq. (65).
We employ two methods to determine m s : In method I we perform a chiral and continuum extrapolation of the A Ω and B Ω parameters separately. Namely, we expand to leading order in ChPT and include O(a 2 ) cut-off effects as follows B Ω (a, m 2 and we limit ourselves to ensembles with m π < 260MeV so that these leading order expression are reliable. In Fig. 17 we illustrate the chiral and continuum extrapolation for the parameters A Ω and B Ω using the value of Z P from method M1a (see Table VII). We note that the values of A Ω and B Ω using the cB211.025.32 and cB211.025.48 ensembles are compatible, demonstrating that finite size effects are small. Using the values of the parameters A Ω and B Ω at the physical pion mass and continuum limit we can extract the strange quark mass in the continuum limit and at the physical pion mass from In method II we adopt an iterative strategy: Namely, we start by fixing a value of the renormalized strange quark mass m s in physical units for all the ensembles. We use Eq. (65) to interpolate to the given m s . We then extrapolate to the continuum limit and physical point using the ChPT result We iterate this procedure changing the value of m s until the resulting value of m Ω given in Eq. (70) at the physical point and continuum limit matches the physical value m (phys.) Ω . In Fig. 18 we illustrate the analysis. The results for the renormalized strange quark mass in the MS scheme at 2 GeV are provided in Table XX using the values of Z P given in Table VII. We compare the different values by plotting them in Fig. 19. As can be seen, despite the different values of Z P at finite lattice spacing, in the continuum limit we obtain very good agreement among different estimates of m s . A similar agreement is also obtained between methods I and II discussed in this Section for the determination of m s . Since the error on the lattice spacing cannot be taken into account in a jackknife    analysis because we used different statistics, we estimate the change in the value of m s by varying the lattice spacing by a standard deviation. As can be seen in Fig. 19 this gives a very small change compared to the statistical error. By increasing t low of the one-state fit by one lattice unit gives an estimate of the systematic error in the extraction of the mass of the Ω. The change in the value of m s is well within the statistical error, as can be seen in Fig. 19. We, thus, average over all the values obtained using the four different determination of Z P and analysis methods I and II (values are given in Table XX and Table XXI of the fitting range by letting t low /a → t low /a + 1 and systematics due to the chiral extrapolation by using ensembles with m π < 190 MeV. Since the error in the lattice spacing cannot be included in the jackknife analysis due to using different statistics for the Ω, we include an additional term in Eq. (42)  give with open symbols the determination using method I to extract ms and with filled symbols using method II. Circles show the results using the mass of the Ω from the one-state fit given in Table XVII. Triangles quantify systematic errors due to the selection of t low (right pointing triangles when t low /a + 1), errors on the lattice spacing (down pointing triangles increasing by a standards deviation the lattice spacings set by the nucleon mass) and errors due to chiral extrapolation ( up pointing triangles obtained using only ensembles with mπ < 190 MeV). The dashed line is the average over the values from method I and II and at different ZP given in Table XX and Table XXI.

Charm quark mass
We employ the same procedure for the determination of m c as described in the previous Section for m s using methods I and II. The mass of Λ c is interpolated linearly in µ c using Eq. (66) within the range spanned by the two  We note that to leading one-loop in ChPT a m 3 π -term with an unknown coefficient is present. Including such term in the fit results in a coefficient consistent with zero and a χ 2 /d.o.f. = 7.2. That such a term is not supported by lattice QCD data was also found in our previous analysis using a larger set of pion masses [19] where this coefficient was found to be consistent with zero. We check that including it does not change the extracted value for m c . Thus, given the larger χ 2 we drop it from our analysis.
Results from method I are reported in Table XXII for all the values of Z P listed in Table VII. We also illustrate in Fig. 20 the chiral and continuum extrapolations of A Λc and B Λc according to Eq. (67) and Eq. (68) respectively, with Z P determined using method M1a. The determination of the values of the parameters of Eq. 72 for method II is carried out as for the case of m s and the results are reported in Table XXIII and, for the M1a case, also in Fig. 21.
The values for m c from methods I and II as well as how they change by varying the lattice spacings by a standard deviation and by the change in m Λc by increasing t low by one lattice spacing in the one-state fit are presented in Fig. 21 Again, we observe a very good agreement between the results obtained via method I and method II and among different determinations of the Z P renormalization constants. We thus average over these values and compute the systematic error in the same way as for µ s . We obtain for the charm quark mass m c in the MS scheme at 3 GeV and the ratio m c /m s the following values where the errors on the ratio are combined in quadrature.

VII. CONCLUSIONS
The focus of this work is the determination of the light, strange and charm quark masses. We perform an analysis of ten N f = 2 + 1 + 1 ensembles simulated at three lattice spacings smaller than 0.1 fm and pion masses in the range from about 350 MeV to 135 MeV. Having two ensembles simulated with the physical value of the pion mass at the two smallest lattice spacings enables us to extrapolate reliably to the physical and continuum limit.
The extraction of the quark masses is done using observables from both the meson sector and the baryon sector. The iso-symmetric values of the pion, kaon and D-meson masses as well as of the pion decay constant are used for the determination of the lattice spacings and the quark masses in the meson analysis. In the baryon sector, we use as inputs the nucleon and pion masses to obtain the lattice spacing and the average light-quark mass, while the masses of the Ω − and Λ c baryons determine the strange and charm quark masses.
In Table XXIV we collect the values of the quark masses obtained in Sections V and VI for the light and strange quark masses in the MS scheme at 2 GeV and for the charm quark mass at 3 GeV. Since the isospin and electromagnetic  (7) 27.23 (10) 11.82 (16) TABLE XXIV. The renormalized quark masses determined in the meson sector (first row) and baryon sector (second row) in the MS scheme. In the third row we give the average over the values obtained in the the meson and baryon sectors, while in the last row we give the latest FLAG averages [13] for N f = 2 + 1 + 1. The light quark mass, m ud (second column), and the strange quark mass, ms (third column), are given at 2 GeV, while the charm quark mass, mc (fourth column), is given at 3 GeV. The second error of the quark masses includes a 0.5% uncertainty (added in quadrature) due to the uncertainty of the conversion of the RCs ZP to the MS scheme. In the fifth and sixth columns we give the ratios ms/m ud and mc/ms, respectively. In the meson sector the error on the ratio is determined in a jackknife analysis. In the baryon sector, since different ensembles are involved in the determination of the quark masses, the error on the ratio is propagated quadratically using the errors on each of the quark masses.
corrections to the nucleon mass are only known for the mass difference between the neutron and proton [37], in our analysis we average over the mass of the proton and neutron. This defines a QCD prescription different from that used in the meson sector. Using the values of the lattice spacing extracted in the meson sector we obtain a nucleon mass in the continuum limit a few MeV smaller than the input value m N,phys = 0.9389 GeV adopted in the baryon sector (see Section III). This results in less than a percent change in the values given in the Table XXIV, which is much smaller than our statistical errors. It is thus justifiable to average over the values obtained in the meson and baryon sectors to produce our final values. In order to perform the above average we adopt the weighted approach given in Eqs. (40)(41)(42). We assume the following weights for the quantities coming from the mesonic and the baryonic sectors, where σ stat M (B) is the first error given in the corresponding rows of Table XXIV. In this way we obtain x ± σ stat ( +σ syst,+ −σ syst,− ) , where The results are given in the last row of Table XXIV and are compared in Fig. 23 with those of the ETM analysis of Ref. [12] and the ones entering the N f = 2 + 1 + 1 averages in the latest FLAG report [13]. The latter ones are based on the results of Refs. [12,38] for the light-quark mass, Refs. [12,[38][39][40] for the strange mass, Refs. [12,19,[38][39][40] for the charm mass, Refs. [12,41,42] for the m s /m ud ratio and Refs. [12,38,40] for the m c /m s ratio. It can be seen that our results are larger by ∼ 2.5 standard deviations in the case of m ud and by ∼ 2 standard deviations in the case of m c with respect to the corresponding FLAG values. Although for the strange quark mass our result coming from the meson sector is larger by ∼ 3 standard deviations, our averaged result is consistent with the FLAG one within our final uncertainty. A good agreement is observed for the mass ratios m s /m ud and m c /m s . We do not believe that these differences can be ascribed to possible uncontrolled effects on the mass renormalization constant 1/Z P . Indeed, the detailed analysis carried out in this work concerning the pion pole subtraction and the residual hadronic contaminations in the RI-MOM determination of the renormalization constant Z P leaves little room for any significant leftover contribution from these terms. Our findings point to the fact that hadronic contaminations are controlled at the level of few per mil. Therefore, we do not consider plausible that the observed tension with the FLAG values may be related to uncontrolled hadronic contaminations on the mass renormalization constant. In this respect, we are considering the possibility of repeating the determination of the quark masses using the same ETM gauge ensembles adopted in this work, but evaluating the mass renormalization in a different scheme, like RI-SMOM, while keeping the same level of control of the hadronic contaminations achieved in this work.
Our final results for the light, strange and charm quark masses as well for the mass ratios m s /m ud and m c /m s are consistent with our previous analysis of Ref. [12] (see also Fig. 23), which was based on Wilson twisted-mass fermions  [12] (blue squares) and the N f = 2 + 1 + 1 averages given in the last FLAG report [13] (black circles). The shorter error bars take into account the statistical error only, while the larger represent the total error, obtained by summing in quadrature the statistical and the systematic errors.
far from the physical pion point. The overall uncertainties for the light and charm quark masses are reduced by a factor of ∼ 1.7−2.0, while in the case of the strange quark mass the uncertainty is almost unchanged, partly due to the difference between the mean values obtained in the meson and baryon sector, which is added to the systematic error. This is also reflected in the two ratios. With respect to the quark mass analysis of Ref. [12] the main improvements are: i) a better control of the chiral extrapolation thanks to gauge ensembles produced close to the physical pion point; ii) a better control of hadronic contaminations in the calculations of the mass renormalization constant; iii) the use of both mesonic and baryonic quantities, which requires simulations of different correlation functions. For all the three masses the contribution from lattice systematics is important, in particular in the case of m s . Our plan is to add at least one further gauge ensemble at a fourth finer value of the lattice spacing at the physical point. This will allow a tightly controlled chiral and continuum extrapolations in both the meson and baryon sectors.