Nonperturbative quark mass renormalisation and running in \(N_{\scriptstyle \mathrm{f}}=3\) QCD
Eur. Phys. J. C
Nonperturbative 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 UAMCSIC & 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 IFCACSIC , 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 nonperturbative 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 finitesize 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 treelevel 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 highenergy 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 lowenergy 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 massindependent 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 nonperturbative 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 longdistance strong interaction effects from first
principles. Lattice QCD (LQCD) is therefore the bestsuited
framework for a highprecision 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 nonperturbatively.
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 highorder perturbative
estimates of the anomalous dimension of quark masses in various
renormalisation schemes exist [
23–25
], a nonperturbative
determination is mandatory to match the current
percentlevel precision of the relevant hadronic observables.
In this work we present a highprecision 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
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 wellestablished finitesize
recursion techniques [
34–45
] to compute the RG running
nonperturbatively. Our main result is a highprecision
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 higherorder coefficients bn≥2, dn≥1 are instead
renormalisation schemedependent. 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 scaleindependent. Furthermore,
the ratios Mi / mi (μ) are independent of the quark flavour
i , due to the massindependence 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 schemedependent, but the ratio
QCD/ QCD between its values in two different schemes can
be calculated exactly using oneloop 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, nonsinglet (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
massindependent 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 nonperturbatively τ (g¯)
from Eq. (2.9b). Note that in order to do so, we need the RG
function β(g¯). As we will note later, the nonperturbative
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 wellcontrolled
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
nonperturbatively 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
nonperturbative 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 spacetime 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
nontrivial 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 nonperturbative 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 nonuniversal
perturbative coefficient of the β function, b2, is known [
62
],
0.483−0.275Nf + 0.0361Nf2 −0.00175Nf3 .
1
b2 = (4π )3
Moreover, the nexttoleading 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 twoloop computation
of the finitevolume 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 highenergy
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 schemeswitching scale μ0/2.
It will be computed nonperturbatively in the GF scheme.
• mi (μpt)/ mi (μ0/2) is the running between the
schemeswitching scale μ0/2 and some highenergy scale μpt. It
will be computed nonperturbatively 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 nonperturbative 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 largevolume 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 highenergy region
3.1 Determination of ZP and
P
Our simulations in the highenergy 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 nonperturbative value of the improvement coefficient
csw [
68
], and the oneloop [
69
] and twoloop [
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
oneloop 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 Oneloop 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 oneloop 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 electroweak 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:ubyu. 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 oneloop 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 continuumextrapolated 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
nontrivial 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.
Higherorder 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
oneloop, 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:ubyu.
The third procedure, labelled as τ :ubyu, starts off just
like σP:ubyu: 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
renormalisedmass 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 oneloop
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 τ :ubyu procedure, having obtained
an estimate for the anomalous dimension τ (u), we again
arrive at an expression for the renormalisedmass ratios using
Eq. (3.10).
The main advantage of the two ubyu 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 uvalues, three datapoints with the twoparameter
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:ubyu 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 ubyu 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 oneloop 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 datacovered range – e.g., the one defined by
g S2F(μpt) = 1.11 – to compute M/ m(μpt), and multiply it
with the nonperturbative 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 highenergy region,
bearing a 0.5% error. As stressed above, our error estimates
can be deemed conservative; an extraconservative 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 lowenergy 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 highenergy
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 treelevel Symanzik improved (LüscherWeisz)
gauge action [
73
], and an O(a)improved fermion action [
67
]
with the nonperturbative value of the improvement
coefficient csw [
74
] and oneloop 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 highenergy 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 highenergy 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 oneloop functional form with coefficient b0.
This in turn means that also the anomalous dimension τ is
poorly approximated by oneloop 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 highenergy 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 oneloop
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 largevolume 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 treelevel
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 twistedmass 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 nonperturbatively
O(a)improved fermion action and a treelevel 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 nonperturbative,
highprecision determination of the quark mass
anomalous dimension in Nf = 3 QCD, spanning from the
electroweak scale to typical hadronic energies. Alongside the
companion nonperturbative determination of the β
function in [
29, 30
], this completes the firstprinciples
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 nonperturbative running at low
energies; another lesson from the present work is that, at this
level of precision, use of perturbation theory in the fewGeV
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 subpercent
renormalisationrelated uncertainties in firstprinciples 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 FPA201568541P and the Centro de
Excelencia Severo Ochoa Programme SEV20120249 and
SEV20160597. 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 FinisTerraeII 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 20072013 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 renormalisationrelated 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
ons.org/licenses/by/4.0/), 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 highenergy
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 continuumextrapolated σ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 continuumlimit 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 highenergy region, where the plaquette gauge
action is used, we are able to use the corresponding twoloop
value of ct [
62
]. In the lowenergy region, where the gauge
action is treelevel Symanzik improved, the twoloop
coefficient is not known, and we take the oneloop value. In the
case of c˜t, we employ the oneloop value [
69
] throughout.
Contributions from these boundary counterterms to the step
scaling function σP start at twoloop 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 treelevel 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 treelevel 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 lowenergy 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 highenergy part is different in the respect that the
computation of the SF coupling uSF, defining the lines of
constant physics, is done with nonvanishing 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 moleculardynamics (MD)
units, the measured integrated autocorrelation times for ZP,
values of the dimensionless bare current quark mass Lm1,
and the boundarytoboundary 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 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)
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 [hepph]
2. S. Heinemeyer et al. [LHC Higgs Cross Section Working Group]. arXiv:1307 .1347 [hepph]
3. L.G. Almeida , S.J. Lee , S. Pokorski , J.D. Wells , Phys. Rev. D 89 ( 3 ), 033006 ( 2014 ). arXiv: 1311 .6721 [hepph]
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 [hepph]
6. S. Aoki et al., Eur. Phys. J. C 74 , 2890 ( 2014 ). arXiv: 1310 .8555 [heplat]
7. S. Aoki et al., Eur. Phys. J. C 77 ( 2 ), 112 ( 2017 ). arXiv: 1607 .00299 [heplat]
8. A. Bazavov et al. [ MILC Collaboration] . PoS CD 09 , 007 ( 2009 ). arXiv: 0910 .2966 [hepph]
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 [heplat]
10. A. Bazavov et al., PoS LATTICE 2010 , 083 ( 2010 ). arXiv: 1011 .1792 [heplat]
11. S. Dürr et al. [BMW Collaboration]. Phys. Lett. B 701 , 265 ( 2011 ). arXiv: 1011 .2403 [heplat]
12. S. Dürr et al. [BMW Collaboration]. JHEP 1108 , 148 ( 2011 ). arXiv: 1011 .2711 [heplat]
13. B. Blossier et al. [ETM Collaboration]. Phys. Rev. D 82 , 114513 ( 2010 ). arXiv: 1010 .3659 [heplat]
14. P. Fritzsch et al. [ALPHA Collaboration]. Nucl. Phys. B 865 , 397 ( 2012 ). arXiv: 1205 .5380 [heplat]
15. N. Carrasco et al. [ETM Collaboration]. JHEP 1403 , 016 ( 2014 ). arXiv:1308 . 1851 [heplat]
16. F. Bernardoni et al. [ALPHA Collaboration]. Phys. Lett. B 730 , 171 ( 2014 ). arXiv: 1311 .5498 [heplat]
17. N. Carrasco et al. [ETM Collaboration]. Nucl. Phys. B 887 , 19 ( 2014 ). arXiv: 1403 .4504 [heplat]
18. C. Alexandrou , V. Drach , K. Jansen , C. Kallidonis , G. Koutsou, Phys. Rev. D 90 ( 7 ), 074501 ( 2014 ). [arXiv: 1406 .4310 [heplat]]
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 [heplat]
20. B. Chakraborty et al. [HPQCD Collaboration], Phys. Rev. D 91 ( 5 ), 054508 ( 2015 ). arXiv: 1408 .4169 [heplat]
21. Y.B. Yang et al. [χQCD Collaboration], Phys. Rev. D 92 ( 3 ), 034517 ( 2015 ). arXiv: 1410 .3343 [heplat]
22. T. Blum et al. [RBC and UKQCD Collaborations], Phys. Rev. D 93 ( 7 ), 074505 ( 2016 ). arXiv: 1411 .7017 [heplat]
23. K.G. Chetyrkin, Phys. Lett. B 404 , 161 ( 1997 ). arXiv:hepph/9703278
24. J.A.M. Vermaseren , S.A. Larin , T. van Ritbergen , Phys. Lett. B 405 , 327 ( 1997 ). arXiv:hepph/9703284
25. J.A. Gracey , Nucl. Phys. B 662 , 247 ( 2003 ). arXiv:hepph/0304113
26. I. Campos et al. [ ALPHA Collaboration] . PoS LATTICE 2015 , 249 ( 2016 ). arXiv: 1508 .06939 [heplat]
27. I. Campos et al. [ ALPHA Collaboration] . EPJ Web Conf . 137 , 08006 ( 2017 ). arXiv: 1611 .06102 [heplat]
28. I. Campos et al. [ ALPHA Collaboration] . PoS LATTICE 2016 , 201 ( 2016 ). arXiv: 1611 .09711 [heplat]
29. M. Dalla Brida et al. [ALPHA Collaboration], Phys. Rev. Lett . 117 ( 18 ), 182001 ( 2016 ). arXiv: 1604 .06193 [hepph]
30. M. Dalla Brida et al. [ALPHA Collaboration], Phys. Rev. D 95 ( 1 ), 014507 ( 2017 ). arXiv: 1607 .06423 [heplat]
31. M. Bruno et al. [ALPHA Collaboration], Phys. Rev. Lett . 119 ( 10 ), 102001 ( 2017 ). arXiv: 1706 .03821 [heplat]
32. M. Lüscher , R. Narayanan , P. Weisz , U. Wolff, Nucl. Phys. B 384 , 168 ( 1992 ). arXiv:heplat/9207009
33. S. Sint, Nucl. Phys. B 421 , 135 ( 1994 ). arXiv:heplat/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:heplat/9810063
35. M. Della , Morte et al. [ALPHA Collaboration]. Nucl. Phys. B 729 , 117 ( 2005 ). arXiv:heplat/0507035
36. M. Guagnelli et al. [ALPHA Collaboration]. JHEP 0603 , 088 ( 2006 ). arXiv:heplat/0505002
37. F. Palombi , C. Pena , S. Sint [ALPHA Collaboration], JHEP 0603 , 089 ( 2006 ). arXiv:heplat/0505003
38. P. Dimopoulos et al., Phys. Lett. B 641 , 118 ( 2006 ). arXiv:heplat/0607028
39. P. Dimopoulos et al. [ALPHA Collaboration]. JHEP 0805 , 065 ( 2008 ). arXiv: 0712 .2429 [heplat]
40. F. Palombi , M. Papinutto , C. Pena , H. Wittig [ALPHA Collaboration], JHEP 0709 , 062 ( 2007 ). arXiv: 0706 .4153 [heplat]
41. B. Blossier , M. Della Morte , N. Garron , R. Sommer [ALPHA Collaboration], JHEP 1006 , 002 ( 2010 ). arXiv: 1001 .4783 [heplat]
42. F. Bernardoni et al. [ALPHA Collaboration]. Phys. Lett. B 735 , 349 ( 2014 ). arXiv: 1404 .3590 [heplat]
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 [heplat]
44. P. Dimopoulos et al. [ ALPHA Collaboration] . arXiv: 1801 .09455 [heplat]
45. C. Pena , D. Preti [ALPHA Collaboration] . arXiv:1706 .06674 [heplat]
46. M. Bruno et al., JHEP 1502 , 043 ( 2015 ). arXiv: 1411 .3982 [heplat]
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 YangMills 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:heplat/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 ). https://doi.org/10.1016/ 0550  3213 ( 95 ) 00352  S . arXiv:heplat/ 9504005
60. P. Fritzsch , A. Ramos , JHEP 1310 , 008 ( 2013 ). arXiv: 1301 .4388 [heplat]
61. A. Ramos, PoS LATTICE 2014 , 017 ( 2015 ). arXiv: 1506 .00118 [heplat]
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:heplat/9911018
63. S. Sint et al. [ALPHA Collaboration]. Nucl. Phys. B 545 , 529 ( 1999 ). arXiv:heplat/9808013
64. U. Wolff, [ALPHA Collaboration], Comput. Phys. Commun . 156 , 143 ( 2004 ). (Erratum: Comput. Phys. Commun . 176 ( 2007 ) 383) . arxiv:heplat/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 CPPACS Collaborations] . Phys. Rev. D 71 , 054505 ( 2005 ). arXiv:heplat/0406028
69. M. Lüscher , P. Weisz , Nucl. Phys. B 479 , 429 ( 1996 ). arXiv:heplat/9606016
70. A. Bode et al. [ALPHA Collaboration], Nucl. Phys. B 540 , 491 ( 1999 ). arXiv:heplat/9809175
71. M. Lüscher , S. Schaefer , Comput. Phys. Commun . 184 , 519 ( 2013 ). arXiv: 1206 .2809 [heplat]
72. M. Lüscher , S. Schaefer, openQCD simulation program for lattice QCD with open boundary conditions . http://luscher.web.cern.ch/ 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 [heplat]
75. S. Takeda , S. Aoki , K. Ide , Phys. Rev. D 68 , 014505 ( 2003 ). arXiv:heplat/0304013
76. Pol Vilaseca, private communication
77. S. Schaefer et al. [ALPHA Collaboration]. Nucl. Phys. B 845 , 93 ( 2011 ). arXiv: 1009 .5228 [heplat]
78. P. Fritzsch , A. Ramos , F. Stollenwerk , PoS Lattice 2013 , 461 ( 2014 ). arXiv: 1311 .7304 [heplat]
79. J. Bulava et al., ALPHA Collaboration. Nucl. Phys. B 896 , 555 ( 2015 ). arXiv: 1502 .04999 [heplat]
80. J. Bulava , M. Della Morte , J. Heitger , C. Wittemeier [ALPHA Collaboration] , Phys. Rev. D 93 ( 11 ), 114513 ( 2016 ). arXiv: 1604 .05827 [heplat]
81. M. Dalla Brida , S. Sint , P. Vilaseca, JHEP 1608 , 102 ( 2016 ). arXiv: 1603 .00046 [heplat]
82. M. Dalla Brida , T. Korzec , S. Sint , P. Vilaseca, in preparation
83. S. Sint, Nucl. Phys. B 847 , 491 ( 2011 ). arXiv: 1008 .4857 [heplat]
84. S. Sint , B. Leder , PoS LATTICE 2010 , 265 ( 2010 ). arXiv: 1012 .2500 [heplat]
85. M. Dalla Brida , S. Sint, PoS LATTICE 2014 , 280 ( 2014 ). arXiv: 1412 .8022 [heplat]
86. D. Mohler , S. Schaefer , J. Simeth, arXiv: 1712 .04884 [heplat]
87. R. Frezzotti et al. [ALPHA Collaboration]. JHEP 0108 , 058 ( 2001 ). arXiv:heplat/0101001
88. [ALPHA Collaboration], to appear
89. M. Bruno , T. Korzec , S. Schaefer, Phys. Rev. D 95 ( 7 ), 074504 ( 2017 ). arXiv: 1608 .08900 [heplat]
90. G. Herdoíza, C. Pena , D. Preti , J.Á. Romero , J. Ugarrio . arXiv: 1711 .06017 [heplat]
91. F. Knechtli et al., ALPHA Collaboration. Phys. Lett. B 774 , 649 ( 2017 ). arXiv: 1706 .04982 [heplat]
92. S. Calì , F. Knechtli , T. Korzec , H. Panagopoulos . arXiv: 1710 .06221 [heplat]
93. F. Knechtli , T. Korzec , B. Leder , G. Moir [ALPHA Collaboration] . arXiv:1710 .07590 [heplat]
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