Non-perturbative quark mass renormalisation and running in \(N_{\scriptstyle \mathrm{f}}=3\) QCD

The European Physical Journal C, May 2018

We determine from first principles the quark mass anomalous dimension in \(N_{\scriptstyle \mathrm{f}}=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 \(\text{ 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 \(\text{ 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\%\).

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

Non-perturbative quark mass renormalisation and running in \(N_{\scriptstyle \mathrm{f}}=3\) QCD

Eur. Phys. J. C Non-perturbative quark mass renormalisation and running in Nf = 3 QCD ALPHA Collaboration I. Campos 2 P. Fritzsch 1 C. Pena 0 D. Preti 5 A. Ramos 4 A. Vladikas 3 0 Instituto de Física Teórica UAM-CSIC & Dpto. de Física Teórica, Universidad Autónoma de Madrid , Cantoblanco, 28049 Madrid , Spain 1 Theoretical Physics Department , CERN, 1211 Geneva 23 , Switzerland 2 Instituto de Física de Cantabria IFCA-CSIC , Avda. de los Castros s/n, 39005 Santander , Spain 3 INFN, Sezione di Tor Vergata c/o Dipartimento di Fisica, Università di Roma “Tor Vergata” , Via della Ricerca Scientifica 1, 00133 Rome , Italy 4 School of Mathematics, Trinity College Dublin , Dublin 2 , Ireland 5 INFN, Sezione di Torino , Via Pietro Giuria 1, 10125 Turin , Italy 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 nonperturbative 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%. 1 Introduction . . . . . . . . . . . . . . . . . . . . . . 2 Strategy . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Quark running and RGI masses . . . . . . . . . . 2.2 Step scaling functions . . . . . . . . . . . . . . . 2.3 Renormalisation schemes . . . . . . . . . . . . . 2.4 Determination of RGI quark masses . . . . . . . 3 Running in the high-energy region . . . . . . . . . . . 3.1 Determination of ZP and P . . . . . . . . . . . 3.2 Determination of the anomalous dimension . . . 3.3 Connection to RGI masses . . . . . . . . . . . . 4 Running in the low-energy region . . . . . . . . . . . 5 Hadronic matching and total renormalisation factor . . - Contents 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– 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 isation constants needed to match to a hadronic scheme at low energies. Conclusions and outlook come in Sect. 6. Several technical aspects of the work are discussed in appendices. 2 Strategy 2.1 Quark running and RGI masses QCD is a theory with Nf + 1 parameters: the coupling constant g and the Nf quark masses {mi , i = 1, . . . , Nf }. When the theory is defined using some regularisation, g and mi are taken to be the bare parameters appearing in the Lagrangian. Removing the regularisation requires to define renormalised parameters g , { mi , i = 1, . . . , Nf } 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 [ 47,48 ], and their scale evolution is given by renormalisation group equations of the form μ μ d dμ d dμ g (μ) = β(g (μ)), mi (μ) = τ (g (μ)) mi (μ), i = 1, . . . , Nf . The renormalisation group functions β and τ admit perturbative expansions of the form 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–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 high-energy 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–25 ], a non-perturbative determination is mandatory to match the current percentlevel 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 high-precision determination of the β function and the QCD parameter in Nf = 3 QCD by the ALPHA Collaboration [ 29–31 ]. We will employ the Schrödinger Functional [ 32,33 ] 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 [ 34–45 ] to compute the RG running nonperturbatively. 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 [ 46 ] can be achieved. The paper is structured as follows. In Sect. 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 renormal1 Preliminary results have appeared as conference proceedings in [ 26– 28 ]. (2.1) (2.2) (2.3) (2.4) (2.5) (2.6) with universal coefficients given by [ 49–55 ] 1 b0 = (4π )2 1 b1 = (4π )4 8 d0 = (4π )2 . 2 11 − 3 Nf , 38 102 − 3 Nf , The higher-order coefficients bn≥2, dn≥1 are instead renormalisation scheme-dependent. It is trivial to integrate Eqs. (2.1, 2.2) formally, to obtain the renormalisation group invariants (RGI) b 1 QCD = μ b0g 2(μ) − 2b102 e− 2b0g 2(μ) × exp − 0 g (μ) dg 1 1 b1 β(g) + b0g3 − b2g 0 Mi = mi (μ) 2b0g 2(μ) . 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 Mi are nonperturbatively defined via the previous expressions. It is easy to check, furthermore, that they are Nf -dependent but μindependent. They can be interpreted as the integration constants of the renormalisation group equations. Also the ratios mi (μ)/ m j (μ), i = j are scale-independent. Furthermore, the ratios Mi / mi (μ) are independent of the quark flavour i , due to the mass-independence of the quark mass anomalous dimension τ . Finally, the values of Mi 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. 2.2 Step scaling functions In our computation, we will access the renormalisation group functions β and τ through the quantities σ and σP, defined as ln 2 = − σP(u) = exp − √u 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. σ (u) = g 2(μ/2) u=g 2(μ) σP(u) = mi (μ) mi (μ/2) u=g 2(μ) . In order to compute the SSF σP, we define renormalised quark masses through the partially conserved axial current (PCAC) relation, ∂μ( AR)μ = ( mi + m j ) PRij , i j where the renormalised, non-singlet (i = j ) axial current and pseudoscalar density operators are given by ( AR)iμj (x ) = ZAψ¯ i (x )γμγ5ψ j (x ), ( PR)i j (x ) = ZPψ¯ i (x )γ5ψ j (x ). In these expressions ZP is the renormalisation constant of the pseudoscalar density in the regularised theory, and ZA is the finite axial current normalisation, required when the QCD regularisation breaks chiral symmetry, as with lattice Wilson fermions. Equation (2.11) implies that, up to the finite current normalisation, current quark masses renormalise with ZP−1. Therefore, the SSF σP of Eq. (2.9b) can be obtained by computing ZP at fixed bare gauge coupling g02 – 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. σP(u) = al→im0 P(g02, aμ) ≡ P(g02, aμ) ZP(g02, aμ/2) ZP(g02, aμ) , g 2(μ)=u , where g02 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 τ (g¯) from Eq. (2.9b). Note that in order to do so, we need the RG function β(g¯). As we will note later, the non-perturbative determination of the β function has already been done in our schemes of choice [ 29,30 ]. In practice, given the β function (and hence σ (u)), and the lattice results for P(g02, 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 ). (2.11) (2.12) (2.13) (2.14) (2.10) 2.3 Renormalisation schemes 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 [ 34 ]. 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 control the connection between hadronic observables and RGI quantities, we will use intermediate finitevolume renormalisation schemes that allow to define fully non-perturbative renormalised gauge coupling and quark 2 Note that the assumption of a lattice regularisation in Eqs. (2.11, 2.14) is inessential: the construction can be applied to any convenient regularisation, provided currents are correctly normalised, and σP is obtained by removing the regulator at constant physics. masses. For that purpose, ZP will be defined by a renormalisation condition imposed using the Schrödinger Functional (SF) [ 32,33 ]. In the following, we will adopt the conventions and notations for the lattice SF setup introduced in [56]. 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 x0 = 0 and x0 = T . The boundary condition for gauge fields has the form 0 1 Uk (x )|x0=0 = P exp a dt Ck (x + (1 − t )akˆ ) , where kˆ is a unit vector in the direction k, P exp is a pathordered exponential, and Ck is some smooth gauge field. A similar expression applies at x0 = T in terms of another field Ck . Fermion fields obey the boundary conditions P+ψ (x )|x0=0 = ρ(x), P−ψ (x )|x0=T = ρ (x), ψ¯ (x ) P− x0=0 = ρ¯(x), ψ¯ (x ) P+ x0=T = ρ¯ (x), (2.17) (2.16) with P± = 21 (1 ± γ0). Gauge fields are periodic in spatial directions, whereas fermion fields are periodic up to a global phase ψ (x + Lkˆ ) = eiθk ψ (x ), ψ¯ (x + Lkˆ ) = ψ¯ (x )e−iθk . The SF is the generating functional Z[C, ρ¯, ρ; C , ρ¯ , ρ ] = D[U ]D[ψ ]D[ψ¯ ] e−S[U,ψ¯ ,ψ], (2.18) (2.19) (2.21) 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 1 Z O = D[U ]D[ψ ]D[ψ¯ ] Oe−S[U,ψ¯ ,ψ] , ρ¯=ρ=ρ¯ =ρ =0 (2.20) where O can involve, in particular, the “boundary fields” δ δ ζ (x) = δρ(x) , ζ¯ (x) = − δρ¯(x) ; δ δ ζ (x) = δρ (x) , ζ¯ (x) = − δρ¯ (x) . 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 [ 32,57,58 ]. In the case of the SF, it has been shown in [ 59 ] that no additional counterterms are needed with respect to the periodic case, except for one boundary term that amounts to d3x Pi j (x ) O ji , 1 f1 = − 3 O i j O ji , 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 finitevolume renormalisation schemes, where the renormalisation scale is identified with the inverse box length, 1 μ = L . (2.15) The renormalisation of the pseudoscalar density, and hence of quark masses, is treated in the same way as in [ 34 ]. We introduce the SF correlation functions (2.22) (2.23) (2.24) (2.25) (2.26) (2.27) where Pi j is the unrenormalised pseudoscalar density, and O, O are operators with pseudoscalar quantum numbers made up of boundary quark fields O O i j i j 1 = L3 1 = L3 d3 y d3 y d3z ζ¯i (y)γ5ζ j (z), d3z ζ¯i (y)γ5ζ j (z). The pseudoscalar renormalisation constant is then defined as fP(L/2) ZP √3 f1 = fP(L/2) √3 f1 tree level , 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: 1 T = L , C = C = 0, θk ≡ θ = 2 . 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,30 ], 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 [32]. Below that scale, we will use the gradient flow (GF) coupling defined in [ 60 ]. As discussed in [ 61 ], 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 ZP is always given by Eq. (2.26). The value of μ0 is implicitly fixed by g S2F(μ0) = 2.0120, where one has [ 29 ] g S2F(μ0/2) = σ (2.0120) = 2.452(11). 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 [ 30 ] g 2GF(μ0/2) = 2.6723(64). When expressed in physical units through the ratio μ0/ QCD, one finds μ0 ≈ 4 GeV [ 31 ] – 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 mi (μ), 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 ZP at all scales. This observation also provides the matching relation for the anomalous dimensions: for any fixed μ we have τSF(g S2F(μ)) = τGF(g 2GF(μ)). 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, b2, is known [ 62 ], 0.483−0.275Nf + 0.0361Nf2 −0.00175Nf3 . 1 b2 = (4π )3 Moreover, the next-to-leading order (NLO) mass anomalous dimension in the SF scheme was computed in [ 63 ], 1 d1 = (4π )2 (0.2168 + 0.084Nf ) . 3 Also here, both couplings are computed from correlation functions projected to the Q = 0 sector of the theory. (2.28) (2.29) (2.30) (2.31) (2.32) (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 bn g 2n, g 2 ∈ [0, 2.45], (2.34) 3 n=0 with b0,1 given by Eq. (2.5), b2 given by Eq. (2.32), and b3 = b3eff a fit parameter with value (4π )4b3eff = 4(3). (2.35) 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 [ 30 ] g 2 ∈ [2.1, 11.3], (2.36a) The reader should note that these values are not exactly the same as those quoted as final results in [ 30 ]. 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 ZP and g 2GF 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 [ 30 ] 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 [ 64 ] 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). 2.4 Determination of RGI quark masses In order to determine RGI quark masses, we will factor Eq. (2.8) as Mi Mi = mi (μpt) mi (μ0/2) mi (μhad) mi (μpt) mi (μ0/2) mi (μhad). (2.37) The three ratios appearing in this expression are flavourindependent running factors:4 • mi (μ0/2)/ mi (μhad) is the running between some lowenergy scale μhad and the scheme-switching scale μ0/2. It will be computed non-perturbatively in the GF scheme. • mi (μpt)/ mi (μ0/2) is the running between the schemeswitching scale μ0/2 and some high-energy scale μpt. It will be computed non-perturbatively in the SF scheme. • Mi / mi (μ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 MW . As discussed in Sect. 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 mi (μhad) at the lowenergy scale is to be computed independently from the running factors, by determining bare PCAC quark masses mi from large-volume lattice simulations at a number of values of the lattice spacing a – or, equivalently, of the bare lattice coupling g02 – and combining them with suitable GF scheme renormalisation factors Zm as mi (μhad) = ali→m0 Zm(g02, aμhad) mi (g02). 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 Zm for the regularisation eventually employed in the computation of mi (g02), within the appropriate range of values of g02. (2.38) 4 The relevant quark masses mi are always renormalised as in Eq. (2.26). This is usually called the SF renormalisation scheme but, as previously explained, in the present work we employ SF and GF renormalisation conditions for the gauge coupling. We use SF and GF to also label our quark mass anomalous dimensions. 3 Running in the high-energy region 3.1 Determination of ZP and P Our simulations in the high-energy range above μ0 have been performed at eight different values of the renormalised Schrödinger Functional coupling uSF ∈ {1.1100, 1.1844, 1.2656, 1.3627, 1.4808, for which we have determined the step scaling function P of Eq. (2.14) at three different values of the lattice spacing L/a = 6, 8, 12. At the strongest coupling uSF = 2.012 we have 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 , vanishes5 at the corresponding value of β = 6/g02. Simulations have been carried out using the plaquette gauge action [ 66 ], and an O(a)-improved fermion action [ 67 ] with the non-perturbative value of the improvement coefficient csw [ 68 ], and the one-loop [ 69 ] and two-loop [ 70 ] values, respectively, of the boundary improvement counterterms c˜t and ct. All the simulations in this paper were performed using a variant of the openQCD code [ 71,72 ]. The data for P can be corrected by subtracting the cutoff effects to all orders in a and leading order in g02, using the one-loop computation of [ 63 ], viz. P(u, a/L) PI(u, a/L) = 1 + uδP(a/L) , δP(a/L) = −d0 ln(2)c(L/a), (3.3) 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 PI 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 ZP at each simulation point, we quote the corresponding values of PI. The first uncertainty 5 Details can be found in [ 65 ]; a discussion of the systematic impact of this procedure on our data is provided in Appendix A1. Table 1 One-loop cutoff effects in P in the SF scheme L/a is statistical, while the second one is an estimate of the systematic uncertainty due to O(an>2) cutoff effects, given by the difference of the one-loop corrected and uncorrected values of the SSF, PI − P. 3.2 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 PI(u, a/L) at fixed u, using the ansatz PI(u, a/L) = σP(u) + ρP(u) L . With u held constant, σP and ρP are fit parameters. A detailed study shows that when the data for PI(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 PI, 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 2030%. Table 6 quotes σP(u) results at the eight values of the coupling, as well as the respective slopes ρP, from this conservative analysis. 0.96 0.955 so as to have a continuous expression for σP(u). The leading non-trivial coefficient is always set to the perturbative universal prediction s1 = −d0 ln(2). The O(u2) parameter can either be left as a free fit parameter or held fixed to the perturbative value c2 = s2 = −d1 ln(2) + ( 1 d2 2 0 − b0d0) ln(2)2. Higher-order coefficients cn>2 are free fit parameters. We label as FITA the one with a free c2 and as FITB the one with c2 = s2. The series expansion of Eq. (3.5) is truncated either at ns = 4 or ns = 5. Finally, the resulting continuous function for σP(u) is readily calculated for the coupling values provided by the recursion (3.5) g S2F(μ0) = 2.012, uk = g S2F(2k μ0) , and the RG evolution can be determined by the renormalisedmass ratios at different scales R(k) m(2k μ0) ≡ m(μ0/2) = k n=0 σP(un), cf. Eq. (2.10). Note that this procedure is the one previously employed for the determination of the running in the Nf = 0 [ 34 ] and Nf = 2 [ 35 ] cases. Our second procedure, labelled as σP:global, is a global analysis of our data, in which the PI(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 ρP(u) = ρnun. nρ n=2 Since our data have been corrected for cutoff effects up to one-loop, we consistently drop terms u0 and u1 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 ns = 4 or ns = 5. Again, the choice of c2 is labelled as FITA (if it is left as a free fit parameter) or FITB (if c2 = s2). 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 PI(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 τ (x ) = −x 2 tn x 2n ns n=0 for the quark mass anomalous dimension. We fix the two leading coefficients to the perturbative asymptotic predictions t0 = d0 ≈ 0.05066 and t1 = d1 ≈ 0.002969 (this is labelled as FITB). The coefficients tn>1 are free fit parameters and the series is truncated at ns = 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 (3.6) (3.7) (3.8) (3.9) ln PI(u, a/L) − ρP(u) = − L (3.11) 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 u0 and u1 from the ρP(u)-polynomial.6 For the τ -series, we always fix the leading universal coefficient to the perturbative asymptotic prediction t0 = d0 ≈ 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 t1 = d1 ≈ 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 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 t1 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 ns = 5 6 We have also fit the unsubtracted P(u, a/L) with a similar global fit approach, obtaining compatible results. We note in passing that in this case consistency requires that only the u0 term is dropped in Eq. (3.8). μ0/2 1 Fig. 2 RG evolution factor towards high energies from the lowest energy scale μ0/2 reached with the SF scheme, given by the ratio m(μ)/ m(μ0/2). The green band is the result from our preferred determination of the anomalous dimension, Eq. (3.12), in the region covered by our data, while the red points are the values of R(k) obtained from the step scaling in factors of 2 based on σP. For the latter we have used the u-by-u FITB results with ns = 5 from Table 7. The scale setting to obtain μ in physical units uses μ0 ≈ 4 GeV, obtained from QCD in [ 31 ]. Perturbative predictions at different orders are also shown for comparison 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 ns = 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 τ (g) = −g2 tn g2n ; 2 n=0 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. 3.3 Connection to RGI masses In order to make the connection with RGI masses, as spelled out in our strategy in Sect. 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 S2F(μ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 nonperturbatively 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 M m(μ0/2) = [2b0 g S2F(μ0/2)]−d0/2b0 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 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 lownr = 2 m(μ0/2)/ m(μhad) 0.5245(41) 0.5222(43) 0.5201(45) 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. (2.14) at three different values of the lattice spacing, now using lattices of size T /a = L/a = 8, 12, 16, and double lattices of size L/a = 16, 24, 32. The bare parameters are chosen so that uGF is approximately equal to one of the seven values {2.12, 2.39, 2.73, 3.20, 3.86, 4.49, 5.29}. (4.1) 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 uGF with respect to uSF [ 29,30,60 ]. Simulations have been carried out using a tree-level Symanzik improved (Lüscher-Weisz) gauge action [ 73 ], and an O(a)-improved fermion action [ 67 ] with the non-perturbative value of the improvement coefficient csw [ 74 ] and one-loop values of the coefficients ct, c˜t for boundary improvement counterterms [ 65,75,76 ]. The chiral point is set using the same strategy as before, cf. Sect. 3. In contrast to the high-energy region, here the coupling constant g 2GF and ZP are measured on the same ensembles, and hereafter we take the resulting correlations into account in our analysis. As discussed in Sect. 2, computations are carried out at fixed topological charge Q = 0. The main motivation for this is the onset of topological freezing [ 77 ] 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 [ 30,78 ]. In practice, this is only an issue for the finest ensembles at the largest values of uGF – for uGF 4 no configurations with nonzero charge have been observed in simulations where the projection is not implemented.7 7 It is noteworthy in this context that in [ 79,80 ] the improvement coefficient cA and the normalisation parameter ZA of the axial current have been measured both for a freely varying Q and in the Q = 0 sector, for several values of the bare gauge coupling in the range listed in Table 4. Results from both definitions were found to be compatible. The results of our simulations are summarised in Table 8. Note that, contrary to the high-energy region, here the value of uGF 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 Sect. 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 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 τ (x ) f (x ) = β(x ) = x 1 nr n=0 fn x 2n , and determine the fit parameters fn via a fit to the usual relation log P(u, a/L) − ρ(u)(a/L)2 = − dx f (x ). Once more, ρ(u) describes the cutoff effects in σP(u). When the fit parameters fn are determined, we can reconstruct the anomalous dimension thanks to the relation τ (g ) f (g ) = β(g ) ⇒ τ (g ) = −g 2 Recall that the parameters pn define our β function in Eq. (2.36). We have tried different ansätze, and as long as nr > 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 nr = 3 and ρ(u) = ρ0 + ρ1u + ρ2u2, which has the best χ 2/dof, and yields m(μ0/2)/ m(μhad) = 0.5226(43). The corresponding τ (x ) function is obtained from Eq. (4.4) by using the coefficients together with the known coefficients for the β function in Eq. (2.36). The covariance between the fn parameters reads (4.5) (4.6) perturbative prediction d0/b0. This is quite surprising, especially taking into account that the β function is not compatible with the one-loop functional form with coefficient b0. 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 f0 with LO perturbation theory is only apparent. Finally, the covariance between the pi and the f j parameters is given by cov( fi , p j ) 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 Sect. 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 f0 = 1.2(2) is compatible within 1.5σ with the one-loop 0 −0.1 −0.2 −0.3 ) −0.4 −0.06 (τg −0.5 −0.08 −0.1 −0.6 −0.12 −0.7 −0.14 −0.16 −0.8 −0.18 1.5 −0.9 0 1 gS2F(μ0/2) g2GF(μ0/2) 5 Hadronic matching and total renormalisation factor The last step in our strategy requires the computation of the PCAC quark mass renormalisation factors Zm at hadronic scales, cf. Eq. (2.38). The latter can be written as Since the values of the axial current normalisation ZA are available from a separate computation [ 80–82 ] (in particular we use the precise values obtained thanks to the chirally rotated SF setup [ 83–85 ]), in order to obtain Zm we just need to determine ZP(g02, aμhad) at a fixed value of μhad for changing bare coupling g02. The values of g02 have to be in the range used in large-volume simulations; for practical purposes, we (4.7) (5.1) 1 0.9 0.8 L/a will thus target the interval β = 6/g02 ∈ [3.40, 3.85] currently covered by CLS ensembles [ 46,86 ]. Our strategy proceeds as follows. We first choose a value of uhad = g 2GF(μhad) such that the relevant g02 range is covered using accessible values of L/a. The precise value of uhad 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 ZP(g02, aμhad) as a function of g2 0 along the line of constant physics fixed by uhad. We have set uhad = 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 [ 31 ], this value of uhad corresponds to an energy scale μhad = 233(8) MeV. Our simulation results are summarised in Table 4. Deviations from the target values of uhad induce a small but visible effect on ZP. 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 g02 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 ZP that enter the fit function are never further away than two standard deviations from the value on the lattice that defines the line of constant physics. The fitted dependences on g02, uGF, and Lm are thus very mild. The measured values of ZP are fitted to a function of the form ZP(g02, uGF, Lm) = ZPhad(g02) + t10 (uGF − uhad) + t20 (uGF − uhad)2 + t01 Lm + t11 (uGF − uhad)Lm, ZPhad(g02) = z0 + z1(β − β0) + z2(β − β0)2, (5.2) where uhad = 9.25 and β0 = 3.79. The terms with coefficients ti j parameterise the small deviations from the intended line of constant physics described above, while Z had is the P interpolating function we are interested in. Note that the ensembles are fully uncorrelated among them, but within each ensemble the values of uGF, Lm, and ZP are correlated, and we take this into account in the fit procedure. We have performed fits setting to zero various subsets of ti j coefficients; we quote as our preferred result the one with t20 = t11 = 0, for which χ 2/dof = 4.43/6, and the coefficients for the interpolating function ZPhad read z0 = 0.348629, z1 = 0.020921, z2 = 0.070613, with covariance matrix cov(zi , z j ) = ⎝ ⎛ The typical precision of the values of ZP extracted from the interpolating function is thus at the few permille level. The fit is illustrated in Fig. 5. Now we can finally assemble the various factors entering Eqs. (2.37, 2.38). Using the results in Sect. 4 we can compute the running m(μ0/2)/ m(μhad) between the schemeswitching scale and hadronic matching scale, and multiply it by the value of M/ m(μ0/2) given in Eq. (3.14). Combining the errors in quadrature, we obtain which has a 0.96% precision; recall μhad = 233(8) MeV. It is important to stress that Eq. (5.4) holds in the continuum, 0.360 . Fig. 5 Projection on the (β, ZP) plane of the fit (grey band) to the results for ZP at the hadronic matching point uhad = 9.25 (filled red points). The data points have been shifted to the target value of uGF and Lm = 0 using the fit function. Horizontal error bars have been assigned to reflect the uncertainties coming from the value of uGF for each point, by defining g02/g02 = uGF/uGF. Vertical dashed lines correspond to the β values of CLS ensembles i.e., it is independent of any detail of the lattice computation. We can then take our interpolating formula for ZP and the known results for ZA, and build the total factor ZM that relates RGI masses to bare PCAC masses computed with a nonperturbatively O(a)-improved fermion action and a tree-level Symanzik improved gauge action, ZM(g02) = M ZA(g02) m(μhad) ZP(g02, aμhad) Recall that the dependence of ZM on μhad has to cancel exactly, up to residual terms contributing to the cutoff effects of ZM(g02)m(g02). Using as input the values of ZA from the chirally rotated SF setup [ 83–85 ], we quote the interpolating function ZM(g02) = Z M(0) + Z M(1)(β − 3.79) M m(μhad) × +Z M(2)(β − 3.79)2 ; Z M(1) = 0.121644, Z M(0) = 2.270073, Z M(2) = −0.464575, with covariance matrix cov(Z M(i), Z M(j)) = ⎝ ⎛ The quoted errors only contain the uncertainties from the determination of ZA and ZP at the hadronic scale. As remarked in [ 34 ], 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 (5.5) (5.6a) (5.6b) Z tm M 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 [ 87 ] 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 mtm(g02, aμhad) = 1 ZP(g02, aμhad) . Z Mtm(g02) = M 1 m(μhad) ZP(g02, aμhad) , The total renormalisation factor then acquires the form and values can be obtained by directly using our interpolating formula for ZPhad. The same comment about the combination of uncertainties as above applies. The values of ZM and Z tm M at the β values of CLS ensembles are provided in Table 5. Equations (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 ZM and Z tm, on the other hand, depend on the action used M in the computation, and hold for a non-perturbatively 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 ZA and ZP in the appropriate interval of values of β, at small computational cost. 6 Conclusions In this paper we have performed a fully non-perturbative, high-precision determination of the quark mass anomalous dimension in Nf = 3 QCD, spanning from the electroweak scale to typical hadronic energies. Alongside the companion non-perturbative determination of the β function in [ 29, 30 ], 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 [31] and the forthcoming publication of renormalised quark masses [ 88 ], a full renormalisation programme of Nf = 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 Nf = 2 + 1 ensembles [ 89, 90 ]. 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 [ 91–93 ]. Acknowledgements P.F., C.P. and D.P. acknowledge support through the Spanish MINECO project FPA2015-68541-P and the Centro de Excelencia Severo Ochoa Programme SEV-2012-0249 and SEV-20160597. C.P. and D.P. acknowledge the kind hospitality offered by CERNTH 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 Sect. 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 Sect. 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. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP3. Appendix A: Systematic uncertainties in the determination of step scaling functions 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 [ 65 ]. 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 |L m| ≤ 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, uSF = 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 1 ∂ P ρκc ≡ L ∂ m u,L , δκc P = |ρκc |tol(L m), which can then be used to assign a systematic uncertainty to P as (A.1) (A.2) where tol(L m) = 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, 8 Our value for the slope ρκc is in the same ballpark as the ones obtained in [ 94 ], where a similar study using Nf = 2 simulations was performed at values of the SF coupling uSF = 0.9793 and uSF = 2.4792, finding ρκc = −0.0755(10) and ρκc = −0.1130(27), respectively. 0.914 ) 6 = a / ,L 0.912 0 2 120. = F S (u 0.910 P Σ (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. Fig. 6 Variation of the SSF P(uSF = 2.0120, L/a = 6) with the value of the quark mass. The shadowed band is a linear fit to the data 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 uSF or uGF are computed with finite precision. (ii) The SF computations of the coupling and the corresponding tuning of β and κ are performed using independent ensembles with respect to the ones employed for the computation of ZP, for reasons explained in Sect. 2. The resulting lines of constant physics are expected to differ by O(a2) effects. 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 nonperturbative β functions (cf. Sect. 2), with little practical differences. In our computation, we use perturbative values for the coefficients ct and c˜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 ct [ 62 ]. 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 of c˜t, we employ the one-loop value [ 69 ] throughout. Contributions from these boundary counterterms to the step scaling function σP start at two-loop order in perturbation theory [ 63 ], and the effects of perturbative truncation in the values of ct, c˜t are therefore expected to be small. A careful study of these systematic uncertainties in the SF scheme was carried out in the Nf = 0 [ 34 ] and Nf = 2 [ 35 ] 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 ZP, considered as a function of, e.g., ct, in a power series of the form ∂ ZP ZP + ∂ct a L ct + · · · where ct is the deviation with respect to the value of ct at which ZP 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 ct, keeping all other simulation parameters fixed, it is possible to estimate the slope (∂ ZP/∂ct). A systematic uncertainty can then be assigned to ZP as (A.3) 0.693 ) 6 = 0.692 a / L 0, 12 0.691 02. = SFu 0.690 ( P Z where δct is some conservative estimate of the perturbative truncation error. Linear error propagation then yields the corresponding systematic uncertainty on P = ZP(2L )/ZP(L ) as δct P P 1 1 2 ZP(2L ) − ZP(L ) a L δct. (A.5) The systematic uncertainty due to the truncation in the value of c˜t can be estimated in exactly the same way. We have performed dedicated simulations to estimate (∂ ZP/∂ ct) and (∂ ZP/∂ c˜t) in the SF scheme at u = 2.0120 and in the GF scheme at u = 4.4901, considering several values (A.4) of ct 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 ZP linearly in the value of the improvement coefficient we find ∂ ZP ∂ct GF;L/a=8 If now we use values of ct and c˜t given by 1 − ctpert and 1 − c˜tpert, respectively – i.e., we assign a 100% uncertainty to the perturbative correction to the tree-level value 1 – we finally obtain 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. (A.8) Appendix B: Simulation details Our data is partly based on ensembles that have been produced in conjunction with the running coupling project [ 29,30 ]. The low-energy ensembles are common to both projects, such that we had to take into account correlations between ZP and uGF as explained in Sect. 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 uSF, defining the lines of constant physics, is done with non-vanishing background gauge field, while ZP 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,65 ]. 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 Nms and their separation τms in molecular-dynamics (MD) units, the measured integrated autocorrelation times for ZP, values of the dimensionless bare current quark mass Lm1, and the boundary-to-boundary correlator f1. Appendix C: Tables Table 6 Results for ZP, PI, and σP in the SF scheme. The last three columns quote the value of σP obtained from a continuum-limit 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) uSF uSF σP(uSF) L a min Type ns nρ FITA FITA FITB FITB FITA FITA FITB FITB FITA FITA FITB FITB FITB FITB FITB FITB FITA FITB FITB FITB FITB FITA FITB* FITB FITB FITB s1, s2 uSF s τms τ Table 8 Results for ZP and P in the GF scheme L/a β ZP(g02, 2L/a) P(g02, L/a) f1 Page 20 of 23 uSF L/a 12 12 12 6 8 6 8 6 8 6 8 6 8 12 10 12 16 β τms τ uGF Nms uGF s 1 2 1 2 1 2 1 2 1 2 1 2 Nms 5001 2000 5001 2400 4602 1404 5001 2000 5001 3000 4600 1205 1. A. Denner , S. Heinemeyer , I. Puljak , D. Rebuzzi , M. Spira , Eur. Phys. J. C 71 , 1753 ( 2011 ). arXiv: 1107 .5909 [hep-ph] 2. S. Heinemeyer et al. [LHC Higgs Cross Section Working Group]. arXiv:1307 .1347 [hep-ph] 3. L.G. Almeida , S.J. Lee , S. Pokorski , J.D. Wells , Phys. Rev. D 89 ( 3 ), 033006 ( 2014 ). arXiv: 1311 .6721 [hep-ph] 4. G.P. Lepage , P.B. Mackenzie , M.E. Peskin . arXiv: 1404 .0319 [hepph] 5. A.A. Petrov , S. Pokorski , J.D. Wells , Z. Zhang , Phys. Rev. D 91 ( 7 ), 073001 ( 2015 ). arXiv: 1501 .02803 [hep-ph] 6. S. Aoki et al., Eur. Phys. J. C 74 , 2890 ( 2014 ). arXiv: 1310 .8555 [hep-lat] 7. S. Aoki et al., Eur. Phys. J. C 77 ( 2 ), 112 ( 2017 ). arXiv: 1607 .00299 [hep-lat] 8. A. Bazavov et al. [ MILC Collaboration] . PoS CD 09 , 007 ( 2009 ). arXiv: 0910 .2966 [hep-ph] 9. C. McNeile , C.T.H. Davies , E. Follana , K. Hornbostel , G.P. Lepage [HPQCD Collaboration ] . Phys. Rev. D 82 , 034512 ( 2010 ). arXiv: 1004 .4285 [hep-lat] 10. A. Bazavov et al., PoS LATTICE 2010 , 083 ( 2010 ). arXiv: 1011 .1792 [hep-lat] 11. S. Dürr et al. [BMW Collaboration]. Phys. Lett. B 701 , 265 ( 2011 ). arXiv: 1011 .2403 [hep-lat] 12. S. Dürr et al. [BMW Collaboration]. JHEP 1108 , 148 ( 2011 ). arXiv: 1011 .2711 [hep-lat] 13. B. Blossier et al. [ETM Collaboration]. Phys. Rev. D 82 , 114513 ( 2010 ). arXiv: 1010 .3659 [hep-lat] 14. P. Fritzsch et al. [ALPHA Collaboration]. Nucl. Phys. B 865 , 397 ( 2012 ). arXiv: 1205 .5380 [hep-lat] 15. N. Carrasco et al. [ETM Collaboration]. JHEP 1403 , 016 ( 2014 ). arXiv:1308 . 1851 [hep-lat] 16. F. Bernardoni et al. [ALPHA Collaboration]. Phys. Lett. B 730 , 171 ( 2014 ). arXiv: 1311 .5498 [hep-lat] 17. N. Carrasco et al. [ETM Collaboration]. Nucl. Phys. B 887 , 19 ( 2014 ). arXiv: 1403 .4504 [hep-lat] 18. C. Alexandrou , V. Drach , K. Jansen , C. Kallidonis , G. Koutsou, Phys. Rev. D 90 ( 7 ), 074501 ( 2014 ). [arXiv: 1406 .4310 [hep-lat]] 19. B. Colquhoun , R.J. Dowdall , C.T.H. Davies , K. Hornbostel , G.P. Lepage [HPQCD Collaboration], Phys. Rev. D 91 ( 7 ), 074514 ( 2015 ). arXiv: 1408 .5768 [hep-lat] 20. B. Chakraborty et al. [HPQCD Collaboration], Phys. Rev. D 91 ( 5 ), 054508 ( 2015 ). arXiv: 1408 .4169 [hep-lat] 21. Y.B. Yang et al. [χQCD Collaboration], Phys. Rev. D 92 ( 3 ), 034517 ( 2015 ). arXiv: 1410 .3343 [hep-lat] 22. T. Blum et al. [RBC and UKQCD Collaborations], Phys. Rev. D 93 ( 7 ), 074505 ( 2016 ). arXiv: 1411 .7017 [hep-lat] 23. K.G. Chetyrkin, Phys. Lett. B 404 , 161 ( 1997 ). arXiv:hep-ph/9703278 24. J.A.M. Vermaseren , S.A. Larin , T. van Ritbergen , Phys. Lett. B 405 , 327 ( 1997 ). arXiv:hep-ph/9703284 25. J.A. Gracey , Nucl. Phys. B 662 , 247 ( 2003 ). arXiv:hep-ph/0304113 26. I. Campos et al. [ ALPHA Collaboration] . PoS LATTICE 2015 , 249 ( 2016 ). arXiv: 1508 .06939 [hep-lat] 27. I. Campos et al. [ ALPHA Collaboration] . EPJ Web Conf . 137 , 08006 ( 2017 ). arXiv: 1611 .06102 [hep-lat] 28. I. Campos et al. [ ALPHA Collaboration] . PoS LATTICE 2016 , 201 ( 2016 ). arXiv: 1611 .09711 [hep-lat] 29. M. Dalla Brida et al. [ALPHA Collaboration], Phys. Rev. Lett . 117 ( 18 ), 182001 ( 2016 ). arXiv: 1604 .06193 [hep-ph] 30. M. Dalla Brida et al. [ALPHA Collaboration], Phys. Rev. D 95 ( 1 ), 014507 ( 2017 ). arXiv: 1607 .06423 [hep-lat] 31. M. Bruno et al. [ALPHA Collaboration], Phys. Rev. Lett . 119 ( 10 ), 102001 ( 2017 ). arXiv: 1706 .03821 [hep-lat] 32. M. Lüscher , R. Narayanan , P. Weisz , U. Wolff, Nucl. Phys. B 384 , 168 ( 1992 ). arXiv:hep-lat/9207009 33. S. Sint, Nucl. Phys. B 421 , 135 ( 1994 ). arXiv:hep-lat/9312079 34. S. Capitani , M. Lüscher , R. Sommer , H. Wittig [ALPHA Collaboration], Nucl. Phys. B 544 ( 1999 ) 669 Erratum: [Nucl. Phys. B 582 ( 2000 ) 762]. arXiv:hep-lat/9810063 35. M. Della , Morte et al. [ALPHA Collaboration]. Nucl. Phys. B 729 , 117 ( 2005 ). arXiv:hep-lat/0507035 36. M. Guagnelli et al. [ALPHA Collaboration]. JHEP 0603 , 088 ( 2006 ). arXiv:hep-lat/0505002 37. F. Palombi , C. Pena , S. Sint [ALPHA Collaboration], JHEP 0603 , 089 ( 2006 ). arXiv:hep-lat/0505003 38. P. Dimopoulos et al., Phys. Lett. B 641 , 118 ( 2006 ). arXiv:hep-lat/0607028 39. P. Dimopoulos et al. [ALPHA Collaboration]. JHEP 0805 , 065 ( 2008 ). arXiv: 0712 .2429 [hep-lat] 40. F. Palombi , M. Papinutto , C. Pena , H. Wittig [ALPHA Collaboration], JHEP 0709 , 062 ( 2007 ). arXiv: 0706 .4153 [hep-lat] 41. B. Blossier , M. Della Morte , N. Garron , R. Sommer [ALPHA Collaboration], JHEP 1006 , 002 ( 2010 ). arXiv: 1001 .4783 [hep-lat] 42. F. Bernardoni et al. [ALPHA Collaboration]. Phys. Lett. B 735 , 349 ( 2014 ). arXiv: 1404 .3590 [hep-lat] 43. M. Papinutto , C. Pena , D. Preti [ALPHA Collaboration] , Eur. Phys. J. C 77 ( 6 ), 376 ( 2017 ). Erratum: [Eur. Phys. J. C 78 ( 1 ), 21 ( 2018 ). arXiv: 1612 .06461 [hep-lat] 44. P. Dimopoulos et al. [ ALPHA Collaboration] . arXiv: 1801 .09455 [hep-lat] 45. C. Pena , D. Preti [ALPHA Collaboration] . arXiv:1706 .06674 [heplat] 46. M. Bruno et al., JHEP 1502 , 043 ( 2015 ). arXiv: 1411 .3982 [hep-lat] 47. S. Weinberg, Phys. Rev. D 8 , 3497 ( 1973 ) 48. G. ' t Hooft, Nucl . Phys. B 61 , 455 ( 1973 ) 49. V.S. Vanyashin , M.V. Terent 'ev , JETP 21 , 375 ( 1965 ) 50. I.B. Khriplovich , Sov. J. Nucl. Phys . 10 ( 1969 ) 235 . (Yad. Fiz. 10 ( 1969 ) 409) 51. G. ' t Hooft, report at the Colloquium on Renormalization of Yang-Mills Fields and Applications to Particle Physics, Marseille, France, June 1972 (unpublished) 52. D.J. Gross , F. Wilczek , Phys. Rev. Lett . 30 , 1343 ( 1973 ) 53. H.D. Politzer , Phys. Rev. Lett . 30 , 1346 ( 1973 ) 54. W.E. Caswell, Phys. Rev. Lett . 33 , 244 ( 1974 ) 55. D.R .T. Jones, Nucl. Phys. B 75 , 531 ( 1974 ) 56. M. Lüscher , S. Sint , R. Sommer , P. Weisz , Nucl. Phys. B 478 , 365 ( 1996 ). arXiv:hep-lat/9605038 57. K. Symanzik , Nucl. Phys. B 190 , 1 ( 1981 ) 58. M. Lüscher , Nucl. Phys. B 254 , 52 ( 1985 ) 59. S. Sint, Nucl. Phys. B 451 , 416 ( 1995 ). 0550 - 3213 ( 95 ) 00352 - S . arXiv:hep-lat/ 9504005 60. P. Fritzsch , A. Ramos , JHEP 1310 , 008 ( 2013 ). arXiv: 1301 .4388 [hep-lat] 61. A. Ramos, PoS LATTICE 2014 , 017 ( 2015 ). arXiv: 1506 .00118 [hep-lat] 62. A. Bode et al., [ALPHA Collaboration], Nucl. Phys. B 576 , 517 ( 2000 ). (Erratum: [Nucl . Phys. B 608 ( 2001 ) 481). (Erratum: [Nucl. Phys. B 600 ( 2001 ) 453) . arXiv:hep-lat/9911018 63. S. Sint et al. [ALPHA Collaboration]. Nucl. Phys. B 545 , 529 ( 1999 ). arXiv:hep-lat/9808013 64. U. Wolff, [ALPHA Collaboration], Comput. Phys. Commun . 156 , 143 ( 2004 ). (Erratum: Comput. Phys. Commun . 176 ( 2007 ) 383) . arxiv:hep-lat/0306017 65. P. Fritzsch, T. Korzec, in preparation 66. K.G. Wilson, Phys. Rev. D 10 , 2445 ( 1974 ) 67. B. Sheikholeslami , R. Wohlert , Nucl. Phys. B 259 , 572 ( 1985 ) 68. N. Yamada et al. [ JLQCD and CP-PACS Collaborations] . Phys. Rev. D 71 , 054505 ( 2005 ). arXiv:hep-lat/0406028 69. M. Lüscher , P. Weisz , Nucl. Phys. B 479 , 429 ( 1996 ). arXiv:hep-lat/9606016 70. A. Bode et al. [ALPHA Collaboration], Nucl. Phys. B 540 , 491 ( 1999 ). arXiv:hep-lat/9809175 71. M. Lüscher , S. Schaefer , Comput. Phys. Commun . 184 , 519 ( 2013 ). arXiv: 1206 .2809 [hep-lat] 72. M. Lüscher , S. Schaefer, openQCD simulation program for lattice QCD with open boundary conditions . luscher/openQCD/ 73. M. Lüscher , P. Weisz, Phys. Lett. 158B , 250 ( 1985 ) 74. J. Bulava , S. Schaefer , Nucl. Phys. B 874 , 188 ( 2013 ). arXiv: 1304 .7093 [hep-lat] 75. S. Takeda , S. Aoki , K. Ide , Phys. Rev. D 68 , 014505 ( 2003 ). arXiv:hep-lat/0304013 76. Pol Vilaseca, private communication 77. S. Schaefer et al. [ALPHA Collaboration]. Nucl. Phys. B 845 , 93 ( 2011 ). arXiv: 1009 .5228 [hep-lat] 78. P. Fritzsch , A. Ramos , F. Stollenwerk , PoS Lattice 2013 , 461 ( 2014 ). arXiv: 1311 .7304 [hep-lat] 79. J. Bulava et al., ALPHA Collaboration. Nucl. Phys. B 896 , 555 ( 2015 ). arXiv: 1502 .04999 [hep-lat] 80. J. Bulava , M. Della Morte , J. Heitger , C. Wittemeier [ALPHA Collaboration] , Phys. Rev. D 93 ( 11 ), 114513 ( 2016 ). arXiv: 1604 .05827 [hep-lat] 81. M. Dalla Brida , S. Sint , P. Vilaseca, JHEP 1608 , 102 ( 2016 ). arXiv: 1603 .00046 [hep-lat] 82. M. Dalla Brida , T. Korzec , S. Sint , P. Vilaseca, in preparation 83. S. Sint, Nucl. Phys. B 847 , 491 ( 2011 ). arXiv: 1008 .4857 [hep-lat] 84. S. Sint , B. Leder , PoS LATTICE 2010 , 265 ( 2010 ). arXiv: 1012 .2500 [hep-lat] 85. M. Dalla Brida , S. Sint, PoS LATTICE 2014 , 280 ( 2014 ). arXiv: 1412 .8022 [hep-lat] 86. D. Mohler , S. Schaefer , J. Simeth, arXiv: 1712 .04884 [hep-lat] 87. R. Frezzotti et al. [ALPHA Collaboration]. JHEP 0108 , 058 ( 2001 ). arXiv:hep-lat/0101001 88. [ALPHA Collaboration], to appear 89. M. Bruno , T. Korzec , S. Schaefer, Phys. Rev. D 95 ( 7 ), 074504 ( 2017 ). arXiv: 1608 .08900 [hep-lat] 90. G. Herdoíza, C. Pena , D. Preti , J.Á. Romero , J. Ugarrio . arXiv: 1711 .06017 [hep-lat] 91. F. Knechtli et al., ALPHA Collaboration. Phys. Lett. B 774 , 649 ( 2017 ). arXiv: 1706 .04982 [hep-lat] 92. S. Calì , F. Knechtli , T. Korzec , H. Panagopoulos . arXiv: 1710 .06221 [hep-lat] 93. F. Knechtli , T. Korzec , B. Leder , G. Moir [ALPHA Collaboration] . arXiv:1710 .07590 [hep-lat] 94. R. Sommer, “ The critical line for the Nf = 2 O(a) improved theory” , ALPHA Collaboration internal notes ( 2002 ) 95. M. Dalla Brida et al., [ALPHA Collaboration] , in preparation

This is a preview of a remote PDF:

I. Campos, P. Fritzsch, C. Pena, D. Preti, A. Ramos, A. Vladikas. Non-perturbative quark mass renormalisation and running in \(N_{\scriptstyle \mathrm{f}}=3\) QCD, The European Physical Journal C, 2018, 387, DOI: 10.1140/epjc/s10052-018-5870-5