Ratio of kaon and pion leptonic decay constants with $N_f = 2 + 1 + 1$ Wilson-clover twisted-mass fermions

We present a determination of the ratio of kaon and pion leptonic decay constants in isosymmetric QCD (isoQCD), $f_K / f_\pi$, making use of the gauge ensembles produced by the Extended Twisted Mass Collaboration (ETMC) with $N_f = 2 + 1 + 1$ flavors of Wilson-clover twisted-mass quarks, including configurations close to the physical point for all dynamical flavors. The simulations are carried out at three values of the lattice spacing ranging from $\sim 0.068$ to $\sim 0.092$ fm with linear lattice size up to $L \sim 5.5$~fm. The scale is set by the PDG value of the pion decay constant, $f_\pi^{isoQCD} = 130.4~(2)$ MeV, at the isoQCD pion point, $M_\pi^{isoQCD} = 135.0~(2)$ MeV, obtaining for the gradient-flow (GF) scales the values $w_0 = 0.17383~(63)$ fm, $\sqrt{t_0} = 0.14436~(61)$ fm and $t_0 / w_0 = 0.11969~(62)$ fm. The data are analyzed within the framework of SU(2) Chiral Perturbation Theory (ChPT) without resorting to the use of renormalized quark masses. At the isoQCD kaon point $M_K^{isoQCD} = 494.2~(4)$ MeV we get $(f_K / f_\pi)^{isoQCD} = 1.1995~(44)$, where the error includes both statistical and systematic uncertainties. Implications for the Cabibbo-Kobayashi-Maskawa (CKM) matrix element $|V_{us}|$ and for the first-row CKM unitarity are discussed.


I. INTRODUCTION
The leptonic decay constants of charged pseudoscalar (P) mesons are the crucial hadronic ingredients necessary for obtaining precise information on the Cabibbo- where [γ] stands for the contribution of virtual and real photons. On the theoretical side, within the SM the above ratio is given by where V ud and V us are the relevant CKM entries, M π(K) is the charged pion(kaon) mass, m µ is the muon mass and δR Kπ represents the isospin breaking (IB) corrections due both to the mass difference (m d − m u ) between the light u-and d-quarks and to the quark electric charges. Finally, in Eq.
(2) f K /f π is the ratio of kaon and pion leptonic decay constants defined in isosymmetric QCD (isoQCD), i.e. with m u = m d and zero quark electric charges.
Recently [4,5] the IB correction δR Kπ has been determined using a non-perturbative approach, based on first principles, through QCD+QED simulations on the lattice, obtaining δR Kπ = −0.0126 (14). From Eq. (1) one has V us V ud f K f π = 0.27683 (29) exp (20) th = 0.27683 (35) , which corresponds to an accuracy of 0.13%. As well known [6], the IB correction δR Kπ and the isoQCD ratio f K /f π separately depend on the prescription used to define what is meant by isoQCD, while only the product (f K /f π ) √ 1 + δR Kπ is independent on such prescription. The hadronic prescription adopted in Refs. [4,5] corresponds to while the quantity (m d − m u ) is obtained from the difference between the experimental charged and neutral kaon masses. The physical pion and kaon masses (4-5) are consistent with those recommended by FLAG-3 [7], and the pion decay constant (6), derived according to Ref. [8] adopting the value of the CKM entry |V ud | from Ref. [9], is used to set the lattice scale 1 .
In this work we present our determination of the leptonic decay constant ratio f K /f π at the physical isoQCD point given by Eqs. (4-6), evaluated using the ETMC gauge ensembles produced with N f = 2 + 1 + 1 flavors of Wilson Clover twisted-mass quarks, including configurations close to the physical point for all dynamical flavors [11,12]. The lattice data will be analyzed within the framework of SU(2) Chiral Perturbation Theory (ChPT) without making use of renormalized quark masses 2 . By means of the pion data we determine the gradient-flow (GF) scales w 0 [13], √ t 0 [14] and t 0 /w 0 adopting the physical value (6) at the pion point (4) to set the lattice scale, obtaining w 0 = 0.17383 (63) fm , √ t 0 = 0.14436 (61) fm , 1 In Ref.
[5] it has been shown that within the precision of the lattice simulations the prescription given by Eqs. (4-6) is equivalent to the Gasser-Rusetsky-Scimemi (GRS) scheme [10], where the renormalized quark masses and coupling constant in a given short-distance scheme (viz. the MS scheme) and at a given scale (viz. 2 GeV) are equal in the full QCD+QED and isoQCD theories. where the error includes both statistical and systematic uncertainties. Our findings (7)(8) are a little larger than the MILC results [15] w 0 = 0.1714 +15 −12 fm and √ t 0 = 0.1416 +8 −5 fm as well as the HPQCD results [16] w 0 = 0.1715 (9) fm and √ t 0 = 0.1420 (8) fm, both obtained using the hadronic value (6) to set the lattice scale. Within 1.5 standard deviations our result (7) is consistent with the recent, precise BMW determination w 0 = 0.17236 (70) fm, obtained in Ref. [17] using the Ω − -baryon mass to set the lattice scale.
As for the ratio f K /f π we determine its value at the physical isoQCD point (4-6) and in the continuum and infinite volume limits, obtaining where again the error includes both statistical and systematic uncertainties.
The IB correction δR Kπ = −0.0126 (14), determined in Refs. [4, 5] and adopted in Eqs. (2-3), stems from the sum of a QED and a strong IB terms, which are both prescription dependent as well as their sum and the isoQCD value (10). Within the GRS prescription (see footnote 1) they are equal respectively to −0.0062 (12) and −0.0064 (7).
Thus, for the ratio of kaon and pion leptonic decay constant including strong IB effects (which we remind is prescription dependent) we get f K + f π + = 1.1957 (44) .
The plan of the paper is as follows.
In Section II some details of the ETMC gauge ensembles and of the simulations are illustrated, while a more complete description is provided in Appendix A. For each gauge ensemble the pion mass and decay constant are extracted from the relevant twopoint correlation functions using a single exponential fit in the appropriate regions of large time distances. Alternatively, in Appendix B the extraction of the ground-state properties is performed through the multiple exponential procedure of Ref. [23]. For one gauge ensemble (cA211.12.48), because of a small deviation from maximal twist, the mass and the decay constant are corrected as described in Appendix C. In Section III the SU(2) ChPT predictions at next-to-leading order (NLO) for the pion decay constant f π , including finite volume effects (FVEs), are presented. For the ensembles cB211. 25.XX, sharing the same light-quark mass and lattice spacing and differing only for the lattice size L, the FVEs are investigated using both the NLO and the resummed NNLO formulae of Ref. [24]. In Section IV, adopting the physical value (6) at the pion point (4), we perform two determinations of the GF scale w 0 using the data for either f π or the quantity X π ≡ (f π M 4 π ) 1/5 , which is found to be less affected by statistical and systematic errors. The two determinations of w 0 agree very nicely, but the one based on the quantity X π turns out to be more precise by a factor of ≈ 2.5. In the same way the other two GF scales √ t 0 and t 0 /w 0 are determined in Appendix D, where our calculations of the relative GF scales w 0 /a, √ t 0 /a and t 0 /(w 0 a) at the physical pion point are also described. In Section V we analyze the data for the decay constant ratio f K /f π using SU(2) ChPT.
In Section VI the implications for V us and the first-row CKM unitarity are discussed.
Finally, our conclusions are collected in Section VII .

