Non-perturbative quark mass renormalisation and running in $N_f=3$ QCD

We determine from first principles the quark mass anomalous dimension in Nf=3 QCD between the electroweak and hadronic scales. This allows for a fully non-perturbative connection of the perturbative and non-perturbative regimes of the Standard Model in the hadronic sector. The computation is carried out to high accuracy, employing massless O(a)-improved Wilson quarks and finite-size scaling techniques. We also provide the matching factors required in the renormalisation of light quark masses from lattice computations with O(a)-improved Wilson fermions and a tree-level Symanzik improved gauge action. The total uncertainty due to renormalisation and running in the determination of light quark masses in the SM is thus reduced to about 1%.


Introduction
In the paradigm provided by the Standard Model (SM) of Particle Physics, quark masses are fundamental constants of Nature. More specifically, Quantum Chromodynamics (QCD), the part of the SM that describes the fundamental strong interaction, is uniquely defined by the values of the quark masses and the strong coupling constant. Apart from this intrinsic interest, precise knowledge of the values of quark masses is crucial for the advancement of frontier research in particle physics -one good illustration being the fact that the values of the bottom and charm quark masses are major sources of uncertainty in several important Higgs branching fractions, e.g., Γ(H → bb) and Γ(H → cc) [1][2][3][4][5].
Quark masses are couplings in the QCD Lagrangian, and have to be treated within a consistent definition of the renormalised theory. A meaningful determination can only be achieved by computing physical observables as a function of quark masses, and matching the result to the experimental values. A non-perturbative treatment of QCD is mandatory to avoid the presence of unquantified systematic uncertainties in such a computation: the asymptotic nature of the perturbative series, and the strongly coupled nature of the interaction at typical hadronic energy scales, implies the presence of an irreducible uncertainty in any determination that does not treat long-distance strong interaction effects from first principles. Lattice QCD (LQCD) is therefore the best-suited framework for a high-precision determination of quark masses. Indeed, following the onset of the precision era in LQCD, the uncertainties on the values of both light and heavy quark masses have dramatically decreased in recent years [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22].
The natural observables employed in a LQCD computation of quark masses are hadronic quantities, considered at energy scales around or below 1 GeV. This requires, in particular, to work out the renormalisation non-perturbatively. Then, in order to make contact with the electroweak scale, where the masses are used to compute the QCD contribution to highenergy observables, the masses have to be run with the Renormalisation Group (RG) across more than two orders of magnitude in energy. While high-order perturbative estimates of the anomalous dimension of quark masses in various renormalisation schemes exist [23][24][25], a nonperturbative determination is mandatory to match the current percent-level precision of the relevant hadronic observables.
In this work we present a high-precision determination of the anomalous dimension of quark masses in QCD with three light quark flavours, as well as of the renormalisation constants required to match bare quark masses. 1 This is a companion project of the recent highprecision determination of the β function and the Λ QCD parameter in N f = 3 QCD by the ALPHA Collaboration [29][30][31][32]. We will employ the Schrödinger Functional [33,34] as an intermediate renormalisation scheme that allows to make contact between the hadronic scheme used in the computation of bare quark masses and the perturbative schemes used at high energies, and employ well-established finite-size recursion techniques [35][36][37][38][39][40][41][42][43][44][45][46] to compute the RG running non-perturbatively. Our main result is a high-precision determination of the mass anomalous dimension between the electroweak scale and hadronic scales at around 200 MeV, where contact with hadronic observables obtained from simulations by the CLS effort [47] can be achieved.
The paper is structured as follows. In Section 2 we describe our strategy, which (similar to the determination of Λ QCD ) involves using two different definitions of the renormalised coupling at energies above and below an energy scale around 2 GeV. Sections 3 and 4 deal with the determination of the anomalous dimension above and below that scale, respectively. Section 5 discusses the determination of the renormalisation constants needed to match to a hadronic scheme at low energies. Conclusions and outlook come in Section 6. Several technical aspects of the work are discussed in appendices.

Quark running and RGI masses
QCD is a theory with N f + 1 parameters: the coupling constant g and the N f quark masses {m i , i = 1, . . . , N f }. When the theory is defined using some regularisation, g and m i are taken to be the bare parameters appearing in the Lagrangian. Removing the regularisation requires to define renormalised parameters g , {m i , i = 1, . . . , N f } at some energy scale µ. In the following we will assume the use of renormalisation conditions that are independent of the values of the quark masses, which leads to mass-independent renormalisation schemes. The renormalised parameters are then functions of the renormalisation scale µ alone [48,49], and their scale evolution is given by renormalisation group equations of the form The renormalisation group functions β and τ admit perturbative expansions of the form −g 2 (d 0 + d 1 g 2 + d 2 g 4 + . . .) , (2.4) with universal coefficients given by [50][51][52][53][54][55][56] b 0 = 1 (4π) 2 (2.6) The higher-order coefficients b n≥2 , d n≥1 are instead renormalisation scheme-dependent. It is trivial to integrate Eqs. (2.1,2.2) formally, to obtain the renormalisation group invariants (RGI) Note that the integrands are finite at g = 0, making the integrals well defined (and zero at universal order in perturbation theory). Note also that Λ QCD and M i are non-perturbatively defined via the previous expressions. It is easy to check, furthermore, that they are N f -dependent but µ-independent. They can be interpreted as the integration constants of the renormalisation group equations. Also the ratios m i (µ)/m j (µ), i = j are scale-independent. Furthermore, the ratios M i /m i (µ) are independent of the quark flavour i, due to the mass-independence of the quark mass anomalous dimension τ . Finally, the values of M i can be easily checked to be independent of the renormalisation scheme. The value of Λ QCD is instead scheme-dependent, but the ratio Λ QCD /Λ QCD between its values in two different schemes can be calculated exactly using one-loop perturbation theory.

