\section{Introduction}
\label{sec:intro}

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.~\cite{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.~\cite{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
$\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 $\alphas$
(see e.g.~refs.~\cite{Aoki:2019cca,dEnterria:2019its,Zyla:2020zbs}).%
\footnote{We here adopt the common notation
		  $\alphas(\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.~\cite{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
$\alphas(M_Z)=0.1179(10)$~\cite{Zyla:2020zbs}, and has a precision
of about 0.8\%. The lattice determinations alone give
$\alphas(M_Z)=0.1182(8)$~\cite{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~\cite{Aoki:2019cca}. It is fair to say that such a
variety of approaches within a PDG subcategory is in fact
unique~\cite{Zyla:2020zbs}.

Despite the tremendous efforts on and off the lattice, however,
the current uncertainty on $\alphas$ 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.~\cite{Salam:2017qdl,dEnterria:2019its}).
An uncertainty on $\alphas(M_Z)$ comfortably below the percent
level is desired for precision  applications. For these reasons,
there are  plans for future phenomenological determinations of
$\alphas(M_Z)$ aiming at reaching the extremely competitive
accuracy of 0.2\% using high-luminosity high-energy data
(see
e.g.~refs.~\cite{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
$\alphas$ by such an important factor is not easy.
Similarly to several phenomenological determinations most lattice
determinations of $\alphas$ are currently limited by systematic
uncertainties related to the use of perturbation theory
at relatively low scales~\cite{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~\cite{Luscher:1991wu,Jansen:1995ck}.
The methods have been recently applied for obtaining one of the
most precise determinations of $\alphas$~\cite{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 $\alphas$ by an
important factor is yet a computationally expensive task,
even employing a step-scaling strategy
(cf.~ref.~\cite{Bruno:2017gxd}).

In this contribution we want to review the recent proposal made
in ref.~\cite{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 $\Nf$-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 $\Nf$-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
$\Nf$ 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^{(\Nf=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 $\Nf$-flavor and pure Yang-Mills theory, that is,
$\Lambda^{(\Nf=0)}_{\overline{\rm MS}}/
\Lambda^{(\Nf)}_{\overline{\rm MS}}$.
Putting these steps together, one obtains
$\Lambda^{(\Nf)}_{\overline{\rm MS}}/\mu_{\rm dec}$,
and given the physical value of $\mu_{\rm dec}$ finds
$\Lambda^{(\Nf)}_{\overline{\rm MS}}$. Considering
$\Nf=3$ or 4, once $\Lambda^{(\Nf)}_{\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^{(\Nf=5)}_{\overline{\rm MS}}$
and from it $\alphas(M_Z)$.

The strategy has already been proven successful in the
determination of
$\Lambda^{(\Nf=3)}_{\overline{\rm MS}}$~\cite{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.~\ref{sec:PrecisionDeterminations} 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.~\ref{sec:Decoupling}, 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 $\Nf$-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 $\alphas$ can be obtained from
results in $\Nf=3$ QCD.

In Sect.~\ref{sec:RenByDec}, the application of heavy-quark
decoupling to the determination of the $\Nf$-flavor QCD coupling
is described in detail and the results of
ref.~\cite{DallaBrida:2019mqg} for
$\Lambda^{(\Nf=3)}_{\overline{\rm MS}}$ are presented.
We conclude in Sect.~\ref{sec:Conclusions} with some
comments on the future prospects for $\alphas$ 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~\cite{Aoki:2019cca} and to other
interesting recent reviews (see
e.g.~refs.~\cite{DelDebbio:2021ryq,Komijani:2020kst}).

\section{Precision determinations: the case of $\alphas$}
\label{sec:PrecisionDeterminations}

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 $\alphas$. 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.

\subsection{Determinations of $\alphas$ 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 $\alphas$ (see
refs.~\cite{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.~\cite{Sommer:2015kza,Salam:2017qdl,DelDebbio:2021ryq} for
some reviews).

Any determination of $\alphas$, whether on the
lattice or not, relies on the following basic strategy.
One considers a short-distance observable $\Obs(q)$
which depends on a characteristic energy scale $q$.
In the limit where $q\to\infty$, this observable is
compared with its theoretical prediction, $\Obs_{\rm th}(q)$,
in terms of a perturbative expansion%
\footnote{Note that in general the coupling $\alphas$
		  to be considered here should be the QCD coupling of
		  the relevant $\Nf$-flavor theory, i.e.
		  $\alphas(\mu)\equiv
		  \alpha_{\overline{\rm MS}}^{(\Nf)}(\mu)$, from
		  which $\alpha_{\overline{\rm MS}}^{(5)}(M_Z)$
		  can eventually be extracted
		  (cf.~e.g.~ref.~\cite{Aoki:2019cca} and
		   Sect.~\ref{subsec:EffectiveCouplings}). At this stage,
		  however, we prefer to keep the discussion simple.
		  We shall return to this point later in detail.}
\begin{equation}
	\label{eq:Oth}
	\Obs_{\rm th}(q)=\sum_{n=0}^N k_n(s)\alphas^n(\mu)
	+{\rm O}(\alphas^{N+1})+{\rm O}\bigg({\Lambda^p\over q^p}\bigg),
	\
	\mu={q\over s}\,.
\end{equation}
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 $\alphas$ is extracted with the scale $q$.
The basic difference between phenomenological and lattice
determinations of $\alphas$ is the choice of the observable
$\Obs(q)$.

Requiring $\Obs_{\rm th}(q)=\Obs(q)$ for some finite $q$, clearly
fixes the value of $\alphas(\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 \Obs(q)$ to which the observable
$\Obs(q)$ is known. This of course depends on the relevant
statistical and systematic uncertainties associated with the
determination of $\Obs(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}(\alphas^{N+1})$ terms in
eq.~(\ref{eq:Oth}). In addition to these there are contaminations
from ``non-perturbative contributions''. These are represented in
eq.~(\ref{eq:Oth}) 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.%
\footnote{Our 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.~\cite{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 $\alphas$ needs to have, \emph{at least},
these general sources of error under control. Note that for
the most part these are systematic in nature.

Lattice determinations of $\alphas$ 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 $\Obs(q)$ as the getaway to extract $\alphas$.
One can therefore \emph{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 $\alphas$, instead, rely on experimental data for the
observable $\Obs(q)$. It is the typical situation that when
$q$ becomes large, and therefore many sources of systematic
uncertainty in eq.~(\ref{eq:Oth}) become small, the experimental
errors $\delta\Obs(q)$ become large. It is thus difficult to find
in general a single experimental quantity $\Obs(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
\emph{carefully} chosen, $\Obs(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.~(\ref{eq:Oth}).

Another advantage for the lattice theorist is that $\Obs$ is
defined within QCD alone. Consequently, the theoretical
description $\Obs_{\rm th}$ of eq.~(\ref{eq:Oth}) does
not need to include contributions besides those from QCD.
In addition, no modeling of hadronization is needed when
comparing the observable $\Obs$ with its perturbative prediction
$\Obs_{\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
$\alphas$ 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 $\alphas$ can
therefore be limited by these factors
\footnote{Experimental observables are in principle sensitive also
	to any New Physics. How this affects the extracted value
	of $\alphas$, 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 $\Obs$ itself. As mentioned earlier,
the uncertainties on the hadronic quantities have, at present,
limited impact on the results for $\alphas$.

All current lattice QCD determinations of $\alphas$,
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.~\ref{sec:Decoupling}).

\subsection{The challenge of reaching high energy}
\label{subsec:ChallengesLattice}

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

As any lattice QCD observable, besides the statistical
uncertainties, $\Obs(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 $\Obs$ 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~\cite{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 $\alphas$ is extracted,
then one is necessarily limited by the simultaneous constraints:
\begin{equation}
	\label{eq:Window}
	L^{-1}\ll \Lambda
	\quad
	\textrm{and}
	\quad
	\Lambda\ll q\ll a^{-1}\,.
\end{equation}
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 $\Obs(q)$
are kept under control, and $\Obs(q)$ can be obtained in
the continuum limit with controlled errors.%
\footnote{Here and in the following we assume that
		  $\Obs(q)$ is a properly renormalized lattice quantity
		  with a well-defined continuum limit and free from
		  infrared divergences. Alternative strategies
		  extract $\alphas$ from bare lattice quantities by
		  taking $q\propto a^{-1}$ and expressing their
		  expansion in lattice perturbation theory in terms of
		  $\alphas$ (see e.g.~ref.~\cite{Mason:2005zx}). In these
		  cases, discretization errors and scale dependence of
		  $\Obs(q)$ are entangled. Addressing systematic
		  uncertainties becomes more subtle and requires a separate
		  discussion (see ref.~\cite{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 $\alphas$ by comparing
$\Obs(q)$ with eq.~(\ref{eq:Oth}).

The typical lattices that are simulated today have sizes
$L/a\lesssim100$. Taking for definiteness $m_\pi L=4$, with
$m_\pi=140\,\MeV$, the first condition in eq.~(\ref{eq:Window})
implies that for such ensembles $q \ll a^{-1}\approx 3.5\,\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 $\alphas$ of ref.~\cite{Cali:2020hrj}.
For their computation the authors employ state-of-the-art
large-volume simulations by the CLS
initiative~\cite{Bruno:2014jqa,Bali:2016umi}.
The two-point functions of both axial and vector currents at
short-distances are used to extract $\alphas$. Given the range of
lattice spacings at their disposal, $a\approx0.04-0.08\,\fm$, the
accessible energies for which continuum limit extrapolations
could be performed in a controlled way are limited to
$q\lesssim 2\,\GeV$.

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

\subsubsection{The finite volume is your friend}
\label{subsubsec:FV}

In order to reach high precision we must tailor the
lattice simulations to the problem at hand. Going back to
eq.~(\ref{eq:Window}), 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 $\Obs(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 \emph{finite-volume
renormalization schemes}~\cite{Luscher:1991wu}. In this case,
the scale $q$ at which the observable $\Obs(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 $\Obs$
considered is a \emph{finite-volume effect}. With this choice,
one computes the non-perturbative RG running of $\Obs(L^{-1})$
by simulating lattices with different physical extent $L$.
This strategy goes under the name of \emph{finite-size
(or step-)scaling}~\cite{Luscher:1991wu,Jansen:1995ck}
(see ref.~\cite{Sommer:2015kza} for a recent account).

More precisely, having fixed the bare QCD parameters through
some hadronic quantities, one computes $\Obs(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\to0}(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\gg1$
and $am_{\rm had}\ll1$, as well as having finite-volume effects
in $m_{\rm had}$ under control.

Secondly, one computes in the continuum limit the change in
$\Obs(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 $\alphas(\mu_{\rm PT}=q_n/s)$ from the value of
$\Obs(q_n)$ (cf.~eq.~(\ref{eq:Oth})).

It is important to emphasize that, if carefully chosen,
the only source of systematic errors that affect the
determination of $\Obs(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\gg1$. This approach elegantly exploits the freedom
that we have in lattice QCD in choosing the observable $\Obs(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 $\Obs(q)$. For the strategy to be
successful in practice, however, this should have a number of
desirable properties (see e.g.~ref.~\cite{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 $\alphas$ through eq.~(\ref{eq:Oth}).
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.

\subsubsection{$\Lambda$-parameters, $\beta$-functions, and all that}
\label{subsubsec:LambdaParams}

As we leave the general discussion for entering a more
quantitative analysis of the challenges of extracting
$\alphas$, we find convenient to reformulate the problem
in slightly different terms. First of all, it is useful
to associate with the observable $\Obs$ a renormalized
coupling $\bar{g}_\Obs$ through the relation:
\begin{equation}
	\label{eq:alphaO}
	\alpha_\Obs(\mu)\equiv{\bar{g}^2_\Obs(\mu)\over4\pi} \equiv
	{\Obs(\mu) -k_0\over k_1} \overset{\mu\to\infty}{=}
	\alphas(\mu) + {k_2\over k_1}\alphas^{2}(\mu)+\ldots\,,
\end{equation}
where the coefficients $k_i\equiv k_i(1)$ are those appearing
in the perturbative expansion, eq.~(\ref{eq:Oth}). This simple
procedure defines a non-perturbative, regularization-independent,
QCD coupling. In terms of these couplings, the extraction of
$\alphas$  is interpreted as the perturbative matching between
$\alpha_{\Obs}$ and $\alphas$, where different observables $\Obs$
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_{\Obs}$ 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, $\Nf$. A more proper notation to use
for the couplings is therefore $\alpha_{\Obs}^{(\Nf)}(\mu)\equiv
[\bar{g}^{(\Nf)}_\Obs(\mu)]^2/(4\pi)$, which emphasizes the fact
that $\bar{g}^{(\Nf)}_\Obs(\mu)$ must be considered as a coupling
within the $\Nf$-flavor theory. For ease of notation, however,
in this section we will often take the liberty of dropping the
superscripts $\Nf$ 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 \emph{mass-independent} (or simply massless)
renormalization schemes~\cite{Weinberg:1951ss}. These are
defined in terms of observables $\Obs$ 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~\cite{Bernreuther:1981sg,Bernreuther:1983zp}.
Hence, in the $\Nf$-flavor theory the latter is characterized by
a fixed number of active flavors corresponding to $\Nf$.%
\footnote{We shall return on these points in more detail
		  in Sect.~\ref{sec:Decoupling}.}

To the coupling $\bar{g}^{(\Nf)}_\Obs$ in a given
mass-independent scheme we can associate a quark-mass independent
$\Lambda$-parameter, $\Lambda^{(\Nf)}_\Obs$, defined as,
\begin{equation}
	\label{eq:LambdaParam}
	\begin{split}
	\Lambda^{(\Nf)}_\Obs&=\mu\,\varphi^{(\Nf)}_{\rm g,\Obs}
	(\bar{g}^{(\Nf)}_\Obs(\mu))\,, \\[1ex]
	\varphi^{(\Nf)}_{\rm g,\Obs}(\bar{g})&=
	(b_0(\Nf)\bar{g}^2)^{-{b_1(\Nf)\over 2b_0(\Nf)^2}}\,
	e^{-{1\over 2b_0(\Nf)\bar{g}^2}}\times\\[1ex]
	&\hspace*{-10.5mm}\times
	\exp\bigg\{
	-\int_{0}^{\bar{g}}
	\rmd g \bigg[{1\over\beta^{(\Nf)}_\Obs(g) } +
	{1\over b_0(\Nf) g^3} - {b_1(\Nf)\over b_0(\Nf)^2 g}\bigg]\bigg\}\,.
	\end{split}
\end{equation}
Through this relation, the value of the coupling
$\bar{g}^{(\Nf)}_\Obs(\mu)$ at any renormalization scale
$\mu$ is in one-to-one correspondence with $\Lambda^{(\Nf)}_\Obs$,
provided that the $\beta$-function,
\begin{equation}
	\label{eq:BetaFunDef}
	\beta^{(\Nf)}_\Obs(\bar{g})
	\equiv\mu{\rmd\bar{g}^{(\Nf)}_\Obs(\mu)\over\rmd \mu}\bigg|_{\bar{g}}\,,
\end{equation}
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:
\begin{equation}
	\label{eq:BetaFunctionPT}
	\beta_\Obs^{(\Nf)\,{\rm PT}}(\bar{g})\equiv
	-\bar{g}^3\sum_{k=0}^{N-1} b_k(\Nf)\,\bar{g}^{2k}\,,
\end{equation}
with
\begin{gather}
	\nonumber
	b_0(\Nf)={1\over(4\pi)^{2}}\bigg(11-{2\over 3}\Nf\bigg)\,,\\
	\label{eq:b0b1}
	b_1(\Nf)={1\over(4\pi)^{4}}\bigg(102-{38\over 3}\Nf\bigg)\,.
\end{gather}
The coefficients $b_0(\Nf),b_1(\Nf)$ are universal and shared
by all mass-independent renormalization schemes. The scheme
dependence only enters through the higher-order coefficients,
$b_i(\Nf)\equiv b^\Obs_i(\Nf)$, 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
$\Nf$-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_{\Obs}$ is \emph{exactly} related to
$\Lambda_{\overline{\rm MS}}$ by the relation:
\begin{equation}
	\label{eq:MStoOLambda}
	{\Lambda_{\overline{\rm MS}}\over\Lambda_{\Obs}}=
	\exp\bigg\{{c_1(1)\over 2b_0}\bigg\}\,.
\end{equation}
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,
\begin{equation}
	\label{eq:MStoOCoupling}
	\bar{g}^2_{\overline{\rm MS}}(s\mu)=
	\bar{g}^{2}_{\Obs}(\mu)
	\Bigg(1+
	\sum_{k=1}^{M}
	c_k(s) \bar{g}^{2k}_{\Obs}(\mu)\Bigg)\,,
	\quad
	s>0\,.
\end{equation}
It is clear that the $\Lambda$-parameters are non-perturbatively
defined once the corresponding couplings and $\beta$-functions are.%
\footnote{Interestingly, eq.~(\ref{eq:MStoOLambda}) provides an
		  indirect non-perturbative definition of
		  $\Lambda_{\overline{\rm MS}}$ through the
		  $\Lambda$-parameter of any non-perturbative
  	  	  scheme.}
In addition, they are \emph{exact} solutions of the Callan-Symanzik
equations~\cite{Callan:1970yg,Symanzik:1970rt,Symanzik:1971vw},
and therefore RG invariants, i.e. $\rmd\Lambda_\Obs/\rmd\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, $\Nf$ (cf.~eq.~(\ref{eq:LambdaParam})).
Each $\Nf$-flavor theory therefore has its own $\Lambda$-parameters.

\subsubsection{Systematic uncertainties in extracting
			  $\Lambda^{(\Nf)}_{\overline{\rm MS}}$}
\label{subsec:SystematicErrorsPT}

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

\paragraph{Strategy 1.} Using eq.~(\ref{eq:LambdaParam})
we first get an \emph{asymptotic} estimate for
$\Lambda_{\Obs}/\mu_{\rm {PT}}$ as,
\begin{equation}
	\label{eq:LambdaMuPT}
	{\Lambda_{\Obs}\over\mu_{\rm PT}}
	=
	\varphi_{\rm g,\Obs}(\bar{g}_{\rm PT})
	\overset{\bar{g}_{\rm PT}\to0}{\approx}
	\varphi^{\rm PT}_{\rm g,\Obs}
	(\bar{g}_{\rm PT})+
	{\rm O}\big(\bar{g}_{\rm PT}^{2N-2}\big)\,.
\end{equation}
In this relation, $\varphi^{\rm PT}_{\rm g,\Obs}$ is defined
analogously to $\varphi_{\rm g,\Obs}$ of eq.~(\ref{eq:LambdaParam})
but with the replacement $\beta_\Obs\to\beta_\Obs^{\rm PT}$,
with $\beta_\Obs^{\rm PT}$ the perturbative $\beta$-function
of eq.~(\ref{eq:BetaFunctionPT}).  From the estimate of
${\Lambda_{\Obs}/\mu_{\rm PT}}$,
${\Lambda_{\overline{\rm MS}}/\mu_{\rm PT}}$
is obtained, with no further approximation, using
eq.~(\ref{eq:MStoOLambda}). Finally, the knowledge of
$\mu_{\rm {PT}}$ in physical units gives us
$\Lambda_{\overline{\rm MS}}$.

As anticipated by eq.~(\ref{eq:LambdaMuPT}), due to the
truncation of $\beta_\Obs$ to $N$-loop order in the
evaluation of ${\Lambda_{\Obs}/\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.~(\ref{eq:LambdaMuPT})
is in fact only asymptotic. Similarly to what we discussed
about eq.~(\ref{eq:Oth}), 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.

\paragraph{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.~(\ref{eq:MStoOCoupling}).
Given this, we can estimate
${\Lambda_{\overline{\rm MS}}/(s\mu_{\rm PT}})$
using eq.~(\ref{eq:LambdaParam}) and the perturbative
expression for the $\beta$-function, eq.~(\ref{eq:BetaFunctionPT}),
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~\cite{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.~(\ref{eq:LambdaMuPT}) and (\ref{eq:MStoOCoupling})).
Secondly, using eq.~(\ref{eq:LambdaParam}) 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}_{\Obs}$ used in lattice QCD
determinations of $\Lambda_{\overline{\rm MS}}$ we have
that $M\leq 3$ (see e.g.~ref.~\cite{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}_\Obs$,
is known (cf.~eq.~(\ref{eq:BetaFunctionPT})).%
\footnote{This is the case because for all these schemes
		  the $N$-loop $\beta$-function, $\beta^{\rm PT}_\Obs$,
		  has been inferred from $\beta^{\rm PT}_{\overline{\rm MS}}$
		  at $N$-loops using the matching relation
		  eq.~(\ref{eq:MStoOCoupling}) and
		  $\beta^{\rm PT}_\Obs$ at $(N-1)$-loop order.}
This means that, for the schemes $\bar{g}_\Obs$ commonly used,
\textbf{Strategy 1.} and \textbf{2.} result in the \emph{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.~(\ref{eq:MStoOCoupling})). 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}\to0$.
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}\to0$.

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_\Obs$. Reducing
the size of the perturbative truncation errors by a given factor
hence requires an \emph{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. \textbf{Strategy 1.} vs \textbf{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}\to0$.%
\footnote{We 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\to0}{\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}\to0$. 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.

\subsection{The accuracy of perturbation theory at high energy}
\label{subsec:AccuracyPT}

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}}$~\cite{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.

\subsubsection{The high-energy regime of $\Nf=3$ QCD}
\label{subsubsec:AccuracyNf3}

\begin{figure}
	\centering

	\resizebox{0.5\textwidth}{!}{%
	\includegraphics{plots/accuracy.pdf}}
	\caption{Determination of $\Lambda_{\overline{\rm MS}}^{(3)}$
	in units of $\mu_{\rm ref}=L_0^{-1}$ for different values
	of $\alpha_{\rm PT}$~\cite{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)$~\cite{DallaBrida:2018rfy}.  As
	the reader can see, when the extraction is performed at
	high-enough energies $(\alpha_{\rm PT}\approx0.1)$, all
	schemes nicely agree.}
	\label{fig:AccuracyNf3}
\end{figure}

We begin with the high-energy studies of
refs.~\cite{Brida:2016flw,DallaBrida:2018rfy} in $\Nf=3$ QCD.
In Fig.~\ref{fig:AccuracyNf3}, we show the results from
these references for $L_0\Lambda^{(\Nf=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\,\GeV$,
is a convenient high-energy reference scale
and, as in previous sections, $\bar{g}_{\rm PT}\equiv\bar{g}_{\Obs}(\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 \textbf{Strategy 1.} of
Sect.~\ref{subsec:SystematicErrorsPT}.
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\,\GeV$ to $140\,\GeV$.
The couplings considered in this study,
$\bar{g}_\Obs(\mu)=\bar{g}_{{\rm SF}_\nu}(L^{-1})$,
belong to a family of finite-volume renormalization schemes
based on the QCD Schr\"odinger functional
(SF)~\cite{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.~\cite{Brida:2016flw,DallaBrida:2018rfy,Sint:2012ae}.%
\footnote{Traditionally only the $\nu=0$ scheme has been
		  considered in applications, see
		  e.g.~refs.~\cite{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}\to0$. 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 $\Nf=3$
by~\cite{Bode:1998hd,Bode:1999sm,DallaBrida:2018rfy}:%
\footnote{For comparison the 3-loop coefficient of
	      $\beta_{\overline{\rm MS}}$ for $\Nf=3$
	      is $(4\pi)^3 b_2^{\overline{\rm MS}}= 0.324$.}
\begin{equation}
	(4\pi)^3 b_{2}^{{\rm SF}_\nu}=
 	-(0.064(27) +\nu\times 1.259(10))\,.
\end{equation}
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.~\ref{fig:AccuracyNf3},
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})$~\cite{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}\approx0.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.~\cite{DallaBrida:2018rfy}, which
corresponds to the gray band in the plot.

\begin{figure*}[!h]

	\centering
	\resizebox{0.4\textwidth}{!}{%
	\includegraphics{plots/msbar_nu+00.pdf}}
	\hspace*{10mm}
	\resizebox{0.4\textwidth}{!}{%
	\includegraphics{plots/msbar_nu-05.pdf}}
	\caption{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~\cite{DallaBrida:2018rfy}.
	The left (right) panel uses the ${\rm SF}_\nu$-scheme with
	$\nu= 0$ ($\nu=-0.5$), cf. text.}
	\label{fig:ScaleChangeNf3}
\end{figure*}

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}\approx0.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}\approx0.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}\approx70\,\GeV$
from the $\nu=0$ scheme (cf.~$\nu=0$ (fit C) in
Fig.~\ref{fig:AccuracyNf3})  \cite{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}\approx0.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 \textbf{Strategy 2}
of Sect.~\ref{subsec:SystematicErrorsPT}.
The corresponding estimates for
$\Lambda^{(3)}_{\overline{\rm MS}}/
\mu_{\rm ref}$ are shown in Fig.~\ref{fig:ScaleChangeNf3}
for the two cases, $\nu=0,-0.5$~\cite{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.~(\ref{eq:MStoOCoupling})). The values of $s$ considered
vary by about a factor $2-3$ around the value of
\emph{fastest apparent convergence}, $s^*$.%
\footnote{The scale factor $s^*$ is defined by
		 (cf.~eq.~(\ref{eq:MStoOLambda})):
		 \[
		 	s^*={\Lambda_{\overline{\rm MS}}\over\Lambda_{\Obs}}=
		 	\exp\bigg\{{c_1(1)\over 2b_0}\bigg\}\,.
		 \]
		 where $b_0$ is the 1-loop coefficient of
		 the $\beta$-function in eq.~(\ref{eq:b0b1}), 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_\Obs(\mu)$
		 (cf.~eq.~(\ref{eq:MStoOCoupling})). With the choice $s=s^*$,
		 the $k=1$ term in eq.~(\ref{eq:MStoOCoupling}) 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.~\cite{Zyla:2020zbs}).
Our intention is to test how this approach works in the
present case.

As one can see from Fig.~\ref{fig:ScaleChangeNf3}, 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}\to0$, and the different determinations come
together once $\alpha_{\rm PT}\lesssim0.1$.

What is clear from the results of
Fig.~\ref{fig:ScaleChangeNf3} 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}\approx0.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}\approx0.1$.

\subsubsection{The case of the pure Yang-Mills theory}
\label{subsubsec:YangMillsLambda}

\begin{figure*}
	\centering

	\resizebox{0.875\textwidth}{!}{%
	\includegraphics{plots/sq8t0xlam.pdf}}
	\caption{The dimensionless product
	$\sqrt{8t_0}\Lambda_{\overline{\rm MS}}^{(0)}$ as a function
	of $\alpha_{\rm PT}$~\cite{Nada:2020jay}. The empty symbols
	represent the data at the given $\alpha_{\rm PT}$, while the
	filled symbols are extrapolations to $\alpha_{\rm PT}\to0$
	(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.~\cite{DallaBrida:2019wur},
	while the data point labeled FlowQCD is the result of
	ref.~\cite{Kitazawa:2016dsl}.}
	\label{fig:PureGaugeLambda}
\end{figure*}

\paragraph{Finite-volume schemes.}

The second example that we consider is taken from the recent
study of ref.~\cite{Nada:2020jay} in the pure Yang-Mills (YM)
theory. This work presents an independent analysis of the results
from a previous study~\cite{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.~\ref{fig:PureGaugeLambda} we show the results for
$\sqrt{8t_0}\Lambda^{(0)}_{\overline{\rm MS}}$ from
ref.~\cite{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$~\cite{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  $\Nf=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)~\cite{Luscher:2010iy,Fodor:2012td,Fritzsch:2013je,DallaBrida:2016kgh}. The interested reader can find more details about the scheme
in the main refs.~\cite{DallaBrida:2019wur,Nada:2020jay}
(see also Sect.~\ref{subsec:CouplingFromDecoupling}).

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, \textbf{Strategy 1.} of
Sect.~\ref{subsec:SystematicErrorsPT} is employed using the
3-loop $\beta$-function in the GF-scheme of
choice~\cite{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 \textbf{Strategy 1.} based on the
${\rm SF}_{\nu=0}$-scheme (SF label in the plot), or by following
\textbf{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^*\approx2$ 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.~\ref{fig:PureGaugeLambda} 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}\approx0.08$ the different strategies give
estimates for $\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm had}$
which vary up to $\approx3\%$. This is about twice as large
as the statistical errors on the points (cf.~Table 3 of
ref.~\cite{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}\to0$. 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.~\cite{Nada:2020jay}  and conclude that only by studying
non-perturbatively the limit $\alpha_{\rm PT}\to0$ 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.

\paragraph{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.~\ref{fig:PureGaugeLambda},
where the very precise results labeled (FlowQCD) from
ref.~\cite{Kitazawa:2016dsl} show a net tension with
the determinations of refs.~\cite{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}\approx0.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.~\cite{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.~\cite{Husung:2020pxg} based on the ${\rm qq}$-coupling,
$\alpha_{\rm qq}$, defined from the static
potential~\cite{Sommer:1993ce}, are particularly interesting in
this respect.%
\footnote{A similar earlier study  on the challenges of
		  extracting $\Lambda_{\overline{\rm MS}}$ in both
		  $\Nf=2$ QCD and the pure-gauge theory using the
		  ${\rm qq}$-coupling can be found in
		  ref.~\cite{Leder:2011pz}. For a recent application
	  	  of this scheme for the computation of $\alphas$
	  	  and a detailed account of the most recent results
	  	  and developments see
	  	  refs.~\cite{Bazavov:2019qoo,Komijani:2020kst}.}
We report them in Fig.~\ref{fig:StaticPotentialC}.
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.~\cite{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
\textbf{Strategy 1.} of Sect.~\ref{subsec:SystematicErrorsPT},
i.e. it corresponds to $\alpha_{\rm PT}$ in our previous discussions.

\begin{figure}
	\centering

	\resizebox{0.5075\textwidth}{!}{%
		\includegraphics{plots/Lambdat0_ntermrho2.pdf}}
	\caption{$\Lambda_{\overline{\rm MS}}^{(0)}$ in units of
		$\sqrt{8t_0}$ determined at various values of
		$\alpha_{\rm qq}\equiv\alpha_{\rm PT}$~\cite{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.~\cite{DallaBrida:2019wur} and
		\cite{Kitazawa:2016dsl}.
		(The reference numbers in the plot are those of
		ref.~\cite{Husung:2020pxg} from which the plot
		is taken.)}
	\label{fig:StaticPotentialC}
\end{figure}

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)}$~\cite{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.~\cite{Husung:2020pxg} lattice spacings
down to $a\approx0.01\,\fm$ are reached while keeping the
lattice extent $L\approx2\,\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\geq3$,
$1\leq m\leq n-2$, which are enhanced at small couplings
(cf.~ref.~\cite{Husung:2020pxg}).

Figure~\ref{fig:StaticPotentialC} 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}\approx0.16-0.35$. As we can see from the
plot, for couplings $\alpha_{\rm qq}\gtrsim0.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}\lesssim0.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~\cite{Husung:2020pxg}. In fact, the errors are too large
to make definite conclusions for the relevant limit
$\alpha_{\rm qq}\to0$.

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$.

\subsection{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 $\Nf=3$ results of ref.~\cite{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.~\cite{Fritzsch:2013je,DallaBrida:2016kgh} and
eq.~(\ref{eq:GFcoupling}) 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~\cite{Luscher:1991wu}. Specifically, having set
$\mu=L^{-1}$,
\begin{equation}
	\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)\,.
\end{equation}
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.~\cite{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.~\cite{DallaBrida:2016kgh} are shown in
Fig.~\ref{fig:SSFGFCoupling}. 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.~\ref{fig:SSFGFCoupling}, 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
\begin{equation}
	\label{eq:FitsGFA}
	\Sigma(u_i,a/L)=\sigma^{(A)}_i+{r}^{(A)}_i\times (a/L)^2
	\quad
	\text{(fit A)}\,,
\end{equation}
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)}\equiv1/\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$.%
\footnote{Perturbation theory predicts:
		$\sigma(u)\overset{u\to0}{\approx} u+s_0u^2+{\rm O}(u^3)$,
		with $s_0=2b_0(\Nf)\ln(2)$, where $b_0(\Nf)$ is the
		universal 1-loop coefficient of the $\beta$-function,
		eq.~(\ref{eq:b0b1})
		(see e.g.~ref.~\cite{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.~\cite{DallaBrida:2016kgh} for a
		detailed discussion about this point.}
This observation suggests us to perform alternative fits to
the data in Fig.~\ref{fig:SSFGFCoupling} considering
the functional form
\begin{equation}
	\label{eq:FitsGFB}
	1/\Sigma(u_i,a/L)=1/\sigma^{(B)}_i+r^{(B)}_i\times (a/L)^2
	\quad
	\text{(fit B)}\,,
\end{equation}
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.~(\ref{eq:FitsGFA}). 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)~\cite{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\to0$
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.~\cite{Balog:2009yj,Balog:2009np,Husung:2019ytz,Husung:2020}).

\begin{figure}
	\centering

	\resizebox{0.495\textwidth}{!}{%
		\includegraphics{plots/SigmaGFcontlim2.pdf}}
	\caption{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$.}
	\label{fig:SSFGFCoupling}
\end{figure}

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.~\ref{fig:SSFGFCoupling} 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.~(\ref{eq:FitsGFA}) and (\ref{eq:FitsGFB})
fit the data well and perfectly agree with each other at finite
$L/a$, the corresponding extrapolations for $a/L\to0$
are clearly affected by some systematics.
In ref.~\cite{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.~(\ref{eq:FitsGFA}) and (\ref{eq:FitsGFB})
(cf.~ref.~\cite{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\to0$ 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.~\cite{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$).%
\footnote{A relevant exception is the case of the SF, for which
		  the leading discretization errors are parametrically
		  of O($a$) (cf.~Sect.~\ref{subsec:FiniteVolumeDecoupling}).
		  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.~\cite{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.~\cite{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.~\cite{DallaBrida:2016kgh,DallaBrida:2018rfy,DallaBrida:2019wur}
		  (cf.~the given references for the details).}
The results of refs.~\cite{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 $\Nf=2-4$.%
\footnote{The results refer to the contributions to
		  discretization effects coming from the lattice
		  action, considering several popular options
		  (cf.~ref.~\cite{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.~\cite{Husung:2019ytz,Husung:2020}).}%
\footnote{In 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.~\cite{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.%
\footnote{It 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.~\cite{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\approx1-3$, in the fit ans\"atze. 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.~\cite{DallaBrida:2016dai,DallaBrida:2017tru}
of the GF-coupling in the pure Yang-Mills theory
using Numerical Stochastic Perturbation Theory.%
\footnote{A recently expanded discussion in ref.~\cite{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.~\cite{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$~\cite{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$.

\begin{figure*}
	\centering

	\resizebox{0.495\textwidth}{!}{%
		\includegraphics{plots/k1c0p3_sp_imp.pdf}}
	\resizebox{0.495\textwidth}{!}{%
		\includegraphics{plots/k1c0p4_sp_imp.pdf}}
	\caption{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.}
	\label{fig:ContinuumLimitsNSPT}
\end{figure*}

In refs.~\cite{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^{(\Nf=0)}_{\overline{\rm MS}}$~\cite{Luscher:1995nr,Luscher:1995np}, one can thus infer
the relation
\begin{equation}
	\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.
\end{equation}
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\to0$. Focusing on the 1-loop coefficient, $k_1(a/L)$,
from SymEFT we expect that
(see refs.~\cite{DallaBrida:2016dai,DallaBrida:2017tru})
\begin{equation}
	\label{eq:k1}
	k_1(a/L)\overset{a/L\to0}{\sim}
	k_1(0)+ \sum_{m=2}^\infty
	\sum_{n=0}^1 r_{mn} (a/L)^m \ln(L/a)^n\,,
\end{equation}
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.~\cite{Husung:2019ytz} and also ref.~\cite{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.~\ref{fig:ContinuumLimitsNSPT} we show the results
from ref.~\cite{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.~\cite{DallaBrida:2016dai,DallaBrida:2017tru}
and eq.~(\ref{eq:GFcoupling})). 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.~\ref{fig:ContinuumLimitsNSPT}), 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.~(\ref{eq:k1}), considering the leading terms $m=2,n=0,1$.
The fits are good, $\chi^2/{\rm dof}\sim1$, 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.~\cite{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.

\section{Heavy-quark decoupling}
\label{sec:Decoupling}

So far we focused on the main challenges
that stand on the way of a precise determination of
$\Lambda^{(\Nf)}_{\overline{\rm MS}}$ and discussed
in detail the cases of $\Nf=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.~\cite{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 $\Nf=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.~\cite{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.~\cite{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 $\Nf=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\,\fm$, we  have that $am_c\gtrsim 0.3$ and
$am_b\gtrsim 1$, where for definiteness we took
$m_c\approx1.27\,\GeV$ and $m_b\approx4.2\,\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.~\cite{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 $\Nf=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 $\Nf=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 $\Nf=3$ flavor theory.

\subsection{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.~\cite{Bruno:2014ufa,Athenodorou:2018wpk}
for a more detailed presentation.

\subsubsection{The effective theory for heavy-quark decoupling}
\label{subsec:HQET}

We begin by considering QCD with $\Nf$ flavors of quarks,
which in short we denote ${\rm QCD}_{\Nf}$. Of these,
$\Nl$ are considered to be light, while the other
$\Nh\equiv\Nf-\Nl$ 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~\cite{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~\cite{Appelquist:1974tg}.

In formulas, the Lagrangian of the effective theory
is of the general form (see e.g.~ref.~\cite{Athenodorou:2018wpk})
\begin{equation}
	\label{eq:EffectiveLagrangian}
	\mathcal{L}_{\rm dec}=\mathcal{L}_0 +
	{1\over M}\mathcal{L}_1 + {1\over M^2}\mathcal{L}_2+\ldots\,,
\end{equation}
where the leading order corresponds to the Lagrangian of
${\rm QCD}$ with $\Nl$ light quarks, i.e.
$\mathcal{L}_0=\mathcal{L}_{\rm QCD_\Nl}$, while the corrections
$\mathcal{L}_k$, $k\geq 1$, consist of linear combinations of
local fields of mass dimension $4+k$, i.e.
\begin{equation}
	\label{eq:CountertermsLagrangians}
	\mathcal{L}_k= \sum_i \omega^{(k)}_i \Phi^{(k)}_i\,,
	\qquad
	[\Phi^{(k)}_i]=4+k\,,
\end{equation}
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}_{\Nf}$ 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}_\Nl$, has a single parameter:
the gauge coupling $\bar{g}^{(\Nl)}(\mu)$. The effective
and fundamental theory are therefore matched at leading order in
$M^{-1}$ once $\bar{g}^{(\Nl)}(\mu)$ is matched. This requires
that $\bar{g}^{(\Nl)}(\mu)$ is properly prescribed at a given
renormalization scale in a given renormalization scheme in terms
of the coupling $\bar{g}^{(\Nf)}(\mu)$ of the fundamental theory
and the heavy-quark masses $M$.%
\footnote{For 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,%
\footnote{This 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.~\ref{subsec:FiniteVolumeDecoupling}.}
it is possible to show that $\mathcal{L}_1=0$ and
therefore O($M^{-1}$) corrections are
absent~\cite{Athenodorou:2018wpk}.%
\footnote{Note that for the sake of argument we
		 exclude the uninteresting case of $\Nl=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_\Nl}$ multiplied by the light-quark
masses $m$~\cite{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.

\subsubsection{The effective theories and their couplings}
\label{subsec:EffectiveCouplings}

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~\cite{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}^{(\Nl)}(\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}^{(\Nf)}(\mu)$.

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

In phenomenological applications of this strategy
(see e.g.~ref.~\cite{Zyla:2020zbs}), the value of
$\bar{g}^{(\Nl)}(\mu_{\rm match})$ is extracted from
the value of the coupling $\bar{g}^{(\Nl)}(\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 $\Obs(q)$ with characteristic energy scale
$q\approx\mu_{\rm low}$, with its experimental
results. The effects of the heavy quarks in $\Obs(q)$ are
expected to be suppressed as O($(q/M)^2$)
(cf.~Sect.~\ref{subsec:HQET}). Hence, assuming that these
effects can be neglected, the perturbative expansion of $\Obs(q)$
can be considered in the $\Nl$-flavor theory, which allows
the coupling $\bar{g}^{(\Nl)}(\mu_{\rm low})$ to be extracted.
As observed above, the determination of the running of the
effective coupling in the $\Nl$-flavor theory does not require
any input from the fundamental theory: one can thus readily
obtain $\bar{g}^{(\Nl)}(\mu_{\rm match})$ from
$\bar{g}^{(\Nl)}(\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.~\cite{Aoki:2019cca}). Firstly, as discussed in
Sect.~\ref{subsec:SystematicErrorsPT}, through the study of the
non-perturbative running of a given massless coupling within
the effective $\Nl$-flavor theory, one determines the ratio
$\Lambda_{\overline{\rm MS}}^{(\Nl)}/\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 $\Nl$-flavor QCD. The physical units of
$\mu_{\rm had}$ and therefore $\Lambda_{\overline{\rm MS}}^{(\Nl)}$
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}}^{(\Nf)}/
\Lambda^{(\Nl)}_{\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}}^{(\Nl)}$,
$\Lambda_{\overline{\rm MS}}^{(\Nf)}$ is obtained.

It is important to note that in extracting
$\Lambda_{\overline{\rm MS}}^{(\Nl)}$ 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.~\cite{Athenodorou:2018wpk}, and start
by reformulating this matching in terms of $\Lambda$-parameters~\cite{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.

\subsection{Perturbative decoupling}

\subsubsection{Definitions}
\label{subsec:DecouplingPT}

As we are interested only in the QCD coupling we
assume for simplicity that the relevant effective theory
is given by massless ${\rm QCD}_{\Nl}$. 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}^{(\Nl)}(\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}^{(\Nl)}$
and of the fundamental theory $\bar{g}^{(\Nf)}$,
reads~\cite{Weinberg:1980wa,Bernreuther:1981sg,Bernreuther:1983zp}
\begin{equation}
	[\bar{g}^{(\Nl)}(\mu/\Lambda^{(\Nl)})]^2=
	F_\Obs\big([\bar{g}^{(\Nf)}(\mu/\Lambda^{(\Nf)})]^2,M/\mu\big)\,,
\end{equation}
where for later convenience we explicitly wrote the dependence
of the couplings on their corresponding $\Lambda$-parameter.
The function $F_\Obs$ depends in principle on the specific
observable $\Obs$ 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}^{(\Nf)}_{\overline{\rm MS}}
(\mu/\Lambda_{\overline{\rm MS}}^{(\Nf)})$ 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:
\begin{gather}
	\nonumber
	\big[\bar{g}^{(\Nl)}_{\overline{\rm MS}}(m_*/\Lambda^{(\Nl)}_{\overline{\rm MS}})\big]^2=
	g_*^2 C(g_*),\,
	\quad
	g_*\equiv \bar{g}^{(\Nf)}_{\overline{\rm MS}}(m_*/\Lambda^{(\Nf)}_{\overline{\rm MS}})\,,\\[1.5ex]
	\label{eq:CouplingMatchingPT}
	C(g_*)=1+h_2 g_*^4 + h_3 g_*^6 + h_4 g_*^8 + \ldots\,,
\end{gather}
where the $h_i$ are pure numbers that depend on both
$\Nh$ and $\Nl$~\cite{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.~(\ref{eq:CouplingMatchingPT}), in terms of RG-invariant
(RGI) quantities. Specifically, this means the
$\Lambda$-parameter of the effective theory,
$\Lambda^{(\Nl)}_{\overline{\rm MS}}$, that of the fundamental
theory $\Lambda^{(\Nf)}_{\overline{\rm MS}}$, and the RGI
quark-mass $M$ of the heavy quarks. The latter can be defined as
\begin{align}
	\nonumber
	M&\equiv\overline{m}^{(\Nf)}_{\overline{\rm MS}}(\mu)
	\varphi^{(\Nf)}_{\rm m, \overline{\rm MS}}(\bar{g}^{(\Nf)}_{\overline{\rm MS}}(\mu))\,,\\
	\nonumber
	\varphi^{(\Nf)}_{\rm m,\overline{\rm MS}}(\bar{g})&\equiv(2b_0(\Nf)\bar{g}^2)^{-{d_0\over 2b_0(\Nf)}}\times\\
	\label{eq:MassRGI}
	&\times\exp \bigg\{-\int_0^{\bar{g}}\rmd g
	\bigg[{\tau_{\overline{\rm MS}}^{(\Nf)}({g})\over
	\beta_{\overline{\rm MS}}^{(\Nf)}({g})} - {d_0\over b_0(\Nf)g} \bigg\}\,,
\end{align}
where $b_0(\Nf)$ is the 1-loop coefficient of the $\beta$-function,
eq.~(\ref{eq:b0b1}), and the function
\begin{equation}
	\tau^{(\Nf)}_{\overline{\rm MS}}(\bar{g})\equiv
	{\mu \over \overline{m}^{(\Nf)}_{\overline{\rm MS}}(\mu)}
	{\rmd \overline{m}^{(\Nf)}_{\overline{\rm MS}}(\mu)\over\rmd \mu}
	\bigg|_{\bar{g}}
\end{equation}
encodes the scale-dependence of the quark masses in the
$\Nf$-flavor theory (in the $\overline{\rm MS}$-scheme).
It has a perturbative expansion
\begin{equation}
	\label{eq:tauMSbar}
	\tau^{(\Nf)}_{\overline{\rm MS}}(\bar{g})\overset{\bar{g}\to0}{\approx}
	-\bar{g}^2\{d_0+\bar{g}^2 d_1+\ldots\}\,,
	\quad
	d_0=8/(4\pi)^2\,,
\end{equation}
known up to
5-loops~\cite{Larin:1993tq,Chetyrkin:1997dh,Vermaseren:1997fq,Baikov:2014qja}, where the actual scheme and $\Nf$-dependence start
from $d_i$, $i\geq 1$. Note that
even though in eq.~(\ref{eq:MassRGI}) 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.~(\ref{eq:LambdaParam}), and the matching relation
eq.~(\ref{eq:CouplingMatchingPT}), it is immediate to conclude
that~\cite{Bruno:2014ufa,Athenodorou:2018wpk}
\begin{equation}
	\label{eq:Plf}
	P_{\ell,\rm f}(M/\Lambda_{\overline{\rm MS}}^{(\Nf)})\equiv
	{\Lambda_{\overline{\rm MS}}^{(\Nl)}\over
	 \Lambda_{\overline{\rm MS}}^{(\Nf)}}\Bigg|_{\rm matched}=
 	{\varphi_{\rm g,\overline{\rm MS}}^{(\Nl)}(g_*\sqrt{C(g_*)})
 	 \over 	\varphi_{\rm g,\overline{\rm MS}}^{(\Nf)}(g_*)}\,.
\end{equation}
As anticipated by our notation,  $P_{\ell,\rm f}$ can be considered
as a function of $M/\Lambda_{\overline{\rm MS}}^{(\Nf)}$, since
the value of the coupling $g_*$ can be expressed in terms of the
RGI parameters through the relation
\begin{equation}
	\label{eq:gstar}
	{M\over\Lambda_{\overline{\rm MS}}^{(\Nf)}}=
	{\varphi_{\rm m,\overline{\rm MS}}^{(\Nf)}(g_*)
	\over
	\varphi_{\rm g,\overline{\rm MS}}^{(\Nf)}(g_*)}\,,
	\qquad
	g_*=\bar{g}^{(\Nf)}_{\overline{\rm MS}}(	{M/\Lambda_{\overline{\rm MS}}^{(\Nf)}}).
\end{equation}
In the limit where $M/\Lambda^{(\Nf)}_{\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$.

\subsubsection{Perturbative uncertainties}
\label{subsubsec:DecouplingPT}

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.~\ref{fig:PerturbativePlf} we present the results
taken from ref.~\cite{Athenodorou:2018wpk} for the relative
deviation
\begin{equation}
	\label{eq:P1Rel}
	{(P_{\ell,\rm f}-P_{\ell,\rm f}^{(1)})\over P^{(1)}_{\ell,\rm f}}
\end{equation}
of $P_{\rm 3,4}$ (left plot) and $P_{\rm 4,5}$ (right plot),
from their \emph{unsystematic} 1-loop approximation
$P^{(1)}_{\ell, \rm f}(M/\Lambda)=(M/\Lambda)^{\eta_0}$,
where $\eta_0=2\Nh/(33-2\Nl)$~\cite{Athenodorou:2018wpk}.%
\footnote{The unsystematic approximation
		  $P^{(1)}_{\ell, \rm f}$ obtained in
		  ref.~\cite{Athenodorou:2018wpk} happens to
		  be quite close to the different orders of
		  the perturbative expansion of $P_{\ell,\rm f}$
		  for the values of $\Nf$ and $\Nl$ 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~\cite{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.~(\ref{eq:Plf})-(\ref{eq:gstar}) are
considered. The $\tau$-function, eq.~(\ref{eq:tauMSbar}),
as well as the matching function $C(g_*)$ of
eq.~(\ref{eq:CouplingMatchingPT}) are considered to
a consistent order in the expansion (see
ref.~\cite{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~\cite{Zyla:2020zbs}.

\begin{figure*}[h]
	\centering

	\resizebox{0.375\textwidth}{!}{%
	\includegraphics{plots/corrections_43.pdf}}
	\hspace*{10mm}
	\resizebox{0.375\textwidth}{!}{%
	\includegraphics{plots/corrections_54.pdf}}
	\caption{Relative differences from the
	``1-loop approximation''
	$P^{(1)}_{\ell, \rm f}(M/\Lambda)=(M/\Lambda)^{\eta_0}$,
	$\eta_0=2\Nh/(33-2\Nl)$, for different orders of the
	perturbative expansion of $P_{\ell, \rm f}(M/\Lambda)$
	as a function of $M/\Lambda$~\cite{Athenodorou:2018wpk}.
	The results for $\Nl=3$, $\Nf=4$
	($\Nl=4$, $\Nf=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.}
	\label{fig:PerturbativePlf}
\end{figure*}

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.

\subsection{Non-perturbative decoupling}

\subsubsection{Definitions}
\label{subsubsec:NonperturbativeDecDef}

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}}^{(\Nf)}$
from $\Lambda_{\overline{\rm MS}}^{(\Nl)}$, let us begin by
recasting the matching of the effective and fundamental theory
in more non-perturbative terms.%
\footnote{For 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^{(\Nl)}$ is a properly chosen function of the scale
$\Lambda^{(\Nf)}$ 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$~\cite{Luscher:2010iy,Sommer:1993ce,Borsanyi:2012zs}.
Matching then means fixing the scale $\Lambda^{(\Nl)}$ through the
condition~\cite{Athenodorou:2018wpk}
\begin{equation}
	\label{eq:MatchingNPT}
	{\Lambda^{(\Nl)}\over \mathcal{S}^{(\Nl)}}=
	P_{\ell,\rm f}^{\mathcal{S}}(M/\Lambda^{(\Nf)})
	{\Lambda^{(\Nf)}\over \mathcal{S}^{(\Nf)}(M)}\,.
\end{equation}
In these equations $\mathcal{S}^{(\Nl)}$ and
$\mathcal{S}^{(\Nf)}(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^{(\Nf)})$ depends
on the scale $\mathcal{S}$ that is considered
in the matching relation. We also stress again that while
the value of $\Lambda^{(\Nf)}$ does not depend on $M$
the one of $\Lambda^{(\Nl)}$ does depend through the
matching condition.

Once $\Lambda^{(\Nl)}$ is properly fixed through
eq.~(\ref{eq:MatchingNPT}) in terms of $\mathcal{S}$,
for any other low-energy quantity $\mathcal{S}'$ we expect that
\begin{equation}
	\mathcal{S}'^{(\Nl)}=\mathcal{S}'^{(\Nf)}(M)
	+{\rm O}\big(\big({\Lambda^{(\Nf)}/ M}\big)^2\big)\,.
\end{equation}
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
\begin{equation}
	\label{eq:RatioS'S}
	{\mathcal{S}^{(\Nl)}\over \mathcal{S}'^{(\Nl)}}=
	{\mathcal{S}^{(\Nf)}(M)\over \mathcal{S}'^{(\Nf)}(M)}
	+{\rm O}\big(\big({\Lambda^{(\Nf)}/ M}\big)^2\big)\,,
\end{equation}
with $\mathcal{S}$ and $\mathcal{S}'$ any two low-energy scales.
Given this observation, from eq.~(\ref{eq:MatchingNPT}) we
conclude that
\begin{equation}
	P_{\ell,\rm f}^{\mathcal{S}}=
	P_{\ell,\rm f}^{\mathcal{S}'}+{\rm O}(({\Lambda^{(\Nf)}/ M})^2)\,.
\end{equation}
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.~(\ref{eq:MatchingNPT}) by
$\mathcal{S}^{(\Nf)}(0)/\Lambda^{(\Nf)}$, where
$\mathcal{S}^{(\Nf)}(0)$ stands for the low-energy
quantity $\mathcal{S}$ computed in the chiral limit
of the $\Nf$-flavor theory. Through this
simple manipulation we find
that~\cite{Bruno:2014ufa,Athenodorou:2018wpk}
\begin{align}
	\nonumber
	{\mathcal{S}^{(\Nf)}(M)\over \mathcal{S}^{(\Nf)}(0)}
	&= Q^{\mathcal{S}}_{\ell,\rm f} \times
	P^{\mathcal{S}}_{\ell,\rm f}(M/\Lambda^{(\Nf)})\\
	\label{eq:FactorizationQP}
	&= Q^{\mathcal{S}}_{\ell,\rm f} \times
	P_{\ell,\rm f}(M/\Lambda^{(\Nf)})
	+{\rm O}\big(\big({\Lambda^{(\Nf)}/ M}\big)^2\big)\,,
\end{align}
where the factor
\begin{equation}
	Q^{\mathcal{S}}_{\ell,\rm f}\equiv
	{\mathcal{S}^{(\Nl)}/\Lambda^{(\Nl)}\over
	\mathcal{S}^{(\Nf)}(0)/\Lambda^{(\Nf)} }\,,
\end{equation}
is defined in terms of the massless $\Nf$- and $\Nl$-flavor
theories. The interesting aspect of eq.~(\ref{eq:FactorizationQP})
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}^{(\Nf)}(M)/\mathcal{S}^{(\Nf)}(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~\cite{Athenodorou:2018wpk}
\begin{equation}
	\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^{(\Nf)}}\,.
\end{equation}
Considering eq.~(\ref{eq:FactorizationQP}),
this can be computed from the mass dependence of
any hadronic quantity as%
\footnote{Note that the number of light quarks on the
		  r.h.s. of this equation is implicitly given
		  by the difference between $\Nf$ and the number
		  $\Nh$ of heavy quarks of which the mass is varied.}
\begin{equation}
	\label{eq:etaM}
	\eta_{\ell,\rm f}^{M,\mathcal{S}}(M) =
	{M\over \mathcal{S}^{(\Nf)}(M)}
	{\partial {\mathcal{S}}^{(\Nf)}(M)\over \partial M}\bigg|_{\Lambda^{(\Nf)}}\,,
\end{equation}
with no need for determining $\mathcal{S}^{(\Nf)}(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^{(\Nf)}/M)^2$)
contributions. In the limit $M/\Lambda^{(\Nf)}\to\infty$, however,
the mass-scaling function becomes universal and can
asymptotically be estimated in perturbation
theory~\cite{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}$.

\subsubsection{Non-perturbative corrections to decoupling}
\label{subsubsec:NonperturbativeDec}

Non-perturbative effects in the decoupling of heavy quarks
have been systematically investigated in a series of recent
papers~\cite{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 $\Nf=2+1$ QCD is to the $\Nf=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
$\Nf=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 $\Nf=2+1$ QCD.

Studying the decoupling of the charm quark through
simulations of $\Nf=2+1+1$ and $\Nf=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.~\cite{Bruno:2014ufa,Athenodorou:2018wpk,Knechtli:2017xgy}
have investigated non-perturbatively a model system
given by QCD with $\Nf=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. $\Nl=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~\cite{Hollwieser:2020qri}. The results collected so far
confirm the findings of the model study.

\begin{figure*}[h]
	\centering

	\resizebox{0.8\textwidth}{!}{%
		\includegraphics{plots/etaM_all_fitc.pdf}}
	\caption{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}$
	\cite{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}}$~\cite{Athenodorou:2018wpk}.
	The vertical dotted lines mark the values of the quark mass
	$M_c,M_c/2$, and $M_c/4$.}
	\label{fig:etaM}
\end{figure*}

\paragraph{Ratios of $\Lambda$-parameters}

\begin{figure*}[h]
	\centering

	\resizebox{0.45\textwidth}{!}{%
	\includegraphics{plots/t0tc_dec_noglob.pdf}}
	\hspace*{10mm}
	\resizebox{0.45\textwidth}{!}{%
	\includegraphics{plots/t0w0_dec_noglob.pdf}}
	\caption{Continuum extrapolated values for
		$\sqrt{t_c/t_0}$ (left) and $\sqrt{t_0}/w_0$ (right)
		as a function of $\Lambda/M$~\cite{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$.}
	\label{fig:LowEnergyRatios}
\end{figure*}

The first set of results that we want to discuss are
from ref.~\cite{Athenodorou:2018wpk} and are shown
in figure~\ref{fig:etaM}. 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.~(\ref{eq:etaM})). 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.~$\Nf=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.%
\footnote{The 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\,\MeV$
		  from ref.~\cite{Fritzsch:2012wq} and
		  $M_c=1510\,\MeV$ from ref.~\cite{Heitger:2013oaa}
	  	  (cf.~ref.~\cite{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~\cite{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 \ref{fig:etaM} and comparing the results for
the different observables $\mathcal{S}$~\cite{Athenodorou:2018wpk}.

Putting all this information together, the authors of
ref.~\cite{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.~\cite{Athenodorou:2018wpk}.%
\footnote{The 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.~(\ref{eq:etaM}) and
	      ref.~\cite{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 \emph{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.

What can we conclude from this about the phenomenologically
interesting case of $P_{3,4}(M_c/\Lambda_{\overline{\rm MS}}^{(4)})$?
We first note that the dependence in perturbation theory of
$\eta^{M}_{\ell,\rm f}$ on $\Nl$ at fixed $\Nh$ is very mild.
At leading order it amounts to about a 20\% effect in going
from $\Nl=0$ to $\Nl=3$~\cite{Athenodorou:2018wpk}. For this
reason, the authors of ref.~\cite{Athenodorou:2018wpk} include
an additional 50\% contribution to their estimate for the
non-perturbative corrections $\Delta\eta^{M}_{0,2}$ to account
for the missing light-quark effects. Secondly, our intuition
from both perturbative and non-perturbative considerations
suggests that most likely the effects of the decoupling of a
single charm quark are about half of those of two quarks. From
these observations, one concludes that one can safely neglect
non-perturbative effects in connecting
$\Lambda_{\overline{\rm MS}}^{(3)}$ and
$\Lambda_{\overline{\rm MS}}^{(4)}$ down to a precision of 1.5\%
or better~\cite{Athenodorou:2018wpk}.

\paragraph{Ratios of hadronic quantities}
The second important category of effects that we
address are non-perturbative contributions from heavy quarks
to dimensionless ratios of low-energy quantities. For
us these are relevant in the context of setting the physical
scale of the theory (see e.g.~ref.~\cite{Sommer:2014mea}).
As well-known, in fact, besides the technical difficulties
of computing the relevant ratios, an important issue that one
faces in setting the scale in lattice QCD is the fact that one
never really simulates the ``real world'', where experiments
are conducted. Hence, when comparing the lattice results with
experiments in order to set the scale of the theory, one must
``correct'' the experimental quantities for effects that are
not taken into account in the lattice simulations, or at least
verify that these are not relevant once compared with the rest
of the uncertainties. The most common examples of effects to
be considered are the difference in the $u$- and $d$-quark
masses, electromagnetic effects, and the unphysical number of
quark flavors.

In the following we focus on the issue of estimating
the effects of the charm quark in low-energy quantities, and
to which extent this can be omitted in lattice simulations.
From the point of view of determining
$\Lambda_{\overline{\rm MS}}^{(4)}$ from the results
of $\Lambda_{\overline{\rm MS}}^{(3)}$, thus, the relevant
question is how accurate the determination of the physical scale
is from simulations in the 3-flavor theory. This amounts to
quantify how accurate it is to compute in $\Nf=2+1$ rather than
$\Nf=2+1+1$ QCD, the ratios of low-energy quantities that are
used to set the physical scale.

As presented in eq.~(\ref{eq:RatioS'S}), the effects of
heavy quarks in dimensionless ratios of low-energy
quantities are expected to be of O($M^{-2}$), provided
that the mass of the heavy quarks is large enough compared
to the energy scales of the observables considered.
In ref.~\cite{Knechtli:2017xgy}
a systematic study was conducted in order to assess the range
of heavy-quark masses for which the O($M^{-2}$) scaling of
the heavy-quark effects predicted by the effective theory
actually sets in. In addition, the authors estimated the
typical size of these contributions when $M\approx M_c$. For
their computations, they considered the very same model of
QCD with two degenerate heavy quarks previously introduced,
and the same range of masses, $M\approx M_c/8-1.2M_c$. The
results for several different ratios of low-energy quantities
obtained in the fundamental $\Nf=2$ QCD theory and in the
effective pure-gauge theory were compared. Note that since
the effective theory is purely gluonic only gluonic quantities
were considered. These, however, are all relevant observables
entering realistic scale-setting determinations.

In figure \ref{fig:LowEnergyRatios} we show two examples
given by the ratios $\sqrt{t_c/t_0}$ (left panel) and
$\sqrt{t_0}/w_0$ (right panel), evaluated in the continuum
limit. The results from $\Nf=2$ QCD for different values of
the heavy-quark masses $M\gtrsim M_c/4$ are shown, as well as
those from the effective pure-gauge theory which correspond
to $M=\infty$. Several fits to the data are proposed. Let us
focus first on the quadratic (blue) and linear (green) fits in
$\Lambda/M$. Among these four fit types, the data seem to favor
the leading-order effective theory predictions: $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$ free parameters, restricted to masses
$M/\Lambda \gtrsim 2.5$, i.e. $M\gtrsim M_c/2$ (blue bands).
Both linear fits in $\Lambda/M$, either excluding (green bands)
or including (green dashed lines) the point at $M\approx M_c/4$
have significantly larger $\chi^2$ per degree of freedom than
the previous fits (cf.~Table 3 in ref.~\cite{Knechtli:2017xgy}).
The pure $(\Lambda/M)^2$ fits which include the $M\approx M_c/4$
results (blue dashed lines), instead, are clearly off.

Although the results do not completely exclude a linear $\Lambda/M$
dependence, the fits obtained by adding to the leading-order
prediction a next-to-leading term $\propto(\Lambda/M)^{4}$ and
including data down to $M\approx M_c/4$ (red dashed lines),
further support the findings from the previous fits. Indeed,
the close agreement for $M\gtrsim M_c/2$ of these fits and the
leading-order predictions restricted to this range, reinforce
the conclusion that the O($M^{-2}$) scaling sets in for masses
$M\gtrsim M_c/2$, while for smaller masses higher-order
contributions become relevant, or the expansion has broken
down entirely.

\begin{figure}[h]
	\centering

	\resizebox{0.515\textwidth}{!}{%
	\includegraphics{plots/ratios.pdf}}
	\caption{Comparison of continuum limit extrapolations
	of $\sqrt{t_0/t_c}$ (bottom) and $\sqrt{t_0}/w_0$ (top)
	of the $\Nf= 3 + 1$ data from ref.~\cite{Hollwieser:2020qri}
	(right) with corresponding $\Nf= 3$ CLS results (left)
	from refs.~\cite{Bruno:2014jqa,Bruno:2016plf} including the
	finest ensemble J500, cf.~ref.~\cite{Bali:2019yiy}.}
	\label{fig:ScaleNf3p1}
\end{figure}

In general the corrections induced by the heavy quarks are
small. They amount to $2-3\%$ at the smallest masses in the plot,
around $M_c/4$, but they are reduced significantly below the
percent level, to about $0.4\%$, once $M\approx M_c$.
As discussed above, we expect that this result provides a reliable
estimate for the magnitude of these effects in the realistic case
of the decoupling of the charm quark in $\Nf=2+1+1$ QCD.
The simultaneous decoupling of two charm-like quarks rather than
just one, likely compensates the missing effects from the
light quarks; this at least in the purely gluonic quantities under
consideration.

These expectations are confirmed by the recent results
of ref.~\cite{Hollwieser:2020qri}. In figure \ref{fig:ScaleNf3p1}
we show their continuum limit extrapolations for different
discretizations of the ratios $\sqrt{t_0/t_c}$ and $\sqrt{t_0}/w_0$.
In this case, the results from $\Nf=3$ and $\Nf=3+1$ QCD are
compared. Both simulations include 3 mass-degenerate quarks with
a mass around the physical average mass of the $u$-, $d$- and
$s$-quarks. The $\Nf=3+1$ simulations include in addition a
fourth quark with a mass set to the physical value of the
charm. By comparing the results of the two set of simulations
one can directly test decoupling in a close-to-real situation.
As one can see from the figure, the differences in these very
precise ratios of low-energy gluonic quantities are far below
the percent level and of the order of magnitude found in the
model. Interestingly, the ensembles generated in
ref.~\cite{Hollwieser:2020qri} open the possibility to study
systematically the effects of the charm quark also in fermionic
low-energy observables.

\paragraph{Conclusions}
Let us summarize what we have learned from the studies above.
1) The ratio $\Lambda_{\overline{\rm MS}}^{(4)}/
\Lambda_{\overline{\rm MS}}^{(3)}$ can be safely computed
in perturbation theory with a precision of \emph{at least}
1.5\%, and realistically much better. In fact, it is important
to stress that the estimate for non-perturbative corrections
to the function $P_{3,4}(M_c/\Lambda_{\overline{\rm MS}}^{(4)})$
is \emph{very} conservative and the actual size of these effects
is much likely quite smaller, i.e. at the level of $0.3\%$ or
so~\cite{Athenodorou:2018wpk}. 2) Power corrections of
O($M^{-2}_c$) in low-energy observables are also found to
be very small, i.e. well-below the percent level
\cite{Knechtli:2017xgy}. All in all, this means that
$\Lambda_{\overline{\rm MS}}^{(5)}$ can be accurately predicted
at the 1-2\% level from $\Lambda_{\overline{\rm MS}}^{(3)}$.
In this respect, we note that the competitive precision of about
0.7\% on the $\alpha_s(M_Z)$ determination of
ref.~\cite{Bruno:2017gxd} (to be discussed below), corresponds to
an uncertainty of 3.5\% on $\Lambda_{\overline{\rm MS}}^{(3)}$.
This uncertainty is also conservative. In conclusions, there is
plenty of room for improvement within $\Nf=3$ QCD.

\section{Renormalization by decoupling}
\label{sec:RenByDec}

\subsection{The QCD coupling from $\Nf=3$ QCD}
\label{subsec:LambdaNf3}

In the previous section we established that by relying
on perturbative decoupling relations for the charm and
bottom quarks, precise determinations of
$\Lambda_{\overline{\rm MS}}^{(5)}$ are possible from
results in the $\Nf=3$ flavor theory. In order to be
able to make this statement the systematic studies on
the non-perturbative effects induced by the charm quark
in low-energy quantities have been instrumental.

Having this settled, as discussed in detail in
Sect.~\ref{sec:PrecisionDeterminations}, the main challenges
for an accurate determination of $\alphas$ on the lattice
are: 1) controlling discretization errors in continuum limit
extrapolations of the chosen non-perturbative coupling(s),
and 2) estimating the uncertainties associated with the use
of perturbation theory in extracting the $\Lambda$-parameter.
In this respect, the combined application of finite-volume
renormalization schemes and finite-size scaling techniques
has proven to be extremely effective in dealing with these
difficulties, paving the way for robust and precise
lattice determinations of the QCD coupling.

The determination of ref.~\cite{Bruno:2017gxd}, in particular,
relies on these techniques to compute
$\Lambda_{\overline{\rm MS}}^{(3)}$. The calculation reaches a
final precision on $\Lambda_{\overline{\rm MS}}^{(3)}$ of
about $3.5\%$, which translates into a $0.7\%$ uncertainty on
$\alphas(M_Z)$. The strength of this result lies in the fact
that all systematic uncertainties are carefully kept under control
at this competitive level of accuracy. The calculation is therefore
a prominent example of the current state-of-the-art determinations
of $\alphas$ from the lattice~\cite{Aoki:2019cca}. Below we want
to briefly recall the main steps that led to this result in order
to understand what the challenges are in improving on it. For a
more detailed presentation we refer the interested reader to the
original
references~\cite{Brida:2016flw,DallaBrida:2016kgh,Bruno:2017gxd,DallaBrida:2018rfy} and
reviews~\cite{Korzec:2017ypb,DallaBrida:2018cmc,Ramos:2019,Sint:2019bmq,DallaBrida:2019zmm}.

\begin{figure}[h]
	\centering

	\resizebox{0.475\textwidth}{!}{%
		\includegraphics{plots/alpha_running.pdf}}
	\caption{Running couplings of $\Nf= 3$ QCD from
		integrating the non-perturbative $\beta$-functions in the
		SF- and GF-schemes~\cite{Brida:2016flw,DallaBrida:2016kgh}.
		They  are  matched non-perturbatively at the scale
		$\mu_0$ defined by  $\bar{g}^2_{\rm SF_{\nu=0}}(\mu_0)  =  2.012$ by computing $\bar{g}^2_{\rm GF}(\mu_0/2)=2.6723(64)$~\cite{DallaBrida:2016kgh}.
		The scales $\mu_{\rm PT}=16\mu_0$ and $\mu_{\rm had}$
		defined by $\bar{g}^2_{\rm GF}(\mu_{\rm had})=11.31$ are
		also shown, as well as the perturbative prediction
		for the ${\rm SF}_{\nu=0}$-coupling for $\mu>\mu_{\rm PT}$
		using the 3-loop $\beta$-function.
		The red curves correspond to the results for
		$\alpha^{(3)}_{\overline{\rm MS}}(\mu)$	obtained from
		$\Lambda_{\overline{\rm MS}}^{(3)}
		=341(12)\,\MeV$~\cite{Bruno:2017gxd}, considering different
		perturbative orders for the $\beta$-function in the
		$\overline{\rm MS}$-scheme.}
	\label{fig:AlphaRunning}
\end{figure}

The $\Lambda_{\overline{\rm MS}}^{(3)}$ determination of
ref.~\cite{Bruno:2017gxd} was obtained from the study of the
non-perturbative running of some convenient finite-volume
schemes from a scale of about $0.2\,\GeV$ to roughly
$70\,\GeV$. The very high energies reached non-perturbatively
allowed for a systematic and robust assessment of the uncertainties
related to the application of perturbation theory. This study
has been presented in Sect.~\ref{subsubsec:AccuracyNf3}, where
the high-energy end of the running in the SF-schemes and the
result for $\Lambda_{\overline{\rm MS}}^{(3)}/\mu_{0}$ with
$\mu_{0}\approx 4.3\,\GeV$, have been discussed in detail.
The rest of the determination is built on the following steps.

Firstly, we have the computation of the lower end of the
running and corresponding determination of the ratio of
finite-volume scales $\mu_{0}/\mu_{\rm had}$
with $\mu_{\rm had}\approx 200\,\MeV$~\cite{DallaBrida:2016kgh}.
For this step, a novel finite-volume coupling defined in terms
of the Yang-Mills gradient flow was employed. This allowed for
reaching much greater precision than otherwise possible using
the schemes considered at high energy~\cite{Brida:2014joa}.
Note that a non-perturbative matching between the finite-volume
schemes at $\mu_0$ was performed in order to continue the running
at lower energies. For illustration, we show in figure
\ref{fig:AlphaRunning} the non-perturbative running of the
finite-volume couplings over the energy range covered,
together with the results for the coupling in the
$\overline{\rm MS}$-scheme obtained from the corresponding
determination of $\Lambda_{\overline{\rm MS}}^{(3)}$ (see below).

In a second step, the relation $\mu_{\rm had}/f_{\pi K}$
was established passing through the intermediate technical
scale $\mu_{\rm ref}^*=1/\sqrt{8t^*_0}$~\cite{Bruno:2017gxd}.
Here, $f_{\pi K}\equiv\frac{1}{3}(2f_K+f_\pi)$ is a
convenient combination of the pion and kaon decay constants,
while the scale $\mu_{\rm ref}^*$ is given in terms of the
flow time $t^*_0$ defined in the SU(3) flavor-symmetric
limit~\cite{Bruno:2016plf}. This step involved a combination of
small-volume and large-volume simulations in the hadronic
regime~\cite{Bruno:2017gxd}. From the experimental value of
$f_{\pi K}$ the precise physical units for $\mu_{\rm had}$ could
be inferred and hence those of $\Lambda_{\overline{\rm MS}}^{(3)}$.
Finally, perturbation theory was used for the functions
$P_{3,4}(M_c/\Lambda_{\overline{\rm MS}}^{(4)})$
and $P_{4,5}(M_b/\Lambda_{\overline{\rm MS}}^{(5)})$
to obtain $\Lambda_{\overline{\rm MS}}^{(5)}$ and from this
$\alphas(M_Z)$. Splitting the determination of
$\Lambda_{\overline{\rm MS}}^{(5)}$ over the above
steps was the key to keep all systematic errors under control.
With the proper choice of observables and techniques
the hard multi-scale problem of relating the low- and
high-energy sectors of QCD could be solved in full confidence.

\begin{figure}[h]
	\centering

	\resizebox{0.475\textwidth}{!}{%
	\includegraphics{plots/budget.pdf}}
	\caption{Contribution in percentage to the relative error
	squared of $\alphas(M_Z)$ from the different steps of the
	determination of ref.~\cite{Bruno:2017gxd}
	(cf. text for more details). As the reader can see, the
	dominant source of uncertainty is the non-perturbative
	running at high energy $\mu\approx4-70\,\GeV$
	($\Lambda_{\overline{\rm MS}}^{(\Nf=3)}/\mu_0$)~\cite{Brida:2016flw,DallaBrida:2018rfy},
	followed by the running at low energy $\mu\approx 0.2-4\,\GeV$
	($\mu_0/\mu_{\rm had}$)~\cite{DallaBrida:2016kgh}, and
	by scale setting
	($\mu_{\rm ref}^*$)~\cite{Bruno:2016plf}.
	Note that the error from decoupling
	($\Lambda_{\overline{\rm MS}}^{(5)}/
	\Lambda_{\overline{\rm MS}}^{(3)}$) is only	perturbative.
	However, even adding the very conservative uncertainty
	estimated in Sect.~\ref{subsubsec:NonperturbativeDec} for
	the non-perturbative contributions to the decoupling of
	the charm quark, the total error is still dominated by
	the one from the running.}
	\label{fig:ErrorBudget}
\end{figure}

It is now instructive to look at the error budget of
this $\alphas$ determination. This is given in figure
\ref{fig:ErrorBudget} which shows the contribution in
percentage to the relative error squared on $\alphas(M_Z)$ from
the different steps described above~\cite{Korzec:2017ypb}.
As it is clear from the figure, the main source of uncertainty
comes from the determination of the non-perturbative running
from $\mu_{\rm had}\approx 0.2\,\GeV$ up to
$\mu_{\rm PT}\approx 70\,\GeV$, where perturbation theory
is applied to extract $\Lambda_{\overline{\rm MS}}^{(3)}/\mu_{\rm PT}$ (cf.~Sect.~\ref{subsubsec:AccuracyNf3}). In particular,
the error accumulated by running from $\mu_{0}$ to $\mu_{\rm PT}$
(labeled as $\Lambda_{\overline{\rm MS}}^{(\Nf=3)}/\mu_0$ in
the plot) contributes roughly 60\% of the total budget.

It is important to recall at this point that the error
coming from the running is completely dominated by
statistical uncertainties. In particular, thanks to
the fact that $\mu_{\rm PT}\approx 70\,\GeV$ was reached
non-perturbatively, the uncertainties due to the use of
perturbation theory are well below the statistical errors
(cf.~Sect.~\ref{subsubsec:AccuracyNf3}). In this respect,
we want to stress the important difference between this
and the majority of other lattice determination of $\Lambda_{\overline{\rm MS}}^{(3)}$,
where perturbation theory is applied at scales
$\mu_{\rm PT}\lesssim 2-3\,\GeV$. In these cases, a large
fraction of the final error comes from the \emph{systematic}
uncertainties related to the truncation of the perturbative
series and possible remnants of non-perturbative contributions
(cf.~ref.~\cite{Aoki:2019cca}). As we have seen, estimating
these systematics reliably is very challenging, particularly so
when precision is desired but the accessible range of
scales is limited to low energy. In this situation,
a reduction of the final uncertainties is highly non-trivial,
and can eventually come only from reaching significantly higher
energy scales. Without a step-scaling approach this is in practice
extremely demanding in QCD given the present computational and
algorithmic capabilities.%
\footnote{We recall that with the current state-of-the-art
		  algorithms, the cost of QCD simulations for reducing
		  the lattice spacing $a$ at fixed physical volume
	  	  and number of independent configurations scales
	  	  $\propto a^{-7}$
	  	  (cf.~e.g.~ref.~\cite{Schaefer:2012tq}). Hence,
	  	  reducing the lattice spacing by a factor $2$ or so
	  	  requires a factor of O($100$) in the cost.}

In the case of the step-scaling method, on the other hand,
reducing the uncertainties on the current $\alphas(M_Z)$
determination is a question of reducing the statistical
uncertainties coming from the running of the coupling(s)
in $\Nf=3$ QCD. This is in principle a straightforward task.
However, reducing the total error by an \emph{important}
factor, say a factor 2 or so, is yet a non-trivial
challenge from the computational point of view.

Rather than following a brute force approach for the reduction
of the error in the computation of the running in $\Nf=3$ QCD,
in the following we shall discuss a novel strategy
which promises the desired error reduction in a
significantly cheaper way~\cite{DallaBrida:2019mqg}.
It is based on the ideas presented in the previous section
on heavy-quark decoupling. The distinct feature of the approach
is that one can replace the non-perturbative computation
of the running in $\Nf=3$ QCD needed to determine
$\Lambda_{\overline{\rm MS}}^{(3)}$ with the
running in the pure Yang-Mills theory.%
\footnote{We take the opportunity to point out the study of
		  ref.~\cite{Brower:2014dfa} where ideas based on
		  heavy-quark decoupling are considered in order to
		  address other RG-related problems in the context
		  of Beyond the Standard Model Physics.}

\subsection{The coupling from decoupling}

\subsubsection{General strategy and master formula}

Let us begin by considering QCD with $\Nf$ flavors of
heavy quarks of RGI mass $M$. In this theory all quarks
are massive and there are no light quarks. As we have
seen in Sect.~\ref{subsec:HQET}, as the mass $M$ becomes
larger and larger this theory is expected to be approximated
better and better by an effective theory given
by the pure Yang-Mills theory.
In particular, once the $\Lambda$-parameters of the fundamental
and effective theory are properly matched, dimensionless
low-energy observables can be computed in the effective theory
up to corrections of O($M^{-n}$), where, in general, $n=2$
if the theories are matched at leading order.

The decoupling of heavy quarks is also valid for
the interesting case of couplings defined in massive
renormalization
schemes~\cite{Appelquist:1974tg,Bernreuther:1981sg}. This
is a direct consequence of the fact that such couplings
are defined in terms of dimensionless observables in the
massive theory. Following the notation of Sect.~\ref{subsubsec:LambdaParams}, we
indicate the generic renormalized  massive coupling in
the $\Nf$-flavor theory as $\bar{g}^{(\Nf)}_{\Obs}(\mu,M)$,
where $\Obs$ denotes the short-distance observable
used to define the coupling (cf.~eq.~(\ref{eq:alphaO})),
while $\mu$ is the renormalization scale. From the decoupling
of heavy quarks applied to the observable $\Obs$ it follows
that,
\begin{equation}
	\label{eq:CouplingDecouplingNP}
	\bar{g}^{(\Nf)}_\Obs(\mu,M)=
	\bar{g}^{(0)}_\Obs(\mu)+{\rm O}(M^{-2})\,,
\end{equation}
where $\bar{g}^{(0)}_\Obs(\mu)$ refers to the corresponding
coupling in the pure Yang-Mills theory evaluated at the
\emph{same} physical scale $\mu$. Note that here and below
we shall  loosely denote as ${\rm O}(M^{-k})$ terms,
contributions that contain terms of
${\rm O}((\Lambda/M)^{k})$ as well as ${\rm O}((\mu/M)^{k})$.

The decoupling relation, eq.~(\ref{eq:CouplingDecouplingNP}),
can be equivalently recast in terms of the renormalization scales
$\mu$ implicitly defined by the couplings. Specifically, given a
numerical value for the coupling $g_M$, we define the scales
$\mu_{\rm dec}^{(\Nf)}$ and $\mu_{\rm dec}^{(0)}$ through
\begin{equation}
	\label{eq:mudec}
	\bar{g}^{(\Nf)}_\Obs(\mu_{\rm dec}^{(\Nf)},M)=
	g_M
	=\bar{g}^{(0)}_\Obs(\mu_{\rm dec}^{(0)})\,.
\end{equation}
From the theory of decoupling it then follows that
\begin{equation}
	\mu_{\rm dec}^{(0)}=
	\mu_{\rm dec}^{(\Nf)}+{\rm O}(M^{-2})\,.
\end{equation}
We stress again that we assume that
the $\Lambda$-parameters in the two theories are
properly matched.

From this basic observation the master formula proposed in
ref.~\cite{DallaBrida:2019mqg} follows. We start by
considering the  relation in eq.~(\ref{eq:MatchingNPT}),
and take for the low-energy scale $\mathcal{S}$ the
renormalization scale $\mu_{\rm dec}$ defined above
in terms of the given couplings. In formulas,
\begin{equation}
	\label{eq:PreMasterEquation0}
	{\Lambda^{(0)}_{\overline{\rm MS}}\over \mu_{\rm dec}^{(0)}}=
	P_{0,\rm f}(M/\Lambda^{(\Nf)}_{\overline{\rm MS}})
	{\Lambda^{(\Nf)}_{\overline{\rm MS}}\over \mu_{\rm dec}^{(\Nf)}}+{\rm O}(M^{-2})\,,
\end{equation}
where for later convenience we took the $\Lambda$-parameters
in the $\overline{\rm MS}$-scheme.%
\footnote{Note that the relation is non-perturbatively valid
		  at this point as $\Lambda_{\overline{\rm MS}}$ can
	  	  be expressed exactly in terms of any non-perturbative
  	  	  scheme (cf.~Sect.~\ref{subsubsec:LambdaParams}).}
Now, rather than interpreting the above equation as a matching
relation for the $\Lambda$-parameters that defines the function
$P_{0,\rm f}$, we shall turn tables and use it to \emph{predict}
the ratio ${\Lambda^{(\Nf)}_{\overline{\rm MS}}/\mu_{\rm dec}^{(\Nf)}}$ in terms of ${\Lambda^{(0)}_{\overline{\rm MS}}/ \mu_{\rm dec}^{(0)}}$.

To this end, we first replace the function $P_{0,\rm f}$ with
its perturbative expansion $P_{0,\rm f}^{\rm PT}$ in the
${\overline{\rm MS}}$-scheme to some order
$n$ (cf.~eqs.~(\ref{eq:Plf}),(\ref{eq:LambdaMuPT})
and Sects.~\ref{subsubsec:DecouplingPT},
\ref{subsubsec:NonperturbativeDecDef}),
\begin{equation}
	\label{eq:PreMasterEquation1}
	P_{0,\rm f}(M/\Lambda^{(\Nf)}_{\overline{\rm MS}})=
	P^{\rm PT}_{0,\rm f}(M/\Lambda^{(\Nf)}_{\overline{\rm MS}})+
	{\rm O}(g_*^{2n-2},M^{-2})\,.
\end{equation}
We therefore have that,
\begin{equation}
	\label{eq:PreMasterEquation2}
	{\Lambda^{(\Nf)}_{\overline{\rm MS}}\over \mu_{\rm dec}^{(\Nf)}}
	P^{\rm PT}_{0,\rm f}(M/\Lambda^{(\Nf)}_{\overline{\rm MS}})=
	{\Lambda^{(0)}_{\overline{\rm MS}}\over \mu_{\rm dec}^{(0)}}
	+{\rm O}(g_*^{2n-2},M^{-2}).
\end{equation}
As a second step, using eq.~(\ref{eq:LambdaParam}) and
the definition in eq.~(\ref{eq:mudec}) we write,
\begin{equation}
	\label{eq:PureGaugeLambda}
	{\Lambda^{(0)}_{\overline{\rm MS}}\over \mu_{\rm dec}^{(0)}}
	= {\Lambda^{(0)}_{\overline{\rm MS}}\over
	   \Lambda^{(0)}_\Obs}\,\varphi^{(0)}_{\rm g,\Obs}(g_M)\,,
\end{equation}
which only involves quantities in the pure-gauge theory.
We recall that $\mathcal{O}$ appearing here is the observable
used to calculate the coupling in
eq.~(\ref{eq:mudec}). Moreover, as discussed in
Sect.~\ref{subsubsec:LambdaParams}, the change of scheme given
by the ratio ${\Lambda^{(0)}_{\overline{\rm MS}}/
\Lambda^{(0)}_\Obs}$ can be computed exactly through a
1-loop calculation (cf.~eq.~(\ref{eq:MStoOLambda})).

Finally, introducing the dimensionless variables,
\begin{equation}
	\rho\equiv
	{\Lambda^{(\Nf)}_{\overline{\rm MS}}/
	\mu_{\rm dec}^{(\Nf)}}\,,
	\qquad
	z\equiv {M/\mu_{\rm dec}^{(\Nf)}}\,,
\end{equation}
using eqs.~(\ref{eq:PreMasterEquation2}) and
(\ref{eq:PureGaugeLambda}) we arrive at the master
equation,
\begin{equation}
	\label{eq:MasterEquation}
	\rho\,P^{\rm PT}_{0,\rm f}(z/\rho) =
	{\Lambda^{(0)}_{\overline{\rm MS}}\over
	\Lambda^{(0)}_\Obs}\,\varphi^{(0)}_{\rm g,\Obs}(g_M)
	+{\rm O}(g_*^{2n-2},M^{-2}),
\end{equation}
which can be solved for $\rho$ once the pure-gauge function
$\varphi^{(0)}_{\rm g,\Obs}(g_M)$ is known. As promised,
the master formula allows us to replace the non-perturbative
computation of the running of the coupling from the low-energy
scale $\mu_{\rm dec}$ up to infinite energy in $\Nf$-flavor
QCD with the corresponding running in the pure-gauge theory.
A few remarks are in order at this point.

First of all, it is important to stress the fact that
eq.~(\ref{eq:MasterEquation}) is \emph{exact} in the limit
where $M\to\infty$. In this limit both the perturbative
O($g_*^{2n-2}$) corrections and the non-perturbative O($M^{-2}$)
contributions vanish. The basic idea that is applied in this
strategy is in fact similar to when we extract
$\Lambda_{\overline{\rm MS}}^{(4)}$ from
$\Lambda_{\overline{\rm MS}}^{(3)}$
by replacing the non-perturbative function
$P_{3,4}(M_c/\Lambda^{(4)}_{\overline{\rm MS}})$ with its
perturbative approximation to some order, neglecting both
higher-order terms and non-perturbative corrections. The
crucial difference in the present case is that the approximation
can be made systematically better by considering larger values
of $M$, which, at least in principle, is a free parameter of
the strategy. In this respect, note that the perturbative
corrections to $P_{\rm 0, f}(M/\Lambda^{(\Nf)}_{\overline{\rm MS}})$
only depend on the value of $M/\Lambda^{(\Nf)}_{\overline{\rm MS}}$,
while the choice of scale $\mu_{\rm dec}^{(\Nf)}$ does not
matter. From the results presented in
Sect.~\ref{subsubsec:DecouplingPT}, we infer that
the perturbative errors due to the truncation of
the perturbative series for
$P_{\rm 0, f}(M/\Lambda^{(\Nf)}_{\overline{\rm MS}})$
are already small for
$M/\Lambda_{\overline{\rm MS}}^{(\Nf)}\approx 5$,
for the relevant values of $\Nf$. On the
other hand, controlling the non-perturbative O($M^{-2}$)
terms requires to have both
O($(\Lambda^{(\Nf)}_{\overline{\rm MS}}/M)^2$) and
O($(\mu_{\rm dec}^{(\Nf)}/M)^2$) terms under control.
Whether this is possible in practice must be assessed
by carefully studying the limit $M\to\infty$ with
the accessible values of $M$.

In comparing  eqs.~(\ref{eq:MatchingNPT}) and
(\ref{eq:PreMasterEquation0}), the attentive reader might
have noticed that we omitted the $M$-dependence
on the scale $\mu^{(\Nf)}_{\rm dec}$. This was intentional.
As we shall see in the later subsection, in practice
it is convenient in fact to define a single scale,
$\mu^{(\Nf)}_{\rm dec}$, common to all the $\Nf$-flavor
theories defined by the different values of $M$.%
\footnote{The fact that we can set a unique scale for all
	theories with different $M$ values is a direct
	consequence of the fact that they all share the very
	same $\Lambda$-parameter and this is $M$ independent.}
This means that, as anticipated by our notation in
eq.~(\ref{eq:mudec}), we will have different values for $g_M$
depending on the value of $M$ considered. The value of $g_M$
is found by computing
$\bar{g}^{(\Nf)}_\Obs(\mu^{(\Nf)}_{\rm dec},M)$
for the given $M$ at the common scale $\mu^{(\Nf)}_{\rm dec}$.
In lattice QCD, setting a common scale among the theories with
different $M$ values can be achieved via the bare
parameters. This amounts to consider the very same lattice
discretization and establish a line of constant physics
along which $\mu_{\rm dec}^{(\Nf)}$ is  kept fixed. The
massive couplings are then evaluated at matching values
of the bare coupling along this line of constant physics
for the different bare quark masses corresponding to the
target RGI masses. From their continuum limit extrapolations
we find the values $g_M$.

Finally, in order to extract $\Lambda_{\overline{\rm MS}}^{(\Nf)}$
from the results for $\rho$ determined from
eq.~(\ref{eq:MasterEquation}) it is necessary to know the
value of $\mu^{(\Nf)}_{\rm dec}$ in physical units. This is
obtained by establishing the relation
$\mu^{(\Nf)}_{\rm dec}/\mu^{(\Nf)}_{\rm phys}$,
where $\mu^{(\Nf)}_{\rm phys}$ denotes a convenient low-energy
scale computed in $\Nf$-flavor QCD at \emph{physical} values of the
quark masses. The scale $\mu^{(\Nf)}_{\rm phys}$ can thus
be related to its experimental counterpart. Of course, it goes
without saying that in order to be able to set the scale
accurately in terms of experimentally measurable quantities,
as well as to perturbatively match
$\Lambda_{\overline{\rm MS}}^{(\Nf)}$ and
$\Lambda_{\overline{\rm MS}}^{(5)}$, we must consider $\Nf=3$
or $4$.

\subsubsection{Another hard multi-scale problem?}
\label{subsubsec:MultiScaleProblem2}

The general strategy presented above is certainly very compelling.
As any other strategy, though, in practical implementations
it comes with its challenges. In particular, in
order to have all systematic effects under control, it is necessary
to carefully address how to accommodate in lattice simulations
the different scales that enter the problem. As we shall see, a
naive approach can easily end up facing severe limitations.

In the general situation, first of all, the space-time volume
has to be large enough for finite-volume effects to be under
control in all relevant observables. This means that the
infrared cutoff set by the linear extent $L$ of the lattice must
be much smaller than all other scales. Secondly, in order to have
small decoupling corrections in eq.~(\ref{eq:MasterEquation}) we
must have that the heavy-quark mass $M$ is larger than all
other physical scales, as in particular  $\mu^{(\Nf)}_{\rm dec}$
and $\Lambda_{\overline{\rm MS}}^{(\Nf)}$. Note that although
$\mu^{(\Nf)}_{\rm dec}$ is in principle arbitrary, in practice
it is not convenient to take this scale to be much larger than
$\Lambda_{\overline{\rm MS}}^{(\Nf)}\sim \mu^{(\Nf)}_{\rm phys}$.
Last but not least, all scales have to be
well below the ultraviolet cutoff set by the lattice spacing.
Putting all these constraints together we find
(cf.~eq.~(\ref{eq:Window})),
\begin{equation}
	\label{eq:WindowProblem2}
	L^{-1}\ll
	\mu^{(\Nf)}_{\rm phys} \sim
	\Lambda_{\overline{\rm MS}}^{(\Nf)}\sim
	\mu^{(\Nf)}_{\rm dec}
	\ll
	M
	\ll
	a^{-1}\,.
\end{equation}
It is clear from this series of inequalities that
having all these scales comfortably resolved on a
\emph{single} lattice is challenging and requires
a very large $L/a$. Just to give an
example, if one attempts the calculations using
a state-of-the-art large-volume ensemble with, say,
$L/a=100$, $m_\pi L=4$, $m_\pi=140\,\MeV$,
which results in $a\approx 0.056\,\fm$,
eq.~(\ref{eq:WindowProblem2})
translates into $M\ll3.5\,\GeV$.

\subsubsection{Finite-volume couplings rescue us again}
\label{subsec:FiniteVolumeDecoupling}

Some of the constraints encoded in eq.~(\ref{eq:WindowProblem2})
can be lifted by considering for the coupling in
eq.~(\ref{eq:CouplingDecouplingNP}) a finite-volume scheme,
i.e. $\bar{g}_\Obs(\mu,M)\equiv\bar{g}_\Obs(L^{-1},M)$
(cf.~Sect.~\ref{subsubsec:FV}). With this choice, the
determination of $\bar{g}_\Obs(\mu^{(\Nf)}_{\rm dec},M)$ does
not require the physical volume to be large. Furthermore, if
the decoupling scale $\mu^{(\Nf)}_{\rm dec}$ is also
defined in terms of a finite-volume coupling, i.e.
$\mu^{(\Nf)}_{\rm dec}\equiv L^{-1}_{\rm dec}$, one has some
additional freedom in the choice of the value of
$\mu^{(\Nf)}_{\rm dec}$. Any sizable scale separation
between $\mu^{(\Nf)}_{\rm dec}$ and the hadronic scale
$\mu^{(\Nf)}_{\rm phys}$ can in fact be bridged through
step-scaling within the $\Nf$-flavor theory
(cf.~Sect.~\ref{subsubsec:FV}). Taking $\mu^{(\Nf)}_{\rm dec}$
larger at fixed $z=M/\mu^{(\Nf)}_{\rm dec}$, allows us to reach
larger values of $M/\Lambda_{\overline{\rm MS}}^{(\Nf)}$,
as well as to profit from simulating at smaller values of
$a$. On the other hand, the larger $\mu^{(\Nf)}_{\rm dec}$
is, the smaller is the range of energy scales for which,
through eq.~(\ref{eq:MasterEquation}), the running of
the coupling in the $\Nf$-flavor theory is replaced
by that in pure-gauge. In practice, choosing
$\mu^{(\Nf)}_{\rm dec}={\rm O}(1\,\GeV)$ is
a good compromise.

By employing finite-volume couplings and scales
the only constraints that we have to meet are to
have small decoupling corrections in eq.~(\ref{eq:MasterEquation}),
and to keep discretization errors under control.
The first condition requires $z=ML_{\rm dec}\gg1$,
while discretization effects are small once
$aM\ll1$ and $a/L_{\rm dec}\ll1$. Putting these
inequalities together we find the conditions
\begin{equation}
	L_{\rm dec}/a\gg z\gg1.
\end{equation}
If we take, say,
$\mu^{(\Nf)}_{\rm dec}=1\,\GeV$ and
$L_{\rm dec}/a=50$, then $M\ll50\,\GeV$.
In summary, using finite-volume schemes, for a given set of
lattice sizes $L/a$ we can reach much finer lattice spacings
$a$, as the extent of the lattice $L$ does not have to be large
in physical units. Smaller $a$ values allow us to consider
larger $M$ values, while having $aM$ reasonably small.
Larger values of $M$ make for a more precise approximation
$P_{\rm 0,f}^{\rm PT}(M/\Lambda)$. Still, large $z$ values
have to be reached in order to control non-perturbative
decoupling corrections.

\paragraph{Heavy-quark decoupling in a finite volume}

As pointed out for the case of computing the running
through a step-scaling procedure (cf.~Sect.~\ref{subsubsec:FV}),
the choice of finite-volume coupling is dictated by several
technical aspects, as for instance, statistical precision and
discretization errors. For the strategy based on decoupling
an additional factor becomes relevant which is the size of
non-perturbative contributions in the decoupling of heavy quarks
(cf.~eq.~(\ref{eq:CouplingDecouplingNP})). In a finite volume,
the situation can be quite different from one coupling
definition to another, as even the leading power in $M^{-1}$
may be different.

Most finite-volume couplings that are used in practice are
based on the QCD Schr\"odinger functional (SF)
\cite{Luscher:1992an,Sint:1993un,Sint:1995rb}. In the
SF the quark and gluon fields satisfy Dirichlet boundary
conditions at the space-time boundaries located at $x_0=0$ and
$T$, where $T$ is the temporal extent of the space-time volume
(cf.~eqs.~(\ref{eq:SFbc})-(\ref{eq:SFbc2})). These boundary
conditions guarantee many compelling
features~\cite{Luscher:1992an,Sint:1993un}. However, they
come with the price of having, for instance, additional
discretization effects of
O($a$)~\cite{Luscher:1992an,Luscher:1996sc}.
Using Symanzik effective theory these can be understood as
dimension 4 counterterms localized at the space-time
boundaries~\cite{Luscher:1992an,Luscher:1996sc}.
In close analogy with Symanzik effective theory an analysis of
the effective Lagrangian for heavy quarks in the presence of SF
boundary conditions shows that the same boundary fields
appear as O($M^{-1}$) counterterms. More precisely, considering
the relevant case $\Nl=0$, the Lagrangian $\mathcal{L}_1$ in
eq.~(\ref{eq:EffectiveLagrangian}) has the form
(cf.~eq.~(\ref{eq:CountertermsLagrangians}))
\begin{gather}
	\nonumber
	\mathcal{L}_1= \omega_b(g)
	\big[
	\mathcal{B}(0)+\mathcal{B}(T)\big]\,,\\[1ex]
	\label{eq:OMTermsSF}
	\mathcal{B}(x_0)
	=-{1\over g^2}\int \rmd \boldsymbol{x}\,
	  \tr\{F_{0k}(x)F_{0k}(x)\}|_{x_0}\,,
\end{gather}
where $F_{\mu\nu}(x)$ is the gluon-field strength tensor,
while $\omega_b$ is a coefficient function that can be fixed
by matching with the fundamental $\Nf$-flavor theory.
Note that for simplicity we listed the only gluonic
operator that is relevant for the class of the SF boundary
conditions normally employed (cf.~eqs.~(\ref{eq:SFbc})).
As a result, if couplings based on the SF are considered, the
decoupling relations
eqs.~(\ref{eq:CouplingDecouplingNP}),(\ref{eq:MasterEquation})
must be corrected to have leading O($M^{-1}$) corrections rather
than O($M^{-2}$)~\cite{Sint:1995ch,Sint:2007ug,DallaBrida:2019mqg}.%
\footnote{The list of O($M^{-1}$) counterterms for the more
	   	  general case, including the situation where the
		  effective theory also contains light quarks,
		  can be inferred from the O($a$) counterterms
		  of the Symanzik effective theory for the
		  SF~\cite{Luscher:1996sc}.
		  It should be clear that the analogy only refers to
		  the fields entering the counterterms and, in general,
		  it does not extend to their coefficient functions
		  $\omega_b$. In particular, one expects different
		  bases of counterterm fields to appear at higher orders
		  in the effective theories, as the two theories have
		  different defining symmetries. From these observations,
		  one concludes that also finite-volume couplings based
		  on open-SF~\cite{Luscher:2014kea} or
		  open~\cite{Luscher:2011kk} boundary conditions are
		  affected by O($M^{-1}$) corrections.}
On the other hand, finite-volume couplings defined through some
variant of periodic boundary conditions, as for instance
regular periodic~\cite{Fodor:2012td} or
twisted~\cite{deDivitiis:1994yz,Ramos:2014kla} boundary
conditions, have leading decoupling corrections of O($M^{-2}$),
as observables in infinite space-time.

There are several possibilities to deal with the issue
of O($M^{-1}$) corrections in SF-based couplings. A
straightforward one is to consider a sufficiently large time
extent $T$, and take for the observable $\Obs$ that defines
the couplings fields which are localized in the middle of the
space-time manifold, i.e. at $x_0=T/2$. This maximizes the
distance from the boundaries and therefore the correlation
between $\Obs$ and the fields responsible for O($M^{-1}$) effects.
In fact, at low energy (i.e. relatively large physical $L$ and $T$)
the O($M^{-1}$) contaminations  are expected to be exponentially
suppressed with the distance of $\Obs$ from the boundaries.

More elegant solutions have been proposed which eliminate
entirely the issue. For example, one could consider for the
heavy quarks a twisted rather than a standard mass~\cite{Frezzotti:2000nk,Frezzotti:2001ea,DellaMorte:2001ys,Pena:2004gb}.
In this case, one can show that $\mathcal{L}_1=0$ and the
decoupling in the SF is realized with O($M_{\rm tw}^{-2}$)
corrections, where $M_{\rm tw}$ is the heavy twisted mass of
the quarks~\cite{Sint:2007ug}. Equivalent in the
continuum is  the situation where the heavy quarks have a
standard mass but the SF boundary conditions are chirally
rotated~\cite{Sint:2005qz,Sint:2007ug,Sint:2010eh}
(see also
refs.~\cite{Sint:2007zz,Sint:2010xy,Lopez:2012as,Lopez:2012mc,Brida:2016rmy,DallaBrida:2018tpn}). The issue with these solutions is that they
require an even number of flavors $\Nf$. They are therefore a
promising approach for the case of $\Nf=4$ QCD. For $\Nf=3$,
one may consider having a doublet of twisted-mass quarks
and a regular massive quark (or equivalently a doublet of
chirally rotated quarks and a regular SF quark
(see e.g.~ref.~\cite{DallaBrida:2018tpn})).
This would reduce the O($M^{-1}$) contributions
to those of only a single flavor.

Another possibility that we want to mention, which is valid
for any choice of $\Nf$, is to match the effective and
fundamental theory at O($M^{-1}$)~\cite{Sommer:2020}. In other
words, by equating the results for some convenient observable
in the effective and fundamental theory one can determine the
coefficient $\omega_b$ appearing in eq.~(\ref{eq:OMTermsSF}).
Once this is determined, the O($M^{-1}$) terms can be taken
into account in the effective theory by computing the insertion
of the counterterm in eq.~(\ref{eq:OMTermsSF}) in the relevant
observables. This guarantees that the decoupling corrections
are of O($M^{-2}$).

\subsection{$\Lambda_{\overline{\rm MS}}^{(3)}$ from
			the decoupling of heavy quarks}
\label{subsec:CouplingFromDecoupling}

\subsubsection{Definitions}

In this subsection we present the results of ref.~\cite{DallaBrida:2019mqg}
where the master formula, eq.~(\ref{eq:MasterEquation}),
was first applied for the computation of
$\Lambda_{\overline{\rm MS}}^{(3)}$. The study considers
$\Nf=3$ QCD which is set on the lattice in terms of
non-perturbatively O($a$)-improved Wilson quarks and the
tree-level Symanzik O($a^2$)-improved gauge
action~\cite{Bulava:2013cta}.%
\footnote{For ease of presentation in the following we
		  use continuum notation to introduce the relevant
		  definitions while referring to the literature
		  for the corresponding lattice expressions.}

The theory is defined in a finite volume with time extent
$T$ and spatial size $L$. The quark fields $\psi,\psibar$
and the gauge field $A_\mu$ satisfy Dirichlet boundary
conditions in the time direction, specifically
($P_\pm\equiv\frac{1}{2}(1\pm\gamma_0)$)~\cite{Luscher:1992an,Sint:1993un,Luscher:1996sc},
\begin{gather}
	\nonumber
	P_+\psi(x)|_{x_0=0}=0
	=
	P_-\psi(x)|_{x_0=T}\,,\\[1ex]
	\nonumber
	\psibar(x)P_-|_{x_0=0}=0=
	\psibar(x)P_+|_{x_0=T}\,,\\[1ex]
	\label{eq:SFbc}
	A_{k}(x)|_{x_0=0,T}=0\,,
\end{gather}
$k=1,2,3$, while in the spatial directions we have
\begin{gather}
	\nonumber
	\psi(x+\hat{k}L)=e^{\frac{i}{2}}\psi(x)\,,
	\quad
	\psibar(x+\hat{k}L)=\psibar(x)e^{-\frac{i}{2}}\,,\\
	\label{eq:SFbc2}
	A_\mu(x+\hat{k}L)=A_\mu(x)\,.
\end{gather}
The finite-volume couplings that we consider in the following
are constructed in terms of the gradient flow field
$B_\mu(t,x)$ which is defined by the
equations~\cite{Narayanan:2006rf,Luscher:2010iy},
\begin{equation}
	\label{eq:GradienFlowEq}
	\partial_t B_\mu(t,x)=D_\nu G_{\nu\mu}(t,x)\,,
	\quad
	B_\mu(0,x)=A_\mu(x)\,,
\end{equation}
where $t>0$ is the flow time and
\begin{equation}
	G_{\mu\nu}=\partial_\mu B_\nu - \partial_\nu B_\mu
	+[B_\mu,B_\nu]\,,
\end{equation}
is the flow-field strength tensor. On the lattice several
discretizations of the GF equations have been considered,
here we employ the Symanzik O($a^2$)-improved definition
proposed in ref.~\cite{Ramos:2015baa}, also known as
Zeuthen flow.

Gauge-invariant composite fields made out of the flow field
$B_\mu(t,x)$ are automatically finite~\cite{Luscher:2011bx},
and are thus ideal quantities to define renormalized
couplings~\cite{Luscher:2010iy,Fodor:2012td,Fritzsch:2013je}.
The definition to which we apply decoupling
(cf.~eq.~(\ref{eq:CouplingDecouplingNP})) is the massive
finite-volume scheme given by
\begin{equation}
	\label{eq:GFTCoupling}
	[\bar{g}^{(3)}_{\rm GFT}(\mu,M)]^2\equiv
	{t^2\over\mathcal{N}}
	\langle E_{\rm mag}(t,x)\rangle_{{\rm SF},Q=0}\big|_{T=2L}^{x_0=L,\mu=L^{-1},\sqrt{8t}=cL}\,,
\end{equation}
where $\mathcal{N}$ is a normalization
constant~\cite{Fritzsch:2013je} and
\begin{equation}
	\label{eq:Emag}
	E_{\rm mag}(t,x)=-\frac{1}{2}\tr\{G_{kl}(t,x)G_{kl}(t,x)\}
\end{equation}
is the magnetic component of the energy density of the
flow field $B_\mu(t,x)$. On the lattice, we define $E_{\rm mag}$
in terms of the O($a^2$)-improved definition given in
ref.~\cite{Ramos:2015baa}.

Note how in eq.~(\ref{eq:GFTCoupling}) the renormalization
scale $\mu$ is set in terms of the finite spatial extent
$L$, as appropriate for a finite-volume coupling. In order
for the coupling to depend on a single scale (apart from $M$),
the flow time $t$ is also linked to $L$ through the constant
$c$~\cite{Fodor:2012td,Fritzsch:2013je}.
The value of this constant is in principle arbitrary, but
experience suggests that $c=0.3$ is a good compromise
between statistical precision and discretization
errors~\cite{Fritzsch:2013je}. From the point of view of
decoupling, we note that the larger the value of $c$ is,
the more sensitive the coupling is to the O($M^{-1}$)
counterterms located at $x_0=0,T$. This is so because for
larger values of $t$ the footprint of the flow field
$B_\mu(t,x)$ extends closer to the boundaries.

In order to attenuate the sensitivity to the O($M^{-1}$)
terms, in eq.~(\ref{eq:GFTCoupling}) we consider a
space-time volume with $T=2L$, and place the energy
density $E_{\rm mag}(t,x)$  at $x_0=T/2$, in order to
maximize the distance from the boundaries
(cf.~Sect.~\ref{subsec:FiniteVolumeDecoupling}).
Taking only the magnetic part of the flow energy
density also helps in reducing the O($M^{-1}$) contaminations.%
\footnote{These choices made for reducing
		  O($M^{-1}$) effects are also effective in
		  reducing the sensitivity to the O($a$)
		  effects stemming from the space-time
		  boundaries
		  (cf.~Sect.~\ref{subsec:FiniteVolumeDecoupling}).}
Lastly, we note that the expectation value
$\langle\cdots\rangle_{\rm SF,Q=0}$ in eq.~(\ref{eq:GFTCoupling})
is meant to be considered in the presence of the SF boundary
conditions, eqs.~(\ref{eq:SFbc})-(\ref{eq:SFbc2}),
and restricted to gauge fields in the trivial
topological sector~\cite{Fritzsch:2013yxa,DallaBrida:2016kgh}.
The latter constraint is imposed in order to circumvent
issues related to topology freezing at small lattice
spacings~\cite{DelDebbio:2002xa,Schaefer:2010hu}.

Having introduced the massive scheme of choice, we now
move to the definition of the decoupling scale
$\mu^{(3)}_{\rm dec}$. We define this in terms of a
\emph{massless} finite-volume coupling. Its definition slightly
differs from that of eq.~(\ref{eq:GFTCoupling}). Specifically,
we take~\cite{Fritzsch:2013je,DallaBrida:2016kgh},
\begin{equation}
	\label{eq:GFcoupling}
	[\bar{g}^{(3)}_{\rm GF}(\mu)]^2\equiv
	{t^2\over\mathcal{N}'}
	\langle E_{\rm mag}(t,x)\rangle_{{\rm SF},Q=0}\big|_{T=L,M=0}^{x_0=L/2,\mu=L^{-1},\sqrt{8t}=cL}\,,
\end{equation}
where as before we set $c=0.3$. The main differences
with respect to the definition in eq.~(\ref{eq:GFTCoupling})
are that the coupling is evaluated at vanishing (renormalized)
quark masses and the temporal extent is shorter, i.e. $T=L$.
The reason for considering this specific definition is because
its non-perturbative running in the $\Nf=3$ theory is known
very precisely in the range of scales $\mu\approx0.2-4\,\GeV$
(see ref.~\cite{DallaBrida:2016kgh} and
Sect.~\ref{subsec:LambdaNf3}). This gives us the freedom to
choose for $\mu^{(3)}_{\rm dec}$ the most convenient value
in this range. Its physical units can in fact be inferred
from combining the knowledge of the non-perturbative
$\beta$-function of
$\bar{g}^{(3)}_{\rm GF}(\mu)$~\cite{DallaBrida:2016kgh} and
the physical scales $\mu_{\rm phys}^{(3)}$ determined in
large-volume hadronic
simulations~\cite{Bruno:2016plf,Bruno:2017gxd}.%
\footnote{The latter are given by a combination of pion and
		  kaon decay constants
		  (cf.~refs.~\cite{Bruno:2016plf,Bruno:2017gxd}
		  and Sect.~\ref{subsec:LambdaNf3} for more
		  details).}

In this study, the decoupling scale is specified
by the condition
\begin{equation}
	\label{eq:LCP}
	[\bar{g}^{(3)}_{\rm GF}(\mu^{(3)}_{\rm dec})]^2=
	3.95\equiv u_0\,,
\end{equation}
which using the information mentioned above is found
to correspond to
\begin{equation}
	\label{eq:mudecPhysUnits}
	\mu_{\rm dec}^{(3)}=789(15)\,\MeV\equiv L_{\rm dec}^{-1}.
\end{equation}
This scale is convenient in practice as, given our choice
of lattice resolutions (see below), it allows us to simulate
at values of the lattice spacing sensibly smaller than
those typically accessible to large-volume simulations. As we
shall see, this enables us to simulate quark masses up to
$M\approx 4M_c$, while having $aM$ effects under good
control. Additionally, we can profit from some perturbative
information in the treatment of O($aM$) effects, which is
expected to be accurate enough at the values of the bare
coupling corresponding to the relevant lattice spacings.
Lastly, the value of $\mu_{\rm dec}^{(3)}$ is low enough in
energy that only a very limited part of the running in $\Nf=3$
QCD is needed in order to connect it to the scales
$\mu^{(3)}_{\rm phys}$ and set its physical units.

\subsubsection{Determinations of the massive coupling}

The next step in the strategy is to determine the value of the
coupling $\bar{g}^{(3)}_{\rm GFT}(\mu^{(3)}_{\rm dec},M)$
for some large quark masses. To this end, we must evaluate the
coupling $\bar{g}^{(3)}_{\rm GFT}(\mu^{(3)}_{\rm dec},M)$ for
several values of the lattice spacing at fixed $\mu^{(3)}_{\rm dec}$
and given $M$, and extrapolate it to the continuum limit.

The condition in eq.~(\ref{eq:LCP}) defines the line of
constant physics along which $\mu_{\rm dec}^{(3)}$ is
constant. Given a set of lattice sizes $L/a$, by
tuning the bare coupling $g_0$ so that the massless GF-coupling
has the prescribed value $u_0$, we can identify the values
of $a$ for which $L=L_{\rm dec}$ is fixed in physical units
as $a/L\to0$. In this respect, note that in order to compute
the massive coupling $\bar{g}^{(3)}_{\rm GFT}(\mu^{(3)}_{\rm dec},M)$ and the massless coupling
$\bar{g}^{(3)}_{\rm GF}(\mu^{(3)}_{\rm dec})$ at matching values
of the lattice spacing up to O($(aM)^2$) corrections, the two
must be evaluated at the same value of the O($a$)-improved bare
coupling $\tilde{g}_0$~\cite{Luscher:1996sc}.
The latter, we recall, is defined as~\cite{Luscher:1996sc},
\begin{equation}
	\label{eq:g0tilde}
	\tilde{g}_0^2\equiv g_0^2(1+b_{\rm g}(\tilde{g}_0)am_{\rm q})\,,
	\quad
	m_{\rm q}=m_0-m_{\rm cr}(\tilde{g}_0)\,,
\end{equation}
where $m_0$ is the bare quark mass, $m_{\rm cr}$ is its critical
value at which the (renormalized) quark masses vanish, and
$b_{\rm g}(\tilde{g}_0)$ is a function of the bare coupling
to be determined. In the massless theory $m_{\rm q}=0$, and the
improved bare coupling coincides with $g_0$. The values of
$g_0$ determined from the condition (\ref{eq:LCP}) in terms of
the massless coupling therefore specify the values of $\tilde{g}_0$
at which the massive couplings should be evaluated. According to
the Symanzik improvement programme, the coefficient
$b_{\rm g}(\tilde{g}_0)$ can be tuned in order to remove
O($am_{\rm q}$) effects in the matching between the massless and
massive renormalization schemes. At present, however, this is only
known to 1-loop order in lattice perturbation theory, where
$b_{\rm g}(g_0)=0.012\,g_0^2\times \Nf+{\rm O}(g_0^4)$~\cite{Sint:1995ch,Luscher:1996sc}.

\begin{figure*}[h]
	\centering

	\resizebox{0.875\textwidth}{!}{%
		\includegraphics{plots/plot_z.pdf}}
	\caption{Continuum limit extrapolations of the massive
	couplings $\Psi^M(u_0,z)$ for $z=1.972,4.0,6.0,8.0$
	(cf.~eq~(\ref{eq:PsiM}))~\cite{DallaBrida:2019mqg}.
	Two cuts $(aM)^2 < 1/8,1/4$ are applied in order to estimate
	the systematic uncertainties in the extrapolations.}
	\label{fig:MassiveCouplings}
\end{figure*}

Together with the bare coupling, the bare quark masses
must be set in order to guarantee a given value for the RGI mass
$M$ in the continuum limit. This is achieved by considering
a value of $z=ML_{\rm dec}$, and by solving for a given
lattice size $L_{\rm dec}/a$ the following
equation for the bare quark masses,
\begin{equation}
	\label{eq:m0}
	z={L_{\rm dec}\over a} {M\over
	\overline{m}^{(3)}_{\rm SF}(\mu^{(3)}_{\rm dec})}
	Z_{\rm m}(\tilde{g}_0,a/L_{\rm dec})a\tilde{m}_{\rm q}\,.
\end{equation}
In this equation,
\begin{equation}
	\tilde{m}_{\rm q}=m_{\rm q}
	(1+b_{\rm m}(\tilde{g}_0)am_{\rm q})\,,
	\quad
	m_{\rm q}=m_0-m_{\rm cr}(\tilde{g}_0)\,,
\end{equation}
is the O($a$)-improved definition for the bare quark mass,
which replaces the regular bare mass $m_{\rm q}$ in order to
eliminate O($am_{\rm q}$) effects in massive
schemes~\cite{Luscher:1996sc}. To this end, the function
$b_{\rm m}(\tilde{g}_0)$ must be properly chosen. The function
$Z_m(\tilde{g}_0,a/L)$, instead, refers to the renormalization
constant that relates the bare quark mass to the
renormalized quark mass $\overline{m}_{\rm SF}(\mu)$ in
the SF-scheme of
refs.~\cite{Sint:1998iq,Capitani:1998mq,DellaMorte:2005kg,Campos:2018ahf}. This, together with the improvement coefficient $b_{\rm m}$,
and the critical mass $m_{\rm cr}$, are known non-perturbatively
for the relevant parameters (see ref.~\cite{DallaBrida:2019mqg}).
The matching factor
$M/\overline{m}^{(3)}_{\rm SF}(\mu^{(3)}_{\rm dec})$ then
allows us to convert the renormalized quark mass
in the SF-scheme at the scale $\mu_{\rm dec}^{(3)}$ to the RGI
mass $M$. It can be obtained from the results of
ref.~\cite{Campos:2018ahf}. Once $am_{\rm q}$ for the given $z$
is known from eq.~(\ref{eq:m0}) at the values of $\tilde{g}_0$
given by the condition eq.~(\ref{eq:LCP}), using the 1-loop
results for $b_{\rm g}$ we can infer from eq.~(\ref{eq:g0tilde})
the values of $g_0$ at which the massive couplings must be
computed in simulations~\cite{DallaBrida:2019mqg}.

Having set the bare parameters we can finally evaluate the
functions,
\begin{equation}
	\label{eq:PsiM}
	\Psi^M(u_0,z)=\lim_{a/L_{\rm dec}\to0}
	\big[\bar{g}^{(3)}_{\rm GFT}(\mu^{(3)}_{\rm dec},M)\big]^2\big|_{[\bar{g}^{(3)}_{\rm GF}(\mu^{(3)}_{\rm dec})]^2=u_0}\,.
\end{equation}
In figure \ref{fig:MassiveCouplings} the results for
the extrapolations in eq.~(\ref{eq:PsiM}) are shown.
Several values of $z$ have been considered, ranging
from $z\approx2-8$, which correspond to RGI masses
$M\approx 1.6-6.3\,\GeV$. The different values of $z$
will allow us to assess the size of the non-perturbative
corrections to decoupling in eq.~(\ref{eq:MasterEquation}).
The range of lattice sizes considered is $L/a=12-32$.

As one can see from the figure, as expected, the continuum limit
extrapolations become more challenging as $z$ becomes larger.
However, at the smaller values of $a/L$, the data seem to
be well described by O($a^2$) discretization errors.
In order to assess systematic effects in the extrapolations,
fits with different cuts in $aM$ have been considered,
specifically $(aM)^2<1/8,1/4$. The results from the different
fits are compatible, with the results for $(aM)^2 < 1/8$
having significantly larger errors at large values of $z$,
where fewer points are left after the cut is imposed.
We take as final results those with cut $(aM)^2<1/8$.

\begin{figure*}[h]
	\centering

	\resizebox{0.875\textwidth}{!}{%
	\includegraphics{plots/final_rho.pdf}}
	\caption{Values for $\rho=\Lambda^{(3)}_{\overline{\rm MS}}/
	\mu_{\rm dec}^{(3)}$ determined from the decoupling relation,
	eq.~(\ref{eq:MasterEquation})~\cite{DallaBrida:2019mqg}.
	As $z=M/\mu^{(3)}_{\rm dec}$ gets larger, the
	results for $\rho$ approach the FLAG value
	$\Lambda^{(3)}_{\overline{\rm MS}}=343(12)\,\MeV$~\cite{Aoki:2019cca}
	in units of $\mu^{(3)}_{\rm dec}=789(15)\,\MeV$
	(cf.~eq.~(\ref{eq:mudecPhysUnits})). The filled symbols illustrate
	possible extrapolations for $z\to\infty$. The results of the
	extrapolations show significantly smaller statistical
	uncertainties than the FLAG result. (Note that
	$\mu_{\rm dec}^{(3)}$ gives a negligible contribution to
	the uncertainty of $\Lambda_{\overline{\rm MS}}^{(3)}/
	\mu_{\rm dec}^{(3)}$ from FLAG.) The uncertainties on
	$\rho$ may be further reduced with modest computational
	effort by improving the determination of
	$\Lambda^{(0)}_{\overline{\rm MS}}/\mu_{\rm dec}^{(0)}$
	from ref.~\cite{DallaBrida:2019wur}.
	The data both at finite $z$ and the extrapolations for
	$z\to\infty$ have been slightly shifted horizontally for
	better clarity.}
	\label{fig:LambdaFromDecoupling}
\end{figure*}

\subsubsection{Results for $\Lambda_{\overline{\rm MS}}^{(3)}$}

A last step separates us from applying eq.~(\ref{eq:MasterEquation})
to compute $\Lambda_{\overline{\rm MS}}^{(3)}$. In order
to use eq.~(\ref{eq:MasterEquation}) directly in the
GFT-scheme of eq.~(\ref{eq:GFTCoupling}) the non-perturbative
running of the corresponding coupling in the pure-gauge
theory, $\bar{g}^{(0)}_{\rm GFT}(\mu)$, should be known. This,
however, has never been computed. On the other hand, the
running of the pure-gauge coupling in the GF-scheme,
$\bar{g}^{(0)}_{\rm GF}(\mu)$, is known very
precisely~\cite{DallaBrida:2019wur}. In other to resolve the
issue, all we have to do is to match non-perturbatively the
GFT- and GF-schemes in the pure Yang-Mills theory. More precisely,
we need to determine the values of the coupling
$g_M=\bar{g}^{(0)}_{\rm GF}(\mu^{(0)}_{\rm dec})$
corresponding to $\bar{g}^{(0)}_{\rm GFT}(\mu^{(0)}_{\rm dec})=
\sqrt{\Psi^M(u_0,z)}$ for the relevant values of $z$.
Given this relation we can compute,
\begin{equation}
	{\Lambda^{(0)}_{\rm GF}\over\mu_{\rm dec}^{(0)}}
	=\varphi_{\rm g,\rm GF}^{(0)}(g_M)
	=\varphi_{\rm g,\rm GFT}^{(0)}(\sqrt{\Psi^M})\,,
	\quad
	g_M=\sqrt{\chi(\Psi^{M})}\,,
\end{equation}
where
$\chi([\bar{g}^{(0)}_{\rm GFT}(\mu)]^2)=
[\bar{g}^{(0)}_{\rm GF}(\mu)]^2$.
The function $\chi$ can easily be obtained in the relevant
range of couplings $[\bar{g}^{(0)}_{\rm GFT}]^2=\Psi^M(u_0,z)$,
by computing $[\bar{g}^{(0)}_{\rm GF}]^2$ and
$[\bar{g}^{(0)}_{\rm GFT}]^2$ at several matching
values of $L/a$ and $g_0$ in this range, and extrapolating
their relation to the continuum limit~\cite{DallaBrida:2019mqg}.

The results for $\rho=\Lambda_{\overline{\rm MS}}^{(3)}/\mu^{(3)}_{\rm dec}$
as obtained from eq.~(\ref{eq:MasterEquation}) for different
values of $z$ are shown in figure
\ref{fig:LambdaFromDecoupling}. For the estimates the function
$P^{\rm PT}_{0,3}(M/\Lambda_{\overline{\rm MS}}^{(3)})$
is evaluated at 5-loop order. In this respect we note that,
as expected from the discussion in
Sect.~\ref{subsubsec:DecouplingPT}, the perturbative
uncertainties in $\rho$  estimated from the effect of the last
known terms of $P^{\rm PT}_{0,3}$ are completely negligible
compared to the other sources of uncertainties (cf. Table 2
of ref.~\cite{DallaBrida:2019mqg}). As one can see from the
figure, excluding the point at $z\approx2$, the non-perturbative
corrections to decoupling are  small. At larger $z$ values they
are compatible with O($z^{-2}$) scaling, indicating that the
O($z^{-1}$) corrections due to the SF boundary conditions are
subdominant. For values of $M\approx 6.3\,\GeV$ (i.e. $z=8$)
the estimated $\rho$ agrees well with the fully $\Nf=3$ flavor
theory results for $\Lambda_{\overline{\rm MS}}^{(3)}/
\mu^{(3)}_{\rm dec}$, where $\Lambda_{\overline{\rm MS}}^{(3)}$
is given by the FLAG average value~\cite{Aoki:2019cca} and
$\mu_{\rm dec}^{(3)}$ is taken from eq.~(\ref{eq:mudecPhysUnits}).
If one attempts a $z\to\infty$ extrapolation of the data the
agreement becomes even better.

\subsection{Summary and miscellaneous remarks}
\label{subsec:Summary}

The results of ref.~\cite{DallaBrida:2019mqg} put on solid
grounds the application of the decoupling relation
eq.~(\ref{eq:MasterEquation}) as a novel strategy to determine
the QCD coupling from lattice QCD. The remarkable feature of
this approach is that the non-perturbative running of the
coupling from the low-energy scale $\mu_{\rm dec}$ up to
high energy is done entirely in the pure-gauge theory. This
opens up the possibility to significantly reduce the current
error on $\alphas$ (cf.~Sects.~\ref{subsec:LambdaNf3} and
\ref{sec:Conclusions}).

In order to translate the results for
$\Lambda_{\overline{\rm MS}}^{(0)}/\mu^{(0)}_{\rm dec}$
into $\Lambda_{\overline{\rm MS}}^{(\Nf)}/\mu^{(\Nf)}_{\rm dec}$,
the strategy relies on two crucial ingredients. The first
ingredient is the computation of a massive coupling
$\bar{g}^{(\Nf)}_{\Obs}(\mu,M)$ at the low-energy
scale $\mu^{(\Nf)}_{\rm dec}$ in an unphysical set-up with
$\Nf$-flavors of degenerate massive quarks of mass
$M\gg\mu^{(\Nf)}_{\rm dec}$. Exploiting the decoupling of the
massive quarks, the scales $\mu^{(\Nf)}_{\rm dec}$ and
$\mu^{(0)}_{\rm dec}$ can be connected through the massive
coupling and so the fundamental $\Nf$-flavor and effective
pure-gauge theories. The second ingredient is the use of
high-order perturbation theory for estimating the ratio of
the $\Lambda$-parameters of the two theories,
$\Lambda_{\overline{\rm MS}}^{(0)}/\Lambda_{\overline{\rm MS}}^{(\Nf)}$, given by the function
$P_{\rm 0,f}(M/\Lambda_{\overline{\rm MS}}^{(\Nf)})$.
Control on the determination of the massive coupling can
be achieved by employing a suitable finite-volume scheme.
Thanks to the fact that the physical volume does not need
to be large, small lattice spacings can be simulated, and
safe continuum limit extrapolations of
$\bar{g}^{(\Nf)}_{\Obs}(\mu_{\rm dec}^{(\Nf)},M)$
with $\mu^{(\Nf)}_{\rm dec}={\rm O}(1\,\GeV)$ can be taken
for quark masses up to a few GeV. At these large masses,
perturbation theory for the function
$P_{\rm 0,f}(M/\Lambda_{\overline{\rm MS}}^{(\Nf)})$
works extremely well and non-perturbative O($M^{-2}$) corrections
to the decoupling relations are found to be small.

From the results of figure \ref{fig:LambdaFromDecoupling}
we can appreciate how the strategy based on decoupling promises
great accuracy. The $z\to\infty$ extrapolations give in fact
results for $\rho$ which are about a factor 2 more precise
than those obtained using the current FLAG estimate for
$\Lambda_{\overline{\rm MS}}^{(3)}$ and $\mu_{\rm dec}^{(3)}$
from eq.~(\ref{eq:mudecPhysUnits}). In order to set this result
on firmer grounds, however, a robust $z\to\infty$ extrapolation
must be performed. To this end, it is important that the
continuum limit extrapolations for the massive couplings
$\Psi^M(u_0,z)$ are made more solid at the largest (most relevant)
$z$ values by investing some additional computational effort.
The lattices used for the computations entering figure
\ref{fig:MassiveCouplings} are in fact rather modest. The
largest simulated lattices have $L/a=32$, $T/a=64$. In addition,
the large quark masses make these simulations significantly
cheaper than the more common massless SF simulations. A factor
two finer lattices are hence  within reach with affordable
computational resources. These lattices will allow us to consider
larger quark masses too, and thus improve even further
the control on the $z\to\infty$ extrapolations. All in all, we
can expect that after these improvements the final determination
for $z\to\infty$ will be at least as precise as the results in
figure \ref{fig:LambdaFromDecoupling} promise, but will include
conservative estimates for all systematics.

It is important to note at this point that a significant
fraction of the error on $\rho$ at finite $z$ comes from
the uncertainty on $\Lambda_{\overline{\rm MS}}^{(0)}/
\mu_{\rm dec}^{(0)}$ from ref.~\cite{DallaBrida:2019wur},
which is about $1.5\%$ (cf.~figure \ref{fig:LambdaFromDecoupling}).
A reduction of this error down to $0.5\%$ or so is desirable
and in principle possible. Given the importance of the result
for the determination of $\Lambda_{\overline{\rm MS}}^{(3)}$,
however, it is crucial that this error reduction is achieved
robustly. As discussed in Sect.~\ref{subsubsec:YangMillsLambda},
even though the determination of $\Lambda$ in the pure Yang-Mills
theory is very much simplified from the computational point
of view compared to QCD, the problem is yet non-trivial and care
must be taken, especially if such a high precision is desired.
For this reason, it is mandatory that
the results for $\Lambda_{\overline{\rm MS}}^{(0)}$ are
corroborated by investigating different strategies where the
estimates of systematic uncertainties are put to a stringent
test. As we have seen, there is currently tension among different
determinations of $\Lambda_{\overline{\rm MS}}^{(0)}$,
some of which quote the desired sub-percent precision
(cf.~Sect.~\ref{subsubsec:YangMillsLambda} and
ref.~\cite{DallaBrida:2019wur}). Studies as the ones of
refs.~\cite{Husung:2020pxg,Nada:2020jay} are hence
encouraged in order to set the actual accuracy at which we
currently know $\Lambda_{\overline{\rm MS}}^{(0)}$.

This corroboration goes hand in hand with the exploration
of new strategies for the determination of
$\Lambda^{(0)}_{\overline{\rm MS}}$. In this respect we point
out the results of ref.~\cite{Nada:2020jay}, where an alternative
way to do step-scaling for GF-based couplings was proposed
and tested. In short, the change of renormalization scale in the
coupling is first achieved by changing the flow time at fixed
physical volume and in a second step the physical volume is
changed at fixed flow time. This has to be compared with the
traditional situation where both flow time and spatial size are
changed at once. One of the interesting features of the approach
is that it amounts to a reanalysis of data gathered from a
traditional step-scaling study. However, the systematics to deal
with are quite independent given the different continuum limit
extrapolations involved. By comparing the two analysis one
can stringently test the assumptions made in one or the other
approach.

Furthermore, it would be interesting to employ
other definitions of finite-volume schemes based on either
different observables and/or set-ups. For instance, the
GF-coupling with twisted boundary conditions explored in
refs.~\cite{Ramos:2014kla,Lin:2015zpa,Ishikawa:2017xam,Bribian:2020xfs} is promising. Differently from the SF case it enjoys
full translational invariance and yet its perturbative
expansion in finite volume appears feasible~\cite{Bribian:2019ybc}.
Of course, standard periodic boundary conditions are also an
option~\cite{Fodor:2012td} despite the difficulties with
perturbation theory in this set-up. In fact, in cases
where the perturbative information is limited or the relevant
perturbative expansion is poorly convergent, a viable option
is to non-perturbatively match the given finite-volume scheme
to some other scheme for which the perturbative $\beta$-function
is known to high-loop order and it is well behaved
(see e.g.~refs.~\cite{DallaBrida:2016kgh,DallaBrida:2019wur}).
This may allow for a significantly more precise determination
of $\Lambda_{\overline{\rm MS}}^{(0)}$
(cf.~Sect.~\ref{subsubsec:YangMillsLambda}). In this respect,
we note that a powerful framework for automated
numerical high-loop calculations in finite volume has been
recently developed and successfully
applied~\cite{DallaBrida:2017tru,DallaBrida:2017pex}.

Another idea that may be interesting to consider is the
determination of $\Lambda^{(0)}_{\overline{\rm MS}}$ based
on the infinite-volume $\beta$-function of GF-based couplings,
following the strategy of
refs.~\cite{Fodor:2017die,Hasenfratz:2019hpg,Fodor:2019ypi}
(see also ref.~\cite{Luscher:2014kea}). In this approach, the
infinite-volume results are obtained by extrapolations from
small-volume simulations, which might already be at hand from
a conventional step-scaling study. If the (non-trivial)
infinite-volume extrapolations can be performed in a controlled
way and the convergence to the perturbative regime of the chosen
scheme is fast enough, this strategy may allow for interesting
crosschecks of the results from step-scaling. In this case, the
framework developed in ref.~\cite{Artz:2019bpr} can be used to
obtain the necessary infinite-volume perturbative information
to high-loop order.

Besides applying different strategies for the calculation of
$\Lambda_{\overline{\rm MS}}^{(0)}$, different techniques
can be considered for the QCD part of the decoupling strategy
as well. A simple extension is to consider different schemes
for the massive finite-volume coupling. In the case of $\Nf=3$
QCD, periodic and twisted boundary conditions would avoid
entirely the issue with O($M^{-1}$) contaminations.
For $\Nf=4$ QCD twisted boundary conditions cannot be
implemented~\cite{Parisi:1984cy}, but twisted-mass fermions
with SF boundary conditions or regular massive quarks with
chirally rotated boundary conditions are available options
(cf.~Sect.~\ref{subsec:FiniteVolumeDecoupling}).

A substantially different approach is to avoid entirely
finite-volume couplings and rely on heavy-quark decoupling
in hadronic quantities. Particularly interesting observables
to consider are the popular gluonic scales
$\mathcal{S}=t_0^{-1/2}, t_c^{-1/2},w_0^{-1},
r_0$~\cite{Luscher:2010iy,Athenodorou:2018wpk,Borsanyi:2012zs,Sommer:1993ce}.
In this case, the decoupling relation is applied more directly in
the form of eq.~(\ref{eq:MatchingNPT}), specifically,
\begin{equation}
	{\Lambda^{(\Nf)}_{\overline{\rm MS}}\over
	\mathcal{S}^{(\Nf)}(M)}
	P^{\rm PT}_{0,\rm f}\bigg({M\over
	{\Lambda^{(\Nf)}_{\overline{\rm MS}}}} \bigg)
	=
	{\Lambda^{(0)}_{\overline{\rm MS}}\over \mathcal{S}^{(0)}}\,
	+{\rm O}(g_*^{2n-2},M^{-2}).
\end{equation}
Once the right hand side is known from computing the running
in pure-gauge theory, the perturbative approximation
$P^{\rm PT}_{0,\rm f}({M/\Lambda^{(\Nf)}_{\overline{\rm MS}}})$
to the ratio of $\Lambda$-parameters can be used to solve
the above equation for
$\Lambda_{\overline{\rm MS}}^{(\Nf)}/\mathcal{S}^{(\Nf)}(M)$.
All that is left to do to determine $\Lambda_{\overline{\rm MS}}^{(\Nf)}$ is then to fix the physical units of the
low-energy quantity $\mathcal{S}^{(\Nf)}(M)$ computed in a
theory with $\Nf$ heavy quarks of RGI mass $M$. This can be
obtained by relating $\mathcal{S}^{(\Nf)}(M)$ to some
convenient physical scale $\mu_{\rm phys}^{(\Nf)}$ evaluated
at physical quark masses via the bare parameters. As discussed
in Sect.~\ref{subsubsec:MultiScaleProblem2}, it might be
difficult to reach large masses $M$ with this approach, while
having discretization errors and finite-volume effects under
control. Some compromises are likely necessary in order to
reach high enough $M$ values to be able to control decoupling
corrections. On the other hand,  the studies of
refs.~\cite{Knechtli:2017xgy,Athenodorou:2018wpk}
show that interesting results may be obtained if masses
close to that of the charm can be reliably reached. This
makes the strategy worth being explored.

\section{Conclusions}
\label{sec:Conclusions}

In this contribution we presented a novel strategy for
the determination of the QCD coupling using lattice QCD.
It exploits the decoupling of heavy quarks at low energy
to connect the pure Yang-Mills theory and
QCD with $\Nf$ flavors of quarks. The main result is that
the computation of the running of the coupling from a
known low-energy scale $\mu_{\rm dec}={\rm O}(1\,\GeV)$
up to high energies can be done entirely in the pure-gauge
theory instead of $\Nf$-flavor QCD. Considering $\Nf=3$ or 4,
this paves the way for unprecedented precision determinations
of $\Lambda_{\overline{\rm MS}}^{(\Nf)}$ from which
$\Lambda_{\overline{\rm MS}}^{(5)}$ and $\alphas$
can be obtained. In ref.~\cite{DallaBrida:2019mqg} the
potential of these methods was successfully established
in the determination of $\Lambda_{\overline{\rm MS}}^{(3)}$.
We now want to put this result into context of a future
precision $\alphas$ extraction.

As presented in Sect.~\ref{subsec:CouplingFromDecoupling},
the results for $\Lambda_{\overline{\rm MS}}^{(3)}/
\mu^{(3)}_{\rm dec}$ from decoupling have an uncertainty
which is about half the one obtained using the FLAG
average $\Lambda_{\overline{\rm MS}}^{(3)}=343(12)\,\MeV$~\cite{Aoki:2019cca}
and $\mu_{\rm dec}^{(3)}=789(15)\,\MeV$ from
eq.~(\ref{eq:mudecPhysUnits})
(cf.~figure	\ref{fig:LambdaFromDecoupling}).
As discussed in Sect.~\ref{subsec:Summary}, by investing some
modest computational effort, this result can be
set on very solid grounds by improving the continuum limits of
the massive couplings (cf.~figure \ref{fig:MassiveCouplings})
and performing a robust $z\to\infty$ extrapolation
(cf.~figure \ref{fig:LambdaFromDecoupling}). Considering lattices
twice as large as the ones simulated in
ref.~\cite{DallaBrida:2019mqg} is in fact affordable.
With such lattices we can expect that the continuum limit
extrapolations of $\Psi^M(u_0,z)$ in eq.~(\ref{eq:PsiM})
can be performed with confidence also at the largest
masses investigated so far ($M\approx6\,\GeV$). In addition,
we will be able to consider some larger $z$ values, too.

The precision on
$\Lambda_{\overline{\rm MS}}^{(3)}/\mu^{(3)}_{\rm dec}$
can be further improved significantly by reducing the
uncertainties coming from the pure-gauge determination of
$\Lambda_{\overline{\rm MS}}^{(0)}/\mu^{(0)}_{\rm dec}$,
which has the (conservative) error of
1.5\%~\cite{DallaBrida:2019wur}. As noticed in
Sect.~\ref{subsec:Summary}, this error can in principle be
reduced by a substantial factor, e.g. down to 0.5\%. However,
while it is certainly possible to reach such a high precision
in a given computation, it is crucial that the
results of different analysis and groups corroborate it.
At present, there is in fact  tension among
determinations of $\Lambda_{\overline{\rm MS}}^{(0)}$
involving results with sub-percent accuracy
(cf.~Sect.~\ref{subsubsec:YangMillsLambda}). It is
important to understand the origin of these
differences. We hope that the renovated interest in this
quantity brought by this new strategy motivates the community
to resolve the issue and contribute to a high-precision
determination.

Once the above steps are achieved the precision on
$\Lambda_{\overline{\rm MS}}^{(3)}$ will be limited
by the present error on $\mu^{(3)}_{\rm dec}$, which is
about 2\%. A reduction of this error down to $1\%$ is
however foreseeable. It requires, first of all, to improve
the results for the running of the GF-coupling
$\bar{g}^{(3)}_{\rm GF}(\mu)$ at energies
$\mu<\mu_{\rm dec}$~\cite{DallaBrida:2016kgh}. We recall
that this is needed in order to connect
$\mu^{(3)}_{\rm dec}$ with the hadronic scales
$\mu^{(3)}_{\rm phys}$ used to set the physical units
of the theory
(cf.~Sect.~\ref{subsec:CouplingFromDecoupling}).
Work in this direction has already started as part of the
HQET efforts by the ALPHA Collaboration
(cf.~ref.~\cite{Fritzsch:2018yag}). Secondly, the scale setting
in terms of the physical scales $\mu_{\rm phys}^{(3)}$ must be
improved as well. In practice, this means to obtain a more precise
determination for a convenient low-energy reference scale in
$\Nf=3$ QCD in physical units,
like for example, $\mu_{\rm ref}^*=1/\sqrt{8t_0^*}$
(cf.~ref.~\cite{Bruno:2016plf}). A precision of $1\%$ or better
on this or similar scales is desirable. This is expected
to be possible by exploiting the new CLS ensembles close to
the physical
point~\cite{Bruno:2014jqa,Bali:2019yiy,Gerardin:2019rua}.
Also in this case, however, corroboration from different
strategies and groups is important.

Through all these steps a determination of
$\Lambda_{\overline{\rm MS}}^{(3)}$ with a final
uncertainty of 1-2\% appears feasible. As discussed in
Sect.~\ref{subsubsec:NonperturbativeDec}, at
this level of precision $\Lambda_{\overline{\rm MS}}^{(5)}$
can yet be obtained from $\Lambda_{\overline{\rm MS}}^{(3)}$
by relying on the perturbative decoupling of the charm quark,
eventually including some conservative estimate for the
unaccounted non-perturbative corrections. A determination of
$\alphas(M_Z)$ at the level of 0.4\% is therefore within reach
thanks to the novel techniques. To further halve the error
on $\alphas(M_Z)$, on the other hand, requires several issues to
be reconsidered. Non-perturbative decoupling effects might
become relevant, and one might need to include electromagnetic
and $m_u\neq m_d$ effects in the lattice computations in order
to set the physical scale of the theory to a greater level of
accuracy (cf.~Sect.~\ref{subsubsec:NonperturbativeDec} and see
also discussion in ref.~\cite{DelDebbio:2021ryq}).

Before concluding we want to note that even though our emphasis
was on the determination of $\Lambda_{\overline{\rm MS}}^{(\Nf)}$,
the ideas presented can be extended to solve other
RG problems. A clear case is that of the quark masses,
where one can replace their running in $\Nf$-flavor QCD
with the one in the quenched approximation. In
ref.~\cite{Cali:2019enm} a similar application was in fact
explored in order to study the non-perturbative charm-quark
effects in the determination of the charm-quark mass itself.
More complicated composite operators, like for instance
four-quark operators, require more thought. First of all,
a study of the quality of their perturbative decoupling relations
is necessary in order to establish whether the strategy has any
chance to be applied in the first place. Then, an investigation
of the non-perturbative decoupling corrections must follow.

In conclusion, we can affirm that the decoupling of
heavy quarks enters at full right in the renormalization
toolkit of the lattice field theorist. Many more
applications of these powerful ideas in lattice QCD and
lattice field theory in general are likely to come.

\section*{Acknowledgments}

I would like to thank the organizers of the Special Issue:
``Lattice Field Theory during the COVID-19 pandemic''
for their kind invitation and the opportunity to share
these results. I am really grateful to all my colleagues for
the rewarding collaboration over the past years.
A special thanks goes to:  Patrick Fritzsch, Roman H\"ollwieser,
Francesco Knechtli, Tomasz Korzec, Alberto Ramos, Stefan Sint,
and Rainer Sommer. In this occasion, I particularly thank
F. Knechtli and R. Sommer for their valuable comments
on this manuscript. Finally, it is with gratitude that I thank
the University of Cyprus, the Cyprus Institute, and the members
of their lattice QCD groups, for their kind hospitality during
the preparation of this work.