11institutetext: Dipartimento di Fisica, Università di Milano-Bicocca, Piazza
della Scienza 3, I-20126 Milano, Italy 22institutetext: INFN, Sezione di
Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy 33institutetext:
Department of Physics, University of Cyprus, P.O. Box 20537, 1678 Nicosia,
Cyprus

# Past, present, and future of precision determinations of the QCD coupling
from lattice QCD

Mattia Dalla Brida Affiliation with institutions 1 and 2 ended on
31/10/2020.112233

(Received: date / Revised version: date)

###### Abstract

Non-perturbative scale-dependent renormalization problems are ubiquitous in
lattice QCD as they enter many relevant phenomenological applications. They
require solving non-perturbatively the renormalization group equations for the
QCD parameters and matrix elements of interest in order to relate their non-
perturbative determinations at low energy to their high-energy counterparts
needed for phenomenology. Bridging the large energy separation between the
hadronic and perturbative regimes of QCD, however, is a notoriously difficult
task. In this contribution we focus on the case of the QCD coupling. We
critically address the common challenges that state-of-the-art lattice
determinations have to face in order to be significantly improved. In
addition, we review a novel strategy that has been recently put forward in
order to solve this non-perturbative renormalization problem and discuss its
implications for future precision determinations. The new ideas exploit the
decoupling of heavy quarks to match ${N_{\rm f}}$-flavor QCD and the pure
Yang-Mills theory. Through this matching the computation of the non-
perturbative running of the coupling in QCD can be shifted to the
computationally much easier to solve pure-gauge theory. We shall present
results for the determination of the $\Lambda$-parameter of ${N_{\rm
f}}=3$-flavor QCD where this strategy has been applied and proven successful.
The results demonstrate that these techniques have the potential to unlock
unprecedented precision determinations of the QCD coupling from the lattice.
The ideas are moreover quite general and can be considered to solve other non-
perturbative renormalization problems.

###### pacs:

12.38.Aw,12.38.Bx,12.38.Gc,11.10.Hi,11.10.Jj
Keywords: QCD, Perturbation Theory, Lattice QCD

##  1 Introduction

Renormalization is a fundamental step in order to extract (meaningful)
phenomenologically relevant results from lattice QCD calculations. For the
lattice theorist it is natural to renormalize the bare parameters of the
lattice QCD Lagrangian and the composite operators of interest in terms of
some hadronic renormalization schemes (cf. refs. Sommer:1997xw ;
Sommer:2015kza ). In order to make the determinations accessible to
phenomenologists, however, it is often necessary to translate the results
obtained in the chosen hadronic schemes to results in the (perturbative)
schemes and at the scales commonly considered in phenomenology. In practice,
this requires the determination of the non-perturbative renormalization group
(RG) running of the renormalized QCD parameters and operators in some
convenient intermediate scheme, from the hadronic scales where they were
originally defined, up to some high-energy scale, where perturbation theory
eventually applies and a matching to phenomenological schemes can be
performed.

Over the last decade or so, lattice QCD has entered a precision era for an
increasingly large set of quantities (cf. ref. Aoki:2019cca ). Renormalization
is a relevant part of many of these computations where it can significantly
impact the quality of the final results. Hence, as we are forced to become
more aware of all possible sources of uncertainties in the determination of
the bare lattice quantities, the same care must be reserved to their
renormalization. In particular, as any other lattice calculation, besides the
statistical errors the determinations of renormalized parameters and operators
have their systematics to deal with, i.e. discretization effects, finite-
volume effects, quark-mass effects, and, when a matching to phenomenological
schemes is necessary, also perturbative uncertainties. It is therefore
important that the development in strategies to compute (bare) lattice
quantities is accompanied with new ideas to improve their renormalization, so
to guarantee a precise and robust end result.

An extreme example of this situation, if we can call it this way, is the
determination of the QCD parameters. In this case we can say that the problem
is entirely a renormalization problem, which, however, has very important
phenomenological applications. On the lattice, the QCD coupling and quark
masses are renormalized in terms of hadronic masses and decay constants, while
in phenomenology the QCD parameters are needed at energies of the order of a
hundred ${\rm GeV}$ and above. One would thus think that lattice QCD is not
the right tool for providing this information given the very high energies
involved. It appears more natural indeed to obtain these parameters directly
from high-energy quantities, rather than from the hadronic spectrum. As we
shall recall later in this contribution, this is actually not the case, as
lattice techniques offer an ideal framework for these computations.

For the last 10-15 years, lattice QCD has consistently delivered some of the
most precise determinations for the QCD parameters, as in particular for the
QCD coupling $\alpha_{\rm s}$ (see e.g. refs. Aoki:2019cca ; dEnterria:2019its
; Zyla:2020zbs ).111We here adopt the common notation $\alpha_{\rm
s}(\mu)\equiv\alpha_{\overline{\rm MS}}^{(5)}(\mu)$, where
$\alpha_{\overline{\rm MS}}^{(5)}(\mu)$ is the QCD coupling of the 5-flavor
theory renormalized in the $\overline{\rm MS}$-scheme (see e.g. ref.
Zyla:2020zbs ). Note that for the ease of notation we often omit to write
explicitly the $\mu$-dependence of the coupling. The current world average for
the QCD coupling evaluated for reference at the $Z$-boson pole mass $M_{Z}$ is
$\alpha_{\rm s}(M_{Z})=0.1179(10)$ Zyla:2020zbs , and has a precision of about
0.8%. The lattice determinations alone give $\alpha_{\rm s}(M_{Z})=0.1182(8)$
Aoki:2019cca , and are the most precise subcategory of those considered by the
PDG. Besides the high precision of the individual state-of-the-art
determinations, it is important to emphasize also their overall consistency.
This is a rather non-trivial result considering the fact that even though all
lattice determinations share some common systematics, these are probed quite
differently by considering very different strategies Aoki:2019cca . It is fair
to say that such a variety of approaches within a PDG subcategory is in fact
unique Zyla:2020zbs .

Despite the tremendous efforts on and off the lattice, however, the current
uncertainty on $\alpha_{\rm s}$ is still large. It is one of the largest
sources of uncertainty in several key processes, particularly so within the
Higgs sector, and it is expected to be a limiting factor in many high-
precision studies at future colliders (see e.g. refs. Salam:2017qdl ;
dEnterria:2019its ). An uncertainty on $\alpha_{\rm s}(M_{Z})$ comfortably
below the percent level is desired for precision applications. For these
reasons, there are plans for future phenomenological determinations of
$\alpha_{\rm s}(M_{Z})$ aiming at reaching the extremely competitive accuracy
of 0.2% using high-luminosity high-energy data (see e.g. refs. Gomez-
Ceballos:2013zzn ; Abada:2019lih ; Agostini:2020fmq ; dEnterria:2020cpv ). The
lattice community needs to meet the challenge.

Reducing the current uncertainties on lattice determinations of $\alpha_{\rm
s}$ by such an important factor is not easy. Similarly to several
phenomenological determinations most lattice determinations of $\alpha_{\rm
s}$ are currently limited by systematic uncertainties related to the use of
perturbation theory at relatively low scales Aoki:2019cca . The issue is due
to the fact that reaching high energy on the lattice requires small lattice
spacings to be simulated and this is in general difficult without a dedicated
strategy.

A way around this has been known since a long time and it is based on the
concepts of finite-volume renormalization schemes and finite-size scaling (or
step-scaling) techniques Luscher:1991wu ; Jansen:1995ck . The methods have
been recently applied for obtaining one of the most precise determinations of
$\alpha_{\rm s}$ Bruno:2017gxd . The key feature of the approach is that it
allows for reaching high energy with all systematics under control. This puts
the lattice determinations in the privileged position of being able to reach
in a clean and controlled way high energies fully non-perturbatively. The
systematics due to the application of perturbation theory, in particular, can
be entirely avoided at the expenses of the statistical errors accumulated in
running from low up to high-enough energy. The net advantage of this situation
is that differently from systematic uncertainties, statistical errors can be
straightforwardly reduced. Nonetheless, a reduction of the current
uncertainties on $\alpha_{\rm s}$ by an important factor is yet a
computationally expensive task, even employing a step-scaling strategy (cf.
ref. Bruno:2017gxd ).

In this contribution we want to review the recent proposal made in ref.
DallaBrida:2019mqg  which may allow for such error reduction in a
substantially cheaper way. The key feature of this proposal is that one can
replace the computation of the RG running of the coupling in ${N_{\rm
f}}$-flavor QCD with that in the pure-gauge theory. It is clear that,
regardless of the chosen strategy, this allows for a substantial
simplification of the problem.

In short, the idea is built on three main steps and exploits the decoupling of
heavy quarks in a couple of ways. In the first step, heavy-quark decoupling is
used to connect a low-energy scale $\mu_{\rm dec}$ in ${N_{\rm f}}$-flavor QCD
with the corresponding scale in the pure-gauge theory. This is achieved
through the computation of a massive renormalized coupling in an (unphysical)
theory with ${N_{\rm f}}$ heavy quarks of mass $M\gg\mu_{\rm dec}$. In a
second step, by computing the non-perturbative RG running in the pure Yang-
Mills theory of a convenient coupling one obtains the pure-gauge
$\Lambda$-parameter in units of $\mu_{\rm dec}$, i.e. $\Lambda^{({N_{\rm
f}}=0)}_{\overline{\rm MS}}/\mu_{\rm dec}$. Finally, perturbative decoupling
relations are invoked at a scale $\mu\approx M$ to estimate the ratio of
$\Lambda$-parameters in the ${N_{\rm f}}$-flavor and pure Yang-Mills theory,
that is, $\Lambda^{({N_{\rm f}}=0)}_{\overline{\rm MS}}/\Lambda^{({N_{\rm
f}})}_{\overline{\rm MS}}$. Putting these steps together, one obtains
$\Lambda^{({N_{\rm f}})}_{\overline{\rm MS}}/\mu_{\rm dec}$, and given the
physical value of $\mu_{\rm dec}$ finds $\Lambda^{({N_{\rm
f}})}_{\overline{\rm MS}}$. Considering ${N_{\rm f}}=3$ or 4, once
$\Lambda^{({N_{\rm f}})}_{\overline{\rm MS}}$ is determined one proceeds as
usual and applies perturbative decoupling relations at the charm and/or bottom
quark-mass scale to estimate $\Lambda^{({N_{\rm f}}=5)}_{\overline{\rm MS}}$
and from it $\alpha_{\rm s}(M_{Z})$.

The strategy has already been proven successful in the determination of
$\Lambda^{({N_{\rm f}}=3)}_{\overline{\rm MS}}$ DallaBrida:2019mqg . The ideas
presented in this reference are however general and may be applied to solve
other non-perturbative scale-dependent renormalization problems that face
analogous challenges.

The outline of this contribution is the following.

We begin in Sect. 2 by recalling the main challenges in solving scale-
dependent renormalization problems on the lattice. The emphasis will be on the
determination of the QCD coupling. Besides introducing important concepts for
later sections, the presentation gives us the opportunity to discuss some
recent interesting determinations. These clearly illustrate the difficulties
that state-of-the-art computations of the coupling have to face in order to be
significantly improved.

In Sect. 3, we introduce the theory of heavy-quark decoupling and present the
results of several recent studies that systematically assessed the size of
non-perturbative effects induced by heavy quarks. More precisely, the accuracy
of using perturbative decoupling relations to match the $\Lambda$-parameters
of different ${N_{\rm f}}$-flavor theories is investigated, as well as the
corrections due to the heavy quarks in low-energy quantities. These studies
not only set the foundation for the renormalization strategy based on
decoupling, but also establish the precision at which $\alpha_{\rm s}$ can be
obtained from results in ${N_{\rm f}}=3$ QCD.

In Sect. 4, the application of heavy-quark decoupling to the determination of
the ${N_{\rm f}}$-flavor QCD coupling is described in detail and the results
of ref. DallaBrida:2019mqg  for $\Lambda^{({N_{\rm f}}=3)}_{\overline{\rm
MS}}$ are presented. We conclude in Sect. 5 with some comments on the future
prospects for $\alpha_{\rm s}$ determinations in view of this new strategy.

We care to note that it is not the aim of the present contribution to discuss
in detail the many different lattice approaches that are currently considered
to determine the QCD coupling. In particular, we do not provide a complete
account of all recent determinations. For such a discussion, we refer the
interested reader to the comprehensive work of FLAG Aoki:2019cca  and to other
interesting recent reviews (see e.g. refs. DelDebbio:2021ryq ;
Komijani:2020kst ).

##  2 Precision determinations: the case of $\alpha_{\rm s}$

Before presenting the renormalization ideas based on decoupling, we believe it
is important to put these into context. The aim of these strategies, in fact,
is not simply that of providing alternative ways to solve non-perturbative
scale-dependent renormalization problems. The goal is to develop a framework
that will allow us to improve significantly our control over the current most
relevant uncertainties. In this section, we thus want to recall what the main
challenges are in solving this class of problems and which are the common
approaches that are used to tackle them. Many of the concepts and observations
that will be presented are for the most part well known. However, these issues
are now more current than ever given the high precision that lattice QCD
calculations have achieved, in particular in the determination of the QCD
parameters. For this reason, we think it is important to address them also
here. This gives us the opportunity to discuss some new insight that has been
gathered from several recent high-precision studies, as well as introducing
key concepts for later sections.

As anticipated, the discussion will focus on the case of the QCD coupling
$\alpha_{\rm s}$. This allow us to analyze in easier terms the main challenges
that we need to face in high-precision non-perturbative determinations of RG
runnings while capturing all the relevant issues. Moreover, lattice
determinations of the strong coupling are a distinct case of competitive
calculations which have the potential to deliver unprecedentedly precise
results for a very relevant and fundamental quantity. Making a significant
progress over the present state-of-the-art determinations by mere brute force,
however, is extremely demanding from the computational point of view. It is
therefore mandatory to develop new strategies with the clear scope of
improving our control on all sources of uncertainty.

###  2.1 Determinations of $\alpha_{\rm s}$ on and off the lattice

As already mentioned, since more than a decade lattice QCD is providing the
Particle Physics community with the most accurate determinations of
$\alpha_{\rm s}$ (see refs. Aoki:2019cca ; dEnterria:2019its ; Zyla:2020zbs ).
The reason behind this is that, as we shall recall, lattice determinations
have some important advantages over their phenomenological counterparts (see
e.g. refs. Sommer:2015kza ; Salam:2017qdl ; DelDebbio:2021ryq  for some
reviews).

Any determination of $\alpha_{\rm s}$, whether on the lattice or not, relies
on the following basic strategy. One considers a short-distance observable
${\mathcal{O}}(q)$ which depends on a characteristic energy scale $q$. In the
limit where $q\to\infty$, this observable is compared with its theoretical
prediction, ${\mathcal{O}}_{\rm th}(q)$, in terms of a perturbative
expansion222Note that in general the coupling $\alpha_{\rm s}$ to be
considered here should be the QCD coupling of the relevant ${N_{\rm
f}}$-flavor theory, i.e. $\alpha_{\rm s}(\mu)\equiv\alpha_{\overline{\rm
MS}}^{({N_{\rm f}})}(\mu)$, from which $\alpha_{\overline{\rm
MS}}^{(5)}(M_{Z})$ can eventually be extracted (cf. e.g. ref. Aoki:2019cca
and Sect. 3.1.2). At this stage, however, we prefer to keep the discussion
simple. We shall return to this point later in detail.

${\mathcal{O}}_{\rm th}(q)=\sum_{n=0}^{N}k_{n}(s)\alpha_{\rm s}^{n}(\mu)+{\rm
O}(\alpha_{\rm s}^{N+1})+{\rm O}\bigg{(}{\Lambda^{p}\over q^{p}}\bigg{)},\
\mu={q\over s}\,.$ (1)

The functions $k_{n}(s)$ appearing in this equation are the coefficient
functions defining the perturbative series. They are known up to some order
$N$ and depend on the scale factor, $s>0$, that relates the renormalization
scale $\mu$ at which $\alpha_{\rm s}$ is extracted with the scale $q$. The
basic difference between phenomenological and lattice determinations of
$\alpha_{\rm s}$ is the choice of the observable ${\mathcal{O}}(q)$.