Step scaling functions
In our computation, we will access the renormalisation group functions β and τ through the quantities σ and σ P , defined as and to which we will refer as coupling and mass step scaling functions (SSFs), respectively. They correspond to the renormalisation group evolution operators for the coupling and quark mass between scales that differ by a factor of two, viz. (2.10) The main advantage of these quantities is that they can be computed accurately on the lattice, with a well-controlled continuum limit for a very wide range of energy scales. This is so thanks to the use of finite size scaling techniques, first introduced for quark masses in [35]. The RG functions can be non-perturbatively computed between the hadronic regime and the electroweak scale, establishing safe contact with the asymptotic perturbative regime In order to compute the SSF σ P , we define renormalised quark masses through the partially conserved axial current (PCAC) relation, where the renormalised, non-singlet (i = j) axial current and pseudoscalar density operators are given by In these expressions Z P is the renormalisation constant of the pseudoscalar density in the regularised theory, and Z A is the finite axial current normalisation, required when the QCD regularisation breaks chiral symmetry, as with lattice Wilson fermions. Eq. (2.11) implies that, up to the finite current normalisation, current quark masses renormalise with Z −1 P . Therefore, the SSF σ P of Eq. (2.9b) can be obtained by computing Z P at fixed bare gauge coupling g 2 0 - i.e., fixed lattice spacing -for scales µ and µ/2. This is repeated for several different values of the lattice spacing a, and the continuum limit of their ratio is taken, viz.
where g 2 0 is the bare coupling, univocally related to a in mass-independent schemes. The condition that the value of the renormalised coupling u = g 2 (µ) is kept fixed in the extrapolation ensures that the latter is taken along a line of constant physics. 2 In this work we will determine non-perturbatively τ (ḡ) from Eq. (2.9b). Note that in order to do so, we need the RG function β(ḡ). As we will note later, the non-perturbative determination of the β function has already been done in our schemes of choice [29,31]. In practice, given the β function (and hence σ(u)), and the lattice results for Σ P (g 2 0 , aµ), one determines the anomalous dimension τ (g ) by extrapolating Σ P (u, aµ) to the continuum (Eq. (2.14)) and then using the relation Eq. (2.9b) to constrain the functional form of τ (g ).

Renormalisation schemes
In order to control the connection between hadronic observables and RGI quantities, we will use intermediate finite-volume renormalisation schemes that allow to define fully non-perturbative renormalised gauge coupling and quark masses. For that purpose, Z P will be defined by a renormalisation condition imposed using the Schrödinger Functional (SF) [33,34]. In the following, we will adopt the conventions and notations for the lattice SF setup introduced in [57].
SF schemes are based on the formulation of QCD in a finite space-time volume of size T 3 × L, with inhomogeneous Dirichlet boundary conditions at Euclidean times x 0 = 0 and x 0 = T . The boundary condition for gauge fields has the form wherek is a unit vector in the direction k, P exp is a path-ordered exponential, and C k is some smooth gauge field. A similar expression applies at x 0 = T in terms of another field C k . Fermion fields obey the boundary conditions with P ± = 1 2 (1 ± γ 0 ). Gauge fields are periodic in spatial directions, whereas fermion fields are periodic up to a global phase The SF is the generating functional where the integral is performed over all fields with the specified boundary values. Expectation values of any product O of fields are then given by where O can involve, in particular, the "boundary fields" . (2.21) The Dirichlet boundary conditions provide an infrared cutoff to the possible wavelengths of quark and gluon fields, which allows to simulate at vanishing quark mass. The presence of non-trivial boundary conditions requires, in general, additional counterterms to renormalise the theory [33,58,59]. In the case of the SF, it has been shown in [60] that no additional counterterms are needed with respect to the periodic case, except for one boundary term that amounts to rescaling the boundary values of quark fields by a logarithmically divergent factor, which is furthermore absent ifρ = ρ =ρ = ρ = 0. It then follows that the SF is finite after the usual QCD renormalisation. The SF naturally allows for the introduction of finite-volume renormalisation schemes, where the renormalisation scale is identified with the inverse box length, The renormalisation of the pseudoscalar density, and hence of quark masses, is treated in the same way as in [35]. We introduce the SF correlation functions where P ij is the unrenormalised pseudoscalar density, and O, O are operators with pseudoscalar quantum numbers made up of boundary quark fields The pseudoscalar renormalisation constant is then defined as with all correlation functions computed at zero quark masses. The renormalisation condition is fully specified by fixing the boundary conditions and the box geometry as follows: Furthermore, for computational convenience (cf. below), all correlation functions will be computed in a fixed topological sector of the theory, chosen to be the one with total topological charge Q = 0. This is just part of the scheme definition, and does not change the ultraviolet structure of the observables. In order to completely fix the renormalisation scheme for quark masses, we still need to provide a definition of the renormalised coupling. This allows to relate the scale µ = 1/L to the bare coupling, and hence to the lattice spacing, in an unambiguous way, so that the continuum limit of Σ P is precisely defined. Following [29,31], we will introduce two different definitions, to be used in qualitatively different regimes. For renormalisation scales larger than some value µ 0 /2, we will employ the non-perturbative SF coupling first introduced in [33]. Below that scale, we will use the gradient flow (GF) coupling defined in [61]. As discussed in [62], this allows to optimally exploit the variance properties of the couplings, so that a very precise computation of the β function, and ultimately of the Λ QCD parameter, is achieved. 3 In our context, the main consequence of this setup is that our quark masses are implicitly defined in two different schemes above and below µ 0 /2; we will refer to them as SF and GF, respectively. Note that the schemes differ only by the choice of renormalized coupling g 2 ; the definition of Z P is always given by Eq. (2.26).
The running of the SF coupling is thus known accurately down to µ 0 /2, and the matching of the two schemes is completely specified by measuring the value of the GF coupling at µ 0 /2, for which one has [31] g 2 GF (µ 0 /2) = 2.6723 (64) .
When expressed in physical units through the ratio µ 0 /Λ QCD , one finds µ 0 ≈ 4 GeV [32] -i.e., the scheme switching takes place at a scale around 2 GeV. It is important to stress that the scheme definition affects different quantities in distinct ways. Obviously, the β function, being a function of the coupling, will be different in the two schemes. The same is true of the mass anomalous dimension τ (g). The renormalised masses m i (µ), on the other hand, are smooth functions across µ 0 /2 by construction, since -unlike the RG functions β and τ , which are functions of g -they are functions of the scale µ, and are fixed by the same definition of Z P at all scales. This observation also provides the matching relation for the anomalous dimensions: for any fixed µ we have Another important motivation for this choice of strategy is the control over the perturbative expansion of the β function and mass anomalous dimension, which becomes relevant at very high energies. In the SF scheme the first non-universal perturbative coefficient of the β function, b 2 , is known [63], Moreover, the next-to-leading order (NLO) mass anomalous dimension in the SF scheme was computed in [64], (2.33) A similar computation in the GF scheme is currently not available, due to the absence of a full two-loop computation of the finite-volume GF coupling in QCD. Let us end this section summarising the results for the β function in our choice of schemes. As discussed above, these results will be essential to the determination of the anomalous dimension τ (g ) in the following sections. On the high-energy side we have [29] β SF (g ) = −g 3 Note that the three leading coefficients are given by the perturbative predictions, which implies that safe contact with the asymptotic perturbative behavior has been made (this is the reason why Eq. (2.34) is accepted as a reliable approximation all the way up to g = 0). On the other hand, on the low energy side, we have [31] β GF (g ) = − g 3 2 n=0 p n g 2n , The reader should note that these values are not exactly the same as those quoted as final results in [31]. There are two reasons for this. First we have added some statistics in some ensembles, where the uncertainty in Σ P was too large. Second, a consistent treatment of the correlations and autocorrelations between Z P and g 2 GF requires knowledge of the joint autocorrelation function in a consistent way. This requirement results in a different covariance matrix between the fit parameters. In any case it is very important to point out that both results, Eqs. (2.36b, 2.36c), and those quoted in [31] are perfectly compatible. The reader can easily check that the differences in the β function are completely negligible within the quoted uncertainties.
All data is analysed using the Γ-method [65] to account for autocorrelations (some integrated autocorrelation times are given in the tables of appendix B). For a full propagation of uncertainties into derived quantities, we subsequently apply standard resampling techniques (boostrap).