II. ETMC ENSEMBLES
In this work we make use of the gauge ensembles produced recently by ETMC in isoQCD with N f = 2 + 1 + 1 flavors of Wilson-clover twisted-mass quarks and described in Refs. [11,12]. The gluon action is the improved Iwasaki one [25], while the fermionic action includes a Clover term [26] with a coefficient fixed by its estimate in one-loop tadpole boosted perturbation theory [27]. Its inclusion turns out to be very beneficial for reducing cutoff effects, in particular on the neutral pion mass, thereby making numerically stable simulations close to the physical pion point [11].
The Wilson mass counterterms of the two degenerate light-quarks as well as of the strange and charm quarks are chosen in order to guarantee automatic O(a)improvement [28,29]. The masses of the strange and charm sea quarks are tuned to their physical values for each ensemble [11,12]. For the valence strange and charm sectors, a mixed action setup employing Osterwalder-Seiler fermions [30], with the same critical mass as determined in the unitary setup, has been adopted in order to avoid any undesired strange-charm quark mixing (through cutoff effects) and to preserve the automatic O(a)-improvement of physical observables [31].
Some properties of the ETMC ensembles, which are relevant for this work, are collected in Table I respectively to a lattice spacing equal to a ≈ 0.082 fm and a ≈ 0.069 fm, the pion mass is simulated quite close to the physical isoQCD value (4).
For each ensemble we compute the pion correlator given by where  (7)). In the last column the number of gauge configurations analyzed for each ensemble is presented.
this way the squared pion mass differs from its continuum counterpart only by terms of O(a 2 µ ) [28,29].
At large time distances one has so that the pion mass M π and the matrix element Z π = | π|q γ 5 q |0 | 2 can be extracted from the exponential fit given in the r.h.s. of Eq. (15).
For maximally twisted fermions the value of Z π determines the pion decay constant f π without the need of the knowledge of any renormalization constant [28,32], namely The time intervals [t min , t max ] adopted for the fit (15) of the pion correlation function (14) as well as the extracted values of the pion mass and decay constant in lattice units are collected in Table II. As anticipated in the Introduction, in this work we will make also use of the data for the quantity X π defined as which turns out to be less affected by lattice artifacts (see below Fig. 1 and later Section III C). The values of X π in lattice units are shown in the last column of Table II.
The statistical errors of the lattice data lie in the range 0.  An alternative way to extract the pion mass and decay constant is the ODE procedure of Ref. [23]. The results obtained by applying this method to the pion correlation function (14) are collected in Appendix B and found to be totally consistent with the findings of the single exponential fit (15) of Table II.
In the case of the ensemble cA211.12.48 due to a small deviation from maximal twist a correction needs to be applied. According to Appendix C the squared pion mass is left uncorrected, while for the pion decay constant f π we use the following formula where m P CAC is the bare untwisted PCAC mass, Z A is the renormalization constant of the axial current and µ is the bare twisted mass of the light valence quarks. For the ensemble cA211.12.48 one has Z A ≈ 0.75 and m P CAC /µ −0.21 (5) [12].
The statistical accuracy of the correlator (14) is significantly improved by using the so-called "one-end" stochastic method [33], which includes spatial stochastic sources at a single time slice randomly chosen. Statistical errors are evaluated using the jackknife procedure.
The results obtained for the pion decay constant w 0 f π and for the quantity w 0 X π (see Eq. (17)), are shown in Fig. 1 for all the gauge ensembles. By comparing the results corresponding to the ensembles cB211.25.XX the FVEs are clearly visible in the case of f π , while they are almost absent in the case of X π . Moreover, also discretization effects in X π turn out to be smaller than those present in f π .
III. THE PION DECAY CONSTANT f π WITHIN SU(2) CHPT Within SU(2) ChPT [34] the pion decay constant f π is given at NLO by where For the squared pion mass one has at NLO where the coefficient C 1 is related to the NLO LEC¯ phys 3 bȳ phys 3 A. Finite volume effects within NLO SU(2) ChPT The structure of FVEs on the pion decay constant can be studied using SU(2) ChPT [34]. At NLO FVEs come entirely from the discretized sum over periodic momenta of the loop contributions. For a finite spatial volume V = L 3 one has where f π (L → ∞) is given by Eq. (20). The correction term ∆ π F V E (L) can be obtained from the chiral log in Eq. (20) via the following replacement where λ ≡ √ 2Bm L = √ ξ 4πf L and Thus, within NLO SU(2) ChPT the quantity ∆ π F V E (L) is given by In the case of the squared pion mass one gets where M 2 π (L → ∞) is given by Eq. (23).