Requiring ${\mathcal{O}}_{\rm th}(q)={\mathcal{O}}(q)$ for some finite $q$,
clearly fixes the value of $\alpha_{\rm s}(\mu)$ only up to some error. This
error comes from several different sources, many of which are common to all
types of determinations. First of all, there is the precision
$\delta{\mathcal{O}}(q)$ to which the observable ${\mathcal{O}}(q)$ is known.
This of course depends on the relevant statistical and systematic
uncertainties associated with the determination of ${\mathcal{O}}(q)$.
Secondly, there is the effect of the truncation of the perturbative series to
a given order, i.e. the size of the ${\rm O}(\alpha_{\rm s}^{N+1})$ terms in
eq. (1). In addition to these there are contaminations from “non-perturbative
contributions”. These are represented in eq. (1) by power corrections to the
perturbative expansion of ${\rm O}(\Lambda^{p}/q^{p})$, where $p>0$ and
$\Lambda$ is some characteristic non-perturbative scale of QCD.333Our
knowledge of the form of non-perturbative effects is in fact rather limited.
In addition, strictly speaking, perturbative and non-perturbative
contributions cannot really be separated due to the asymptotic nature of the
series (see e.g. ref. Martinelli:1996pk ). However, as the discussion is at
this point qualitative we simply adopt the simplistic representation of non-
perturbative effects as power corrections to the perturbative expansion. Thus,
regardless of the chosen strategy, an accurate determination of $\alpha_{\rm
s}$ needs to have, _at least_ , these general sources of error under control.
Note that for the most part these are systematic in nature.

Lattice determinations of $\alpha_{\rm s}$ are in principle favored in several
ways in succeeding at this task. Firstly, on the lattice the QCD parameters
are first renormalized in terms of some precisely measured hadronic quantities
(e.g. hadron masses, decay constants, etc.), for which experimental
uncertainties typically contribute only marginally to the end result. Once
these are fixed, one has lots of freedom in choosing an observable
${\mathcal{O}}(q)$ as the getaway to extract $\alpha_{\rm s}$. One can
therefore _devise_ convenient observables which have small statistical and
systematic uncertainties; in particular there is no need for these quantities
to be accessible in experiments. Phenomenological determinations of
$\alpha_{\rm s}$, instead, rely on experimental data for the observable
${\mathcal{O}}(q)$. It is the typical situation that when $q$ becomes large,
and therefore many sources of systematic uncertainty in eq. (1) become small,
the experimental errors $\delta{\mathcal{O}}(q)$ become large. It is thus
difficult to find in general a single experimental quantity ${\mathcal{O}}(q)$
that allows one to accurately follow its scale dependence over a wide range of
$q$-values. On the lattice, on the other hand, if _carefully_ chosen,
${\mathcal{O}}(q)$ can be computed precisely from low- up to very high-energy
scales. This gives a unique handle on controlling non-perturbative corrections
and the contribution of the missing perturbative orders in eq. (1).

Another advantage for the lattice theorist is that ${\mathcal{O}}$ is defined
within QCD alone. Consequently, the theoretical description
${\mathcal{O}}_{\rm th}$ of eq. (1) does not need to include contributions
besides those from QCD. In addition, no modeling of hadronization is needed
when comparing the observable ${\mathcal{O}}$ with its perturbative prediction
${\mathcal{O}}_{\rm th}$. Different is again the situation for
phenomenological determinations. In these cases, other Standard Model (SM)
contributions may be needed in order to extract $\alpha_{\rm s}$ and some
modeling of hadronization is necessary. Depending on the process, these are
known only up to some accuracy and typically depend on the value of other SM
parameters as well. The precision one can aim for $\alpha_{\rm s}$ can
therefore be limited by these factors.444Experimental observables are in
principle sensitive also to any New Physics. How this affects the extracted
value of $\alpha_{\rm s}$, however, is something hard to assess. Of course,
lattice QCD determinations are not entirely exempted from this issue. In this
case, however, the problem of disentangling the QCD contributions from “the
rest” is confined to the hadronic quantities entering the renormalization of
the theory, rather than to the observable ${\mathcal{O}}$ itself. As mentioned
earlier, the uncertainties on the hadronic quantities have, at present,
limited impact on the results for $\alpha_{\rm s}$.

All current lattice QCD determinations of $\alpha_{\rm s}$, on the other hand,
have to deal with the fact that their calculations are performed with an
unphysical number of quark flavors. The bottom and typically also the charm
quark are in fact not included in the simulations. This brings up the issue of
having to account for their missing contributions. We shall leave this very
important discussion aside for the moment and come back to it in detail in the
following (see Sect. 3).

###  2.2 The challenge of reaching high energy

Having for the most part presented the pros that lattice determinations of
$\alpha_{\rm s}$ in principle have, we now address what the main difficulties
are in practice.

As any lattice QCD observable, besides the statistical uncertainties,
${\mathcal{O}}(q)$ is affected by several systematics that need to be
controlled. These include general ones, i.e. discretization errors, finite-
volume effects, an unphysical number of quarks, and quark-mass effects, as
well as others which depend on the specific choice of ${\mathcal{O}}$ and set-
up that we make (e.g. excited-state contaminations, finite-temperature
effects, Gribov copies, topology freezing, etc.). Finite quark-mass effects
are typically not a relevant issue in determinations of the QCD coupling
Aoki:2019cca . As anticipated, we then leave the problem of having an
unphysical quark-content for later. We also ignore observable specific issues.
Here we focus instead on discretization and finite-volume effects. The
combination of having the two under control, in fact, can severely restrict
the accessible range of $q$-values, if the renormalization strategy is not
carefully chosen.

In particular, if one is determined in resolving within the same lattice
simulation both the hadronic energy scales relevant for the renormalization of
the lattice theory, and the energy scale $q$ at which $\alpha_{\rm s}$ is
extracted, then one is necessarily limited by the simultaneous constraints:

$L^{-1}\ll\Lambda\quad\textrm{and}\quad\Lambda\ll q\ll a^{-1}\,.$ (2)

The first inequality expresses the fact that finite-volume effects must be
under control in hadronic quantities. The infrared cutoff set by the finite
extent $L$ of the lattice must be smaller than the typical non-perturbative
scales of QCD, denoted here by $\Lambda$. The second inequality, instead,
encodes two separate requirements. On the one hand, the scale $q$ must be much
lower than the ultraviolet cutoff set by the lattice spacing $a$. In this way,
discretization errors in ${\mathcal{O}}(q)$ are kept under control, and
${\mathcal{O}}(q)$ can be obtained in the continuum limit with controlled
errors.555Here and in the following we assume that ${\mathcal{O}}(q)$ is a
properly renormalized lattice quantity with a well-defined continuum limit and
free from infrared divergences. Alternative strategies extract $\alpha_{\rm
s}$ from bare lattice quantities by taking $q\propto a^{-1}$ and expressing
their expansion in lattice perturbation theory in terms of $\alpha_{\rm s}$
(see e.g. ref. Mason:2005zx ). In these cases, discretization errors and scale
dependence of ${\mathcal{O}}(q)$ are entangled. Addressing systematic
uncertainties becomes more subtle and requires a separate discussion (see ref.
Aoki:2019cca ). On the other hand, $q$ needs to be much larger than the scales
$\Lambda$. Only in this situation perturbation theory can reliably be applied
to extract $\alpha_{\rm s}$ by comparing ${\mathcal{O}}(q)$ with eq. (1).

The typical lattices that are simulated today have sizes $L/a\lesssim 100$.
Taking for definiteness $m_{\pi}L=4$, with $m_{\pi}=140\,{\rm MeV}$, the first
condition in eq. (2) implies that for such ensembles $q\ll a^{-1}\approx
3.5\,{\rm GeV}$. Of course this is a crude estimate and somewhat higher
energies might be achieved by compromising at different stages of the
calculation (e.g. considering heavier pion masses or smaller volumes in order
to reach smaller values of $a$, or taking $aq\lesssim 1$). Nonetheless, it is
clear that although convenient in practice, considering lattices devised for
studying hadronic physics to compute short-distance quantities necessarily
poses severe challenges on the feasibility of the approach, as the accessible
energies $q$ are quite limited.

As a concrete example of this situation, we can point to the recent
determination of $\alpha_{\rm s}$ of ref. Cali:2020hrj . For their computation
the authors employ state-of-the-art large-volume simulations by the CLS
initiative Bruno:2014jqa ; Bali:2016umi . The two-point functions of both
axial and vector currents at short-distances are used to extract $\alpha_{\rm
s}$. Given the range of lattice spacings at their disposal, $a\approx
0.04-0.08\,{\rm fm}$, the accessible energies for which continuum limit
extrapolations could be performed in a controlled way are limited to
$q\lesssim 2\,{\rm GeV}$.

As we shall see below with explicit examples, a limited range of low
$q$-values unavoidably implies a limited attainable precision for $\alpha_{\rm
s}$ as perturbative truncation errors become a major issue.

####  2.2.1 The finite volume is your friend

In order to reach high precision we must tailor the lattice simulations to the
problem at hand. Going back to eq. (2), there is in fact no reason to try to
satisfy simultaneously the two conditions as these belong to separate
problems. On the one hand, there is the determination of the low-energy
quantities used for the hadronic renormalization of the lattice theory, while,
on the other hand, there is the determination of ${\mathcal{O}}(q)$ for large
$q$.

A more natural strategy is therefore to split the problem over several sets of
lattice simulations, each one covering a different range of energy scales. In
this way, systematic effects can be more easily kept under control, as the
relevant conditions will be milder for each individual set of simulations. The
way to effectively achieve this in practice is to employ what are known as
_finite-volume renormalization schemes_ Luscher:1991wu . In this case, the
scale $q$ at which the observable ${\mathcal{O}}(q)$ is evaluated is
identified with the inverse linear extent of the finite volume, i.e.
$q=L^{-1}$. One may say that, in fact, the observable ${\mathcal{O}}$
considered is a _finite-volume effect_. With this choice, one computes the
non-perturbative RG running of ${\mathcal{O}}(L^{-1})$ by simulating lattices
with different physical extent $L$. This strategy goes under the name of
_finite-size (or step-)scaling_ Luscher:1991wu ; Jansen:1995ck  (see ref.
Sommer:2015kza  for a recent account).

More precisely, having fixed the bare QCD parameters through some hadronic
quantities, one computes ${\mathcal{O}}(L^{-1}_{{\rm had}})$ at a low-energy
scale $q_{\rm had}\equiv L^{-1}_{{\rm had}}\approx\Lambda$, and determines
$L_{{\rm had}}$ in physical units. This is achieved by computing $\lim_{a\to
0}(am_{\rm had})(L_{\rm had}/a)={\rm O}(1)$, where $m_{\rm had}$ is a known
low-energy scale. No large scale separations are involved in this step, and at
common bare parameters one can satisfy the conditions $L_{\rm had}/a\gg 1$ and
$am_{\rm had}\ll 1$, as well as having finite-volume effects in $m_{\rm had}$
under control.

Secondly, one computes in the continuum limit the change in
${\mathcal{O}}(L^{-1})$ as $L$ is varied by a known factor, say, $L\to L/2$.
This step is repeated a number of times $n$, going from each new $L$ to the
next one. Once the energy scale reached, $q_{n}=2^{n}/L_{\rm had}$, is large
compared to the hadronic scales, perturbation theory can safely be applied to
extract $\alpha_{\rm s}(\mu_{\rm PT}=q_{n}/s)$ from the value of
${\mathcal{O}}(q_{n})$ (cf. eq. (1)).

It is important to emphasize that, if carefully chosen, the only source of
systematic errors that affect the determination of ${\mathcal{O}}(L^{-1})$ are
discretization effects. In particular, no matter what the scale $q=L^{-1}$ is,
discretization effects are under control once $L^{-1}=q\ll a^{-1}$, i.e.
$L/a\gg 1$. This approach elegantly exploits the freedom that we have in
lattice QCD in choosing the observable ${\mathcal{O}}(q)$ to completely
circumvent the issue of necessarily having a finite volume. In particular,
within this strategy the computational power is entirely invested into
controlling a single systematic uncertainty.

In principle there is quite some freedom in choosing the finite-volume
observable ${\mathcal{O}}(q)$. For the strategy to be successful in practice,
however, this should have a number of desirable properties (see e.g. ref.
Sommer:2015kza ). First of all, it should be easily and precisely measurable
in Monte Carlo simulations. It should be computable in perturbation theory to
a sufficiently high-loop order in order to guarantee good precision in
extracting $\alpha_{\rm s}$ through eq. (1). It should preferably be gauge
invariant, in order to avoid issues with Gribov copies once studied non-
perturbatively, and also be directly measurable for zero quark masses (see
below). Finally, it should have, in general, small lattice artifacts. In fact,
it is not straightforward to find a single observable that has all these nice
features for any range of $q$-values. This, however, is not a real issue as
different observables may be accurately combined in order to cover all the
relevant range of scales $q$. We shall see explicit examples of good
complementary observables in forthcoming sections.

####  2.2.2 $\Lambda$-parameters, $\beta$-functions, and all that

As we leave the general discussion for entering a more quantitative analysis
of the challenges of extracting $\alpha_{\rm s}$, we find convenient to
reformulate the problem in slightly different terms. First of all, it is
useful to associate with the observable ${\mathcal{O}}$ a renormalized
coupling $\bar{g}_{\mathcal{O}}$ through the relation:

$\alpha_{\mathcal{O}}(\mu)\equiv{\bar{g}^{2}_{\mathcal{O}}(\mu)\over
4\pi}\equiv{{\mathcal{O}}(\mu)-k_{0}\over
k_{1}}\overset{\mu\to\infty}{=}\alpha_{\rm s}(\mu)+{k_{2}\over
k_{1}}\alpha_{\rm s}^{2}(\mu)+\ldots\,,$ (3)

where the coefficients $k_{i}\equiv k_{i}(1)$ are those appearing in the
perturbative expansion, eq. (1). This simple procedure defines a non-
perturbative, regularization-independent, QCD coupling. In terms of these
couplings, the extraction of $\alpha_{\rm s}$ is interpreted as the
perturbative matching between $\alpha_{{\mathcal{O}}}$ and $\alpha_{\rm s}$,
where different observables ${\mathcal{O}}$ define different renormalization
schemes. The common normalization allows us to compare the value of the
couplings in different schemes as they approach the high-energy limit. This is
useful in assessing the regime of applicability of the perturbative matching
as the latter can be characterized by the value that $\alpha_{{\mathcal{O}}}$
should reach.

We recall at this point that the non-perturbative couplings studied within
lattice QCD are implicitly defined for a given number of quark flavors,
${N_{\rm f}}$. A more proper notation to use for the couplings is therefore
$\alpha_{{\mathcal{O}}}^{({N_{\rm f}})}(\mu)\equiv[\bar{g}^{({N_{\rm
f}})}_{\mathcal{O}}(\mu)]^{2}/(4\pi)$, which emphasizes the fact that
$\bar{g}^{({N_{\rm f}})}_{\mathcal{O}}(\mu)$ must be considered as a coupling
within the ${N_{\rm f}}$-flavor theory. For ease of notation, however, in this
section we will often take the liberty of dropping the superscripts ${N_{\rm
f}}$ and leave these understood, unless they are needed to avoid any confusion
or for later reference.

Having this noticed, a particularly convenient class of schemes to consider
are _mass-independent_ (or simply massless) renormalization schemes
Weinberg:1951ss . These are defined in terms of observables ${\mathcal{O}}$
evaluated for vanishing quark masses. As a result, the RG running of these
couplings is decoupled from that of the renormalized quark masses and
therefore simpler to solve. On the other hand, differently from the physical
case of massive schemes, quarks do not decouple in the RG running of massless
schemes Bernreuther:1981sg ; Bernreuther:1983zp . Hence, in the ${N_{\rm
f}}$-flavor theory the latter is characterized by a fixed number of active
flavors corresponding to ${N_{\rm f}}$.666We shall return on these points in
more detail in Sect. 3.

To the coupling $\bar{g}^{({N_{\rm f}})}_{\mathcal{O}}$ in a given mass-
independent scheme we can associate a quark-mass independent
$\Lambda$-parameter, $\Lambda^{({N_{\rm f}})}_{\mathcal{O}}$, defined as,

$\begin{split}\Lambda^{({N_{\rm f}})}_{\mathcal{O}}&=\mu\,\varphi^{({N_{\rm
f}})}_{\rm g,{\mathcal{O}}}(\bar{g}^{({N_{\rm
f}})}_{\mathcal{O}}(\mu))\,,\\\\[4.30554pt] \varphi^{({N_{\rm f}})}_{\rm
g,{\mathcal{O}}}(\bar{g})&=(b_{0}({N_{\rm f}})\bar{g}^{2})^{-{b_{1}({N_{\rm
f}})\over 2b_{0}({N_{\rm f}})^{2}}}\,e^{-{1\over 2b_{0}({N_{\rm
f}})\bar{g}^{2}}}\times\\\\[4.30554pt]
&\times\exp\bigg{\\{}-\int_{0}^{\bar{g}}{\rm d}g\bigg{[}{1\over\beta^{({N_{\rm
f}})}_{\mathcal{O}}(g)}+{1\over b_{0}({N_{\rm f}})g^{3}}-{b_{1}({N_{\rm
f}})\over b_{0}({N_{\rm f}})^{2}g}\bigg{]}\bigg{\\}}\,.\end{split}$ (4)

Through this relation, the value of the coupling $\bar{g}^{({N_{\rm
f}})}_{\mathcal{O}}(\mu)$ at any renormalization scale $\mu$ is in one-to-one
correspondence with $\Lambda^{({N_{\rm f}})}_{\mathcal{O}}$, provided that the
$\beta$-function,