Determination of RGI quark masses
In order to determine RGI quark masses, we will factor Eq. (2.8) as The three ratios appearing in this expression are flavour-independent running factors: 4 • m i (µ 0 /2)/m i (µ had ) is the running between some low-energy scale µ had and the schemeswitching scale µ 0 /2. It will be computed non-perturbatively in the GF scheme.
• m i (µ pt )/m i (µ 0 /2) is the running between the scheme-switching scale µ 0 /2 and some high-energy scale µ pt . It will be computed non-perturbatively in the SF scheme.
• M i /m i (µ pt ) can be computed perturbatively using NLO perturbation theory in the SF scheme. This perturbative matching would be safe, entailing a small systematic uncertainty due to perturbative truncation, provided µ pt is large enough -say, µ pt of order M W . As discussed in Sec. 3, we will actually use our non-perturbative results for the mass anomalous dimension at high energies to constrain the truncation systematics, and obtain a very precise matching to perturbation theory.
Finally, the renormalised quark mass m i (µ had ) at the low-energy scale is to be computed independently from the running factors, by determining bare PCAC quark masses m i from largevolume lattice simulations at a number of values of the lattice spacing a -or, equivalently, of the bare lattice coupling g 2 0 -and combining them with suitable GF scheme renormalisation factors Z m as Therefore, the complete renormalisation programme for light quark masses requires the computation of each of the three running factors to high precision, as well as the determination of Z m for the regularisation eventually employed in the computation of m i (g 2 0 ), within the appropriate range of values of g 2 0 .
3 Running in the high-energy region  also simulated an extra finer lattice with L/a = 16, in order to have a stronger crosscheck of our control over continuum limit extrapolations in the less favourable case. The value of the hopping parameter κ is tuned to its critical value κ c , so that the bare O(a)-improved PCAC mass vanishes 5 at the corresponding value of β = 6/g 2 0 . Simulations have been carried out using the plaquette gauge action [67], and an O(a)-improved fermion action [68] with the nonperturbative value of the improvement coefficient c sw [69], and the one-loop [70] and twoloop [71] values, respectively, of the boundary improvement countertermsc t and c t . All the simulations in this paper were performed using a variant of the openQCD code [72,73].
The data for Σ P can be corrected by subtracting the cutoff effects to all orders in a and leading order in g 2 0 , using the one-loop computation of [64], viz.
where c(a/L) is given in Table 1. This correction is bigger than our statistical uncertainties for L/a = 6, but smaller than the ones for L/a > 6. As will be discussed below, the scaling properties of Σ I P are somewhat better than those of the unsubtracted Σ P -and, more importantly, the study of the impact of the perturbative subtraction will allow us to assign a solid systematic uncertainty related to the continuum limit extrapolation.
The results of our simulations are summarised in Table 6. Alongside the results for Z P at each simulation point, we quote the corresponding values of Σ I P . The first uncertainty is statistical, while the second one is an estimate of the systematic uncertainty due to O(a n>2 ) cutoff effects, given by the difference of the one-loop corrected and uncorrected values of the SSF, Σ I P − Σ P .