C
(1) with i being NLO LECs that have a logarithmic pion mass dependence Finally, in Eqs. Mπ are defined in the Appendix A of Ref. [24], but useful approximate analytic formulae are given by [24] S (4) with The expansion variable ξ π is defined as [24] Different choices of the expansion variable are possible: one can replace f π with the LO LEC f and/or replace M 2 π with 2Bm (and correspondingly M π L with √ 2Bm L in the arguments of the functions g 1 and g 2 ). At NLO (i.e., for the GL formula) the above changes are equivalent, since any difference represents a NNLO effect. Instead, in the CDH formula additional terms appear at NNLO, which can be found in Ref. [35]. Here we consider only the alternative definition which requires the addition to the r.h.s of Eq. (31) of the term f π (∞) 8ξ 2 π 4 g 1 (M π L) and to the r.h.s. of Eq. (32) of the term M π (∞) −2ξ 2 π 4 g 1 (M π L) .
The GL formula corresponds to Eqs.  Table I, evaluated using the GL and CDH formulae and assuming respectively the two definitions (42) and (43) for the expansion variable ξ π .
In the case of CDH formula we adopt the following values of the NLO LECs: It can be seen that the GL formula applied to both the pion mass and decay constant works quite well for M π L 3, particularly in the case of the definition (43) of the expansion variable ξ π . The above condition is satisfied by all ETMC ensembles of Table I except the ensemble cB211.25.24. C. FVEs for the quantity X π The interesting feature of the quantity X π , given by Eq. (17), is the absence of NLO chiral logs in its SU(2) ChPT expansion (see Eqs. (20) and (23) Let's now apply the SU(2) ChPT predictions for interpolating the pion data to the physical pion mass and for extrapolating them to the continuum and infinite volume limits. The goal is to determine the GF scale w 0 adopting the physical value (6) at the pion point (4) without resorting to the use of the renormalized light-quark mass. In the next two subsections we separately analyze the pion decay constant f π and the quantity X π , respectively.
A. Determination of w 0 using the data for f π Using the simulated values aM π and af π in lattice units we evaluate the expansion variable ξ π , defined (from now on) as which depends on neither w 0 nor w 0 /a. Then, for each gauge ensemble we calculate the and we re-express the quantity ξ (see Eq. (21)) in terms of the pion mass in the infinite volume limit (see Eq. (30)) where only the knowledge of w 0 /a is required to calculate the pion mass in units of w 0 and the free parameter becomes w 0 f .
We correct the data of the pion decay constant w 0 f π (L) for FVEs (see Eq. (25)), Analogously, for the pion mass w 0 M π (L) one has The data for w 0 f π (L → ∞) are fitted in terms of the variable ξ (see Eq. (46)) using the following functional form where with respect to a pure NLO ansatz we have added a possible higher-order term quadratic in ξ as well as discretization effects proportional to a 2 and a 2 M 2 π . The free parameters appearing in Eq. (49) are w 0 f , A 1 , A 2 , D 0 , D 1 and their values are obtained from a standard χ 2 -minimization. From the value of w 0 f the GF scale w 0 can be determined as follows. Let us consider the physical value of the variable (44), Using Eq. (49) in the continuum limit the physical value of the variable (46), namely can be obtained by solving the relation In this way the value of the LEC f in physical units is given by and, therefore, w 0 can be determined using the value of w 0 f .
We start by considering a pure NLO fit, i.e. A 2 = 0, including only the discretization effect proportional to a 2 , i.e. D 1 = 0 in Eq. (49), and we apply it to all pion data up to M π 350 MeV. The discretization coefficient D 0 turns out to be quite small, cretization effects proportional to a 2 M 2 π (i.e. D 1 = 0), but limited to pion masses below 190 MeV (4 data points and 3 parameters). In this case one gets w 0 = 0.1736 (15) where () stat+f it incorporates the uncertainties induced by both the statistical errors and the fitting procedure itself, () syst corresponds to the uncertainty related to chiral interpolation, discretization and finite-volume effects, while the last error is their sum in quadrature. More precisely, the various systematic uncertainties are estimated by considering the results obtained with A 2 = 0 or A 2 = 0 in the case of the chiral extrapolation, with D 1 = 0 or D 1 = 0 (but limited to M π < 190 MeV) for the discretization effects and with κ F V E = 1 or κ F V E = 1 for the FVEs.
B. Determination of the GF scale w 0 using the data for X π In this Section we illustrate the results of the analysis of the lattice data for the quantity w 0 X π adopting the following fitting function where the variable ξ is defined by Eq. (46), given in terms of the pion mass corrected for the FVEs using the GL formula (45), and the coefficient A 1 is related to the LEC¯ phys We have performed several fits similar to those adopted in Section IV A and the corresponding results are collected in Table III. The quality of the NLO fit with D 1 = 0 (47)  Results for w 0 obtained by fitting the lattice data for w 0 X π using Eq. (56) and adopting the isoQCD values (4) and (6) for fixing the lattice scale at the physical pion point.
is illustrated in Fig. 6, where it is also clearly visible the presence of discretization effects proportional to a 2 M 2 π , as already observed in the case of w 0 f π (see Fig. 5). We stress that for both quantities, w 0 f π and w 0 X π , the inclusion of a discretization term proportional to a 2 M 2 π leads to higher values of w 0 . This result is reassuringly confirmed also by a NLO fit without such a discretization term (i.e. D 1 = 0), but limited to pion masses below 190 MeV (see the fourth row of Table III).
By averaging the last four results of Table III one has which nicely agree with the corresponding results obtained by the analysis of f π given by Eqs. (53)(54)(55). Note that the determination of w 0 obtained using X π is more precise than the one from f π by a factor equal to ≈ 2.5.
Our result (57) is slightly larger than both the MILC result w 0 = 0.1714 +15 −12 fm from Ref. [15] and the HPQCD result w 0 = 0.1715 (9) fm from Ref. [16], obtained using the value (6) to set the lattice scale. Within 1.5 standard deviations it is consistent with the recent, precise BMW determination w 0 = 0.17236 (70), obtained in Ref. [17] using the Ω − -baryon mass to set the lattice scale. Furthermore, the difference with the recent result w 0 = 0.1709(11) fm, obtained in Ref. [18] using the Ω − -baryon mass to set the lattice scale, is within ∼ 2 standard deviations.
In Appendix D 2 the procedure used in this Section to determine the GF scale w 0 is repeated in the case of the scales √ t 0 and t 0 /w 0 , obtaining and Our finding (60) is larger than the MILC result √ t 0 = 0.1416 +8 −5 fm from Ref. [15] and the HPQCD result √ t 0 = 0.1420 (8) fm from Ref. [16], while within 1.5 standard deviations it is consistent with the recent result √ t 0 = 0.1422 (14) fm from Ref. [18].
The values of the lattice spacing corresponding to the three GF scales are collected in Table XII of Appendix D 2. At large time distances one has which allows the extraction of the kaon mass M K and the matrix element Z K = | K|q s γ 5 q |0 | 2 from the exponential fit given in the r.h.s. of Eq. (67). The kaon decay constant f K is given by and, using the pion data (16) for f π , the ratio f K /f π is evaluated at each simulated strange bare quark mass. The time intervals [t min , t max ] adopted for the fit (67) of the kaon correlation function (66) are the same as those used for the case of the pion correlator, collected in Table II.
As in the case of the pion data (see Section II), due to a small deviation from maximal twist, a correction should be applied to observables of the ensemble cA211.12.48. We use the following formula (see Appendix B) where, we remind, m P CAC is the bare untwisted PCAC mass, Z A is the renormalization constant of the axial current and µ i is the bare twisted mass of the valence quarks. In the degenerate case m s = m one gets K f = K , i.e. Eq. (19), while for m s >> m one has K f 1/cos(θ /2).
Since the LECs of the SU(2) ChPT depend on the value of the (renormalized) strange quark mass m s , we need to interpolate the ratio f K /f π at an approximately fixed value of m s . To this end we take advantage of the fact that the meson mass combination π is proportional to m s at LO in ChPT. Thus, for each gauge ensemble, adopting a simple quadratic spline, the lattice data for f K /f π are interpolated at a reference kaon mass given by The results obtained for the ratio f K /f π interpolated at the kaon reference mass (72) are shown in Fig. 7 for all the ETMC gauge ensembles. The statistical errors of the data lie in the range 0.1 ÷ 0.6%.
We now apply the correction for FVEs using the GL formula and the expansion variable ξ π defined as where f is fixed at the value given by Eq. (58). For the pion and kaon decay constants the NLO FVE corrections are respectively given by [24] ∆ so that the overall FVE correction for f K /f π is given by Finally, in terms of the variable ξ, defined in Eq. (46), the data for (f K /f π )(L → ∞) are fitted using the following ansatz where with respect to the well-known SU(2) ChPT prediction at NLO a quadratic term in ξ as well as discretization effects proportional to a 2 and a 2 M 2 π have been added. The free parameters appearing in Eq. (77) are R 0 , R 1 , R 2 , D 0 , D 1 and their values are obtained by a straightforward χ 2 -minimization procedure. We have performed several fits based on Eq. (77) and the results for the ratio (f K /f π ) isoQCD at the physical pion point (4) are collected in Table IV. The quality of the NLO fit with R 2 = D 1 = 0 is illustrated in Fig. 8. It can be seen that FVEs are properly taken care of and that discretization effects are quite small.  4) and (5), obtained using the fitting function (77). As a check of the impact of FVEs we multiply the GL correction in Eq. (76) by a factor κ F V E , which is treated as a further free parameter in the NLO fit. The factor κ F V E turns out to be consistent with unity, κ F V E = 1.19 (24), and the NLO result Putting together all the various results we obtain where we remind that () stat+f it incorporates the uncertainties induced by both the statistical errors and the fitting procedure itself. Adopting the results of the ODE procedure (see Appendix B) for the extraction of the pion and kaon masses and decay constants the analysis of the ratio f K /f π yields which compares very well with the finding (78).