$\beta^{({N_{\rm f}})}_{\mathcal{O}}(\bar{g})\equiv\mu{{\rm
d}\bar{g}^{({N_{\rm f}})}_{\mathcal{O}}(\mu)\over{\rm
d}\mu}\bigg{|}_{\bar{g}}\,,$ (5)

is known. The $\beta$-function describes the dependence of the coupling on the
renormalization scale. In perturbation theory, it has an expansion which at
the $N$-loop order reads:

$\beta_{\mathcal{O}}^{({N_{\rm f}})\,{\rm
PT}}(\bar{g})\equiv-\bar{g}^{3}\sum_{k=0}^{N-1}b_{k}({N_{\rm
f}})\,\bar{g}^{2k}\,,$ (6)

with

$\displaystyle b_{0}({N_{\rm f}})={1\over(4\pi)^{2}}\bigg{(}11-{2\over
3}{N_{\rm f}}\bigg{)}\,,$ $\displaystyle b_{1}({N_{\rm
f}})={1\over(4\pi)^{4}}\bigg{(}102-{38\over 3}{N_{\rm f}}\bigg{)}\,.$ (7)

The coefficients $b_{0}({N_{\rm f}}),b_{1}({N_{\rm f}})$ are universal and
shared by all mass-independent renormalization schemes. The scheme dependence
only enters through the higher-order coefficients, $b_{i}({N_{\rm f}})\equiv
b^{\mathcal{O}}_{i}({N_{\rm f}})$, with $i\geq 2$.

A first compelling property of $\Lambda$-parameters is that, differently from
the case of couplings, their scheme dependence is in fact trivial. Leaving the
${N_{\rm f}}$-dependence implicit and taking the $\Lambda$-parameter in the
$\overline{\rm MS}$-scheme, $\Lambda_{\overline{\rm MS}}$, as reference, we
have that any other $\Lambda$-parameter $\Lambda_{{\mathcal{O}}}$ is _exactly_
related to $\Lambda_{\overline{\rm MS}}$ by the relation:

${\Lambda_{\overline{\rm
MS}}\over\Lambda_{{\mathcal{O}}}}=\exp\bigg{\\{}{c_{1}(1)\over
2b_{0}}\bigg{\\}}\,.$ (8)

In this equation $c_{1}(1)$ is the 1-loop coefficient of the perturbative
matching relation between the corresponding couplings, which at $M$-loop order
reads,

$\bar{g}^{2}_{\overline{\rm
MS}}(s\mu)=\bar{g}^{2}_{{\mathcal{O}}}(\mu)\Bigg{(}1+\sum_{k=1}^{M}c_{k}(s)\bar{g}^{2k}_{{\mathcal{O}}}(\mu)\Bigg{)}\,,\quad
s>0\,.$ (9)

It is clear that the $\Lambda$-parameters are non-perturbatively defined once
the corresponding couplings and $\beta$-functions are.777Interestingly, eq.
(8) provides an indirect non-perturbative definition of
$\Lambda_{\overline{\rm MS}}$ through the $\Lambda$-parameter of any non-
perturbative scheme. In addition, they are _exact_ solutions of the Callan-
Symanzik equations Callan:1970yg ; Symanzik:1970rt ; Symanzik:1971vw , and
therefore RG invariants, i.e. ${\rm d}\Lambda_{\mathcal{O}}/{\rm d}\mu=0$.
From a non-perturbative perspective this makes the $\Lambda$-parameters
natural quantities to compute, as they provide reference scales for both the
low- and high-energy regimes of QCD. We finally stress that, as their
corresponding couplings, the $\Lambda$-parameters are defined for a given
number of quark flavors, ${N_{\rm f}}$ (cf. eq. (4)). Each ${N_{\rm
f}}$-flavor theory therefore has its own $\Lambda$-parameters.

####  2.2.3 Systematic uncertainties in extracting $\Lambda^{({N_{\rm
f}})}_{\overline{\rm MS}}$

Once the overall energy scale of the given ${N_{\rm f}}$-flavor theory has
been set,888The issue of fixing the scale of a theory with an unphysical value
of ${N_{\rm f}}$ will be addressed in Sect. 3. the non-perturbative value of
the coupling, $\bar{g}_{\rm PT}\equiv\bar{g}_{\mathcal{O}}(\mu_{\rm PT})$, in
any scheme, at some high-energy scale $\mu_{\rm PT}$, allows for estimating
$\Lambda_{\overline{\rm MS}}^{({N_{\rm f}})}$. Below we present two commonly
employed strategies.

##### Strategy 1.

Using eq. (4) we first get an _asymptotic_ estimate for
$\Lambda_{{\mathcal{O}}}/\mu_{\rm{PT}}$ as,

${\Lambda_{{\mathcal{O}}}\over\mu_{\rm PT}}=\varphi_{\rm
g,{\mathcal{O}}}(\bar{g}_{\rm PT})\overset{\bar{g}_{\rm PT}\to
0}{\approx}\varphi^{\rm PT}_{\rm g,{\mathcal{O}}}(\bar{g}_{\rm PT})+{\rm
O}\big{(}\bar{g}_{\rm PT}^{2N-2}\big{)}\,.$ (10)

In this relation, $\varphi^{\rm PT}_{\rm g,{\mathcal{O}}}$ is defined
analogously to $\varphi_{\rm g,{\mathcal{O}}}$ of eq. (4) but with the
replacement $\beta_{\mathcal{O}}\to\beta_{\mathcal{O}}^{\rm PT}$, with
$\beta_{\mathcal{O}}^{\rm PT}$ the perturbative $\beta$-function of eq. (6).
From the estimate of ${\Lambda_{{\mathcal{O}}}/\mu_{\rm PT}}$,
${\Lambda_{\overline{\rm MS}}/\mu_{\rm PT}}$ is obtained, with no further
approximation, using eq. (8). Finally, the knowledge of $\mu_{\rm{PT}}$ in
physical units gives us $\Lambda_{\overline{\rm MS}}$.

As anticipated by eq. (10), due to the truncation of $\beta_{\mathcal{O}}$ to
$N$-loop order in the evaluation of ${\Lambda_{{\mathcal{O}}}/\mu_{\rm PT}}$,
our estimate for ${\Lambda_{\overline{\rm MS}}/\mu_{\rm PT}}$, comes with a
systematic uncertainty of ${\rm O}\big{(}\bar{g}_{\rm PT}^{2N-2}\big{)}$. It
is important to stress that eq. (10) is in fact only asymptotic. Similarly to
what we discussed about eq. (1), our estimate for ${\Lambda_{\overline{\rm
MS}}/\mu_{\rm PT}}$ is in principle also affected by “non-perturbative
corrections”, if the coupling $\bar{g}_{\rm PT}$ is not small enough for
perturbation theory to be in the regime of applicability. We shall come back
to this issue shortly.

##### Strategy 2.

A second possibility to estimate $\Lambda_{\overline{\rm MS}}$ is to first
obtain from the non-perturbative value of $\bar{g}_{\rm PT}$ an estimate for
$\bar{g}_{\overline{\rm MS}}(s\mu_{\rm{PT}})$ using the $M$-loop relation, eq.
(9). Given this, we can estimate ${\Lambda_{\overline{\rm MS}}/(s\mu_{\rm
PT}})$ using eq. (4) and the perturbative expression for the $\beta$-function,
eq. (6), in the $\overline{\rm MS}$-scheme. As before, the knowledge of
$\mu_{\rm{PT}}$ in physical units then gives us $\Lambda_{\overline{\rm MS}}$.

In order to establish the perturbative uncertainties associated with this
second approach, we first recall that the $\beta$-function in the
$\overline{\rm MS}$-scheme is currently known to 5-loop order Baikov:2016tgj ;
Luthe:2016ima ; Herzog:2017ohr . This introduces a systematic uncertainty in
the determination of ${\Lambda_{\overline{\rm MS}}/\mu_{\rm PT}}$ of ${\rm
O}\big{(}\bar{g}^{8}_{\overline{\rm MS}}(s\mu_{\rm{PT}})\big{)}\approx{\rm
O}\big{(}\bar{g}_{\rm PT}^{8}\big{)}$ (cf. eqs. (10) and (9)). Secondly, using
eq. (4) it is easy to show that the perturbative matching at $M$-loop order
between $\bar{g}_{\rm PT}$ and $\bar{g}_{\overline{\rm MS}}(s\mu_{\rm{PT}})$
translates into a systematic uncertainty in ${\Lambda_{\overline{\rm
MS}}/\mu_{\rm PT}}$ of ${\rm O}\big{(}\bar{g}_{\rm{PT}}^{2M}\big{)}$.

For all schemes $\bar{g}_{{\mathcal{O}}}$ used in lattice QCD determinations
of $\Lambda_{\overline{\rm MS}}$ we have that $M\leq 3$ (see e.g. ref.
Aoki:2019cca ). The systematic errors coming from the matching between
couplings is therefore parametrically larger than the one from the truncation
of $\beta_{\overline{\rm MS}}$. Moreover, note that for all these schemes we
have that $M=N-1$, where $N$ is the loop order at which the corresponding
perturbative $\beta$-function, $\beta^{\rm PT}_{\mathcal{O}}$, is known (cf.
eq. (6)).999This is the case because for all these schemes the $N$-loop
$\beta$-function, $\beta^{\rm PT}_{\mathcal{O}}$, has been inferred from
$\beta^{\rm PT}_{\overline{\rm MS}}$ at $N$-loops using the matching relation
eq. (9) and $\beta^{\rm PT}_{\mathcal{O}}$ at $(N-1)$-loop order. This means
that, for the schemes $\bar{g}_{\mathcal{O}}$ commonly used, Strategy 1. and
2. result in the _same_ parametric uncertainties of ${\rm
O}\big{(}\bar{g}_{\rm{PT}}^{2N-2}\big{)}$. Clearly, although parametrically
the same, the actual size of the corrections might be different. In this
second strategy, in particular, when matching the couplings we have the
freedom to choose the parameter $s$ (cf. eq. (9)). Different choices can
result in different perturbative corrections to ${\Lambda_{\overline{\rm
MS}}/\mu_{\rm PT}}$.

Devising different strategies like the ones above and comparing their outcome
can help us assessing the systematic uncertainties in $\Lambda_{\overline{\rm
MS}}$ coming from the use of perturbation theory at $\mu_{\rm{PT}}$. A truly
systematic study, however, requires to compare the determination of
$\Lambda_{\overline{\rm MS}}/\mu_{\rm ref}$, where $\mu_{\rm ref}$ is a common
reference scale, for several different values of $\bar{g}_{\rm PT}$ as
$\bar{g}_{\rm PT}\to 0$. Only if agreement is found among all determinations,
possibly including different strategies, one may be reassured that ${\rm
O}\big{(}\bar{g}_{\rm PT}^{2N-2}\big{)}$ terms, as well as non-analytic terms
in the coupling, are negligible within the statistical uncertainties. In the
case where, instead, the results for $\Lambda_{\overline{\rm MS}}/\mu_{\rm
ref}$ show a clear dependence on $\bar{g}_{\rm PT}$, one should first confirm
that this is actually compatible with the expected ${\rm O}\big{(}\bar{g}_{\rm
PT}^{2N-2}\big{)}$ corrections. If this is the case, one may be confident that
the asymptotic regime of the perturbative expansion is reached, no “non-
perturbative corrections” are relevant, and one can therefore take as final
estimate for $\Lambda_{\overline{\rm MS}}/\mu_{\rm ref}$ the extrapolated
result for $\bar{g}_{\rm PT}\to 0$.

Clearly, the program above is ambitious. The running of the coupling
$\bar{g}_{\rm PT}$ at high energies is only logarithmic in $\mu_{\rm
PT}/\Lambda_{\mathcal{O}}$. Reducing the size of the perturbative truncation
errors by a given factor hence requires an _exponentially_ larger change in
the energy scale $\mu_{\rm{PT}}$. In order to accurately estimate the
systematic uncertainties coming from the use of perturbation theory one
therefore needs to cover, non-perturbatively, a wide range of energies,
reaching up to very high scales.

If the chosen strategy to determine $\Lambda_{\overline{\rm MS}}$ does not
allow for this and the accessible range of $\bar{g}_{\rm PT}$ is quite
limited, one might be tempted to estimate the uncertainties due to the
application of perturbation theory in more simplistic ways. For example, one
might opt for simply adding to the final result an uncertainty
$\delta\Lambda_{\overline{\rm MS}}/\Lambda_{\overline{\rm MS}}=k\,\bar{g}_{\rm
PT}^{2N-2}$, where $k$ is estimated in some way from the available
perturbative information. Alternatively, one might estimate these
uncertainties based on the spread of the results obtained at the smallest
available $\bar{g}_{\rm{PT}}$ from different strategies (e.g. Strategy 1. vs
Strategy 2.). Given the asymptotic nature of the perturbative expansion,
however, these practices cannot be considered reliable in general. From the
very definition of asymptotic series the only reliable way to assess its
accuracy is to compare the series with the full function as $\bar{g}_{\rm
PT}\to 0$.101010We recall that a series is said to be asymptotic to the
function $f(\lambda)$, $\lambda\in\mathbb{R}$, if:
$|f(\lambda)-\sum_{j=0}^{N}a_{j}\lambda^{j}\big{|}\overset{\lambda\to
0}{\rightarrow}{\rm O}(\lambda^{N+1})$, $\forall N$. Note in particular that
at fixed $\lambda$, larger $N$ does not necessarily imply a better
approximation of the series to the function. In order to do so, the coupling
$\bar{g}_{\rm PT}$ must be varied by a sensible amount reaching down to small
values.

For the same reasons, it is not advisable to estimate the size of “non-
perturbative corrections” using some model assumption, or use some model to
extrapolate the results for $\Lambda_{\overline{\rm MS}}/\mu_{\rm ref}$ to
$\bar{g}_{\rm PT}\to 0$. Our knowledge of the form of non-perturbative effects
is rather limited and the separation between what is perturbative and non-
perturbative is all but well defined. Hence, it is always debatable whether
any model that tries to capture non-analytic terms in the coupling is really
adequate to describe the data within the given accuracy. Moreover, if the
coupling $\bar{g}_{\rm PT}$ cannot be varied much, it is difficult to really
distinguish, e.g. a power correction, from some higher-order term in
$\bar{g}_{\rm PT}$, when statistical errors and other uncertainties are
present. A more reliable practice is thus to avoid regions of large
$\bar{g}_{\rm PT}$ where the ${\rm O}\big{(}\bar{g}_{\rm PT}^{2N-2}\big{)}$
behavior has not clearly set in.

###  2.3 The accuracy of perturbation theory at high energy

In this section we want to review some recent determinations of
$\Lambda_{\overline{\rm MS}}$ which paid particular attention to the issue of
the accuracy of perturbation theory in extracting $\Lambda_{\overline{\rm
MS}}$ Brida:2016flw ; DallaBrida:2018rfy ; DallaBrida:2019wur ; Husung:2020pxg
; Nada:2020jay . As we shall see, the concerns exposed in the previous
sections are legitimate once the precision goals become competitive. A robust
analysis of perturbative truncation errors is essential to reach high
accuracy.

####  2.3.1 The high-energy regime of ${N_{\rm f}}=3$ QCD

Figure 1: Determination of $\Lambda_{\overline{\rm MS}}^{(3)}$ in units of
$\mu_{\rm ref}=L_{0}^{-1}$ for different values of $\alpha_{\rm PT}$
DallaBrida:2018rfy . The extraction in different ${\rm SF}_{\nu}$-schemes
$(\nu=-0.5,0,0.3)$ is shown, as well as a comparison with the final result
$\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}=0.0791(19)$
DallaBrida:2018rfy . As the reader can see, when the extraction is performed
at high-enough energies $(\alpha_{\rm PT}\approx 0.1)$, all schemes nicely
agree.