Determination of the anomalous dimension
Once the lattice step scaling function Σ P (u, a/L) is known, we are in a position to determine the RG running of the light quark masses between the hadronic and electro-weak energy scales.
This we do using four methods; though equivalent in theory, they are distinct numerical procedures. Thus they give us insight into the magnitude of the systematic errors of our results. Two of these procedures consist in extrapolations of Σ P (u, a/L) to the continuum limit. Knowledge of the continuum SSF σ P (u) is adequate for controlling the RG evolution between energy scales. The other two procedures essentially extract the quark mass anomalous dimension τ from σ P , using Eq. (2.9b). In this way we have multiple crosschecks on the final result.
The first procedure is labelled as σ P :u-by-u. It starts with the continuum extrapolation of Σ I P (u, a/L) at fixed u, using the ansatz With u held constant, σ P and ρ P are fit parameters. A detailed study shows that when the data for Σ I P (u, a/L) are extrapolated to the continuum linearly in (a/L) 2 , the effect of the subtraction of cutoff effects at one-loop becomes noticeable, cf. Fig. 1. The fits to the unsubtracted values of Σ P have a total χ 2 /dof = 12.9/9, while the fits to the subtracted data have χ 2 /dof = 8.6/9 (with "total" above meaning χ 2 /dof, summed over all fits). Based on this observation, we add the systematic uncertainty quoted in Table 6 in quadrature to the statistical uncertainty of Σ I P , and use the result as input for our fits. This procedure increases the size of the uncertainties of our continuum-extrapolated values by about a 20-30%. Table 6 quotes σ P (u) results at the eight values of the coupling, as well as the respective slopes ρ P , from this conservative analysis.
The next step in this procedure consists in fitting the eight extrapolated σ P (u) results to a polynomial of the form so as to have a continuous expression for σ P (u). The leading non-trivial coefficient is always set to the perturbative universal prediction s 1 = −d 0 ln(2). The O(u 2 ) parameter can either be left as a free fit parameter or held fixed to the perturbative value Higher-order coefficients c n>2 are free fit parameters. We label as FITA the one with a free c 2 and as FITB the one with c 2 = s 2 . The series expansion of Eq. (3.5) is truncated either at n s = 4 or n s = 5.
Finally, the resulting continuous function for σ P (u) is readily calculated for the coupling values provided by the recursion and the RG evolution can be determined by the renormalised-mass ratios at different scales cf. Eq. (2.10). Note that this procedure is the one previously employed for the determination of the running in the N f = 0 [35] and N f = 2 [36] cases. Our second procedure, labelled as σ P :global, is a global analysis of our data, in which the Σ I P (u, a/L)-extrapolation is performed with respect to both variables u and a/L. We extrapolate our dataset using Eq. (3.4), with σ P (u) expanded according to Eq. (3.5), and ρ P (u) expanded according to Since our data have been corrected for cutoff effects up to one-loop, we consistently drop terms u 0 and u 1 from the above polynomial. This series expansion is truncated at either n ρ = 2 or n ρ = 3. The series expansion of Eq. (3.5), is truncated either at n s = 4 or n s = 5. Again, the choice of c 2 is labelled as FITA (if it is left as a free fit parameter) or FITB (if c 2 = s 2 ). Once σ P (u) has been obtained from the global extrapolation, the RG running is determined from Eq. (3.7), just like in procedure σ P :u-by-u.
The third procedure, labelled as τ :u-by-u, starts off just like σ P :u-by-u: at constant u, we fit the datapoints Σ I P (u, a/L) with Eq. (3.4), obtaining σ P (u). Then the continuum values σ P (u) are fit with Eq. (2.9b), where in the integrand we use the results of Sect. 2 for the β function (cf. Eqs. (2.34-2.36)) and the polynomial ansatz for the quark mass anomalous dimension. We fix the two leading coefficients to the perturbative asymptotic predictions t 0 = d 0 ≈ 0.05066 and t 1 = d 1 ≈ 0.002969 (this is labelled as FITB).
The coefficients t n>1 are free fit parameters and the series is truncated at n s = 2, . . . , 5.
Having thus obtained an estimate for the anomalous dimension τ (u), we arrive at another determination of the renormalised-mass ratios, using the expression with the couplings determined through Eq. (3.6).
Our fourth procedure, labelled as τ :global, consists in performing the continuum extrapolation of Σ P (u, a/L) and the determination of the anomalous dimension τ (u) simultaneously. Once more, Σ P (u, a/L) is treated as a function of two variables. This approach is based on the relation with ρ P (u) and τ (u) parameterised by the polynomial ansätze (3.8) and (3.9) respectively, and the β function being again provided by Eqs. (2.34-2.36). In practice the ρ P (u)-series is truncated at n ρ = 2, 3. Again we account for the one-loop correction of our data for cutoff effects by consistently dropping terms u 0 and u 1 from the ρ P (u)-polynomial. 6 For the τ -series, we always fix the leading universal coefficient to the perturbative asymptotic prediction t 0 = d 0 ≈ 0.05066, while we either leave the rest of the coefficients to be determined by the fit (labelled FITA), or fix the O(x 4 ) coefficient to the perturbative prediction t 1 = d 1 ≈ 0.002969 (labelled FITB). Like in the previous τ :u-by-u procedure, having obtained an estimate for the anomalous dimension τ (u), we again arrive at an expression for the renormalised-mass ratios using Eq. (3.10).
The main advantage of the two u-by-u analyses is that one has full control over the continuum extrapolations, since they are performed independently of the determination of the anomalous dimension. The disadvantage is that having to fit, for most u-values, three datapoints with the two-parameter function (3.4), we are forced to include our L/a = 6 results, which are affected by the largest discretisation errors. As far as the two global analyses are concerned, they have two advantages. First, one can choose not to include the coarser lattice 10 100  Table 7.
The scale setting to obtain µ in physical units uses µ 0 ≈ 4 GeV, obtained from Λ QCD in [32]. Perturbative predictions at different orders are also shown for comparison.
data points, i.e., one can fit only to the data with L/a > 6. This provides an extra handle on the control of cutoff effects. Second, slight mistunings in the value of the coupling u can be incorporated by the fit. In order to have a meaningful quantitative comparison of the four procedures, we display in Table 7 results for R (k) , obtained from all four methods and for a variety of fit ansätze. In general, the agreement is good, though it is clear that the fit quality improves when the data with L/a = 6 are discarded. Fits for τ that do not use the known value for t 1 tend to have larger errors, as expected, since the asymptotic perturbative behavior is not constrained by the known analytic results. We therefore focus on FITB. Concerning fits from procedures σ P :u-by-u and σ P :global, the parameterisation with n s = 5 tends to provide a better description of our data. Anyway, the key point is that the analysis coming from the recursion of the step scaling functions is in good agreement with that coming from a direct determination of the anomalous dimension.
Based on this discussion, we quote as our preferred result the determination of τ from a global FITB without the L/a = 6 data, and n s = n ρ = 2. We refer to this determination as FITB* in the following. The result for the anomalous dimension at high energies is thus given  by t n g 2n ; 12) and is illustrated in Fig. 3. We stress that this result comes from a conservative approach, since we drop the L/a = 6 data, and the statistical error of the points is increased by the value of the one-loop prediction for cutoff effects. To illustrate the good agreement of the various determinations of the running mass, Fig. 2 illustrates the comparison between our preferred fit for the anomalous dimension and the values obtained from the recursion based on σ P , demonstrating the excellent level of agreement between otherwise fairly different analyses, as well as the comparison with perturbative predictions.