VI. IMPLICATIONS FOR V us AND THE FIRST-ROW CKM UNITARITY
Inserting our isoQCD result (78) into Eq. (3) the ratio of the CKM entries V us and V ud is given by Using the value |V ud | = 0.97370 (14) which is in good agreement with the latest estimate |V us | = 0.2252 (5) from leptonic modes provided by the PDG [3].
Using the values |V ub | = 0.00382 (24) [3] and |V ud | = 0.97370 (14) [3, 22] our result (81) implies for the unitarity of the first-row of the CKM matrix the value which in turn would imply a 3σ tension with unitarity from leptonic modes. Had we used the result V ud = 0.97420 (21) from Ref. [9] the first-row CKM unitarity would be fulfilled within one standard deviation, i.e. within a precision of 0.5 permil.
Another source of information on V us is represented by the semileptonic K 3 decay.
In this case the relevant hadronic quantity is the vector form factor at zero momentum transfer f + (0). From the high-precision experimental data on K 3 decays one has V us f + (0) = 0.2165 (4) [37].
Using the ETMC determination f + (0) = 0.9709 (46) obtained with Wilson twistedmass quarks in Ref. [38], one gets the semileptonic result V us = 0.2230 (11) to be compared with the leptonic one given in Eq. (81). The above finding is combined with Eq. (80) to obtain the red ellipse in Fig. 9, which represents a 68% likelihood contour. For comparison the blue ellipse corresponds to the FLAG-4 contour for N f = 2 + 1 + 1 [19], defined by the bands corresponding to V us = 0.2231 (7) and V us /V ud = 0.2313 (5). The two determinations of V ud obtained in Refs. [9] and [22] are also shown. Finally, the dotted line represents the correlation between V us and V ud when the CKM matrix is taken to be unitary.

VII. CONCLUSIONS
We have presented a determination of the ratio of kaon and pion leptonic decay constants in isoQCD, f K /f π , adopting the gauge ensembles produced by ETMC with for N f = 2 + 1 + 1 [19] and in this work by using Eq. (80) and the semileptonic result of Ref. [38]. The determinations of V ud obtained from superallowed nuclear β transitions obtained in Refs. [9,22] are also shown as green and orange bands, labelled respectively HT and SGP R.
The dotted line indicates the correlation between V ud and V us that follows assuming the unitarity of the CKM-matrix. The ellipses represent 68% likelihood contours.
The simulations are carried out at three values of the lattice spacing ranging from ∼ 0.068 to ∼ 0.092 fm with linear lattice size up to L ∼ 5.5 fm. The scale is set using the value the pion decay constant f isoQCD π = 130.4 (2) MeV taken from Ref. [8]. Two observables, f π and (f π M 4 π ) 1/5 , have been analyzed within the framework of SU(2) ChPT without making use of renormalized quark masses. The latter quantity is found to be marginally affected by lattice artifacts and provides a precise determination of the GF scales, namely: w 0 = 0.17383 (63) fm, √ t 0 = 0.14436 (61) fm and t 0 /w 0 = 0.11969 (62) fm.
As for the decay constant ratio f K /f π we get at the physical isoQCD point, defined by Eqs. (4-6), the result where the error includes both statistical and systematic uncertainties in quadrature.

ACKNOWLEDGMENTS
We thank all the ETMC members for a very productive collaboration. We warmly thank G.C. Rossi for his continuous support and for a careful reading of the manuscript.
We are grateful to B. Joó for his kind support with our refactoring and extension of the Booster [39] and Juwels [40] at the Jülich Supercomputing Centre (JSC) and we grate-

Integrator Setups
In the generation of gauge ensembles via the Hybrid Monte Carlo algorithm, the effective lattice action can be represented by a sum over monomials corresponding to different contributions to the partition function as defined below. In the integration of the equations of motion, the forces contributed by the different monomials differ by orders of magnitude, allowing them to be integrated on different time scales accordingly, as detailed in Table VI below.

a. Monomial Types
We define below the different types of monomials that we employ in our effective lattice action to simulate QCD using N f = 2 + 1 + 1 twisted mass clover fermions.