We begin with the high-energy studies of refs. Brida:2016flw ;
DallaBrida:2018rfy  in ${N_{\rm f}}=3$ QCD. In Fig. 1, we show the results
from these references for $L_{0}\Lambda^{({N_{\rm f}}=3)}_{\overline{\rm MS}}$
as a function of $\alpha^{2}_{\rm PT}\equiv\bar{g}^{4}_{\rm PT}/(4\pi)^{2}$.
In this plot, $\mu_{\rm ref}\equiv L_{0}^{-1}\approx 4.3\,{\rm GeV}$, is a
convenient high-energy reference scale and, as in previous sections,
$\bar{g}_{\rm PT}\equiv\bar{g}_{{\mathcal{O}}}(\mu_{\rm{PT}})$ is the value of
the coupling in the given scheme at which perturbation theory is applied to
extract $\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}$. The ratio
$\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}$ is obtained following
Strategy 1. of Sect. 2.2.3. The scales at which perturbation theory is used
correspond to $\mu_{\rm{PT}}=2^{n}\mu_{\rm ref}$, with $n=0,\ldots,5$, and
range from about $4\,{\rm GeV}$ to $140\,{\rm GeV}$. The couplings considered
in this study, $\bar{g}_{\mathcal{O}}(\mu)=\bar{g}_{{\rm SF}_{\nu}}(L^{-1})$,
belong to a family of finite-volume renormalization schemes based on the QCD
Schrödinger functional (SF) Luscher:1992an ; Sint:1993un ; Sint:1995rb .
Different schemes within the family are identified by different values of the
parameter $\nu$. The precise definition of the schemes is not important and
can be found in refs. Brida:2016flw ; DallaBrida:2018rfy ; Sint:2012ae
.111111Traditionally only the $\nu=0$ scheme has been considered in
applications, see e.g. refs. Luscher:1993gh ; Sint:1995ch ; Bode:1998hd ;
Bode:1999sm ; DellaMorte:2004bc ; Aoki:2009tf  for some important examples.

In order to estimate $\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}$ the
3-loop approximation to the relevant $\beta$-functions, $\beta_{{\rm
SF}_{\nu}}$, is used. The results for $\Lambda^{(3)}_{\overline{\rm
MS}}/\mu_{\rm ref}$ are therefore expected to show corrections of ${\rm
O}(\alpha_{\rm PT}^{2})$ as $\alpha_{\rm PT}\to 0$. It is important to note at
this point that the 3-loop coefficients of the $\beta$-functions in the
different ${\rm SF}_{\nu}$-schemes are given for ${N_{\rm f}}=3$ by
Bode:1998hd ; Bode:1999sm ; DallaBrida:2018rfy :121212For comparison the
3-loop coefficient of $\beta_{\overline{\rm MS}}$ for ${N_{\rm f}}=3$ is
$(4\pi)^{3}b_{2}^{\overline{\rm MS}}=0.324$.

$(4\pi)^{3}b_{2}^{{\rm SF}_{\nu}}=-(0.064(27)+\nu\times 1.259(10))\,.$ (11)

Hence, from the perturbative point of view all schemes with $|\nu|\lesssim 1$
appear to be on similar footing and the perturbative expansion of their
$\beta$-functions is well behaved.

From this observation, one might naively expect that the ${\rm O}(\alpha_{\rm
PT}^{2})$ corrections to $\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}$
obtained from different intermediate schemes are similar too. Going back to
Fig. 1, we see that in all cases the results are well described by a ${\rm
O}(\alpha_{\rm PT}^{2})$ dependence over the whole range of investigated
couplings. This is compatible with the expectation from the known leading non-
analytic term in the expansion which is expected to be quite small at these
couplings, i.e. ${\rm O}(e^{-2.6/\alpha})$ DallaBrida:2018rfy . However, we
clearly see a substantial difference in the size of the ${\rm O}(\alpha_{\rm
PT}^{2})$ corrections depending on the ${\rm SF}_{\nu}$-scheme that is
considered. While one can find cases ($\nu=0.3$) where the ${\rm
O}(\alpha_{\rm PT}^{2})$ corrections are insignificant within errors, other
schemes ($\nu=-0.5$) show significant corrections. The results for $\nu=-0.5$
at $\alpha_{\rm PT}\approx 0.12$, for instance, show a $7-8\%$ deviation from
the final result $\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}=0.0791(19)$
quoted in ref. DallaBrida:2018rfy , which corresponds to the gray band in the
plot.

Figure 2: Determination of $L_{0}\Lambda_{\overline{\rm
MS}}^{(3)}=\Lambda_{\overline{\rm MS}}^{(3)}/\mu_{\rm ref}$ at different
values of $\alpha\equiv\alpha_{\rm PT}$, and using different renormalization
scales (values of $s$) to match with the $\overline{\rm MS}$-scheme
DallaBrida:2018rfy . The left (right) panel uses the ${\rm SF}_{\nu}$-scheme
with $\nu=0$ ($\nu=-0.5$), cf. text.

As the value of the coupling at which perturbation theory is applied becomes
smaller, any significant difference among the different determinations of
$\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}$ steadily fades away. In
particular, once $\alpha_{\rm PT}\approx 0.1$ is reached, any difference is
well below the statistical uncertainties on $\Lambda^{(3)}_{\overline{\rm
MS}}/\mu_{\rm ref}$, which at these couplings are at the level of $2-3\%$. A
robust estimate for $\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}$ can
therefore be obtained by taking its value at $\alpha_{\rm PT}\approx 0.1$ from
one of the schemes that show milder perturbative corrections. The result
$\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}=0.0791(19)$ given above, for
instance, corresponds to $\mu_{\rm PT}=2^{4}\mu_{\rm ref}\approx 70\,{\rm
GeV}$ from the $\nu=0$ scheme (cf. $\nu=0$ (fit C) in Fig. 1)
DallaBrida:2018rfy .

The first important message from this study is that it is in fact impossible
to predict the actual size of perturbative truncation errors only from the
available perturbative information. To reliably assess these errors,
perturbation theory must be tested against non-perturbative data over a wide
range of energy scales. From the study we presented, in particular, we
conclude that in order to be able to quote in full confidence the competitive
precision of $2-3\%$ on $\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}$, one
must reach non-perturbatively $\alpha_{\rm PT}\approx 0.1$. At these couplings
perturbative truncation errors are fully under control and the error on
$\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}$ is entirely dominated by the
statistical uncertainties coming from the non-perturbative running of the
coupling.

It is now instructive to look at the result of the analysis of the same data
according to Strategy 2 of Sect. 2.2.3. The corresponding estimates for
$\Lambda^{(3)}_{\overline{\rm MS}}/\mu_{\rm ref}$ are shown in Fig. 2 for the
two cases, $\nu=0,-0.5$ DallaBrida:2018rfy . The different determinations in
each plot are obtained by varying the parameter $s$ entering the perturbative
matching between the $\overline{\rm MS}$-coupling and the ${\rm
SF}_{\nu}$-couplings (cf. eq. (9)). The values of $s$ considered vary by about
a factor $2-3$ around the value of _fastest apparent convergence_ ,
$s^{*}$.131313The scale factor $s^{*}$ is defined by (cf. eq. (8)):
$s^{*}={\Lambda_{\overline{\rm
MS}}\over\Lambda_{{\mathcal{O}}}}=\exp\bigg{\\{}{c_{1}(1)\over
2b_{0}}\bigg{\\}}\,.$ where $b_{0}$ is the 1-loop coefficient of the
$\beta$-function in eq. (7), and $c_{1}(1)$ is the 1-loop coefficient of the
matching relation between the $\overline{\rm MS}$-coupling
$\bar{g}^{2}_{\overline{\rm MS}}(\mu)$ and the coupling of interest
$\bar{g}^{2}_{\mathcal{O}}(\mu)$ (cf. eq. (9)). With the choice $s=s^{*}$, the
$k=1$ term in eq. (9) vanishes. In phenomenological determinations of the QCD
coupling the spread of the results obtained by varying $s$ around some
“optimal” value, typically by a factor $2$ or so, is commonly used to get an
estimate of perturbative truncation errors (see e.g. ref. Zyla:2020zbs ). Our
intention is to test how this approach works in the present case.

As one can see from Fig. 2, for all choices of $s$ the data show the expected
${\rm O}(\alpha_{\rm PT}^{2})$ scaling. The slope of the data, however, can
vary significantly depending on the choice of the parameter $s$. As expected,
the significance of these differences is reduced as $\alpha_{\rm PT}\to 0$,
and the different determinations come together once $\alpha_{\rm PT}\lesssim
0.1$.

What is clear from the results of Fig. 2 is that the procedure of assigning a
systematic error based on the spread of the results with $s$ at some fixed
coupling is not always reliable. In the case of the $\nu=0$ scheme (left
panel), the spread in the results between, say, $s^{*}/2$ and $2s^{*}$,
encloses the final estimate (gray band in the plot) for all coupling values in
the range. If this uncertainly was added to the statistical errors, it would
give a conservative estimate for the total uncertainly. On the other hand, in
the case of the $\nu=-0.5$ scheme (right panel), the procedure significantly
underestimates the actual size of the ${\rm O}(\alpha_{\rm PT}^{2})$
corrections. Again $\alpha_{\rm PT}\approx 0.1$ has to be reached for the
perturbative uncertainties to be small compared to the statistical ones.

From this second analysis we reaffirm the conclusion that it is very difficult
to reliably estimate perturbative truncation errors if the coupling
$\alpha_{\rm PT}$ cannot be varied much, and if this is confined to values
significantly larger than $\alpha_{\rm PT}\approx 0.1$.

####  2.3.2 The case of the pure Yang-Mills theory

Figure 3: The dimensionless product $\sqrt{8t_{0}}\Lambda_{\overline{\rm
MS}}^{(0)}$ as a function of $\alpha_{\rm PT}$ Nada:2020jay . The empty
symbols represent the data at the given $\alpha_{\rm PT}$, while the filled
symbols are extrapolations to $\alpha_{\rm PT}\to 0$ (shifted for better
readability) of the different approaches to the perturbative matching (see
text for more details). The gray band is the result of ref. DallaBrida:2019wur
, while the data point labeled FlowQCD is the result of ref. Kitazawa:2016dsl
.

##### Finite-volume schemes.

The second example that we consider is taken from the recent study of ref.
Nada:2020jay  in the pure Yang-Mills (YM) theory. This work presents an
independent analysis of the results from a previous study DallaBrida:2019wur ,
using novel techniques. Before entering the discussion, we care to stress that
the case of the pure Yang-Mills theory is not just a curious example. As we
shall see in the following section, through the strategy of renormalization by
decoupling precise results for $\Lambda^{(5)}_{\overline{\rm MS}}$ can be
obtained from the accurate knowledge of $\Lambda^{(0)}_{\overline{\rm MS}}$.
From this perspective, a robust determination of the $\Lambda$-parameter of
the pure YM theory is very relevant.

In Fig. 3 we show the results for $\sqrt{8t_{0}}\Lambda^{(0)}_{\overline{\rm
MS}}$ from ref. Nada:2020jay  as a function of $\alpha_{\rm PT}^{2}$. The
scale $\mu_{\rm had}=1/\sqrt{8t_{0}}$ is defined in terms of the flow time
$t_{0}$ Luscher:2010iy , while $\alpha_{\rm PT}$ is once again the value of
the relevant coupling at the renormalization scale $\mu_{\rm{PT}}$ where
perturbation theory is applied. Similarly to the case of ${N_{\rm f}}=3$ QCD
discussed above, different schemes and strategies have been considered in
order to extract $\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm had}$ and study
the perturbative truncation errors. In all cases, the non-perturbative RG
running from $\mu_{\rm had}$ up to $\mu_{\rm{PT}}$ is obtained using a finite-
volume scheme based on the YM gradient flow (GF) Luscher:2010iy ; Fodor:2012td
; Fritzsch:2013je ; DallaBrida:2016kgh . The interested reader can find more
details about the scheme in the main refs. DallaBrida:2019wur ; Nada:2020jay
(see also Sect. 4.3).

Once $\mu_{\rm{PT}}$ is reached, $\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm
had}$ is estimated in a number of ways. For the case labeled as (GF) in the
plot, Strategy 1. of Sect. 2.2.3 is employed using the 3-loop $\beta$-function
in the GF-scheme of choice DallaBrida:2017tru . In the other cases, the GF-
scheme is first non-perturbatively matched to the ${\rm SF}_{\nu=0}$-scheme
introduced in the previous section. Perturbation theory is then applied either
following Strategy 1. based on the ${\rm SF}_{\nu=0}$-scheme (SF label in the
plot), or by following Strategy 2. and matching the ${\rm SF}_{\nu=0}$\- and
$\overline{\rm MS}$-coupling ($\overline{\rm MS}(s=1,2)$ in the figure). In
the latter case, two values of the $s$-parameter, $s=1,2$, are studied; note
that $s^{*}\approx 2$ in this case. In all cases, the leading parametric
uncertainties in $\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm had}$ from the
truncation of the perturbative expansion are of ${\rm O}(\alpha_{\rm
PT}^{2})$.

Going back to Fig. 3 we see how two out of the four strategies ((SF) and
$\overline{\rm MS}(s=2)$)) give results which are essentially independent on
$\alpha_{\rm PT}$ over the whole range of couplings considered for the
extraction of $\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm had}$. Note that in
going from the largest to the smallest couplings the energy scale varies by a
factor $32$ while $\alpha_{\rm PT}$ changes by about a factor 2. On the other
hand, the other two types of determinations ((GF) and $\overline{\rm
MS}(s=1)$)) show a significant $\alpha_{\rm PT}$ dependence, roughly
compatible with the expected ${\rm O}(\alpha_{\rm PT}^{2})$ scaling. What is
remarkable is that even considering values of $\alpha_{\rm PT}\approx 0.08$
the different strategies give estimates for $\Lambda^{(0)}_{\overline{\rm
MS}}/\mu_{\rm had}$ which vary up to $\approx 3\%$. This is about twice as
large as the statistical errors on the points (cf. Table 3 of ref.
Nada:2020jay ). In the case of the (GF) and ($\overline{\rm MS}(s=1)$)
determinations, it is clear that a trustworthy estimate for
$\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm had}$ can be quoted only by
extrapolating the results for $\alpha_{\rm PT}\to 0$. In general, perturbative
truncation errors are large also in the pure YM theory given the precision one
can reach.

The results above show us once again the importance of an explicit non-
perturbative calculation of the running of the coupling over a significant
range of values, reaching down to small couplings, in order to assess the
actual size of the perturbative corrections. We join the authors of ref.
Nada:2020jay  and conclude that only by studying non-perturbatively the limit
$\alpha_{\rm PT}\to 0$ one can avoid the dangerous game of estimating
perturbative uncertainties at some finite (potentially large) value of
$\alpha_{\rm PT}$. Without studying this limit, the determinations can easily
be affected by perturbative truncation errors, even at surprisingly small
values of the coupling.

##### Large-volume schemes.

A precise determination of the $\Lambda$-parameter in the pure YM theory is
certainly very much facilitated from the computational point of view with
respect to the case of QCD. However, as we have seen in the previous example,
it is yet a non-trivial challenge to control perturbative truncation errors
once a $1-2\%$ precision in $\Lambda_{\overline{\rm MS}}$ is reached.

The disagreement among some recent determinations of
$\Lambda^{(0)}_{\overline{\rm MS}}$ is a clear signal that these difficulties
should not be underestimated. The issue is well illustrated in Fig. 3, where
the very precise results labeled (FlowQCD) from ref. Kitazawa:2016dsl  show a
net tension with the determinations of refs. DallaBrida:2019wur ; Nada:2020jay
. We recall that the former result is based on extracting
$\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm had}$ from the plaquette
expectation value calculated in large-volume simulations. Bare lattice
perturbation theory at couplings $\alpha_{\rm PT}\approx 0.095-0.12$ is used,
with parametric uncertainties of O($\alpha_{\rm PT}^{2}$). We refer the reader
to the given reference for the details. Here we just note that all the above
determinations satisfy the most stringent criteria set by FLAG (cf. ref.
Aoki:2019cca ). Yet, one or more of these results have underestimated
uncertainties.

Other groups have recently engaged in a precision determination of
$\Lambda^{(0)}_{\overline{\rm MS}}$, also with the intent of resolving the
disagreement above. The recent results of ref. Husung:2020pxg  based on the
${\rm qq}$-coupling, $\alpha_{\rm qq}$, defined from the static potential
Sommer:1993ce , are particularly interesting in this respect.141414A similar
earlier study on the challenges of extracting $\Lambda_{\overline{\rm MS}}$ in
both ${N_{\rm f}}=2$ QCD and the pure-gauge theory using the ${\rm
qq}$-coupling can be found in ref. Leder:2011pz . For a recent application of
this scheme for the computation of $\alpha_{\rm s}$ and a detailed account of
the most recent results and developments see refs. Bazavov:2019qoo ;
Komijani:2020kst . We report them in Fig. 4. In the case of $\alpha_{\rm qq}$
the corresponding $\beta$-function, $\beta_{\rm qq}$, is known up to 4-loop
order, and some partial information is available also at 5-loops (see e.g.
ref. Tormo:2013tha ). Determinations of $\Lambda^{(0)}_{\overline{\rm
MS}}/\mu_{\rm had}$ from $\alpha_{\rm qq}$ are hence expected to have
asymptotically ${\rm O}(\alpha_{\rm qq}^{3})$ corrections. In the plot,
$\alpha_{\rm qq}$ refers to the coupling at which perturbation theory is used
according to Strategy 1. of Sect. 2.2.3, i.e. it corresponds to $\alpha_{\rm
PT}$ in our previous discussions.