Connection to RGI masses
In order to make the connection with RGI masses, as spelled out in our strategy in Sec. 2, we could apply NLO perturbation theory directly at an energy scale µ pt at the higher end of our data-covered range -e.g., the one defined by g 2 SF (µ pt ) = 1.11 -to compute M/m(µ pt ), and multiply it with the non-perturbative result for m(µ pt )/m(µ 0 /2). We however observe that the perturbative description of the running above µ pt can be constrained by employing our non-perturbatively determined form of τ , which at very small values of the coupling agrees with the asymptotic perturbative expression by construction. It is thus possible to directly compute the quantity using as input the τ function from different fits, in order to assess the systematic uncertainty of the procedure. The result of this exercise for a selection of global fits for τ , both with and without the L/a = 6 data, is provided in Table 2. The agreement among different parameterisations of the anomalous dimension is very good. The value coming from our preferred fit is M m(µ 0 /2) = 1.7505(89) . (3.14) This is our final result coming from the high-energy region, bearing a 0.5% error. As stressed above, our error estimates can be deemed conservative; an extra-conservative error estimate could be obtained by adding in quadrature the maximum spread of central values in Table 2, which would yield a 0.7% final uncertainty. We however consider the latter an overestimate, and stick to Eq. (3.14) as our preferred value.
4 Running in the low-energy region As already explained, at energies µ < µ 0 /2 it is convenient to change to the GF scheme. The objective of the low-energy running is to compute the ratio m(µ 0 /2)/m(µ had ), that, together with the ratio M/m(µ 0 /2) of eq. (3.13), will provide the total running factor. We have again determined the step scaling function Σ P of Eq. Note that the lattice sizes are larger than in the high-energy region. This allows to better tackle cutoff effects, which are expected to be larger because of the stronger coupling, and the larger scaling violations exhibited by u GF with respect to u SF [29,31,61]. Simulations have been carried out using a tree-level Symanzik improved (Lüscher-Weisz) gauge action [74], and an O(a)-improved fermion action [68] with the non-perturbative value of the improvement coefficient c sw [75] and one-loop values of the coefficients c t ,c t for boundary improvement counterterms [66,76,77]. The chiral point is set using the same strategy as before, cf. Section 3.
In contrast to the high-energy region, here the coupling constant g 2 GF and Z P are measured on the same ensembles, and hereafter we take the resulting correlations into account in our analysis. As discussed in Section 2, computations are carried out at fixed topological charge Q = 0. The main motivation for this is the onset of topological freezing [78] within the range of bare coupling values covered by our simulations; setting Q = 0 allows to circumvent the large computational cost required to allow the charge to fluctuate in the ensembles involved. The projection is implemented as proposed in [31,79]. In practice, this is only an issue for the finest ensembles at the largest values of u GF -for u GF 4 no configurations with nonzero charge have been observed in simulations where the projection is not implemented. 7 The results of our simulations are summarised in Table 8. Note that, contrary to the highenergy region, here the value of u GF is not exactly tuned to a constant value for different L/a. In practice, the slight mistuning is not visible when extrapolating our data to the continuum, but our data set simply begs to be analysed using the global approach described in Section 3. This approach only requires to have pairs of lattices of sizes L/a and 2L/a simulated at the same values of the bare parameters.
Moreover, there is no guarantee that when computing m(µ 0 /2)/m(µ had ) the scale factor µ had /µ 0 is an integer power of two. This speaks in favour of performing the analysis using the  Table 3: Results for the running factor m(µ 0 /2)/m(µ had ) anomalous dimension. As in the previous section, we will use the available information on the β function, eq. (2.36). We will parameterise the ratio of RG functions as and determine the fit parameters f n via a fit to the usual relation Once more, ρ(u) describes the cutoff effects in σ P (u). When the fit parameters f n are determined, we can reconstruct the anomalous dimension thanks to the relation Recall that the parameters p n define our β function in Eq. (2.36).
We have tried different ansätze, and as long as n r > 1 one gets a good description of the data. The largest systematic dependence in our results comes from the functional form of ρ(u). Various simple polynomials have been used, as described in Table 3. As the reader can check, all fit ansätze produce results for m(µ 0 /2)/m(µ had ) that agree within one sigma. As final result we choose the fit with n r = 3 and ρ(u) = ρ 0 + ρ 1 u + ρ 2 u 2 , which has the best χ 2 /dof, and yields m(µ 0 /2)/m(µ had ) = 0.5226 (43) .   The functional form of the anomalous dimension with its uncertainty can thus be easily reproduced, and is shown in Fig. 3. Together with the result for τ in the high-energy region discussed in Section 3, and the result in Eq. (3.14), one has then all the ingredients needed to reconstruct the scale dependence of light quark masses in the full energy range relevant to SM physics, as shown in Fig. 4. Let us end this section by pointing out that the coefficient f 0 = 1.2(2) is compatible within 1.5σ with the one-loop perturbative prediction d 0 /b 0 . This is quite surprising, especially taking into account that the β function is not compatible with the one-loop functional form with coefficient b 0 . This in turn means that also the anomalous dimension τ is poorly approximated by one-loop perturbation theory. One is thus driven to conclude that the agreement of f 0 with LO perturbation theory is only apparent.