Degenerate Determinant: [det(ρ)]
The action contribution of a degenerate doublet of clover-improved twisted mass quarks is given by in the twisted basis and in the hopping parameter normalisation, where T is the clover term. In our simulations we use the conventional value r = 1.
For convenience, we defineμ ≡ 2κµ and absorb 2κc SW into T , defining the twoflavour operator and the Hermitian operator Q sw = γ 5 D sw , where in turn D sw is the cloverimproved Wilson Dirac operator. We then have Q ± = Q sw ± iμ, such that Q † + = Q − and Q + Q − = Q 2 sw +μ 2 . The contribution to the partition function of the mass-degenerate (light) twisted mass quark doublet is thus given by det (Q + Q − ) = det Q 2 sw +μ 2 .
An even-odd Schur decomposition of the sub-matrices Q ± then gives from which we obtainQ ± defined only on the odd sites of the latticê The light quark determinant can then be reexpressed as In order to implement mass preconditioning, theQ ± can be shifted by a constant through the addition of a further twisted mass:Ŵ ± (ρ) =Q ± ± iρ, such that W +Ŵ− =Q +Q− + ρ 2 . It should be noted that this shift is applied to the evenodd-preconditioned operator, such that the factor M ± ee remains independent of ρ since its inverse is non-trivial.
In terms of pseudofermion fields, one thus obtains a contribution to the partition function which we refer to as the degenerate determinant and a corresponding contribution which we refer to as a determinant ratio.
Eq. (A8) generalises to the introduction of multiple shifts ρ 1 , ρ 2 , . . . , ρ n to contributions of the form: The pseudofermion fields φ i are defined only on the odd sites of the lattice and are generated from a random spinor field R i , sampled from a normalized Gaussian distribution at the beginning of each molecular dynamics trajectory. In the case of the determinant, we have φ i =Q + R i , while in the case of the determinant ratio The complete mass-preconditioned contribution with n shifts is thus given by where the last factor has the form of Eq. (A8) with the target twisted quark mass inQ ± . In general, the different contributions are integrated on multiple time scales because their contributions to the force differ by orders of magnitude.

Rational Approximation Partial Fraction: [rat(n , n k )]
The Dirac operator for the non-degenerate flavour doublet employed in the strange-charm sector reads with the property Equivalently, as used (without the clover term) in Ref. [41], one may write As before, we define Q h = γ 5 D h and the implementation of even-odd preconditioning translates straightforwardly from the mass-degenerate case, although the construction of M h ee has to take into account the additional (off-diagonal) flavour structure.
The operatorQ h , defined only on the odd sites, has the propertyQ † h = τ 1 fQ h τ 1 f and the non-degenerate quark doublet contributes a factor to the partition function, which we simulate via where we made use of the shorthand notationQ 2 We use a rational approximation of order N (see Refs. [42][43][44]) For this, we employ the Zolotarev solution [45] for the optimal approximation to 1/ √ y, where the coefficients a i satisfy the property The amplitude A, the coefficients a i and the maximal deviation of the rational approximation δ = max y |1 − √ yR(y)| are computed analytically at given order N and lower bound < y < 1. These are where sn(u, k) and cs(u, k) = cn(u, k)/sn(u, k) are Jacobi elliptic functions and K(k) is the complete elliptic integral. In all our simulations we use N = 10 and = λ min /λ max where λ min and λ max are respectively the lower and upper bound of the eigenvalues ofQ 2 h . In order to have all the eigenvalues λ in the range < λ < 1, we re-scaleQ 2 h with λ max . The values of λ min and λ max per ensemble are given in Table V. In the simulation, we explicitly check that the eigenvalues ofQ 2 h remain within these bounds.
The factors in the approximation can be grouped where We perform partial fraction expansions of the terms such that the necessary matrix inverse can be calculated efficiently using a multishift solver. The coefficients q i are given by We can further define µ i = √ a 2i and ν i = √ a 2i−1 and express q i as At the beginning of each trajectory, pseudofermion fields are generated as follows: again a random spinor field R is sampled from a Gaussian distribution. Now, we need to compute pseudofermion fields φ from and, therefore, we need operators C † and C with the property C is given by (inspired by twisted mass) which can again be written as a partial fraction The rational approximation R can be applied to a vector using a multi-mass solver and the partial fraction representation. The same works for C: after solving N equations simultaneously for (Q 2 h +ν 2 i ) −1 , i = 1, ..., N , we have to multiply every term with (Q h − iν i ). The hermitian conjugate of C is given by For the acceptance step only the application of R is needed.

Rational Approximation Correction Factor: [ratcor(n)]
The rational approximation R only has a finite precision. This finite precision can be accounted for during the acceptance step in the HMC by estimating [44] 1 − |Q h |R, if the rational approximation is precise enough. This can be achieved by including a monomial det(|Q h |R) in the simulation, for which one needs an Following Ref. [44], B can be written as The series converges rapidly and can, thus, be truncated after a few terms, m + 1. The convergence can actually be controlled during the simulation and the truncation does not need to be fixed. We choose to sum the series until the contribution of the given term to the acceptance Hamiltonian is below the residual precision squared, r 2 , that we employ for the solution of the linear systems involved in the approximation of R in the acceptance step, such that |c m Z m φ| 2 < r 2 , which we typically choose to be at the limit of double precision arithmetic.
For this monomial the pseudofermion field is computed from where R is again a Gaussian random vector, see above.

b. Simulation parameters
In Table VI  for λ is larger than 1/6. Namely assuming unity of the second order commutators and neglecting any correlations leads to the optimal value of λ ≈ 0.1931833275. In the usage of the 2MN integrator with multiple time scales, experience suggests that further deviations from this optimal value improve the acceptance rate, such that we often use schemes with increasing values of λ from the innermost to the outermost time scale, as shown in Table VI.

Software details
The simulations presented in this study have been generated using the tmLQCD [47][48][49] software suite, which provides all the necessary components to perform N f = 2+1+1 simulations of twisted mass clover fermions, including implementations of the polynomial and rational HMC algorithms for the non-degenerate determinant. To enable multigrid solvers to be used in simulations [50], tmLQCD provides an interface to DDαAMG [51], a multigrid solver library optimized for twisted mass (clover) fermions [52]. The force calculation of some monomials in the light quark sector is accelerated by a 3-level multi-grid approach. Moreover, we extended the DDαAMG method for the mass non-degenerate twisted mass operator. The multi-grid solver used in the rational approximation [53] is particularly helpful for the lowest terms of the rational approximation, as well as for 3 2MN 14 0.193 detrat(0, 0.0012), rat (7,9)  does not benefit at all. To improve efficiency, tmLQCD also provides an interface to the QPhiX [54] lattice QCD library, which we have refactored and extended [55] to support twisted mass clover fermions, including the non-degenerate doublet. For solves related to the degenerate determinant and determinant ratios, this allows us to efficiently and flexibly combine mixed-precision CG and SIMD vector lengths of 8 or 16 as required by AVX512. On KNL, single-precision QPhiX kernels are up to a factor of 5 more efficient than their tmLQCD-native equivalents. Also in the MMS-CG solves in the non-degenerate sector, the double-precision kernels in QPhiX are up to a factor of 2 more efficient than the tmLQCD-native equivalents on KNL. Combined, these efficiency improvements lead to overall speedup factors of 2-3 in the HMC on this architecture with smaller overall gains on Skylake and EPYC.
Appendix B: Extraction of aM π and af π using the ODE procedure The spectral decomposition of the pion correlation function (14) can be investigated adopting the ODE procedure of Ref. [23]. This method is able to extract exponential signals from the temporal dependence of a lattice correlator without any a priori knowledge of the multiplicity of each signal and it does not require any starting input for the masses and the amplitudes of the signals.
The ODE approach is sensitive to the noise of the correlator, so that pure oscillating signals (conjugate pairs of imaginary masses) may typically appear in the ODE spectral decomposition. Therefore, we adopt an improved version of the ODE procedure, in which a subsequent χ 2 -minimization procedure is applied to the non-noisy part of the ODE spectral decomposition [23]. In this way the accuracy of the physical (i.e. non-noisy) part of the ODE spectral decomposition is improved.
The time intervals [t min , t max ] adopted for the analysis and the extracted values of the pion mass and decay constant in lattice units are collected in Table VII.
Within the ODE procedure we searched for 8 exponential signals in the time intervals of Table VII and in all cases at least two physical (non-noisy) exponential signals were found. Then, a χ 2 -minimization procedure was applied using the physical ODE solution as the starting point. The minimized values of the χ 2 variable turned out to be always less than 1.
The extracted values as well as their statistical errors of the ground-state mass and decay constant, collected in Table VII, are nicely consistent with the corresponding ones obtained by the direct single exponential fit (15) shown in Table II  Using the above pion data for f π the NLO analysis of Section IV A, including the discretization term proportional to a 2 M 2 π , yields for the GF scale w 0 the value in agreement with the result (52). Analogously, the use of the data for X π and of the NLO fit (56) with A 2 = F F V E = 0 (see Section IV B) leads to in agreement with the corresponding result shown in the second row of Table III.