Figure 4: $\Lambda_{\overline{\rm MS}}^{(0)}$ in units of $\sqrt{8t_{0}}$
determined at various values of $\alpha_{\rm qq}\equiv\alpha_{\rm PT}$
Husung:2020pxg . The results using different orders of perturbation theory for
$\beta_{\rm qq}$ are shown, as well as a comparison with the determinations of
refs. DallaBrida:2019wur  and Kitazawa:2016dsl . (The reference numbers in the
plot are those of ref. Husung:2020pxg  from which the plot is taken.)

Despite the accurate perturbative knowledge there are a few challenges when
using the ${\rm qq}$-scheme for precision determinations of
$\Lambda_{\overline{\rm MS}}^{(0)}$ Husung:2020pxg . The most relevant ones
for our discussion are, first of all, that the scheme is conventionally
defined in an infinite space-time volume. In order to measure the coupling at
small lattice spacings one therefore needs large lattice sizes to maintain the
physical extent of the lattice large. In the computation of ref.
Husung:2020pxg  lattice spacings down to $a\approx 0.01\,{\rm fm}$ are reached
while keeping the lattice extent $L\approx 2\,{\rm fm}$. This means simulating
lattices with up to $L/a\approx 200$. Secondly, the perturbative expansion of
$\alpha_{\rm qq}$ displays some infrared divergences starting at 3-loop order
in $\beta_{\rm qq}$. When resummed these give rise to terms of the form,
$\alpha_{\rm qq}^{n}\log(\alpha_{\rm qq})^{m}$, $n\geq 3$, $1\leq m\leq n-2$,
which are enhanced at small couplings (cf. ref. Husung:2020pxg ).

Figure 4 shows the results for $\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm
had}$ as a function of $\alpha_{\rm qq}^{3}$. The range of couplings covered
by the data is $\alpha_{\rm qq}\approx 0.16-0.35$. As we can see from the
plot, for couplings $\alpha_{\rm qq}\gtrsim 0.21$ the results for
$\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm had}$ have good precision, but
perturbative uncertainties are large. This can be seen by looking at the
difference between the 3-loop and 4-loop results (or analogously between the
4-loop and 4-loop + 5-loop log-terms results). At these large couplings, the
perturbative expansion seems to have reached its limit of applicability. This
severely limits the precision one can aim at for $\Lambda^{(0)}_{\overline{\rm
MS}}/\mu_{\rm had}$ if one is restricted to this range of couplings. For
couplings $\alpha_{\rm qq}\lesssim 0.21$, the different orders of perturbation
theory seem to start converging. On the other hand, the errors on the data
become large. This is due to the difficulties in extrapolating the results to
the continuum limit Husung:2020pxg . In fact, the errors are too large to make
definite conclusions for the relevant limit $\alpha_{\rm qq}\to 0$.

All in all, we see from this last example that a precise determination of
$\Lambda_{\overline{\rm MS}}^{(0)}$ is a challenge. Finite-volume
renormalization schemes allow us to cover a wide range of couplings, reaching
down to rather small values. Yet, having control on perturbative truncation
errors requires care. When using large-volume schemes the situation is further
complicated by controlling continuum limit extrapolations at the smallest
(most relevant) couplings. Small couplings require small lattice spacings,
which require large lattice sizes in order to keep the physical volume large.
As a result, even in the computationally simpler case of the pure YM theory,
one might have precise data confined for the most part to a region of
couplings too large to have perturbative uncertainties fully under control,
while at smaller couplings the data is not precise enough for a competitive
determination of $\Lambda$.

###  2.4 The tricky business of continuum extrapolations

Having discussed the difficulties of estimating perturbative truncation errors
in precision determinations of the $\Lambda$-parameters, we now want to touch
on the issue of systematic uncertainties related to the continuum limit
extrapolations of the relevant couplings.

To give the reader a feeling of the pitfalls that these continuum
extrapolations can conceal, we first consider the ${N_{\rm f}}=3$ results of
ref. DallaBrida:2016kgh . The relevant quantity to look at in this case is the
step-scaling function (SSF) of the finite-volume GF-coupling $\bar{g}^{2}_{\rm
GF}(\mu)$ with SF boundary conditions (see refs. Fritzsch:2013je ;
DallaBrida:2016kgh  and eq. (53) for the definition of this scheme). We recall
that the SSF, $\sigma(u)$, encodes the change in the coupling $u$ when the
renormalization scale is varied by a factor of 2 Luscher:1991wu .
Specifically, having set $\mu=L^{-1}$,

$\sigma(u)\equiv\bar{g}^{2}_{\rm GF}(\mu/2)|_{\bar{g}^{2}_{\rm
GF}(\mu)=u},\quad\sigma(u)=\lim_{a/L\to 0}\Sigma(a/L,u)\,.$ (12)

It is clear from its definition that the SSF is a discrete version of the
$\beta$-function. The latter can in fact be obtained once the SSF is known in
a range of couplings (cf. ref. DallaBrida:2016kgh ).

On the lattice, the SSF is determined by extrapolating to the continuum limit
its discrete approximations, $\Sigma(a/L,u)$. In order to compute the latter
one must first identify a set of lattice sizes $L/a$ and corresponding values
of the bare coupling $g_{0}$ for which $\bar{g}^{2}_{\rm GF}(L^{-1})=u$, with
$u$ a specific value. The lattice SSFs $\Sigma(a/L,u)$ are then given by the
couplings $\bar{g}^{2}_{\rm GF}((2L)^{-1})$ measured at the values of $g_{0}$
previously determined but on lattices with sizes $2L/a$.

The results for the lattice SSFs of the GF-coupling of ref. DallaBrida:2016kgh
are shown in Fig. 5. They correspond to 9 values of the coupling
$u_{i}\in[2.1,6.5]$, $i=1,\ldots,9$. As one can see from the figure, the
lattice data are very precise. On the other hand, discretization effects are
in general large, particularly so at the largest couplings. The results for
$\Sigma(u,a/L)$ vary in fact by up to 20% in the range of $L/a$ considered,
which is quite a significant change compared to the statistical errors on the
points.

Given the results in Fig. 5, we may expect that a simple fit of the data
linear in $(a/L)^{2}$ is all that is needed to extrapolate these to the
continuum limit. In particular, we may consider individual continuum
extrapolations for each $u_{i}$ value using the functional form

$\Sigma(u_{i},a/L)=\sigma^{(A)}_{i}+{r}^{(A)}_{i}\times(a/L)^{2}\quad\text{(fit
A)}\,,$ (13)

where $\sigma^{(A)}_{i}$, ${r}^{(A)}_{i}$ are fit parameters. Within the
uncertainties, linearly in $(a/L)^{2}$ is in fact excellent and the above fits
are very good ($\chi^{2}/{\rm dof}\approx 0.7$). One is thus tempted to take
the precise values for $\sigma^{(A)}_{i}$ as estimates for the continuum SSF.

The continuum results so obtained are well described by the simple relation:
$\Delta\sigma_{i}^{(A)}\equiv 1/\sigma^{(A)}_{i}-1/u_{i}\approx-0.082$. Note
that this is the functional form expected from the perturbative expansion of
$\sigma(u)$ at 1-loop order, although the coefficient predicted by
perturbation theory is slightly different, i.e.
$\approx-0.079$.151515Perturbation theory predicts: $\sigma(u)\overset{u\to
0}{\approx}u+s_{0}u^{2}+{\rm O}(u^{3})$, with $s_{0}=2b_{0}({N_{\rm
f}})\ln(2)$, where $b_{0}({N_{\rm f}})$ is the universal 1-loop coefficient of
the $\beta$-function, eq. (7) (see e.g. ref. DallaBrida:2018rfy ). The close
agreement between the non-perturbative data for $\sigma(u)$ and 1-loop
perturbation theory is quite peculiar, considering the fact that it holds up
to $\sigma(u)={\rm O}(10)$. We refer the interested reader to ref.
DallaBrida:2016kgh  for a detailed discussion about this point. This
observation suggests us to perform alternative fits to the data in Fig. 5
considering the functional form

$1/\Sigma(u_{i},a/L)=1/\sigma^{(B)}_{i}+r^{(B)}_{i}\times(a/L)^{2}\quad\text{(fit
B)}\,,$ (14)

with $\sigma^{(B)}_{i}$, ${r}^{(B)}_{i}$ new fit parameters. The quality of
these fits is as good as for the fits A of eq. (13). Distinguishing between
the two fit forms would require significantly higher statistical precision
than the present one.

It is important to note at this point that any functional form that we
consider for the continuum extrapolations, necessarily comes with assumptions.
Both fits A and B above, for instance, assume discretization errors of
O($a^{n}$), $n\geq 3$, to be negligible. Moreover, even focusing only on the
leading O($a^{2}$) effects, we know from Symanzik effective theory (SymEFT)
Symanzik:1981hc ; Symanzik:1983dc ; Symanzik:1983gh  that these are not simply
given by a “classical” term $\propto(a/L)^{2}$. They are in fact a non-trivial
combination of different terms which in the limit $a/L\to 0$ are
asymptotically $\propto(a/L)^{2}\ln(L/a)^{-\Gamma_{i}}$, or rather
$\propto(a/L)^{2}[\bar{g}^{2}(a^{-1})]^{\Gamma_{i}}$, where
$\Gamma_{i}\in\mathbb{R}$, $i=1,2,\ldots,$ and $\bar{g}^{2}(a^{-1})$ is the
given renormalized coupling of the effective theory evaluated at a scale
$\mu=a^{-1}$ (cf. refs. Balog:2009yj ; Balog:2009np ; Husung:2019ytz ;
Husung:2020 ).

Figure 5: Continuum extrapolations of $\Sigma(a/L,u)/u$ for 9 different values
of $u\in[2.1,6.5]$ and 3 lattice sizes $L/a=8,12,16$.

If terms of higher order than $a^{2}$ as well as logarithmic corrections to
pure $a^{2}$ scaling were completely negligible in the data, the fit
parameters $\sigma^{(A)}_{i}$ and $\sigma^{(B)}_{i}$ should perfectly agree.
From the results in Fig. 5 we see that there is in fact agreement within one
standard deviation. However, the difference between the results from the two
fits is clearly systematic, with the results from fit B being always larger
than those from fit A.

The issue becomes more evident if one tries to obtain a smooth
parameterization for the continuum SSF from the fitted continuum values
$\sigma_{i}$. As noticed earlier, a fit of
$\Delta\sigma_{i}=1/\sigma_{i}-1/u_{i}$ to a constant $\Delta\sigma$ provides
a good description of the continuum data $\sigma_{i}$ in the whole range of
$u\in[2.1,6.5]$; this is the case for both $\Delta\sigma^{(A)}_{i}$ and
$\Delta\sigma^{(B)}_{i}$ ($\chi^{2}/{\rm dof}<1$). Given the 9 independent
values of $\Delta\sigma_{i}$ for each fit, the results for the corresponding
$\Delta\sigma$ are 3 times more precise than the individual
$\Delta\sigma_{i}$. The systematic effect then becomes clearly noticeable as
one finds: $\Delta\sigma^{(A)}=-0.0823(4)$ and
$\Delta\sigma^{(B)}=-0.0832(4)$, for the constant fits to
$\Delta\sigma^{(A)}_{i}$ and $\Delta\sigma^{(B)}_{i}$, respectively.

The previous considerations show that the description of discretization errors
as pure $(a/L)^{2}$ effects is in this case not accurate enough for the level
of precision claimed in the continuum limit. Even though the different
functional forms in eqs. (13) and (14) fit the data well and perfectly agree
with each other at finite $L/a$, the corresponding extrapolations for $a/L\to
0$ are clearly affected by some systematics. In ref. DallaBrida:2016kgh , more
conservative error estimates and robust central values for the continuum
results are eventually obtained by carefully accounting as systematic
uncertainties in the data the not entirely negligible effects of the higher-
order terms neglected in eqs. (13) and (14) (cf. ref. DallaBrida:2016kgh  for
the details).

The example above might not seem too pessimistic. However, it should come as a
warning for the more general situation. Estimating properly the systematic
uncertainties related to continuum limit extrapolations of high-precision data
can easily become a challenge, particularly so if discretization errors are
not small.

As recalled earlier, the leading asymptotic dependence of renormalized lattice
quantities on the lattice spacing as $a\to 0$ is given by a combination of
terms $\propto a^{n}[\bar{g}^{2}(a^{-1})]^{\Gamma_{i}}$, where the number and
values of the $\Gamma_{i}$, as well as $n$, depend on the chosen
discretization and set-up. The $\Gamma_{i}$ are in fact inferred from the
anomalous dimensions of the fields defining the O($a^{n}$) counterterms in the
SymEFT, and the order of perturbative improvement that has been possibly
implemented (cf. ref. Husung:2019ytz ). Hence, when one considers a pure
$a^{n}$ dependence for the discretization errors, one is implicitly assuming
that all $\Gamma_{i}\approx 0$. This, however, cannot be taken for granted.

For most cases of interest, the leading discretization effects have $n=2$,
i.e. they are of O($a^{2}$).161616A relevant exception is the case of the SF,
for which the leading discretization errors are parametrically of O($a$) (cf.
Sect. 4.2.3). In applications, however, the O($a$) effects are subdominant
with respect to the O($a^{2}$) effects, and often also compared to the
statistical errors. The precision studies of refs. DallaBrida:2016kgh ;
DallaBrida:2018rfy ; DallaBrida:2019wur  thus opt for treating the O($a$)
effects as (small) systematic uncertainties in the data, and perform continuum
extrapolations assuming leading O($a^{2}$) effects. In this respect, we note
that in refs. Husung:2019ytz ; Husung:2020  the $\Gamma_{i}$ relevant for the
O($a$) effects in the (pure-gauge) SF have been computed. The results support
the treatment of O($a$) effects pursued in refs. DallaBrida:2016kgh ;
DallaBrida:2018rfy ; DallaBrida:2019wur  (cf. the given references for the
details). The results of refs. Husung:2019ytz ; Husung:2020  then show that in
the case of QCD we have O(10) different terms that contribute in general, and
$\Gamma_{i}\approx[-0.1,3]$ for several common discretizations and values of
${N_{\rm f}}=2-4$.171717The results refer to the contributions to
discretization effects coming from the lattice action, considering several
popular options (cf. ref. Husung:2020 ). If the relevant observable is not a
spectral quantity, additional effects originating from the lattice fields that
define it are present. These depend on the specific observable and choice of
discretization (see, e.g. refs. Husung:2019ytz ; Husung:2020 ).181818In the
case of the pure-gauge theory only two terms from the lattice action
contribute to the O($a^{2}$) effects. The $\Gamma_{i}$ for different options
can be found in ref. Husung:2019ytz . In all cases, $\Gamma_{i}\gtrsim 0.6$.
Having all $\Gamma_{i}\gtrsim 0$ is certainly positive. In particular, the
contributions relevant in the massless theory all have $\Gamma_{i}>0$, which
implies a faster approach to the continuum limit with respect to pure $a^{2}$
terms. However, the large number of terms contributing makes for a complicated
pattern of discretization errors in the general case, with no clear
contribution(s) dominating. As a result, it may be difficult in practice to
identify the terms that are actually relevant. Moreover, terms of the form
$a^{2}[\bar{g}^{2}(a^{-1})]^{\Gamma_{i}}$ with $\Gamma_{i}\approx 2-3$ can be
hard to distinguish from $a^{3}$ or $a^{4}$ terms in a limited range of
lattice spacings when statistical uncertainties are present.191919It is clear
that even though the SymEFT can predict the form of the leading asymptotic
discretization errors, it cannot predict the region where these dominate over
formally suppressed contributions. In practice, it may thus be difficult to
establish the regime of applicability of the results from SymEFT. The
continuum estimates obtained by including different contributions, on the
other hand, may vary appreciably. In this situation, precise and robust final
estimates are not easily achieved.

We stress that it is particularly important to take these considerations into
account when aiming for precise determinations of short-distance quantities
like the couplings. As discussed in previous sections, in the most interesting
region of high energy, $\mu\gg\Lambda$, $a\mu$ may not be so small. Continuum
extrapolations are thus likely to be difficult and require special attention.
Following the lines of refs. Husung:2019ytz ; Husung:2020  one should take the
non-trivial $a$-dependence predicted by SymEFT into account, provided the
information is available. If this is not the case, one should try at least to
estimate the uncertainties associated with neglecting logarithmic corrections
to classical scaling, e.g. by considering terms $\propto
a^{2}[\bar{g}^{2}(a^{-1})]^{\Gamma_{i}}$, with $\Gamma_{i}\approx 1-3$, in the
fit ansätze. Ideally, one would like to be in the situation where within the
uncertainties the continuum estimates do not sensibly depend on whether these
terms are considered or not.