Hadronic matching and total renormalisation factor
The last step in our strategy requires the computation of the PCAC quark mass renormalisation factors Z m at hadronic scales, cf. Eq. (2.38). The latter can be written as Since the values of the axial current normalisation Z A are available from a separate computation [81][82][83] [32]. Perturbative predictions at different orders are also shown for comparison.
setup [84][85][86]), in order to obtain Z m we just need to determine Z P (g 2 0 , aµ had ) at a fixed value of µ had for changing bare coupling g 2 0 . The values of g 2 0 have to be in the range used in large-volume simulations; for practical purposes, we will thus target the interval β = 6/g 2 0 ∈ [3.40, 3.85] currently covered by CLS ensembles [47,87].
Our strategy proceeds as follows. We first choose a value of u had = g 2 GF (µ had ) such that the relevant g 2 0 range is covered using accessible values of L/a. The precise value of u had is fixed by simulating at one of the finest lattices. Finally, other lattice sizes are simulated, such that we can obtain an interpolating formula for Z P (g 2 0 , aµ had ) as a function of g 2 0 along the line of constant physics fixed by u had . We have set u had = 9.25, fixed by simulating on an L/a = 20 lattice at β = 3.79. Lattice sizes L/a = 16, 12, 10 have then been used at smaller β, and two L/a = 24 lattices have been added so that the finest β = 3.85 point can be safely interpolated. Using the results in [32], this value of u had corresponds to an energy scale µ had = 233(8) MeV.
Our simulation results are summarised in Table 4. Deviations from the target values of u had induce a small but visible effect on Z P . The measured values of the PCAC mass are often beyond our tolerance |Lm| 0.001 (cf. App. A), especially for the smaller lattice sizes. This is a consequence of the fact that the interpolating formula for κ c as a function of g 2 0 loses precision for L/a = 10, and of the small cutoff effect induced by computing at zero topological charge (which is part of our renormalisation condition). Note that the values of Z P that enter the fit function are never further away than two standard deviations from the value on the lattice L/a β κ u GF Lm Z P (g 2 0 , L/a) # cfg 10 3.4000 0.1368040500 9.282(39) −0.0221 (31) Table 4: Results for Z P in the GF scheme, used to determine quark mass renormalization constants at u had = 9.25. Alongside the values of u GF and Z P , we also quote the value of the PCAC mass m in units of the physical lattice length, and the statistics for each ensemble.
that defines the line of constant physics. The fitted dependences on g 2 0 , u GF , and Lm are thus very mild.
The measured values of Z P are fitted to a function of the form where u had = 9.25 and β 0 = 3.79. The terms with coefficients t ij parameterise the small deviations from the intended line of constant physics described above, while Z had P is the interpolating function we are interested in. Note that the ensembles are fully uncorrelated among them, but within each ensemble the values of u GF , Lm, and Z P are correlated, and we take this into account in the fit procedure. We have performed fits setting to zero various subsets of t ij coefficients; we quote as our preferred result the one with t 20 = t 11 = 0, for which χ 2 /dof = 4.43/6, and the coefficients for the interpolating function Z had which has a 0.96% precision; recall µ had = 233(8) MeV. It is important to stress that Eq. (5.4) holds in the continuum, i.e., it is independent of any detail of the lattice computation. We can then take our interpolating formula for Z P and the known results for Z A , and build the total factor Z M that relates RGI masses to bare PCAC masses computed with a non-perturbatively O(a)-improved fermion action and a tree-level Symanzik improved gauge action, .