Appendix C: Maximal twist corrections for masses and decay constants
We follow the general approach of Ref. [31] to the mixed action formulation of twisted mass lattice QCD, which ensures an unitary continuum limit (provided sea and valence renormalized quark masses are matched). Here however we allow for small deviations (due e.g. to numerical errors) from the maximal twist case, i.e. for m 0 = m cr .
1. N f =2+1+1 isosymmetric QCD with twisted clover Wilson quarks The lattice action can be conveniently written in terms of gauge, sea quark and valence quark plus valence ghost fields. If the sea quarks are arranged in two-flavour fields χ and χ h and the valence quarks are described by one-flavour fields χ f , with f = u, d, ..., we have with the valence sector given by where ellipses stand for possible replica (χ f ) of the valence quarks with µ f = −µ f and the ghost terms exactly cancel the valence fermion contributions to the effective action.
Here we find it convenient to express all fermion fields in the canonical quark basis for untwisted Wilson fermions and denote by D W clov the well-known clover improved (gauge covariant) Dirac matrix: We start by discussing the light valence quark sector, the extension to heavier flavours is straightforward. Following Refs. [28,32], we define (as customary) the twist angle ω in terms of the bare mass parameters of the light valence quark (u, with Z A the renormalization constant ofX γ µ γ 5 (τ 1,2,3 /2)X , which, being independent of quark mass parameters, is defined in the chiral limit µ f → 0, m 0 → m cr . Maximal twist corresponds to |ω | = π/2, i.e. to angle θ ≡ π/2 − ω equal to zero or π. We thus have cos θ = sin ω = 1 .

(C4)
Here µ is the bare twisted mass parameters for the (u, d) doublet and m P CAC denotes the untwisted bare quark mass of the (u, d) doublet as obtained from the non-singlet WTI's -hence a function of m 0 plus the other bare parameters. We recall that m P CAC ∝ m 0 − m cr . The renormalized twisted and untwisted quark mass parameters that appear in the chiral WTI's read (up to discretization effects) where Z P and Z S 0 are the renormalization constants of the pseudoscalar non-singlet and the scalar singlet densities (in the canonical basis for untwisted Wilson quarks).
Defining the twist angle and hence formulating the maximal twist condition in terms of m P CAC , as measured on the ensembles with 2+1+1 dynamical flavours, effectively takes care of (compensates for) all the UV cutoff effects related to the breaking of chiral symmetry, including those coming from the 2+1+1 sea quark flavours.

Pion mass and decay constant
We argue here that in the case of small enough numerical deviations from maximal twist the lattice charged pion quantities with P 1 =X γ 5 (τ 1 /2)X and X = (χ u , χ d ) T , approach M π and f π as a → 0 with lattice artifacts having numerically small, and (we shall see) within errors immaterial, differences as compared to the O(a 2 ) cutoff effects occurring at maximal twist. Of course these values of M π and f π correspond to the light quark renormalized mass The numerical information on M π and f π comes from the simple correlator C 11 P P (x 0 ) = a 3 x P 1 (x)P 1 (0) (and C 22 P P (x 0 )). The large-x 0 behaviour of C 11 P P (x 0 ) determines M π and an exact lattice WTI relates the operator P 1 to the four-divergence of a conserved lattice (backward one-point split) current, which we denote byV 2 χ,µ , v.i.z.
implying that the pion-to-vacuum matrix element ofV 2 χ,µ gives information on f π (barring the case of cos θ = 0) . In Eq. (C7) P 1 R = Z P P 1 and the equalities hold at operator level for finite lattice spacing (a > 0). Hence the l.h.s. of Eq. (C7) is a renormalized operator and information on the approach of its matrix elements to the continuum limit can be obtained by studying the behaviour as a → 0 of the corresponding matrix elements of 2µ R P 1 R . Taking the matrix element of Eq. (C7) between the vacuum and a one-π 1 state of zero three-momentum and noting that in the continuum limit where ψ = (u, d) T obeys the (continuum) e.o.m. (γ · D + M R )ψ = 0, for a > 0 one obtains 2[µ π 1 (0)|P 1 |0 ]| L = [cos θ M 2 π f π + sin θ π 1 (0)|∂ 0 (ψγ 0 This relation implies that as a → 0 the ratio at generic θ = ±π/2. Hence at a > 0 the ratio (C10) represents a bona fide lattice estimator of f π , while its discretization errors depend on the lattice artifacts in M 2 π , θ (or equivalently m R , µ R ) and the renormalized quantity 2µ G π 1 = 2µ π 1 (0)|P 1 |0 .
We are interested here in situations where Z A m P CAC /µ < 1 but not negligibly small, say slightly above 0.1. This situation indeed occurs in our gauge configuration ensemble cA211.12.48, at a ∼ 0.095 fm and aµ = 0.0012, where we find Z A m P CAC /µ ∼ −0.15.
In this case an analysisà la Symanzik of M 2 π , m R and 2µ G π 1 shows (see below) that the change in the lattice artifacts of our lattice estimator of f π , with respect to those purely O(a 2n ) (with n integer) that occur at maximal twist, is smaller than 0.001f π , therefore numerically immaterial within statistical errors that are typically of order 0.005f π .
a. The change in the lattice artifacts for M 2 π and f π If |am P CAC | is non-zero, though definitely smaller than |aµ |, the same holds for |a(m 0 − m cr )| and one expects that the appropriate lattice estimators of M π , f π and any other physical quantity will be altered already at O(a) as compared to their counterparts at maximal twist. Correcting analytically the lattice estimators for the deviation from maximal twist at order a 0 is hence not enough and one must also check that the out- where X val = (χ u , χ d ) T describes the valence light quarks in the same basis as in Eq. (C2) while ellipses stands for d ≤ 4 terms involving heavier valence quarks and ghost terms.
Upon taking the continuum limit in isosymmetric QCD with 2+1+1 dynamical flavours, we must have coinciding sea and valence renormalized masses for each flavour and keep constant as a → 0 the renormalized parameters g 2 R , µ R , µ R h , R h . The LEL terms L n are suitable linear combinations of the d = n > 4 operator terms allowed by the symmetries of the lattice theory (C1). In particular it turns out that where O 5 (O 4 ) is a linear combination of the m-independent terms allowed in L 5 (L 4 ).
Among the latter, since we employ in our correlators flavour diagonal OS valence quark fields χ val f with twisted mass µ f > 0, or χ val f fields with twisted mass µ f = −µ f , for the purposes of this Section only the terms bilinear in the X val andX val , or X val and X val , fields are relevant, which we may call mO val 5, and m 2 O val 4, . We note that exact or spurionic lattice symmetries rule out 4 the L 5 terms of the form µ f i F ·F , µ fXf τ 0,1,2,3 iγ 5 γ ·DX f , µ fXf τ 0,1,2,3 γ ·DX f , as well as the analogous L 6 terms of the form mµ f i F ·F , mµ fXf τ 0,1,2,3 iγ 5 γ ·DX f , mµ fXf τ 0,1,2,3 γ ·DX f .

The undesired O(a) modification in M 2
π is estimated to be smaller than 0.001 M 2 π , hence immaterial within our statistical errors of a few permil. The O(a 2 ) change in M 2 π due to the non-zero value of am is also of order 0.001 M 2 π or smaller, because of the form (C13) of L 6 , with |am| 0.0002 and a 2 Λ 2 QCD sin θ 0.001.
d. The discretization effects on f π |L Based on Eq. (C10) the lattice artifacts on f π can be estimated in terms of the cutoff effects in (µ / cos θ ) π 1 (0)|P 1 |0 | L and in M 2 π | L . We discussed above the lattice artifacts of M 2 π | L . Concerning (µ / cos θ ) π 1 (0)|P 1 |0 | L , we can see it as the product of the renormalized quantities Z −1 P (µ / cos θ ) = M R and Z P π 1 (0)|P 1 |0 ]| L = G R π 1 , and then discuss separately the lattice artifacts in each of these two factors 5 .
As for Z −1 P (µ / cos θ ) = (m R ) 2 + (µ R ) 2 , the form of the m-dependent terms in L 5 , i.e. those with coefficients b m and b µ (which are all in modulus 1), implies that the O(a) corrections to m R and µ R , and hence to M R are only of relative size < |am| 0.0002.
Even smaller are the corrections to M R of order a 2 m.
Let us then consider the out-of-maximal-twist induced cutoff artifacts in the matrix element Z P π 1 (0)|P 1 |0 ]| L = G R π 1 . They clearly arise from from the lattice two-point correlator C 11 P (x 0 ). At order a, since for the operator P 1 =ψ γ 5 (τ 1 /2)ψ it is known that c P = 0,b P = 0 and b P = 1 + O(g 2 0 ), implying |b P am| ∼ 0.0002, the numerically dominant lattice artifacts stem from the term a d 4 y P 1 (x)P 1 (0)L 5 (y) | cont in the e. Analysis of M 2 π and f π in terms of the renormalized light quark mass The discussion above is valid also in case the observables M 2 π and f π are analysed as functions of the renormalized light quark mass As already noted in Section C 2 d, for the ensemble cA211.12.48 the observed small deviation from maximal twist leads to an undesired O(am) relative change in M R that is only of order |am| 0.0002 and thus fully negligible in comparison to other errors.
A rather obvious but practically important and general caveat follows in the case one insists in analysing the lattice data, e.g. for M π or f π , obtained on gauge ensembles with non-zero am-values in terms of µ R rather than M R . Since for am = 0 a generic observable Q obs actually refers to M R > µ R , one should, besides possibly applying an analytic correction to the datum for Q obs (as requested e.g. for f π , but not for M π ), also shift the value of the observable itself according to where in practice, since typically |am| < 0.001, only the terms linear in µ R − M R are numerically important.

Mass and decay constant of the kaon and heavier PS-mesons
Generalizing the arguments of Section C 2 for the mass and decay constant of the pion to the case of the kaon or heavier pseudoscalar (PS) mesons is rather straightforward.
Indeed, within the N f = 2 + 1 + 1 lattice QCD framework of Section C 1, as far as the control of the effects of a small but non-zero value of am = aZ A m P CAC is concerned, there is not much difference between a PS meson made out of two light valence quarks (pion, with |µ u | = |µ d |) and a PS meson made out of a light (mass µ ) and a heavier (mass µ x ) valence quark. This is a consequence of the fact that the extraction of the PS meson mass and decay constant relies on general positivity properties of two-point correlation functions (close enough to the continuum limit) and on exact chiral WTI's, which hold valid irrespectively of the value of valence quark masses. One might thus deal in a fully analogous way with PS mesons made out of two non-light valence quarks.
For definiteness, however, we shall here focus on the case of the kaon, i.e. µ x = µ s µ and µ d = µ > 0.
It turns out that in the case of small enough numerical deviations from maximal twist, say 0.1 aµ < |am| aµ aµ s , the lattice charged kaon quantities with P s,d =X s,d γ 5 (τ 1 /2)X s,d and X s,d = (χ s , χ d ) T , approach M K and f K as a → 0 with lattice artifacts having numerically small, and within errors immaterial, differences as compared to the O(a 2 ) cutoff effects occurring at maximal twist. These values of M K and f K correspond to the light (d) quark renormalized mass M R = (m R ) 2 + (µ R ) 2 and to the strange (s) quark renormalized mass M R s = (m R ) 2 + (µ R s ) 2 . In full analogy to the definitions adopted for light valence quarks (see Section C 1) we define µ R s = µ s /Z P , where Z P is (can be conveniently taken as) the same massindependent renormalization constant of the PS non-singlet quark bilinear operator as above, and cos θ s = sin ω s = 1 1 + (Z 2 A m 2 P CAC )/µ 2 s , sin θ s = cos ω s = 1 . (C20) a. On the lattice estimators (C19) of M K and f K at O(a 0 ) The numerical information on M K and f K comes in fact from the simple correlator C s,d KK (x 0 ) = a 3 x P s,d (x)P d,s (0) . The large-x 0 behaviour of C s,d KK (x 0 ) determines M K as the kaon mass that, owing to renormalizability (and unitarity in the continuum limit) of the lattice theory of Section C 1, corresponds to the light and strange quark masses s . An exact lattice WTI relates the operator P s,d to the four-divergence of a conserved lattice (backward one-point split) current, which we denote byV s,d χ,µ , viz.
implying that the kaon-to-vacuum matrix element ofV s,d χ,µ gives information on f K , barring the case of cos((θ + θ s )/2) = 0 . In Eq. (C21) P s,d R = Z P P s,d and the equalities hold at operator level for a > 0. Hence the l.h.s. of Eq. (C21) is a renormalized operator and information on the approach of its matrix elements to the continuum limit can be obtained by studying the behaviour as a → 0 of the corresponding matrix elements of (µ R + µ R s )P S,d R .
Taking the matrix element of Eq. (C21) between the vacuum and a one-K state of zero three-momentum and noting that in the continuum limit for the light valence quark and gluonic terms in L 5 , one finds that the discussion for the kaon case closely follows the one for the pion (see Section C 2 a), provided one replaces the valence quark field pair X val = (χ u , χ d ) T , used for the latter, with the valence quark field pair X val s,d = (χ s , χ d ) T relevant for the kaon, as well as θ with (θ s + θ )/2. This implies that for M 2 K and f K the numerically dominant changes in the lattice artifacts as compared to the maximal twist case, which for M 2 π and f π were proportional to sin θ , turn out to be proportional to thereby getting reduced by a factor of about two with respect to the pion case.
In conclusion, if |am| ∼ 0.0002 (as it happens for our ensemble cA211.12.48), we estimate a lattice artifact modification, with respect to the case of maximal twist, that does not exceed 0.0005 M 2 K for M 2 K and 0.0002 f K for f K and is hence safely negligible as compared to our current statistical errors.
From the arguments above it should also be clear that the same quantitative estimates hold also for the lattice artifact changes induced by am = 0 in the mass and the decay constant of heavy-light PS mesons with a charm or even heavier non-light valence quark having a mass µ x µ s . In fact the heavier the valence quark, the smaller |θ x | |m/µ x | and the more numerically irrelevant the effect of a small non-zero value of am.
absolute scales √ t 0 and t 0 /w 0 using the SU(2) ChPT analysis of the data for X π , carried out in Section IV B in the case of the GF scale w 0 .

Determination of the relative GF scales
In this section we provide the details of the determinations of the gradient flow (GF) scales at the physical point. Our analysis is based on the values of the gradient flow scales in Table VIII as obtained on the ensembles in Table I. We calculate the scales following the definitions in [13,14,58] using the standard Wilson action for the gradient flow evolution and the symmetrized discretization of the action density as described in [14].
Apart from the usual scales s 0 /a ≡ √ t 0 /a and w 0 /a, we also consider the derived scale (t 0 /w 0 )/a as well as the dimensionless ratio (s 0 /w 0 ). The former is interesting, because it exhibits very mild quark-mass dependence and reduced autocorrelations, while the dimensionless ratio (s 0 /w 0 ) can be used for assessing lattice artefacts and crosschecking the consistency of the various analysis procedures and determinations of the scales in physical units.
The errors are calculated by taking into account the autocorrelations of the various quantities, which are also listed in Table VIII. It is well known that the autocorrelations can be sizable for the GF scales. In particular on coarse lattices it is known that they can be larger than the ones for the topological charge [59]. More importantly, they are also expected to grow towards the chiral limit. The data shown in Table VIII phys. , i.e., X(M π )/a = (X/a) phys. + c · ∆ 2 . Note that for the ensembles of the scale (t 0 /w 0 )/a is very small, i.e., the corrections are less than 0.5% at our largest pion mass ensemble cA211.53.24, to be compared to 2.5% for √ t 0 /a and 4.5% for w 0 /a.
We note that the quark-mass dependence exhibits clearly visible lattice artefacts, but the dependence seems to become weaker towards the continuum limit. One peculiar feature is the fact that for the scale (t 0 /w 0 )/a the slope of the quark-mass dependence changes sign when going from the coarser to the finer lattice spacing. This certainly warrants further investigation, once more data is available, however, one should keep in mind that for this quantity the slope is consistent with zero within less than 3σ.
In order to examine the lattice artefacts of the GF scales further, we now turn to the dimensionless ratio s 0 /w 0 . In Figure 11 we  which includes a light quark-mass dependence proportional to ∆ 2 described by B 0 and lattice artefacts proportional to a 2 /w 2 0 described by A 1 and B 1 . The latter coefficient describes the lattice artefacts on the quark-mass dependence. The global fit suggests that B 0 , describing the quark-mass dependence in the continuum, is well consistent with zero, i.e., B 0 = 0.001 (7). That is, in the continuum the ratio √ t 0 /w 0 appears to have no dependence on the pion mass at all and the observed pion-mass dependence at finite lattice spacings is apparently just a lattice artefact. However, given the fact that we do not have data for the ratio at the lattice spacing cC211 away from the physical point, and hence no information on the quark-mass dependence at the finest lattice spacing, it is not clear how solid this conclusion is. Nevertheless, we may attempt to fit our data with B 0 = 0 fixed, and in figure 12 we show the results for this global fit. The We note that the ratio is determined with a precision in the sub-permille region, i.e., better than 0.8 permille. As such, it provides an interesting consistency crosscheck on any other, independent determination of the scales, e.g., through hadronic quantities. and D 1 obtained in the case of the three GF scales w 0 , √ t 0 and t 0 /w 0 using the data on X π and adopting fitting functions similar to Eq. (56) with A 2 = F F V E = 0.
Finally, the values of the lattice spacing a corresponding to the three GF scales (D3-D5) and to the relative scales given in Table X are shown in Table XII. The three determinations of a differ by O(a 2 ) effects, as shown in Fig. 14. In particular, we get: a( √ t 0 )/a(w 0 ) 1 − 0.09 (2) a 2 (w 0 )/w 2 0 and a(t 0 /w 0 )/a(w 0 ) 1 − 0.18 (2) a 2 (w 0 )/w 2 0 .
[1] N. Cabibbo, Phys. Rev. Lett. 10, 531 (1963 FIG. 14. Ratio of the lattice spacing a obtained from the GF scales √ t 0 (red circles) and t 0 /w 0 (blue squares) with the one determined from the GF scale w 0 (see Table XII). The solid and dashed lines are linear fits.