Given the observations above, we want to bring the reader’s attention to a
recent study where the non-trivial $a$-dependence of discretization effects
was found to be a relevant issue. Specifically, we consider the computations
of refs. DallaBrida:2016dai ; DallaBrida:2017tru  of the GF-coupling in the
pure Yang-Mills theory using Numerical Stochastic Perturbation Theory.202020A
recently expanded discussion in ref. Nada:2020jay  provides another clear
illustration of the difficulties of continuum limit extrapolations of precise
coupling data using results from the pure-gauge theory (cf. Figs. 6 and 7 of
this reference and related discussion). We strongly recommend the interest
reader to consult this reference. We moreover refer to the important
pioneering studies of refs. Balog:2009yj ; Balog:2009np  in the non-linear
$\sigma$-model in two dimensions. In this framework, the lattice theory is
numerically solved through a Monte Carlo simulation up to a finite order in
the bare coupling $g_{0}$ DiRenzo:1994sy ; DiRenzo:2004hhl . From expectation
values in this “truncated theory” one can obtain the perturbative coefficients
of the expansion of lattice quantities in $g_{0}$.

Figure 6: Continuum extrapolations for $k_{1}(a/L)$. Results for $c=0.3$ (left
panel) and $c=0.4$ (right panel) are shown, for two different discretizations
of the relevant observable. Different fits to the data are considered, cf.
text.

In refs. DallaBrida:2016dai ; DallaBrida:2017tru , the GF-coupling with SF
boundary conditions $\bar{g}^{2}_{\rm GF}(\mu)$ has been computed up to two-
loop order in $g_{0}^{2}$. Using the relation between $\alpha_{0}\equiv
g_{0}^{2}/(4\pi)$ and $\alpha_{\overline{\rm MS}}\equiv\alpha^{({N_{\rm
f}}=0)}_{\overline{\rm MS}}$ Luscher:1995nr ; Luscher:1995np , one can thus
infer the relation

$\alpha_{\rm GF}(\mu)=\alpha_{\overline{\rm
MS}}(\mu)+k_{1}(a/L)\alpha^{2}_{\overline{\rm
MS}}(\mu)+k_{2}(a/L)\alpha^{3}_{\overline{\rm MS}}(\mu)+\ldots.$ (15)

The coefficients $k_{1}(a/L)$, $k_{2}(a/L)$ are functions of the resolution
$a/L$ considered for the lattice. In order to obtain the matching relation
between the couplings in the continuum limit the coefficients must be
extrapolated for $a/L\to 0$. Focusing on the 1-loop coefficient, $k_{1}(a/L)$,
from SymEFT we expect that (see refs. DallaBrida:2016dai ; DallaBrida:2017tru
)

$k_{1}(a/L)\overset{a/L\to
0}{\sim}k_{1}(0)+\sum_{m=2}^{\infty}\sum_{n=0}^{1}r_{mn}(a/L)^{m}\ln(L/a)^{n}\,,$
(16)

with $r_{mn}$ some constants. Note that the coefficient $r_{21}$ of the
leading term $\propto\ln(L/a)$ implicitly depends on the $\Gamma_{i}$
predicted by SymEFT (cf. Sect. 5.2. of ref. Husung:2019ytz  and also ref.
Husung:2020 ). Compared to the case of the full theory, the results from the
truncated theory have a simpler (yet non-trivial) cutoff dependence. Given the
high precision reached in these calculations, this allows for a clean
illustration of the difficulties in continuum extrapolations.

In Fig. 6 we show the results from ref. DallaBrida:2016dai  for $k_{1}(a/L)$
for two different values of the parameter $c$, $0.3$ and $0.4$, that specifies
the GF-scheme (cf. refs. DallaBrida:2016dai ; DallaBrida:2017tru  and eq.
(53)). Two different discretizations of the observable defining the coupling
($E^{\rm clov}_{\rm mag}$, $E^{\rm plaq}_{\rm mag}$) are also considered. The
simulated lattices have sizes $L/a=10-32$.

Starting from the results for $c=0.3$ (left panel of Fig. 6), we see how the
data is very precise but discretization effects are sizable. In the plot we
then show two types of extrapolations to the continuum limit. For the first
type (solid lines), lattices with $L/a=12-32$ are fitted to the asymptotic
form, eq. (16), considering the leading terms $m=2,n=0,1$. The fits are good,
$\chi^{2}/{\rm dof}\sim 1$, and the extrapolated results for the two
discretizations agree well. The $m=2,n=1$, term is in fact crucial to obtain
good fits. For the second set of extrapolations (dashed lines), we consider
instead lattices with $L/a=12-24$. In this case the data can be very well
described by a pure $(a/L)^{2}$ term ($m=2,n=0$) over the whole range of
lattice sizes. The continuum extrapolated values obtained from these fits have
significantly smaller statistical errors than the ones from the previous fits,
and yet there is perfect agreement between the two discretizations. On the
other hand, the results deviate from the previous estimates by several of
their standard deviations.

The results for $c=0.4$ exhibit qualitatively the same features, although the
statistical errors on $k_{1}(a/L)$ are about a factor 2 larger and the two
discretizations now show rather different lattice artifacts. On the other
hand, cutoff effects are generally smaller than for $c=0.3$, and we thus
include $L/a=10$ in the fits. It is clear that in both cases, $c=0.3,0.4$, a
reliable continuum extrapolation for $k_{1}(a/L)$ is challenging due to the
non-trivial $a$-dependence of the data. In particular, larger lattices than
the ones considered here are clearly needed in order to obtain accurate
continuum results (cf. ref. DallaBrida:2017tru  for the final determination).

In conclusion, through these examples we saw how assessing the systematics
related to the continuum limit extrapolations of couplings can be challenging.
This is especially true when one wants to maintain the high precision reached
on the lattice data also in the continuum limit, but discretization errors are
large. It then becomes hard to avoid systematic biases in the final
determinations. To this end, it is crucial to test all the assumptions that
enter the functional forms chosen for the extrapolations. In particular, we
must keep in mind that good fits do not necessarily mean good results for
parameters, especially for extrapolations outside the range covered by the
data.

##  3 Heavy-quark decoupling

So far we focused on the main challenges that stand on the way of a precise
determination of $\Lambda^{({N_{\rm f}})}_{\overline{\rm MS}}$ and discussed
in detail the cases of ${N_{\rm f}}=0,3$. The interesting quantity for
phenomenology, however, is $\Lambda^{(5)}_{\overline{\rm MS}}$. At present,
lattice estimates of $\Lambda^{(5)}_{\overline{\rm MS}}$ are for the most part
based on determinations of $\Lambda^{(3)}_{\overline{\rm MS}}$, while just a
handful are obtained from $\Lambda^{(4)}_{\overline{\rm MS}}$ (cf. ref.
Aoki:2019cca ). As we shall recall in the next subsection, the most common
strategy to obtain $\Lambda^{(5)}_{\overline{\rm MS}}$ is in fact to non-
perturbatively compute $\Lambda^{(3)}_{\overline{\rm MS}}$ through simulations
of the ${N_{\rm f}}=3$ theory and then rely on perturbative decoupling
relations for the heavy quarks to estimate the ratios
$\Lambda^{(4)}_{\overline{\rm MS}}/\Lambda^{(3)}_{\overline{\rm MS}}$ and
$\Lambda^{(5)}_{\overline{\rm MS}}/\Lambda^{(4)}_{\overline{\rm MS}}$ (see
e.g. ref. Aoki:2019cca ).

The main reason for this is because, as is well-known, simulating the charm
quark dynamically is at present challenging, let alone the case of the bottom
quark. While the inclusion of the charm quark in the computation of the
running of the QCD coupling may be only moderately challenging with a suitable
strategy (see e.g. ref. Tekin:2010mm ), it does pose important difficulties in
large-volume hadronic simulations. Besides the increased computational cost in
simulating an additional quark with respect to ${N_{\rm f}}=2+1$ simulations,
and the more complicated tuning of the bare QCD parameters necessary to define
proper lines of constant physics, discretization effects are a serious source
of concern. Given the currently most accessible lattice spacings in hadronic
simulations, say $a\gtrsim 0.05\,{\rm fm}$, we have that $am_{c}\gtrsim 0.3$
and $am_{b}\gtrsim 1$, where for definiteness we took $m_{c}\approx 1.27\,{\rm
GeV}$ and $m_{b}\approx 4.2\,{\rm GeV}$. In the hadronic regime it is
therefore a real challenge to control the discretization effects induced by
including the charm quark in simulations, and unrealistic for the case of the
bottom quark. This is particularly true for the case of Wilson quarks where
the charm quark can potentially introduce large O($am_{c}$) effects, unless a
complete Symanzik O($a$) improvement programme is carried out, which is
certainly no simple task (see e.g. ref. Fritzsch:2018kjg ).

In this situation, it is mandatory to assess the reliability of the strategy
presented above for the determination of $\Lambda^{(5)}_{\overline{\rm MS}}$.
To this end, in the following we shall recall the general theory of decoupling
of heavy quarks and critically address its application in lattice
determinations of $\Lambda^{(5)}_{\overline{\rm MS}}$. This includes both the
usage of perturbation theory for the inclusion of heavy-quark loops in the
running of the QCD coupling, that is to estimate the ratios
$\Lambda_{\overline{\rm MS}}^{(4)}/\Lambda_{\overline{\rm MS}}^{(3)}$,
$\Lambda_{\overline{\rm MS}}^{(5)}/\Lambda_{\overline{\rm MS}}^{(4)}$, as well
as the determination of the physical units of $\Lambda^{(5)}_{\overline{\rm
MS}}$ from scale setting in the ${N_{\rm f}}=2+1$ theory. As we shall see,
given the current precision on $\Lambda^{(3)}_{\overline{\rm MS}}$, accounting
for heavy-quark effects by means of perturbation theory in the running of the
QCD coupling is remarkably accurate, even for the case of the charm quark. In
addition, charm-quark effects in (dimensionless) low-energy quantities are
found to be quite small, supporting the fact that ${N_{\rm f}}=2+1$ QCD is
accurate enough for establishing the physical scale. As a result, competitive
determinations of $\Lambda^{(5)}_{\overline{\rm MS}}$ are possible from
results in the ${N_{\rm f}}=3$ flavor theory.

###  3.1 The effective theory for heavy-quark decoupling and the QCD couplings

In this subsection we introduce the effective theory of heavy quarks and
recall how this is conventionally applied in the determination of
$\Lambda_{\overline{\rm MS}}^{(5)}$. We refer the reader to refs.
Bruno:2014ufa ; Athenodorou:2018wpk  for a more detailed presentation.

####  3.1.1 The effective theory for heavy-quark decoupling

We begin by considering QCD with ${N_{\rm f}}$ flavors of quarks, which in
short we denote ${\rm QCD}_{{N_{\rm f}}}$. Of these, ${N_{\ell}}$ are
considered to be light, while the other ${N_{\rm h}}\equiv{N_{\rm
f}}-{N_{\ell}}$ are heavy. For simplicity, we assume that the light quarks are
degenerate with mass $m$, while the heavy quarks are also degenerate but with
a mass $M\gg\Lambda$. The effective theory associated with the decoupling of
the heavy quarks is formally obtained by integrating out in the functional
integral the fields associated with the heavy quarks Weinberg:1980wa . The
field theory that results is characterized by having an infinite number of
non-renormalizable interactions, which are suppressed at low energies by
negative powers of the heavy-quark masses $M$. The couplings of the effective
theory can be fixed order by order in $M^{-1}$ by requiring that, at each
given order, a finite number of observables is equal to the corresponding ones
in the fundamental theory. Once the couplings are fixed up to a certain order
$M^{-n}$, the effective theory is said to be matched to the fundamental one at
this order, and can be used to describe the effects of the heavy quarks at low
energies up to corrections of O($M^{-n-1}$). In this sense, we say that as
$M\to\infty$ the heavy quarks decouple from low-energy physics as their
effects eventually fade away Appelquist:1974tg .

In formulas, the Lagrangian of the effective theory is of the general form
(see e.g. ref. Athenodorou:2018wpk )

$\mathcal{L}_{\rm dec}=\mathcal{L}_{0}+{1\over M}\mathcal{L}_{1}+{1\over
M^{2}}\mathcal{L}_{2}+\ldots\,,$ (17)

where the leading order corresponds to the Lagrangian of ${\rm QCD}$ with
${N_{\ell}}$ light quarks, i.e. $\mathcal{L}_{0}=\mathcal{L}_{\rm
QCD_{N_{\ell}}}$, while the corrections $\mathcal{L}_{k}$, $k\geq 1$, consist
of linear combinations of local fields of mass dimension $4+k$, i.e.

$\mathcal{L}_{k}=\sum_{i}\omega^{(k)}_{i}\Phi^{(k)}_{i}\,,\qquad[\Phi^{(k)}_{i}]=4+k\,,$
(18)

with $\omega^{(k)}_{i}$ dimensionless couplings. The fields $\Phi^{(k)}_{i}$
are built from the light-quark and gluon fields, and include possible powers
of the light-quark masses. They must respect the symmetries of the fundamental
${\rm QCD}_{{N_{\rm f}}}$ theory, as in particular gauge invariance, Euclidean
(or Lorentz) symmetry, and chiral symmetry.

In the case where the light quarks are massless, the leading-order effective
theory, ${\rm QCD}_{N_{\ell}}$, has a single parameter: the gauge coupling
$\bar{g}^{({N_{\ell}})}(\mu)$. The effective and fundamental theory are
therefore matched at leading order in $M^{-1}$ once
$\bar{g}^{({N_{\ell}})}(\mu)$ is matched. This requires that
$\bar{g}^{({N_{\ell}})}(\mu)$ is properly prescribed at a given
renormalization scale in a given renormalization scheme in terms of the
coupling $\bar{g}^{({N_{\rm f}})}(\mu)$ of the fundamental theory and the
heavy-quark masses $M$.212121For ease of notation we do not use any symbol to
indicate the generic scheme of renormalized couplings. We however assume that,
unless otherwise stated, the couplings are defined in a mass-independent
scheme. In addition, provided that the fundamental theory is defined on a
manifold without boundaries,222222This means that the theory lives in infinite
space-time or in a finite volume with (some variant of) periodic boundary
conditions. The special but relevant case of a finite volume with boundaries
will be considered in Sect. 4.2.3. it is possible to show that
$\mathcal{L}_{1}=0$ and therefore O($M^{-1}$) corrections are absent
Athenodorou:2018wpk .232323Note that for the sake of argument we exclude the
uninteresting case of ${N_{\ell}}=1$, for which a proof of this result is to
our knowledge missing. In this situation, the leading-order corrections
induced by the heavy quarks are suppressed as $M^{-2}$ at low energy.

In the case where the light quarks have a non-vanishing mass, the only mass-
dimension 5 fields allowed in $\mathcal{L}_{1}$ are given by the fields of the
leading-order Lagrangian $\mathcal{L}_{\rm QCD_{N_{\ell}}}$ multiplied by the
light-quark masses $m$ Athenodorou:2018wpk . Their effect can be reabsorbed
into a redefinition of the gauge coupling and light-quark masses of O($m/M$).
Of course, in the case of massive light quarks the matching of the effective
and fundamental theory at leading order also requires the matching of the
light-quark masses. In addition, the couplings $\omega_{i}^{(k)}$ of the
corrections $\mathcal{L}_{k\geq 2}$ will depend in general on the light-quark
masses, too.

####  3.1.2 The effective theories and their couplings

The application of the effective theory for heavy quarks in the determination
of the QCD couplings was first advocated by Weinberg in his seminal paper on
effective field theories Weinberg:1980wa . The idea is based on the
observation that for mass-independent renormalization schemes the RG equations
of the renormalizable couplings of the effective theory completely decouple
from the others. This means that the couplings of the non-renormalizable
interactions can be completely ignored when determining the variation of the
running coupling $\bar{g}^{({N_{\ell}})}(\mu)$ of the effective theory with
the energy scale $\mu$. The heavy quarks affect the value of the coupling of
the effective theory only through the matching with the coupling of the
fundamental theory, $\bar{g}^{({N_{\rm f}})}(\mu)$.

The matching relation between $\bar{g}^{({N_{\ell}})}(\mu)$ and
$\bar{g}^{({N_{\rm f}})}(\mu)$ can in principle be established in perturbation
theory. This is best done at a scale $\mu_{\rm match}\approx M$
Weinberg:1980wa ; Bernreuther:1981sg . Assuming the validity of perturbation
theory at this scale, if the running coupling $\bar{g}^{({N_{\ell}})}(\mu_{\rm
match})$ in the effective theory is known, one can turn tables and obtain
$\bar{g}^{({N_{\rm f}})}(\mu_{\rm match})$ from inverting the matching
conditions.