(5.5)
Recall that the dependence of Z M on µ had has to cancel exactly, up to residual terms contributing to the cutoff effects of Z M (g 2 0 )m(g 2 0 ). Using as input the values of Z A from the chirally rotated SF setup [84][85][86], we quote the interpolating function   Recall that the quoted errors do not contain the contribution from the running factor M/m(µ had ). The correlations between the errors at different β can be readily obtained from the covariance matrices provided in Eqs. (5.6b) and (5.3b), respectively.
The quoted errors only contain the uncertainties from the determination of Z A and Z P at the hadronic scale. As remarked in [35], the error of the total running factor M/m(µ had ) in Eq. (5.4) has to be added in quadrature to the error in the final result for the RGI mass, since it only affects the continuum limit, and should not be included in the extrapolation of data for current quark masses to the continuum limit.
Finally, we note that our results can be also used to obtain renormalised quark masses when a twisted-mass QCD Wilson fermion regularisation [88] is employed in the computation. In that case the bare PCAC mass coincides with the bare twisted mass parameter, which can be renormalised with Z tm m (g 2 0 , aµ had ) = 1 Z P (g 2 0 , aµ had ) .

(5.7)
The total renormalisation factor then acquires the form , (5.8) and values can be obtained by directly using our interpolating formula for Z had P . The same comment about the combination of uncertainties as above applies. The values of Z M and Z tm M at the β values of CLS ensembles are provided in Table 5.
Eqs. (5.4,5.3,5.6) are the final, and most important, results of this work. We stress once more that the result for M/m(µ had ), which is by far the most computationally demanding one, holds in the continuum, and is independent of the lattice regularisation employed in its determination, as well as any other computational detail. The expressions for Z M and Z tm M , on the other hand, depend on the action used in the computation, and hold for a nonperturbatively O(a)-improved fermion action and a tree-level Symanzik improved gauge action. Repeating the computation for a different lattice action would only require obtaining the values of Z A and Z P in the appropriate interval of values of β, at small computational cost.

Conclusions
In this paper we have performed a fully non-perturbative, high-precision determination of the quark mass anomalous dimension in N f = 3 QCD, spanning from the electroweak scale to typical hadronic energies. Alongside the companion non-perturbative determination of the β function in [29,31], this completes the first-principles determination of the RG functions of fundamental parameters for light hadron physics. Together with the determination of the Λ QCD parameter [32] and the forthcoming publication of renormalised quark masses [89], a full renormalisation programme of N f = 3 QCD will have been achieved. For the purpose of the latter computation, we have also provided a precise computation of the matching factors required to obtain renormalised quark masses from PCAC bare quark masses obtained from simulations based on CLS N f = 2 + 1 ensembles [90,91]. The total uncertainty introduced by the matching factor to RGI quark masses is at the level of 1.1%.
A slight increase in this precision is achievable within the same framework. This would require however a significantly larger numerical effort, adding larger lattices and hence finer lattice spacings to the continuum limit extrapolation, and augmenting the precision for the tuning of bare parameters. One lesson from the present work is that the 1% ballpark is not much above the irreducible systematic uncertainty achievable with the methodology employed. Improvements of the latter will thus be necessary to reduce the uncertainty purely due to renormalisation to the few permille level. It is important to stress that the uncertainty is completely dominated by the contribution from the non-perturbative running at low energies; another lesson from the present work is that, at this level of precision, use of perturbation theory in the few-GeV region is not necessarily satisfactory.
Apart from increasing the precision of the computation, one obvious next step is the inclusion of heavier flavours. This would ideally result in reaching sub-percent renormalisationrelated uncertainties in first-principles computations of the charm and bottom quark masses, whose values play a key role in frontier studies of B-and Higgs physics. First steps in this direction, at the level of the computation of the running coupling, have already been taken [92][93][94]. 2016-0597. C.P. and D.P. acknowledge the kind hospitality offered by CERN-TH at various stages of this work. The simulations reported here were performed on the following HPC systems: Altamira, provided by IFCA at the University of Cantabria, on the FinisTerrae-II machine provided by CESGA (Galicia Supercomputing Centre), and on the Galileo HPC system provided by CINECA. FinisTerrae II was funded by the Xunta de Galicia and the Spanish MINECO under the 2007-2013 Spanish ERDF Programme. Part of the simulations reported in Section 5 were performed on a dedicated HPC cluster at CERN. We thankfully acknowledge the computer resources offered and the technical support provided by the staff of these computing centers. We are grateful to J. Koponen for helping us complete some of the simulations discussed in Section 5 in the concluding stages of this project. We are indebted to our fellow members of the ALPHA Collaboration for many valuable discussions and the synergies developed with other renormalisation-related projects; special thanks go to M. Dalla Brida, T. Korzec, R. Sommer, and S. Sint.

A.1 Tuning of the critical mass
The tuning of the chiral point, as part of setting our lines of constant physics, is treated in detail in [66]. Briefly, a set of tuning runs at various values of β and κ are used to compute the PCAC mass m in Eq. (2.11), and interpolate at fixed β values for κ c such that |Lm| ≤ 0.001, with an uncertainty of at most the same order. This implies that the values of the quark mass at which the renormalisation condition in Eq. (2.26) is imposed are not exactly zero.
In order to assess the relevant systematics, we have performed dedicated simulations at the strongest coupling covered in the SF scheme, u SF = 2.0120, for which Σ P (2.0120, L/a = 6) has been computed at the value of β = 6.2735 indicated in Table 6, and four different values of κ, besides the nominal one for κ c given in Table 6. This is expected to be the simulation point within the high-energy regime where systematics may have a stronger impact. The result of this exercise is shown in Fig. 6. A linear fit to the data allows to estimate the slope coefficient which can then be used to assign a systematic uncertainty to Σ P as where tol(Lm) = 0.001 is our tolerance for defining the critical point. We obtain ρ κc = −0.15 (15), which yields for the systematic uncertainty δ κc Σ P = 0.00015 (15). 8 Besides being compatible with zero within 1σ, the central value is more than four times smaller than the statistical uncertainty quoted for Σ P (2.0120, L/a = 6) in Table 6. For larger lattices or smaller renormalised couplings these systematics are expected to decrease further. On the other hand, they will increase as the renormalised coupling increases, but there is no obvious reason why its relative size with respect to the statistical uncertainty will grow as well. We thus conclude that the systematic uncertainty related to the tuning to zero quark mass is negligible at the level of precision we attain in the computation of σ P .

A.2 Tuning of the gauge coupling
Our lines of constant physics are formally defined by a nominal value of either the SF or the GF coupling, that is to be kept fixed in the lattices at which the denominator of Eq. (2.14) is computed. In practice, this is so only within some finite precision, due to two different reasons: (i) Both couplings u SF or u GF are computed with finite precision. The shadowed band is a linear fit to the data.
In order to assess whether the finite precision on the value of the gauge coupling has an impact on the computation of the continuum σ P , we have: (i) Repeated the continuum-limit extrapolations of Σ P at fixed u, introducing a horizontal error on the value of (a/L) 2 propagated from the uncertainty on u at each value of a/L. To that purpose one can use either the perturbative relation between u and L, or the known non-perturbative β functions (cf. Section 2), with little practical differences.
(ii) Repeated fits to the continuum-extrapolated σ P as a function of u, introducing an uncertainty on u that covers the spread of the computed values along that particular line of constant physics.
In either case, we have found that the impact of introducing the additional uncertainties in the final description of σ P are completely negligible within our current level of precision. This source of systematic uncertainty is therefore ignored in our final analyses.

A.3 Perturbative values of boundary improvement coefficients
In our computation, we use perturbative values for the coefficients c t andc t that appear in Schrödinger Functional boundary O(a) improvement counterterms, employing the highest available order for the relevant lattice action. In computations in the high-energy region, where the plaquette gauge action is used, we are able to use the corresponding two-loop value of c t [63]. In the low-energy region, where the gauge action is tree-level Symanzik improved, the two-loop coefficient is not known, and we take the one-loop value. In the case ofc t , we employ the one-loop value [70] throughout. Contributions from these boundary counterterms to the step scaling function σ P start at two-loop order in perturbation theory [64], and the effects of perturbative truncation in the values of c t ,c t are therefore expected to be small. A careful study of these systematic uncertainties in the SF scheme was carried out in the N f = 0 [35] and N f = 2 [36] computations. Especially in the former case, a precise statement could be made that perturbative truncation effects do not change the result for the continuum limit of Σ P even at the largest values of the coupling. In our computation, we have performed a dedicated analysis at two values of the coupling, to check the size of the resulting effects in both the SF and GF schemes. In order to have a quantitative handle on the effect from a shift on boundary improvement coefficients, let us formally expand Z P , considered as a function of, e.g., c t , in a power series of the form where ∆c t is the deviation with respect to the value of c t at which Z P is computed. The factor (a/L) is made explicit to stress that the perturbative truncation is leaving uncancelled O(a) terms, which we are parameterising. By simulating at a number of values of c t , keeping all other simulation parameters fixed, it is possible to estimate the slope (∂Z P /∂c t ). A systematic uncertainty can then be assigned to Z P as where δc t is some conservative estimate of the perturbative truncation error. Linear error propagation then yields the corresponding systematic uncertainty on Σ P = Z P (2L)/Z P (L) as The systematic uncertainty due to the truncation in the value ofc t can be estimated in exactly the same way.
We have performed dedicated simulations to estimate (∂Z P /∂c t ) and (∂Z P /∂c t ) in the SF scheme at u = 2.0120 and in the GF scheme at u = 4.4901, considering several values of c t and c t in an interval given by artificially augmenting the size of the perturbative correction to the tree-level value 1 by factors of up to 4. The simulations are performed in L/a = 6 and L/a = 8 lattices, which should be affected by the largest uncertainty. The results are illustrated in Fig. 7. By fitting the results for Z P linearly in the value of the improvement coefficient we find In all cases, these figures are much smaller than the quoted statistical error of Σ P . This justifies neglecting this source of systematic uncertainty in our analysis.

Appendix B Simulation details
Our data is partly based on ensembles that have been produced in conjunction with the running coupling project [29,31]. The low-energy ensembles are common to both projects, such that we had to take into account correlations between Z P and u GF as explained in section 4. To reach the desired precision in the computation of the mass anomalous dimension, the statistics on some of those ensembles had to be increased. Hence, we give a comprehensive list of the corresponding simulations in Table 10 which enter our data analysis. The high-energy part is different in the respect that the computation of the SF coupling u SF , defining the lines of constant physics, is done with non-vanishing background gauge field, while Z P is better computed with zero background field. As a result we have produced an independent set of ensembles summarised in Table 9. The bare parameters are inherited from the line of constant physics condition [29,66].
Tables 9, 10 list the line of constant physics (L/a, β, κ) for several fixed values of the renormalized coupling u = g 2 (L), corresponding to a fixed scale L in the continuum. The data are complemented by the number of measurements N ms and their separation τ ms in moleculardynamics (MD) units, the measured integrated autocorrelation times for Z P , values of the dimensionless bare current quark mass Lm 1 , and the boundary-to-boundary correlator f 1 .  Table 6: Results for Z P , Σ I P , and σ P in the SF scheme. The last three columns quote the value of σ P obtained from a continuumlimit extrapolation fitting all three points linearly in (a/L) 2 , the corresponding slope parameters, and the values of χ 2 per degree of freedom (note that the numbers of degrees of freedom is always 1, save for the case u = 2.0120 where it is 2).     Table 9: Details for SF simulations of (s × L/a) 4 lattices with vanishing background field and plaquette gauge action, cf. Section 3. Each trajectory has a length of τ = 2 MD units, and lines of constant physics (fixed u SF ) are set with background field as reported in Refs. [29,96] (continues on the next page).