In phenomenological applications of this strategy (see e.g. ref. Zyla:2020zbs
), the value of $\bar{g}^{({N_{\ell}})}(\mu_{\rm match})$ is extracted from
the value of the coupling $\bar{g}^{({N_{\ell}})}(\mu_{\rm low})$ at some
lower energy scale, $\mu_{\rm low}\ll M$. The latter is obtained by comparing
the perturbative expansion for some process ${\mathcal{O}}(q)$ with
characteristic energy scale $q\approx\mu_{\rm low}$, with its experimental
results. The effects of the heavy quarks in ${\mathcal{O}}(q)$ are expected to
be suppressed as O($(q/M)^{2}$) (cf. Sect. 3.1.1). Hence, assuming that these
effects can be neglected, the perturbative expansion of ${\mathcal{O}}(q)$ can
be considered in the ${N_{\ell}}$-flavor theory, which allows the coupling
$\bar{g}^{({N_{\ell}})}(\mu_{\rm low})$ to be extracted. As observed above,
the determination of the running of the effective coupling in the
${N_{\ell}}$-flavor theory does not require any input from the fundamental
theory: one can thus readily obtain $\bar{g}^{({N_{\ell}})}(\mu_{\rm match})$
from $\bar{g}^{({N_{\ell}})}(\mu_{\rm low})$. Clearly, for this strategy to
work in practice the energy scale $\mu_{\rm low}\ll M$ must be yet
sufficiently high for perturbation theory to apply.

In non-perturbative applications on the lattice, the strategy presented above
is realized in the following way (see e.g. ref. Aoki:2019cca ). Firstly, as
discussed in Sect. 2.2.3, through the study of the non-perturbative running of
a given massless coupling within the effective ${N_{\ell}}$-flavor theory, one
determines the ratio $\Lambda_{\overline{\rm MS}}^{({N_{\ell}})}/\mu_{\rm
had}$, where $\mu_{\rm had}$ is a convenient (not necessarily physical) low-
energy scale. Assuming that the effects of the heavy quarks can be neglected
in the ratio of low-energy scales $\mu_{\rm had}/\mu_{\rm phys}$, where
$\mu_{\rm phys}$ is an experimentally accessible hadronic quantity, this can
also be computed within ${N_{\ell}}$-flavor QCD. The physical units of
$\mu_{\rm had}$ and therefore $\Lambda_{\overline{\rm MS}}^{({N_{\ell}})}$ can
then be established by taking $\mu_{\rm phys}$ from experiments. As a second
step, as we shall see in detail in the following subsection, the matching
relation between the couplings of the effective and fundamental theory is
expressed as a relation between their $\Lambda$-parameters. The ratio
$\Lambda_{\overline{\rm MS}}^{({N_{\rm
f}})}/\Lambda^{({N_{\ell}})}_{\overline{\rm MS}}$ is thus estimated using
perturbation theory at a scale $\mu\approx M$. From this estimate and the
results for $\Lambda_{\overline{\rm MS}}^{({N_{\ell}})}$,
$\Lambda_{\overline{\rm MS}}^{({N_{\rm f}})}$ is obtained.

It is important to note that in extracting $\Lambda_{\overline{\rm
MS}}^{({N_{\ell}})}$ from the running of the chosen non-perturbative scheme,
perturbation theory can be applied at arbitrarily large energy scales. The
perturbative matching between the effective and fundamental theory, instead,
is best performed at a scale $\mu\approx M$, where $M$ is set by the mass of
the given heavy quarks that decouple. Whether perturbation theory is accurate
in this step hence depends on how heavy these quarks are. As we shall see, in
the $\overline{\rm MS}$-scheme the perturbative matching is remarkably
accurate already for masses $M$ close to that of the charm quark.

In the next subsection we describe in detail the perturbative matching of the
QCD couplings. We follow the lines of the presentation in ref.
Athenodorou:2018wpk , and start by reformulating this matching in terms of
$\Lambda$-parameters Bernreuther:1981sg . After doing so, we investigate the
accuracy of a perturbative matching within perturbation theory itself, and
later present a discussion on the typical size of non-perturbative corrections
one can expect in these relations.

###  3.2 Perturbative decoupling

####  3.2.1 Definitions

As we are interested only in the QCD coupling we assume for simplicity that
the relevant effective theory is given by massless ${\rm QCD}_{{N_{\ell}}}$.
As discussed above, at leading order in $M^{-1}$ the only parameter of the
effective theory that needs to be fixed is therefore the running coupling
$\bar{g}^{({N_{\ell}})}(\mu)$. The effective theory hence predicts any
observable once the coupling in the chosen scheme is specified at some scale.

The general form of the relation between the couplings of the leading-order
effective theory $\bar{g}^{({N_{\ell}})}$ and of the fundamental theory
$\bar{g}^{({N_{\rm f}})}$, reads Weinberg:1980wa ; Bernreuther:1981sg ;
Bernreuther:1983zp

$[\bar{g}^{({N_{\ell}})}(\mu/\Lambda^{({N_{\ell}})})]^{2}=F_{\mathcal{O}}\big{(}[\bar{g}^{({N_{\rm
f}})}(\mu/\Lambda^{({N_{\rm f}})})]^{2},M/\mu\big{)}\,,$ (19)

where for later convenience we explicitly wrote the dependence of the
couplings on their corresponding $\Lambda$-parameter. The function
$F_{\mathcal{O}}$ depends in principle on the specific observable
${\mathcal{O}}$ that is used to establish the matching between the two
theories. The dependence on the observable, however, is suppressed by powers
of $M^{-1}$. In perturbation theory, these power corrections can be uniquely
isolated from the logarithmic terms in $M$ and can therefore be dropped. This
is consistent with matching the theories at leading order in $M^{-1}$. In
doing so, the relation between the couplings becomes universal, i.e
independent on the specific matching condition. It only depends on the
renormalization schemes chosen for the couplings.

In the $\overline{\rm MS}$-scheme the matching relation (also referred to as
decoupling relation) is known up to 4-loop order. Below, we consider this only
for the convenient choice of matching scale $\mu=m_{*}$, where $m_{*}$ is
implicitly defined by the equation $\overline{m}_{\overline{\rm
MS}}(m_{*})=m_{*}$, where $\overline{m}_{\overline{\rm
MS}}(\mu)\equiv\overline{m}^{({N_{\rm f}})}_{\overline{\rm
MS}}(\mu/\Lambda_{\overline{\rm MS}}^{({N_{\rm f}})})$ are the running masses
of the heavy quarks in the fundamental theory in the $\overline{\rm
MS}$-scheme. Given this choice, the matching relation reads:

$\displaystyle\big{[}\bar{g}^{({N_{\ell}})}_{\overline{\rm
MS}}(m_{*}/\Lambda^{({N_{\ell}})}_{\overline{\rm
MS}})\big{]}^{2}=g_{*}^{2}C(g_{*}),\,\quad g_{*}\equiv\bar{g}^{({N_{\rm
f}})}_{\overline{\rm MS}}(m_{*}/\Lambda^{({N_{\rm f}})}_{\overline{\rm
MS}})\,,$ $\displaystyle
C(g_{*})=1+h_{2}g_{*}^{4}+h_{3}g_{*}^{6}+h_{4}g_{*}^{8}+\ldots\,,$ (20)

where the $h_{i}$ are pure numbers that depend on both ${N_{\rm h}}$ and
${N_{\ell}}$ Grozin:2011nk ; Chetyrkin:2005ia ; Schroder:2005hy ;
Kniehl:2006bg ; Gerlach:2018hen . As explained in these references, the
particular choice of matching scale makes all contributions
$\propto\log(\overline{m}_{\overline{\rm MS}}(\mu)/\mu)$ appearing in the
matching relation vanish, and it is considered to be optimal. This implies in
particular that the $g_{*}^{2}$ term in $C(g_{*})$ is absent. In the general
case, the scales $\overline{m}_{\overline{\rm MS}}(\mu)$ and $\mu$ should
anyway not be chosen too separated. Large coefficients otherwise appear in the
perturbative matching relation which can compromise its applicability.

As anticipated, from the perspective of lattice applications it is compelling
to recast the matching relation, eq. (20), in terms of RG-invariant (RGI)
quantities. Specifically, this means the $\Lambda$-parameter of the effective
theory, $\Lambda^{({N_{\ell}})}_{\overline{\rm MS}}$, that of the fundamental
theory $\Lambda^{({N_{\rm f}})}_{\overline{\rm MS}}$, and the RGI quark-mass
$M$ of the heavy quarks. The latter can be defined as

$\displaystyle M$ $\displaystyle\equiv\overline{m}^{({N_{\rm
f}})}_{\overline{\rm MS}}(\mu)\varphi^{({N_{\rm f}})}_{\rm m,\overline{\rm
MS}}(\bar{g}^{({N_{\rm f}})}_{\overline{\rm MS}}(\mu))\,,$
$\displaystyle\varphi^{({N_{\rm f}})}_{\rm m,\overline{\rm MS}}(\bar{g})$
$\displaystyle\equiv(2b_{0}({N_{\rm f}})\bar{g}^{2})^{-{d_{0}\over
2b_{0}({N_{\rm f}})}}\times$
$\displaystyle\times\exp\bigg{\\{}-\int_{0}^{\bar{g}}{\rm
d}g\bigg{[}{\tau_{\overline{\rm MS}}^{({N_{\rm
f}})}({g})\over\beta_{\overline{\rm MS}}^{({N_{\rm f}})}({g})}-{d_{0}\over
b_{0}({N_{\rm f}})g}\bigg{\\}}\,,$ (21)

where $b_{0}({N_{\rm f}})$ is the 1-loop coefficient of the $\beta$-function,
eq. (7), and the function

$\tau^{({N_{\rm f}})}_{\overline{\rm
MS}}(\bar{g})\equiv{\mu\over\overline{m}^{({N_{\rm f}})}_{\overline{\rm
MS}}(\mu)}{{\rm d}\overline{m}^{({N_{\rm f}})}_{\overline{\rm
MS}}(\mu)\over{\rm d}\mu}\bigg{|}_{\bar{g}}$ (22)

encodes the scale-dependence of the quark masses in the ${N_{\rm f}}$-flavor
theory (in the $\overline{\rm MS}$-scheme). It has a perturbative expansion

$\tau^{({N_{\rm f}})}_{\overline{\rm MS}}(\bar{g})\overset{\bar{g}\to
0}{\approx}-\bar{g}^{2}\\{d_{0}+\bar{g}^{2}d_{1}+\ldots\\}\,,\quad
d_{0}=8/(4\pi)^{2}\,,$ (23)

known up to 5-loops Larin:1993tq ; Chetyrkin:1997dh ; Vermaseren:1997fq ;
Baikov:2014qja , where the actual scheme and ${N_{\rm f}}$-dependence start
from $d_{i}$, $i\geq 1$. Note that even though in eq. (21) we conveniently
defined the RGI mass $M$ through the $\overline{\rm MS}$-scheme, its value is
in fact scheme independent, as long as mass-independent schemes are considered
for the quark masses. This means in particular that $M$ can be non-
perturbatively defined through any non-perturbative (massless) renormalization
scheme.

Using the above definitions, together with the definitions in eq. (4), and the
matching relation eq. (20), it is immediate to conclude that Bruno:2014ufa ;
Athenodorou:2018wpk

$P_{\ell,\rm f}(M/\Lambda_{\overline{\rm MS}}^{({N_{\rm
f}})})\equiv{\Lambda_{\overline{\rm
MS}}^{({N_{\ell}})}\over\Lambda_{\overline{\rm MS}}^{({N_{\rm
f}})}}\Bigg{|}_{\rm matched}={\varphi_{\rm g,\overline{\rm
MS}}^{({N_{\ell}})}(g_{*}\sqrt{C(g_{*})})\over\varphi_{\rm g,\overline{\rm
MS}}^{({N_{\rm f}})}(g_{*})}\,.$ (24)

As anticipated by our notation, $P_{\ell,\rm f}$ can be considered as a
function of $M/\Lambda_{\overline{\rm MS}}^{({N_{\rm f}})}$, since the value
of the coupling $g_{*}$ can be expressed in terms of the RGI parameters
through the relation

${M\over\Lambda_{\overline{\rm MS}}^{({N_{\rm f}})}}={\varphi_{\rm
m,\overline{\rm MS}}^{({N_{\rm f}})}(g_{*})\over\varphi_{\rm g,\overline{\rm
MS}}^{({N_{\rm f}})}(g_{*})}\,,\qquad g_{*}=\bar{g}^{({N_{\rm
f}})}_{\overline{\rm MS}}({M/\Lambda_{\overline{\rm MS}}^{({N_{\rm f}})}}).$
(25)

In the limit where $M/\Lambda^{({N_{\rm f}})}_{\overline{\rm MS}}\to\infty$,
the coupling $g_{*}$ goes to zero and the function $P_{\ell,\rm f}$ admits an
asymptotic perturbative expansion in terms of $g_{*}^{2}$.

####  3.2.2 Perturbative uncertainties

We now want to study the behavior of the perturbative expansion of
$P_{\ell,\rm f}(M/\Lambda)$. We shall focus our attention on the functions
$P_{\rm 3,4}$ and $P_{\rm 4,5}$, which are the relevant ones for estimating
$\Lambda^{(5)}_{\overline{\rm MS}}$ given $\Lambda^{(3)}_{\overline{\rm MS}}$.

In fig. 7 we present the results taken from ref. Athenodorou:2018wpk  for the
relative deviation

${(P_{\ell,\rm f}-P_{\ell,\rm f}^{(1)})\over P^{(1)}_{\ell,\rm f}}$ (26)

of $P_{\rm 3,4}$ (left plot) and $P_{\rm 4,5}$ (right plot), from their
_unsystematic_ 1-loop approximation $P^{(1)}_{\ell,\rm
f}(M/\Lambda)=(M/\Lambda)^{\eta_{0}}$, where $\eta_{0}=2{N_{\rm
h}}/(33-2{N_{\ell}})$ Athenodorou:2018wpk .242424The unsystematic
approximation $P^{(1)}_{\ell,\rm f}$ obtained in ref. Athenodorou:2018wpk
happens to be quite close to the different orders of the perturbative
expansion of $P_{\ell,\rm f}$ for the values of ${N_{\rm f}}$ and ${N_{\ell}}$
considered and $M/\Lambda\lesssim 30$. For this reason it is used here to get
the overall magnitude of the function $P_{\ell,\rm f}$. On the other hand, the
different orders of perturbation theory are not expected to converge to
$P^{(1)}_{\ell,\rm f}$ as $M/\Lambda\to\infty$ since this approximation does
not have the correct asymptotic limit Athenodorou:2018wpk . The results are
shown as a function of $M/\Lambda$ and for different orders of perturbation
theory. The order in perturbation theory refers to the order at which the
$\beta$-functions entering eqs. (24)-(25) are considered. The $\tau$-function,
eq. (23), as well as the matching function $C(g_{*})$ of eq. (20) are
considered to a consistent order in the expansion (see ref.
Athenodorou:2018wpk  for the details). The values of the RGI masses of the
charm and bottom quarks ($M_{c}$ and $M_{b}$, respectively) in units of the
relevant $\Lambda$-parameter are also shown. These are inferred from the
results of the PDG Zyla:2020zbs .

Figure 7: Relative differences from the “1-loop approximation”
$P^{(1)}_{\ell,\rm f}(M/\Lambda)=(M/\Lambda)^{\eta_{0}}$, $\eta_{0}=2{N_{\rm
h}}/(33-2{N_{\ell}})$, for different orders of the perturbative expansion of
$P_{\ell,\rm f}(M/\Lambda)$ as a function of $M/\Lambda$ Athenodorou:2018wpk .
The results for ${N_{\ell}}=3$, ${N_{\rm f}}=4$ (${N_{\ell}}=4$, ${N_{\rm
f}}=5$) are given in the left (right) panel. The values for the RGI charm
($M_{c}$) and bottom ($M_{b}$) quark masses in units of the proper
$\Lambda$-parameters are marked by vertical lines.

Starting from the case of $P_{4,5}$, we see how the difference between the
3-loop and 2-loop results is around 2% at the $b$-quark mass. The 4- and
5-loop results then differ only by very tiny corrections from the 3-loop
approximation. Looking at the behavior of the perturbative series alone, the
perturbative prediction for the decoupling of the $b$-quark appears to be very
reliable and accurate. Similar conclusions can be drawn for the decoupling of
the charm quark. Also in the case of $P_{3,4}$, the difference between the 3-
and 2-loop results is about 2% at the charm-quark mass, and higher-order
corrections are all much smaller.

In conclusion, the perturbative results suggest that a perturbative treatment
of the matching between the relevant effective and fundamental theories
introduces only errors at the sub-percent level in the functions $P_{\ell,\rm
f}$, and therefore in the connection of their $\Lambda$-parameters.

###  3.3 Non-perturbative decoupling

####  3.3.1 Definitions

Judging from perturbation theory alone, the perturbative description of the
decoupling of heavy quarks seems to work remarkably well. The rapid
convergence of the different orders of the expansion of $P_{\ell,\rm f}$ at
the values of both the charm- and bottom-quark masses, seems to suggest that
the series is well within its regime of applicability and higher-order
corrections are small. On the other hand, the perturbative expansion cannot
tell us anything about the size of non-perturbative corrections to the
decoupling relations and whether perturbation theory actually applies at all.
Whether non-perturbative effects are significant within the target precision
at the relevant quark masses can only be established through a non-
perturbative investigation.

In order to set the grounds for estimating the error that one makes when using
a perturbative approximation for $P_{\ell,\rm f}$ to extract
$\Lambda_{\overline{\rm MS}}^{({N_{\rm f}})}$ from $\Lambda_{\overline{\rm
MS}}^{({N_{\ell}})}$, let us begin by recasting the matching of the effective
and fundamental theory in more non-perturbative terms.252525For simplicity we
consider also here the case where the light quarks are massless.

The leading-order effective theory describes the fundamental theory at low
energy only when the corresponding $\Lambda$-parameter
$\Lambda^{({N_{\ell}})}$ is a properly chosen function of the scale
$\Lambda^{({N_{\rm f}})}$ of the fundamental theory and of the RGI mass $M$ of
the heavy quarks. To be more precise, consider a low-energy mass-scale
$\mathcal{S}\ll M$. This could be, for instance, a hadronic mass or any of the
popular technical scales, $t^{-1/2}_{0},r^{-1}_{0},w^{-1}_{0}$ Luscher:2010iy
; Sommer:1993ce ; Borsanyi:2012zs . Matching then means fixing the scale
$\Lambda^{({N_{\ell}})}$ through the condition Athenodorou:2018wpk

${\Lambda^{({N_{\ell}})}\over\mathcal{S}^{({N_{\ell}})}}=P_{\ell,\rm
f}^{\mathcal{S}}(M/\Lambda^{({N_{\rm f}})}){\Lambda^{({N_{\rm
f}})}\over\mathcal{S}^{({N_{\rm f}})}(M)}\,.$ (27)

In these equations $\mathcal{S}^{({N_{\ell}})}$ and $\mathcal{S}^{({N_{\rm
f}})}(M)$ refer to the given low-energy scale computed in the effective and
fundamental theory, respectively. Note in particular that the matching
function $P_{\ell,\rm f}^{\mathcal{S}}(M/\Lambda^{({N_{\rm f}})})$ depends on
the scale $\mathcal{S}$ that is considered in the matching relation. We also
stress again that while the value of $\Lambda^{({N_{\rm f}})}$ does not depend
on $M$ the one of $\Lambda^{({N_{\ell}})}$ does depend through the matching
condition.

Once $\Lambda^{({N_{\ell}})}$ is properly fixed through eq. (27) in terms of
$\mathcal{S}$, for any other low-energy quantity $\mathcal{S}^{\prime}$ we
expect that

$\mathcal{S}^{\prime({N_{\ell}})}=\mathcal{S}^{\prime({N_{\rm f}})}(M)+{\rm
O}\big{(}\big{(}{\Lambda^{({N_{\rm f}})}/M}\big{)}^{2}\big{)}\,.$ (28)

Note that ratios of low-energy scales, instead, do not depend on the value of
the $\Lambda$-parameters and are therefore insensitive to their matching. For
these, it readily holds that

${\mathcal{S}^{({N_{\ell}})}\over\mathcal{S}^{\prime({N_{\ell}})}}={\mathcal{S}^{({N_{\rm
f}})}(M)\over\mathcal{S}^{\prime({N_{\rm f}})}(M)}+{\rm
O}\big{(}\big{(}{\Lambda^{({N_{\rm f}})}/M}\big{)}^{2}\big{)}\,,$ (29)

with $\mathcal{S}$ and $\mathcal{S}^{\prime}$ any two low-energy scales. Given
this observation, from eq. (27) we conclude that

$P_{\ell,\rm f}^{\mathcal{S}}=P_{\ell,\rm f}^{\mathcal{S}^{\prime}}+{\rm
O}(({\Lambda^{({N_{\rm f}})}/M})^{2})\,.$ (30)

In other words, at the non-perturbative level the function $P_{\ell,\rm
f}^{\mathcal{S}}$ intrinsically comes with O($M^{-2}$) ambiguities. For this
reason we shall often simply write it as $P_{\ell,\rm f}$, keeping the
dependence on $\mathcal{S}$ and the related ambiguities understood.

An interesting consequence of the above relations follows from multiplying eq.
(27) by $\mathcal{S}^{({N_{\rm f}})}(0)/\Lambda^{({N_{\rm f}})}$, where
$\mathcal{S}^{({N_{\rm f}})}(0)$ stands for the low-energy quantity
$\mathcal{S}$ computed in the chiral limit of the ${N_{\rm f}}$-flavor theory.
Through this simple manipulation we find that Bruno:2014ufa ;
Athenodorou:2018wpk

$\displaystyle{\mathcal{S}^{({N_{\rm f}})}(M)\over\mathcal{S}^{({N_{\rm
f}})}(0)}$ $\displaystyle=Q^{\mathcal{S}}_{\ell,\rm f}\times
P^{\mathcal{S}}_{\ell,\rm f}(M/\Lambda^{({N_{\rm f}})})$
$\displaystyle=Q^{\mathcal{S}}_{\ell,\rm f}\times P_{\ell,\rm
f}(M/\Lambda^{({N_{\rm f}})})+{\rm O}\big{(}\big{(}{\Lambda^{({N_{\rm
f}})}/M}\big{)}^{2}\big{)}\,,$ (31)

where the factor

$Q^{\mathcal{S}}_{\ell,\rm
f}\equiv{\mathcal{S}^{({N_{\ell}})}/\Lambda^{({N_{\ell}})}\over\mathcal{S}^{({N_{\rm
f}})}(0)/\Lambda^{({N_{\rm f}})}}\,,$ (32)

is defined in terms of the massless ${N_{\rm f}}$\- and ${N_{\ell}}$-flavor
theories. The interesting aspect of eq. (31) is that the ratio on the l.h.s.
can be computed within the fundamental theory, while the r.h.s. is a
consequence of the decoupling of the heavy quarks. In particular, we see how
the mass dependence of the ratio ${\mathcal{S}^{({N_{\rm
f}})}(M)/\mathcal{S}^{({N_{\rm f}})}(0)}$ is expressed in terms of the
function $P^{\mathcal{S}}_{\ell,\rm f}$, while the factor
$Q^{\mathcal{S}}_{\ell,\rm f}$ is just an overall constant. In the limit of
large mass $M$, the mass dependence of this ratio is therefore expected to be
universal and described by perturbation theory. As a result, this relation
allows us to put at test the perturbative expansion of $P_{\ell,\rm f}$ and
estimate the typical size of non-perturbative corrections.

To this end, it is convenient in practice to introduce the mass-scaling
function Athenodorou:2018wpk

$\eta_{\ell,\rm f}^{M,\mathcal{S}}(M)\equiv{M\over P^{\mathcal{S}}_{\ell,\rm
f}}{\partial P^{\mathcal{S}}_{\ell,\rm f}\over\partial
M}\bigg{|}_{\Lambda^{({N_{\rm f}})}}\,.$ (33)

Considering eq. (31), this can be computed from the mass dependence of any
hadronic quantity as262626Note that the number of light quarks on the r.h.s.
of this equation is implicitly given by the difference between ${N_{\rm f}}$
and the number ${N_{\rm h}}$ of heavy quarks of which the mass is varied.

$\eta_{\ell,\rm f}^{M,\mathcal{S}}(M)={M\over\mathcal{S}^{({N_{\rm
f}})}(M)}{\partial{\mathcal{S}}^{({N_{\rm f}})}(M)\over\partial
M}\bigg{|}_{\Lambda^{({N_{\rm f}})}}\,,$ (34)

with no need for determining $\mathcal{S}^{({N_{\rm f}})}(0)$ or $Q_{\ell,\rm
f}$. As for the functions $P^{\mathcal{S}}_{\ell,\rm f}$, the $\eta_{\ell,\rm
f}^{M,\mathcal{S}}$ obtained from different low-energy quantities
$\mathcal{S}$ are expected to differ by O($(\Lambda^{({N_{\rm f}})}/M)^{2}$)
contributions. In the limit $M/\Lambda^{({N_{\rm f}})}\to\infty$, however, the
mass-scaling function becomes universal and can asymptotically be estimated in
perturbation theory Athenodorou:2018wpk . In the following we shall see how by
studying the mass-scaling function we will be able to obtain valuable insight
on the applicability of perturbation theory in computing $P_{\ell,\rm f}$.

####  3.3.2 Non-perturbative corrections to decoupling

Non-perturbative effects in the decoupling of heavy quarks have been
systematically investigated in a series of recent papers Bruno:2014ufa ;
Athenodorou:2018wpk ; Knechtli:2017xgy ; Cali:2017brl ; Cali:2019enm ;
Hollwieser:2020qri , which focus on the relevant case of the charm. The main
question that these works contribute to answer is how good of an approximation
${N_{\rm f}}=2+1$ QCD is to the ${N_{\rm f}}=2+1+1$ flavor theory. From our
perspective, we are particularly interested in understanding how precisely we
can expect to obtain $\Lambda_{\overline{\rm MS}}^{(4)}$ (and consequently
$\Lambda_{\overline{\rm MS}}^{(5)}$) from results in ${N_{\rm f}}=2+1$ QCD. As
already mentioned, the first issue is to understand how accurately we can
estimate $\Lambda_{\overline{\rm MS}}^{(4)}$ from $\Lambda_{\overline{\rm
MS}}^{(3)}$ by relying on a perturbative approximation for
$P_{3,4}(M_{c}/\Lambda_{\overline{\rm MS}}^{(4)})$. Secondly, we must question
how much we can rely on setting the physical units of the theory using
${N_{\rm f}}=2+1$ QCD.

Studying the decoupling of the charm quark through simulations of ${N_{\rm
f}}=2+1+1$ and ${N_{\rm f}}=2+1$ QCD is very challenging from the
computational point of view. As a result, the physical effects one is after
can easily end up being masked by the final uncertainties. For this reason,
the authors of refs. Bruno:2014ufa ; Athenodorou:2018wpk ; Knechtli:2017xgy
have investigated non-perturbatively a model system given by QCD with ${N_{\rm
f}}=2$ degenerate heavy quarks. When the mass of the doublet of quarks becomes
large the heavy quarks eventually decouple, and the theory is expected to be
described at low energy by the pure Yang-Mills theory, i.e. ${N_{\ell}}=0$
flavor QCD. Studying this model rather than the realistic case, avoids the
complications of simulating light quarks. This allows one to reach much finer
lattice spacings than typically possible in large-volume simulations with
light quarks, because the volumes are similar to those in pure-gauge theory.
Fine lattice spacings are essential to have discretization errors under
control in the presence of quarks with masses close to that of the charm.

By studying this model one expects to be able to reliably estimate the typical
size of the effects induced by the charm in low-energy physics. In particular,
the absence of the light quarks is not expected to change the picture very
much and their effect is likely more than compensated by the additional heavy
quark present in the model. In fact, as we shall report below, some first
results have been recently obtained for the more realistic situation where a
single charm quark decouples in the presence of 3 mass-degenerate lighter
quarks Hollwieser:2020qri . The results collected so far confirm the findings
of the model study.

Figure 8: Mass dependence of the mass-scaling functions
$\eta^{M,\mathcal{S}}_{0,2}$ obtained from the hadronic scales
$\mathcal{S}=t_{0}^{-1/2},t_{c}^{-1/2},w_{0}^{-1}$ Athenodorou:2018wpk . The
data for a given mass $M$ are slightly displaced horizontally for better
clarity. The non-perturbative results are compared to the perturbative
estimates at 1- and 4-loop order. The dash-dotted lines are fits to the data
for $\mathcal{S}=t_{c}^{-1/2}$ and $w_{0}^{-1}$ used to estimate the size of
the non-perturbative corrections to $\eta_{0,2}^{M,\mathcal{S}}$
Athenodorou:2018wpk . The vertical dotted lines mark the values of the quark
mass $M_{c},M_{c}/2$, and $M_{c}/4$.

##### Ratios of $\Lambda$-parameters

Figure 9: Continuum extrapolated values for $\sqrt{t_{c}/t_{0}}$ (left) and
$\sqrt{t_{0}}/w_{0}$ (right) as a function of $\Lambda/M$ Knechtli:2017xgy .
The line in the blue band is the leading-order effective theory prediction:
$R(M)=R(\infty)+k_{R}\times\Lambda^{2}/M^{2}$,
$R=\sqrt{t_{c}/t_{0}},\sqrt{t_{0}}/w_{0}$, with $k_{R}$ a free parameter,
fitted through points from $M=\infty$ down to $M/\Lambda=2.5$. The line in the
green band is instead a fit linear in $\Lambda/M$. For comparison the dashed
lines represent the quadratic (blue) and linear (green) fit through points
from $M=\infty$ down to $M/\Lambda=1.28$. Also shown by the dashed-dotted red
line is a fit in this range adding to the leading-order prediction a next-to-
leading term proportional to $\Lambda^{4}/M^{4}$.

The first set of results that we want to discuss are from ref.
Athenodorou:2018wpk  and are shown in figure 8. They correspond to the
determination of the mass-scaling functions $\eta^{M,\mathcal{S}}_{0,2}$ based
on the gluonic low-energy scales
$\mathcal{S}=t_{0}^{-1/2},t_{c}^{-1/2},w_{0}^{-1}$ (cf. eq. (34)). The precise
definition of these scales can be found in the given reference. The results
refer to the continuum limit of the model system introduced above, i.e.
${N_{\rm f}}=2$ QCD with two heavy quarks. The RGI mass $M$ of the heavy
quarks is varied from about $M_{c}/8$ to $1.2M_{c}$, where $M_{c}$ is the RGI
mass of the charm quark.272727The value of the charm-quark mass in simulations
is set by targeting $M_{c}/\Lambda_{\overline{\rm MS}}^{(2)}=4.87$, which is
obtained using $\Lambda_{\overline{\rm MS}}^{(2)}=310\,{\rm MeV}$ from ref.
Fritzsch:2012wq  and $M_{c}=1510\,{\rm MeV}$ from ref. Heitger:2013oaa  (cf.
ref. Athenodorou:2018wpk  for the details). The plot also includes the results
for $\eta^{M}_{0,2}$ at 1- and 4-loop order in perturbation theory
Athenodorou:2018wpk .

As one can see from the figure, the results for the mass-scaling functions
corresponding to different low-energy observables significantly differ at the
smaller values of $M$ in the range. As expected, however, they consistently
approach each other as $M\to\infty$. In particular, for values of the mass
close to the charm-quark mass all determinations well agree within errors,
indicating the smallness of the O($M_{c}^{-2}$) corrections.

The results for $\eta^{M}_{0,2}$ are about $1/10$ for $M\approx M_{c}$, both
in perturbation theory and non-perturbatively. In fact, within the
uncertainties of roughly 10% the perturbative and non-perturbative results
perfectly agree. Note that even though the relative precision on
$\eta^{M}_{0,2}$ might not seem impressive, it corresponds to an absolute
error on $\eta^{M}_{0,2}$ of about $0.01$. This means, in particular, that it
is reasonable to assume that the difference, $\Delta\eta^{M}_{0,2}$, between
the non-perturbative $\eta^{M}_{0,2}$ and its (4-loop) perturbative
approximation is bounded by this error for masses $M\gtrsim M_{c}$. The
scaling of the non-perturbative corrections $\Delta\eta^{M}_{0,2}$ to
$\eta^{M}_{0,2}$ as a function of $M$ can then be assessed by studying the
$M$-dependence seen in figure 8 and comparing the results for the different
observables $\mathcal{S}$ Athenodorou:2018wpk .

Putting all this information together, the authors of ref. Athenodorou:2018wpk
obtain a conservative estimate for the size of the non-perturbative
contributions $\Delta P_{0,2}$ to $P_{0,2}$. We shall not report their
detailed discussion here and refer the interested reader to ref.
Athenodorou:2018wpk .282828The deviation $\Delta\log(P_{0,2})$ of the full
function $\log(P_{0,2})$ from its perturbative approximation is obtained by
integrating $\Delta\eta^{M}_{0,2}$ in $\log(M/\Lambda)$ from the value of
$M/\Lambda$ of interest up to $M/\Lambda=\infty$ (cf. eq. (34) and ref.
Athenodorou:2018wpk ). In order to translate to $P_{0,2}$ the estimate for the
non-perturbative corrections to $\eta^{M}_{0,2}$ made for $M\lesssim M_{c}$,
one therefore needs to make some assumptions on how these approach the limit
$M/\Lambda\to\infty$. Depending on the assumptions, more or less conservative
estimates are obtained. Summarizing their conclusions, given the results for
$\Delta\eta^{M}_{0,2}$ obtained in the model, one can safely state that the
non-perturbative contributions to $P_{0,2}(M_{c}/\Lambda_{\overline{\rm
MS}}^{(2)})$ are _at most_ 2% and quite likely at the level of 0.4%. This
translates into at least a 2% precision of perturbation theory in the
conversion of $\Lambda$-parameters for the investigated case.