We now create a new Core with name b2 by adding the bottom PS mass $m_{b}^{\rm
PS}(2\,\mbox{GeV})=4.483$ GeV to the Core named O4, reusing the parameters
mbPSRef and mufB defined in the previous example as well as the matching scale
factor $2.0$, specified by the option parameter fnQ

      In[]:= AddPSMass["O4", "b2", mbPSRef, mufB, 5, 4.2,
              4.2, fnQ -> 2.0]

Internally, first the 5-flavor running mass is determined and subsequently
used to create the 5-flavor Core named b2. To match precisely the conversion
method used in Ref. [29] we set the running mass scale and the renormalization
scale of the strong coupling for that conversion to $4.2$ GeV here.

Reading out the resulting bottom quark standard running mass results in

      In[]:= MassMS["b2", 5]
      Out[]= 4.1621220500397325

with the expected $1$ MeV difference to the value $\overline{m}_{b}=4.163$ GeV
quoted in the reference paper. This difference arises here again because the
numerical conversion formulae employed in REvolver are exactly invertible,
i.e. they produce the same numerical mass differences regardless in which way
the conversion is carried out.

Note that the default REvolver routines convert from the MSR mass to a low-
scale short-distance mass at the intrinsic scale of the low-scale short-
distance mass (see Sec. 6.2.2) to resum potentially large logarithms involving
the ratio of the quark mass and the renormalization scale. In the example
above, we have, however, explicitly set the input parameter nfConv of the
function MassPS to $5$ to enforce conversion at the scale of the
$\overline{\mathrm{MS}}$ mass (which is the approach used in Ref. [29] and
which does not resum these logarithms). The result including log-resummation
via R-evolution differs by around $15$ MeV and can be extracted by using the
corresponding REvolver commands with default parameter settings. For example,
converting from the standard running mass to the PS mass we obtain

      In[]:= MassPS["b", 5, mufB]
      Out[]= 4.468247788253957

where nfConv is automatically set to $4$ and $R$ as well as the
renormalization scale of the strong coupling are set to mufB.

The corresponding (log-resummed) conversion from the PS mass to the standard
running mass is achieved by adding to the Core named O4 the bottom quark PS
mass $m_{b}^{\rm PS}(2\,\mbox{GeV})=4.483$ GeV, specifying log resummation via
R-evolution. The associated new Core named b3 is created by

      In[]:= AddPSMass["O4", "b3", mbPSRef, mufB, fnQ -> 2.0]

resulting in a bottom quark standard running mass of

      In[]:= MassMS["b3", 5]
      Out[]= 4.176206116625182

At this point it should be mentioned that for bottom quarks large scale
hierarchies cannot arise, such that the log-resummed conversion is not
superior and the difference between the fixed-order and the log-resummed
conversions may be better considered as a scheme variation. To illustrate this
we show, order by order, the bottom quark PS mass computed in the fixed-order
expansion as well as by utilizing R-evolution:

      In[]:= MassPS["b", 5, mufB, 5, mbmb, mbmb, 1.0, #] & /@
              {1, 2, 3, 4}
             MassPS["b", 5, mufB, 4, mufB, mufB, 1.0, #] & /@
              {1, 2, 3, 4}
      Out[]= {4.371486051575019, 4.452182487749768,
              4.484416838176134, 4.484037803646092}
      Out[]= {4.456398324160689, 4.487678320598545,
              4.494296775616467, 4.468247788253957}

Next, we consider the 1S scheme. Comparing the 1S bottom quark mass
$m_{b}^{\rm 1S}=4.670$ GeV quoted Ref. [29] and obtained in REvolver using the
relativistic counting and conversion at the high scale $\overline{m}_{b}$
(which again agrees with the approach used in Ref. [29] and does not account
for the resummation of logarithms involving the bottom mass and inverse Bohr
radius) we obtain

      In[]:= mb1SRef = 4.670;
             mb1SRevo = Mass1S["b", 5, 5, mbmb, "relativistic"]
             mb1SRevo - mb1SRef
      Out[]= 4.671281708897272
      Out[]= 0.0012817088972720825

Again, the $1$ MeV discrepancy emerges from the perturbative difference
between converting from or to the standard running mass. For comparison, we
also present the corresponding 1S mass value accounting for log-resummation
via R-evolution. This can be achieved by calling the routine Mass1S without
any optional parameters. This implies that the default setting is used, namely
that the 1S mass is converted from the MSR mass at the inverse Bohr radius
$M_{q,B}$ (the intrinsic scale of the 1S mass, see Sec. 6.2.2) and that non-
relativistic counting is applied. The bottom 1S mass is, however, sensitive to
the ultra-soft scale $M^{2}_{q,B}/m_{b}$, related to logarithms of the strong
coupling. This entails that for the conversion from the MSR mass (at the
inverse Bohr radius) a scale larger than $M_{q,B}$ may be adopted as the
renormalization scale for the strong coupling to avoid perturbative
instabilities.666This issue does not arise for the top quark due to is large
mass value. Adopting $2M_{q,B}$ for the renormalization scale of the strong
coupling, the 1S mass accounting for log-resummation is obtained by

      In[]:= MbB = MBohr["b", 5]
             Mass1S["b", 5, 4, MbB, "nonrelativistic", 2*MbB]
      Out[]= 1.8792648142678285
      Out[]= 4.6676688462595175

The result differs by $4$ MeV to the one using fixed-order conversion. As for
the bottom PS mass, there is no conceptual improvement using log-resummation.
This can again be seen from comparing the 1S mass values at lower orders.
Order by order, the values of the bottom quark 1S mass, derived in fixed-order
and with R-evolution, respectively, are given by

      In[]:= Mass1S["b", 5, 5, mbmb, "relativistic", mbmb, #] & /@
              {1, 2, 3, 4}
             Mass1S["b", 5, 4, MbB, "nonrelativistic", 2*MbB,
              #] & /@ {1, 2, 3, 4}
      Out[]= {4.516545706062886, 4.64055194539697,
              4.680062686477036, 4.671281708897272}
      Out[]= {4.609340303152403, 4.666301702298438,
              4.6809114692112495, 4.6676688462595175}

Of course one can also convert from the 1S scheme to the standard running mass
using the routine Add1SMass, shortly demonstrated in the following.

Reusing the Core named O4, which contains $4$ massless quarks, we create a new
Core with name b4 by adding the bottom 1S mass using fixed-order conversion

      In[]:= Add1SMass["O4", "b4", mb1SRef, 5, mbmb,
              "relativistic", mbmb, fnQ -> 2.0]

which results in the following bottom quark standard running mass:

      In[]:= MassMS["b4", 5]
      Out[]= 4.161807379452769

again with the expected $1$ MeV difference to the value
$\overline{m}_{b}=4.163$ GeV quoted in Ref. [29], related to invertible
numerical conversion algorithm employed in REvolver. The analogous procedure
for converting the 1S mass to the standard running mass using R-evolution can
be performed by executing

      In[]:= Add1SMass["O4", "b5", mb1SRef, 4, MbB,
              "nonrelativistic", 2*MbB, fnQ -> 2.0]
      In[]:= MassMS["b5", 5]
      Out[]= 4.165169213761735

#### Charm and bottom quark RS masses

In Ref. [30] the charm quark RS mass $m_{c}^{\rm RS}(\nu_{f})=1.202$ GeV at
the scale $\nu_{f}=1$ GeV was extracted from charmonium bound states masses,
given in Eq. (2.15) of that paper. The RS mass was converted to the standard
running charm mass $\overline{m}_{c}=1.217$ GeV with 4-loop fixed-order
formulae for $\alpha_{s}^{(5)}(m_{Z})=0.1184$ and using $\mu_{c}=1.27$ GeV as
the charm threshold matching scale. The renormalization scale of the strong
coupling was set to $\nu_{c}=1.5$ GeV and the pole mass renormalon
normalization constant for $n_{f}=3$ active flavors was quoted to be
$N_{m}=0.5626$. This relation can be easily reproduced with REvolver. First,
we set the relevant variables

      In[]:= {amZ, Nm} = {0.1184, 0.5626};
             {mcmc, mcmcMatch, nufc, nuc} =
              {1.217, 1.27, 1.0, 1.5};
             {mbmb, mbmbMatch} = {4.185, 4.2};

and create the respective Core

      In[]:= CoreDeleteAll[]
             CoreCreate["cb", 5, {5, amZ, mZdef},
              {{4, mcmc, mcmc}, {5, mbmb, mbmb}},
              fMatch -> {mcmcMatch/mcmc, mbmbMatch/mbmb}]

We include the standard running bottom mass $\overline{m}_{b}=4.185$ GeV with
the corresponding threshold matching scale $\mu_{b}=4.2$ GeV to enable
automatic matching of the strong coupling constant which is given at
$\mu=m_{Z}$. The charm quark RS mass computed by REvolver and the difference
to the value quoted in the reference paper are given by

      In[]:= mcRSRef = 1.202;
             mcRSREvo = MassRS["cb", 4, nufc, 4, mcmc, nuc,
               4, 4, Nm]
             mcRSRef - mcRSREvo
      Out[]= 1.201144351209422
      Out[]= 0.0008556487905779786

showing agreement at the sub-MeV level.

In Ref. [30] also the bottom quark RS mass value $m_{b}^{\rm
RS}(\nu_{f})=4.379$ GeV at the scale $\nu_{f}=2$ GeV was extracted from
bottomonium meson masses [ see Eq. (2.3) of that paper ] and subsequently
converted to the standard running bottom mass $\overline{m}_{b}=4.379$ GeV
with 4-loop fixed-order formulae using $\mu_{b}=4.2$ GeV as the bottom
threshold matching scale and $\nu_{b}=2.5$ GeV as the renormalization scale of
the strong coupling. The charm mass corrections in the perturbative relation
between the RS and the $\overline{\mathrm{MS}}$ mass have been implemented in
an effective way by evaluating the conversion series that relates the
$\overline{\mathrm{MS}}$ and the RS masses for the bottom quark and $n_{\ell}$
massless quarks for $n_{\ell}=3$ dynamical flavors, but using the 3-flavor
strong coupling $\alpha_{s}^{(3)}$ computed with the charm mass threshold
properly accounted for. This treats the charm quark as a decoupled flavor. The
REvolver RS mass routines do not directly provide this functionality, but in
the following we show how this evaluation can still be carried out in
REvolver. We first set the necessary parameters not yet defined before

      In[]:= {nufb, nub} = {2.0, 2.5};

and compute the 3-flavor strong coupling $\alpha_{s}^{(3)}(\mu_{c})$ at the
charm matching scale $\mu_{c}=1.27$ GeV, taking into account the bottom as
well as charm quark thresholds in the usual way

      In[]:= as3 = AlphaQCD["cb", mcmcMatch, 3];

Next, we define a Core with decoupled charm quarks by using the precomputed
$n_{f}=3$ strong coupling and specifying a total number of $4$ flavors,
including $3$ massless flavors and the massive bottom quark as the 4-th flavor

      In[]:= CoreCreate["b", 4, {3, as3, mcmcMatch},
              {4, mbmb, mbmb}, fMatch -> {mbmbMatch/mbmb}]

Setting all previously defined scales we get

      In[]:= mbRSRef = 4.379;
             mbRSREvo = MassRS["b", 4, nufb, 4, mbmb, nub,
               4, 4, Nm]
             mbRSREvo - mbRSRef
      Out[]= 4.379013026999629
      Out[]= 0.00001302699962923981

for the bottom quark mass in the RS scheme with a negligible difference to the
value quoted in the reference article.

We can also convert from the bottom RS to the standard running mass using the
REvolver routine AddRSMass. To do that, we first create a Core containing the
$3$ massless non-decoupled flavors with

      In[]:= CoreCreate["O3", {3, as3, mcmcMatch}]

and subsequently create a new Core with name b2 by adding the bottom RS mass
to the Core with name O3 setting the various input parameters in analogy to
the example above, i.e. using fixed-order conversion

      In[]:= AddRSMass["O3", "b2", mbRSRef, nufb, 4, 4.2, nub, 4,
              4, Nm]

Internally, first the 5-flavor running mass is determined and subsequently
used to create the Core named b2. We set the running mass scale for that
conversion to $4.2$ GeV.

Extracting the bottom quark standard running mass from the newly created Core
results in

      In[]:= MassMS["b2", 4]
      Out[]= 4.184987620121892

#### Bottom quark kinetic masses

In Ref. [21] the 3-loop corrections to the perturbative relation between the
pole and kinetic quark mass schemes have been determined for the case of one
massive quark and $n_{\ell}$ massless quarks. Here we show how to employ
REvolver to reproduce the values for the bottom quark kinetic mass $m_{b}^{\rm
kin}(\mu)$ at the scale $\mu=1$ GeV given in Eq. (8) of that paper, obtained
by converting from the standard running mass $\overline{m}_{b}=4.163$ GeV.

The conversion from $\overline{m}_{b}$ to the bottom kinetic mass was
considered using the perturbative series computed with the massive bottom
quark and $n_{\ell}$ massless quarks for $n_{\ell}=4$ (i.e. with the charm
quark treated as massless) with the result $m_{b}^{\rm kin}(\mu)=4.523$ GeV,
as well as for $n_{\ell}=3$ (i.e. with the charm quark treated as decoupled)
with the result $m_{b}^{\rm kin}(\mu)=4.521$ GeV. The first case can be
reproduced in a straightforward way with REvolver as it corresponds to a
realistic physical scenario.

First we define the relevant parameters

      In[]:= {mbmb, amZ, muCutb} = {4.163, 0.1179, 1.0};

and create a Core in which the charm quark is treated as massless

      In[]:= CoreDeleteAll[];
             CoreCreate["b4", 5, {5, amZ, mZdef},
              {{5, mbmb, mbmb}}]

The value obtained by REvolver and the difference to the value quoted in Ref.
[21] is

      In[]:= mKinb4Ref = 4.523;
             mKinb4REvo = MassKin["b4", 5, muCutb, 5, mbmb, mbmb]
             mKinb4REvo - mKinb4Ref
      Out[]= 4.523457226246508
      Out[]= 0.0004572262465085686

i.e. there is perfect agreement.

The second case is in close analogy to the treatment of the RS mass just
discussed above and requires that the conversion series, which relates the
$\overline{\mathrm{MS}}$ and the kinetic mass for the bottom quark and
$n_{\ell}$ massless quarks, is evaluated for $n_{\ell}=3$ dynamical flavors,
but using the 3-flavor strong coupling $\alpha_{s}^{(3)}$ computed with the
charm mass threshold properly accounted for. REvolver does not provide
functionality to carry out this conversion directly, but in the following we
show how this evaluation can still be performed.

To obtain the $n_{f}=3$ strong coupling in the usual way, we first create a
new Core involving a massive bottom as well as charm quark

      In[]:= mcmc = 1.263;
             CoreCreate["bc", 5, {5, amZ, mZdef},
              {{4, mcmc, mcmc}, {5, mbmb, mbmb}}]

and extract $\alpha_{s}^{(3)}(\overline{m}_{b})$

      In[]:= a3mbmb = AlphaQCD["bc", mbmb, 3];

Next, we create the Core to be employed for the mass scheme conversion. We use
the computed $n_{f}=3$ strong coupling and specify a total number of $4$
flavors with $3$ massless flavors and the massive bottom quark

      In[]:= CoreCreate["b3", 4, {3, a3mbmb, mbmb},
              {{4, mbmb, mbmb}}]

The value of the kinetic bottom quark mass computed by REvolver and the
difference to the corresponding result quoted in Ref. [21] can now be obtained
by

      In[]:= mKinb3Ref = 4.521;
             mKinb3REvo = MassKin["b3", 4, muCutb, 4, mbmb, mbmb]
             mKinb3Ref - mKinb3REvo
      Out[]= 4.520784332346221
      Out[]= 0.00021566765377922792

The numbers are in perfect agreement.

In Ref. [22] the lighter flavor mass corrections to the relation between the
pole and kinetic masses were computed explicitly up to ${\cal
O}(\alpha_{s}^{3})$. These corrections are implemented in REvolver. They have
the property that they exclusively come from the flavor number decoupling
relations of the strong coupling, and they are also referred to as “scheme B”
in that reference. In Eq. (79) of that reference the bottom quark kinetic mass
$m_{b}^{\mathrm{kin}}(1\,\mathrm{GeV})=4.526$ GeV was obtained using
$\alpha_{s}^{(5)}(m_{Z})=0.1179$ for the strong coupling, the charm quark
running mass $\overline{m}_{c}(2\,\mathrm{GeV})=0.993$ GeV and the bottom
quark standard running mass $\overline{m}_{b}=4.136$ GeV. For the calculation,
fixed-order conversion was applied using the standard running bottom mass
$\overline{m}_{b}$ as the renormalization scale of the strong coupling. To
reproduce the result with REvolver we define the input values with

      In[]:= {amZ, mc3, mbmb} = {0.1179, 0.993, 4.163};

and the related Core with

      In[]:= CoreDelete["bc"]
             CoreCreate["bc", 5, {5, amZ, mZdef},
               {{4, mc3, 3}, {5, mbmb, mbmb}}]

The bottom quark kinetic mass extracted by REvolver and its difference to the
reference value are returned by

      In[]:= mbKinRef = 4.526;
             mbKinREvo = MassKin["bc", 5, 1.0, 5, mbmb, mbmb]
             mbKinRef - mbKinREvo
      Out[]= 4.527084826155039
      Out[]= -0.0010848261550391314

The small deviation of $1$ MeV originates from the fact that in Ref. [22] the
light massive flavor corrections related to the charm quark are parametrized
in terms of $\overline{m}_{c}(3\,\mathrm{GeV})$, while in REvolver they are
parametrized in terms of the charm standard running mass $\overline{m}_{c}$.

Performing the same conversion with R-evolution to resum logarithms of the
intrinsic physical scales of the mass schemes gives

      In[]:= MassKin["bc", 5, 1.0]
      Out[]= 4.53472077059131

i.e. the bottom quark kinetic mass with log-resummation is larger by about $7$
MeV. We note that, following Ref. [22], the default value for the kinetic mass
intrinsic scale (where the default log-resummed conversion between the running
and the kinetic mass is carried out) is set to be twice its renormalization
scale, that is, $2$ GeV in the example above.

####  7.2.4 Top quark mass at low scales

We return to Ref. [29] and investigate the top quark mass conversion between
the $\overline{\mathrm{MS}}$ and various short-distance schemes with 4-loop
accuracy. This concerns Table I of that article.

Following Ref. [29] we first define the values for the strong QCD coupling
$\alpha_{s}^{(5)}(m_{Z})=0.1185$, the top quark standard running mass
$\overline{m}_{t}=163.643$ GeV and the renormalization scale of the PS mass
$\mu_{f}=20$ GeV

      In[]:= {amZ, mtmt, mufT} = {0.1185, 163.643, 20.0};

and a Core with a top threshold matching scale of twice the top standard
running mass

      In[]:= CoreDeleteAll[]
             CoreCreate["t", 6, {5, amZ, mZdef}, {{6, mtmt, mtmt}},
              fMatch -> {2.0}]

In analogy to Sec. 7.2.3, we evaluate the PS and 1S top quark mass employing
REvolver and determine the difference to the values $m_{t}^{\rm 1S}=172.227$
GeV and $m_{t}^{\rm PS}(20\,\mbox{GeV})=171.792$ GeV quoted in the reference
paper by converting directly from the $\overline{\mathrm{MS}}$ scheme, leading
to

      In[]:= mtPSRef = 171.792;
             mtPSRevo = MassPS["t", 6, mufT, 6, mtmt, mtmt]
             mtPSRevo - mtPSRef
      Out[]= 171.7950829141975
      Out[]= 0.0030829141975061702

for the PS scheme, and

      In[]:= mt1SRef = 172.227;
             mt1SRevo = Mass1S["t", 6, 6, mtmt, "relativistic"]
             mt1SRevo - mt1SRef
      Out[]= 172.2298841816133
      Out[]= 0.0028841816132967324

for the 1S scheme, respectively.

As already described in Sec. 7.2.3, the small difference (of about $3$ MeV)
results from the fact that in REvolver conversion is based on formulae where
the (standard) running mass is taken as the input, while in Ref. [29] the
conversion was carried out starting from the PS and 1S masses.

We now investigate the influence of large logarithms of the ratio between the
top quark mass and the intrinsic scale of the respective short-distance
schemes. For the top quark the impact of the summation of these logarithms can
be significant. To this end we consider the convergence of the perturbative
series relating the standard running top mass $\overline{m}_{t}$ to the PS
mass $m_{t}^{\rm PS}(20\,\mbox{GeV})$, the 1S mass $m_{t}^{\rm 1S}$ as well as
the running mass $m_{t}^{(5)}(2\,\mathrm{GeV})$.

Using fixed-order conversion (without log-resummation) from $\overline{m}_{t}$
and using $\overline{m}_{t}$ as the renormalization scale, the top PS mass at
$4$ loops, and the corresponding corrections of order
$\mathcal{O}(\alpha_{s}^{n})$ with $1\leq n\leq 4$ amount to

      In[]:= MassPS["t", 6, mufT, 6, mtmt, mtmt]
             Table[
              MassPS["t", 6, mufT, 6, mtmt, mtmt, 1.0, n] -
               MassPS["t", 6, mufT, 6, mtmt, mtmt, 1.0, n - 1],
              {n, 1, 4}]
      Out[]= 171.7950829141975
      Out[]= {6.636196514422409, 1.199280138792858,
              0.26779429813998945, 0.048811962842222556}

On the other hand, employing R-evolution to first evolve down to the PS mass
renormalization scale $\mu_{f}=20$ GeV and then converting to the PS mass at
that scale, which is REvolvers default procedure, we obtain

      In[]:= MassPS["t", 6, mufT]
             Table[
              MassPS["t", 6, mufT, 5, mufT, mufT, 1.0, n] -
               MassPS["t", 6, mufT, 5, mufT, mufT, 1.0, n - 1],
              {n, 1, 4}]
      Out[]= 171.79912206134247
      Out[]= {0., 0.06895809427228983,
              0.001838267705210228, -0.01976315025407871}

We observe that the corrections are considerably smaller when R-evolution is
accounted for and that there is a difference of around $4$ MeV in the final
4-loop converted values for the two approaches. The $\mathcal{O}(\alpha_{s})$
correction term is zero in the case of the log-resummed conversion since the
1-loop perturbative coefficients of the PS and MSR masses coincide.

For the conversion to the 1S mass the observation regarding fixed-order versus
log-resummed conversion is similar. In the fixed-order case we get

      In[]:= Mass1S["t", 6, 6, mtmt, "relativistic", mtmt]
             Table[
              Mass1S["t", 6, 6, mtmt, "relativistic", mtmt, n] -
               Mass1S["t", 6, 6, mtmt, "relativistic",
                mtmt, n - 1], {n, 1, 4}]
      Out[]= 172.2298841816133
      Out[]= {7.129280949172028, 1.227243021981991,
              0.21914263057425387, 0.01121757988499894}

while the log-resummed evaluation gives again substantially smaller
corrections:

      In[]:= mBohr = MBohr["t", 6]
             Mass1S["t", 6, 5]
             Table[
              Mass1S["t", 6, 5, mBohr, "nonrelativistic",
                mBohr, n] - Mass1S["t", 6, 5, mBohr,
                "nonrelativistic", mBohr, n - 1], {n, 1, 4}]
      Out[]= 32.11684612826197
      Out[]= 172.20797480544945
      Out[]= {1.1666583987086199, 0.18536848855083576,
              -0.047109547951663444, -0.010351024770614004}

where the inverse Bohr radius, i.e. the intrinsic scale of the 1S mass, is
shown as the first output for completeness.

Finally, we investigate the effect of log-resummation on the computation of
the top quark $5$-flavor running mass at $2$ GeV,
$m_{t}^{(5)}(2\,\mathrm{GeV})$, when converting from
$m_{t}^{(5)}(\overline{m}_{t})$. For simplicity we create a new Core named t,
this time where for the top threshold matching scale the default value, the
top quark standard running mass, is adopted:

      In[]:= CoreDelete["t"]
             CoreCreate["t", 6, {5, amZ, mZdef},
               {{6, mtmt, mtmt}}];

Now we extract the value of the 5-flavor MSR mass
$m_{t}^{(5)}(\overline{m}_{t})$ at the scale of the standard top running mass
from this Core as a reference

      In[]:= m5mt = MassMS["t", 6, mtmt, 5]
      Out[]= 163.67571467667918

REvolver does not provide fixed-order conversions of the running masses
between two different renormalization scales. It is, however, possible to
access these fixed-order corrections through the pole mass routine MassPoleFO.
Care has to be taken that the two calls of the MassPoleFO routines involve the
same renormalization scale (which here is the standard running mass) to ensure
that the pole mass renormalon is properly canceled:

      In[]:= FOTab =
               MassMS["t", 6, 2.0, 5] -
                  MassPoleFO["t", 6, 5, 2.0, mtmt, #] +
                  MassPoleFO["t", 6, 5, mtmt, mtmt, #] & /@
                {0, 1, 2, 3, 4};
             FOTab[[5]]
             Table[FOTab[[n]] - FOTab[[n - 1]], {n, 2, 5}]
      Out[]= 173.2337130005733
      Out[]= {7.46778271952536, 1.5242807051151317,
              0.4246371655552821, 0.14129773369833742}

The difference of the two calls of the MassPoleFO routine removes the log-
resummed corrections from the output of MassMS and replaces it by the
corresponding fixed-order terms. As for the above examples, we first display
the $m_{t}^{(5)}(2\,\mathrm{GeV})$ value at $\mathcal{O}(\alpha_{s}^{4})$ and
the fixed-order perturbative correction terms at order
$\mathcal{O}(\alpha_{s}^{n})$ with $1\leq n\leq 4$.

For the investigation of the correction terms in the case of R-evolution we
create four Cores, setting the respective loop order $n$ of the R-evolution
equation to $1\leq n\leq 4$

      In[]:= CoreCreate["t-R" <> ToString[#], 6, {5, amZ, mZdef},
                 {{6, mtmt, mtmt}}, runMSR -> #] & /@ {1, 2, 3, 4};

The four Cores differ only by the loop order used for R-evolution, but employ
the default 4-loop matching to determine $m_{t}^{(5)}(\overline{m}_{t})$ from
the standard running mass $\overline{m}_{t}$. The value of
$m_{t}^{(5)}(2\,\mathrm{GeV})$ corresponding to $4$-loop running and the
corrections coming from the individual running orders can be extracted with

      In[]:= MassMS["t-R4", 6, 2.0]
             RevoTab = MassMS["t-R" <> ToString[#], 6, 2.0] & /@
                {1, 2, 3, 4};
             RevoTab = Prepend[RevoTab, m5mt];
             Table[RevoTab[[n]] - RevoTab[[n - 1]], {n, 2, 5}]
      Out[]= 173.36186564820926
      Out[]= {8.841888795217045, 0.8531892592812653,
              0.035471795296302844, -0.044398878264530595}

where the first entry in the curly brackets is the difference between
$m_{t}^{(5)}(\overline{m}_{t})$ and $m_{t}^{(5)}(2\,\mathrm{GeV})$ obtained
with 1-loop R-evolution, while the subsequent entries refer to the differences
to the previous order in $m_{t}^{(5)}(2\,\mathrm{GeV})$ obtained when adding
the 2-, 3- and 4-loop terms to the R-evolution anomalous dimension. Once again
we observe that the convergence is much better when using R-evolution, and
that the corresponding results are more precise than using fixed-order
conversion. Using 4-loop R-evolution leads to a value for
$m_{t}^{(5)}(2\,\mathrm{GeV})$ that is about $130$ MeV higher than when using
4-loop fixed-order conversion. This difference is consistent with the size of
the 4-loop fixed-order correction of $141$ MeV, when adopting the latter as an
uncertainty for the 4-loop fixed-order conversion.

We see that the higher order corrections to the R-evolution anomalous
dimension lead to very small effects, so that one could worry that their size
may not reflect the perturbative uncertainty at the corresponding loop order.
A different way to estimate perturbative uncertainties is to perform
$\lambda$-variation in the R-evolution equation, which we shortly demonstrate
in the following.

We create a set of Cores with $\lambda$ values in the range
$1/2\leq\lambda\leq 2$ and for various loop orders in the R-evolution
equation. The Cores are created with

      In[]:= list2 = 2^Range[-1, 1, 2/49];
             Outer[
               CoreCreate["t-l" <> ToString[#2] <> "-R" <>
                  ToString[#1], 6, {5, amZ, mZdef},
                 {{6, mtmt, mtmt}}, lambdaMSR -> #2,
                 runMSR -> #1] &, {1, 2, 3, 4}, list2];

and the table containing the corresponding values of
$m_{t}^{(5)}(2\,\mathrm{GeV})$ is generated by

      In[]:= mListRevo =
               Outer[
                MassMS["t-l" <> ToString[#2] <> "-R" <>
                   ToString[#1], 6, 2.0, 5] &, {1, 2, 3, 4},
                list2];

Order by order, the values of the $\lambda$ variations divided by two are then
given by

      In[]:= Table[(Max[mListRevo[[n]]] - Min[mListRevo[[n]]])/2,
              {n, 1, 4}]
      Out[]= {1.082493387031036, 0.2525991201191715,
              0.041115232629010734, 0.02465406235192802}

The values of $\lambda$ variations are consistent with the size of the
corrections to the R-evolution equation and substantially smaller than the
size of the fixed-order corrections, illustrating that resumming logarithms
via R-evolution leads to a more precise mass conversion than using fixed-order
corrections.

####  7.2.5 Running masses for complex renormalization scales

REvolver provides the functionality to determine running masses at complex
renormalization scales. In analogy to Sec. 7.1.2, we demonstrate the quality
of the running mass evolution in the complex plane by showing numerical
consistency for the high-scale ($\overline{\mathrm{MS}}$) running top quark
mass concerning the Cauchy’s integral formula, where we adopt
$\mathtt{muMS}=350$ GeV for the central scale and $174$ GeV for the radius of
the Cauchy integral. We set up a Core with the top quark standard running mass
$\overline{m}_{t}=163$ GeV and subsequently determine the central running mass
and the corresponding value using the residue theorem

      In[]:= CoreCreate["O", 6, {5, amZdef, mZdef},
              {{6, 163, 163}}]
      In[]:= {muMS, rMS} = {350, 174};
             mMSDirect = MassMS["O", 6, muMS, 6]
             mMSResidue =
              1/(2 Pi)
               NIntegrate[MassMS["O", 6, muMS + rMS E^(I phi), 6],
                {phi, 0, 2 Pi}, PrecisionGoal -> 10]
             mMSResidue - mMSDirect
      Out[]= 154.072755231406
      Out[]= 154.07275523140612 + 4.240739575284689*^-16 I
      Out[]= -1.1368683772161603*^-13 - 4.240739575284689*^-16 I

Both results are in perfect agreement reflecting the high numerical precision
of the evolution routines implemented in REvolver.

This high numerical precision is also maintained at much smaller
renormalization scales (where perturbation theory in general is less
reliable). We demonstrate this for the top quark running (MSR) mass at a
central scale of $2$ GeV and a radius of $1$ GeV for the Cauchy integral:

      In[]:= muR = 2;
             mRDirect = MassMS["O", 6, muR, 5]
             mRResidue =
              1/(2 Pi)
               NIntegrate[MassMS["O", 6, muR + E^(I phi), 5],
                {phi, 0, 2 Pi}, PrecisionGoal -> 10]
             mRResidue - mRDirect
      Out[]= 172.64605095372784
      Out[]= 172.64605095372804 + 4.417437057588218*^-18 I
      Out[]= 1.9895196601282805*^-13 + 4.417437057588218*^-18 I

####  7.2.6 Top quark pole masses

The pole mass scheme suffers from an ${\cal O}(\Lambda_{\rm QCD})$ renormalon,
which entails that its perturbative relation to a short-distance mass involves
a factorially divergent perturbative series and an associated ambiguity in its
value [31, 32, 3, 1]. The pole mass (and its value) can either be treated as
an order-dependent concept or one can assign its value to be in the region
where the perturbative series reaches its minimal correction term (and the
partial sum of the perturbative series increases linearly), sometimes called
the asymptotic region. We call the latter value the “asymptotic pole mass”.
The associated ambiguity arises from the principle ignorance where precisely
to truncate the partial sum within the asymptotic region. Different lines of
reasoning have been proposed concerning the determination of the asymptotic
value and the associated ambiguity, see Refs. [1, 3] for recent analyses.

REvolver provides functionalities to extract the order-dependent as well as
the asymptotic pole mass value from a Core. In case of the asymptotic value
and the associated ambiguity, REvolver provides routines employing a number of
different methods and allowing for various options as explained in Sec. 6.2.3.
In the following we demonstrate some of these functionalities reproducing
results quoted in Refs. [3] and [1] focusing on the scenario of the top quark
with massless bottom and charm quarks. Scenarios with massive bottom and charm
quarks are treated in Sec. 7.3.4.

#### Minimal correction approach for the pole mass ambiguity

First, we consider Ref. [3] where the top quark pole mass renormalon ambiguity
quoted in Eq. (4.7) is determined from the perturbative series between the
pole mass and the standard running mass for the case where the bottom and
charm mass effects are ignored. To reproduce the result we define the top
standard running mass value used in Ref. [3]

      In[]:= mtmt = 163.508;

and create an associated Core

      In[]:= CoreDeleteAll[]
             CoreCreate["t", 6, {5, amZdef, mZdef},
              {{6, mtmt, mtmt}}]

We execute the function MassPoleDetailed to extract the asymptotic pole mass
value, the ambiguity, and the order of the smallest correction term. The
method adopted in Ref. [3] to determine the pole mass ambiguity was associated
to the size of the minimal correction and is accessed in REvolver specifying
the min method:

      In[]:= MassPoleDetailed["t", 6, mtmt, mtmt, "min"]
      Out[]= {173.60839591880665, 0.06443134543422957, 8}

Furthermore the third and fourth arguments are set to the standard running top
mass mtmt to mimic the choices adopted in Ref. [3]. The choice of arguments
when calling the function MassPoleDetailed ensures that the analyzed series is
the one for the difference between pole and the (5-flavor) MSR mass at the
scale mtmt and that the renormalization scale for the strong coupling is mtmt
as well. The values for the asymptotic pole mass and its ambiguity quoted in
Eq. (4.8) in Ref. [3] are $173.608$ GeV and $67$ MeV, respectively. The
asymptotic value is in perfect agreement. The small discrepancy of $3$ MeV in
the ambiguity is related to the fact that REvolver uses the perturbative
series for the relation between the pole and the MSR mass for the routine
MassPoleDetailed while in Ref. [3] the relation between the pole and the
standard running mass was considered.

#### Asymptotic series for the pole-$\overline{\mathrm{MS}}$ mass relation

The method of estimating higher order coefficients in the
pole-$\overline{\mathrm{MS}}$ relation employed by REvolver is described in
Sec. 4.4 of Ref. [1] and relies on an asymptotic formula that can, depending
on the specified options, reproduce the exactly known coefficients, see Eq.
(36) in Sec. A.5.

The values of the asymptotic coefficients can be easily read out using the
function MassPoleFO. For the orders $5$–$9$ REvolver returns

      In[]:= (MassPoleFO["t", 6, 5, mtmt, mtmt, #] -
                  MassPoleFO["t", 6, 5, mtmt, mtmt, # - 1])/
                (mtmt * (AlphaQCD["t", mtmt, 5]/(4 Pi))^#) & /@
              Range[5, 9]
      Out[]= {1.4248584469453057*^7, 1.1661884055877335*^9,
              1.1323781055413771*^11, 1.2729796430260885*^13,
              1.6260702903901682*^15}

for $a_{n}^{\mathrm{MSR\,\prime}}$ (i.e. the series coefficients for the MSR
and pole mass difference) and

      In[]:= (MassPoleFO["t", 6, 6, mtmt, mtmt, #] -
                  MassPoleFO["t", 6, 6, mtmt, mtmt, # - 1])/
                (mtmt * (AlphaQCD["t", mtmt, 6]/(4 Pi))^#) & /@
              Range[5, 9]
      Out[]= {1.429074531848818*^7, 1.1687368433901165*^9,
              1.1344063409072752*^11, 1.2749318929376549‘*^13,
              1.6282469370459092*^15}

for $a_{n}^{\mathrm{\overline{\mathrm{MS}}\,\prime}}$ (i.e. the series
coefficients for the standard running and pole mass difference), agreeing well
with Table 2 of Ref. [1] within errors. The deviation from the central values
quoted in that table arises because in Ref. [1] the central values of
asymmetric uncertainty intervals have been quoted.

#### Asymptotic series for the pole-$\overline{\mathrm{MS}}$ mass relation and
the asymptotic pole mass

Having access to the pole mass at in principle arbitrary order using
conversion from $\overline{\mathrm{MS}}$ and MSR masses at arbitrary
renormalization scales, REvolver allows to easily reproduce the numbers
utilized to produce Fig. 6 of Ref. [2], where the order-dependent top quark
pole mass was shown. We create a new Core to switch to the top standard
running mass value $\overline{m}_{t}=163$ GeV in accordance with Ref. [2]

      In[]:= mtmt = 163;
             CoreDeleteAll[]
             CoreCreate["t", 6, {5, amZdef, mZdef},
              {{6, mtmt, mtmt}}]

Converting directly from the $\overline{\mathrm{MS}}$ mass at the scale
$\overline{m}_{t}$, including scale variation in the range
$\overline{m}_{t}/2\leq\mu\leq 2\overline{m}_{t}$ the pole mass, order by
order $0\leq n\leq 12$, is given by

      In[]:= Table[{{ord, MassPoleFO["t", 6, 6, mtmt, mtmt, ord]},
               {MassPoleFO["t", 6, 6, mtmt, mtmt/2, ord] -
                 MassPoleFO["t", 6, 6, mtmt, mtmt, ord],
                MassPoleFO["t", 6, 6, mtmt, 2*mtmt, ord] -
                 MassPoleFO["t", 6, 6, mtmt, mtmt, ord]}},
              {ord, 0, 12}]
      Out[]= {{{0, 163.}, {0., 0.}},
              {{1, 170.5097731904309},
               {0.7122980457679375, -0.5972354956839183}},
              {{2, 172.11277687226905},
               {0.27618516105792423, -0.3086769323029728}},
              {{3, 172.6079790748946},
               {0.12091854036228256, -0.15116686757511388}},
              {{4, 172.80297651883814},
               {0.06337072738099891, -0.08027361964516899}},
              {{5, 172.9150350500228},
               {0.04822992022724293, -0.05421541398965246}},
              {{6, 172.99420261176041},
               {0.04210709824803871, -0.043768064396573436}},
              {{7, 173.06058298165766},
               {0.042838521236177485, -0.04090645714344987}},
              {{8, 173.1250293570062},
               {0.04963473302441912, -0.04338206485874707}},
              {{9, 173.19612979388887},
               {0.0644919645042421, -0.05147133915082236}},
              {{10, 173.28398617671283},
               {0.09283257854019666, -0.06758882965976909}},
              {{11, 173.40418238677069},
               {0.14657118914180955, -0.09737521610099975}},
              {{12, 173.58454177427788},
               {0.2518136273276639, -0.1527834480848469}}}

Figure 4: Top quark pole mass $m^{\mathrm{pole}}_{t}$ as a function of the
perturbative order for massless bottom and charm quarks, converted from the
$\overline{\mathrm{MS}}$ mass (upper left panel) or the MSR mass (other
panels) at different scales of $R$. The central dots are obtained for the
default renormalization scales for the strong coupling, the error bars
represent the scale variation. This corresponds to the black dots and error
bars in Fig. 6 of Ref. [2].

Analogous commands can be used to extract the relevant values of the pole mass
when converting from the MSR scheme at various scales of $R$. Utilizing these
numbers it is now straightforward to reproduce Fig. 6 of Ref. [2] for
vanishing bottom and charm mass (black dots in Ref. [2]), see Fig. 4.777Note
that the upper left panel of Fig. 4 (corresponding to the numbers given above)
shows the conversion directly from the standard $\overline{\mathrm{MS}}$ mass
instead of the MSR mass at the high scales. The asymptotic value which can be
assigned to the pole mass lies in the region where the series grows linearly,
which is in Fig. 4, roughly in the region around $173$ GeV. This is confirmed
by the output of the command MassPoleDetailed, which gives

      In[]:= MassPoleDetailed["t", 6, mtmt, mtmt, "min"]
      Out[]= {173.07508395890844, 0.06445302708283407, 8}

      In[]:= MassPoleDetailed["t", 6, mtmt, mtmt, "range"]
      Out[]= {173.07976320487433, 0.12605058936964042, 8}

      In[]:= MassPoleDetailed["t", 6, mtmt, mtmt, "drange"]
      Out[]= {173.1081654841447, 0.16284602132711257, 8}

for the three methods to estimate the asymptotic pole mass supported by
REvolver from the relation between the pole and the MSR mass
$m_{t}^{(5)}(\overline{m}_{t})$ at the scale $163$ GeV. While the estimation
of the ambiguity can vary by more than a factor of $2$ depending on the
employed strategy, the estimated asymptotic pole mass is fairly stable.

Figure 5:  Comparison of the pole mass ambiguity estimating strategies
available in REvolver obtained from the series between the pole mass and the
running mass $m_{t}^{(5)}(R)=m_{t}^{\rm MSR}(R)$ for the case that all quarks
lighter than the top quark are considered to be massless. The figure shows the
size of the ambiguity of the top quark pole mass in GeV as estimated by the
min (blue lower curve), drange (orange upper curve) and range (green middle
curve) strategies over the scale $R$.

Figure 5 shows a comparison of the three strategies of estimating the pole
mass ambiguity implemented in REvolver from the perturbative series between
the pole mass and the running (MSR) mass $m_{t}^{(5)}(R)=m_{t}^{\rm MSR}(R)$
at the scale $R$. The plot is produced with the command

      In[]:= LogLinearPlot[
              {MassPoleDetailed["t", 6, R, "min"][[2]],
               MassPoleDetailed["t", 6, R, "drange"][[2]],
               MassPoleDetailed["t", 6, R, "range"][[2]]},
              {R, 2, mtmt}, PlotRange -> {0.05, 0.26}]

omitting all options related to plot styling and annotations for brevity.

The figure shows the size of the ambiguity of the top quark pole mass in GeV
as estimated by the strategies over the scale $R$, specifying the scale in the
perturbative relation between the pole and running mass. All quarks lighter
than the top quark are considered to be massless. Due to its discrete nature,
the drange method produces a discontinuous curve, while the continuous version
range provides a very stable result. The min method leads to an ambiguity
estimate that is logarithmically decreasing with $R$.

###  7.3 Cores with Multiple Massive Quarks

####  7.3.1 Strong coupling and $\Lambda_{\rm QCD}$

#### Strong coupling evolution through multiple flavor thresholds

In Ref. [33] the 3-flavor strong coupling value $\alpha_{s}^{(3)}(m_{\tau})$
was determined from hadronic $e^{+}e^{-}$ $R$-ratio data for c.m. energies
below the charm production threshold. For this determination, fixed-order (FO)
as well as contour-improved (CI) perturbation theory were applied, and the
corresponding values $\alpha_{s}^{(3)}(m_{\tau})_{\rm FO}=0.298$ and
$\alpha_{s}^{(3)}(m_{\tau})_{\rm CI}=0.304$ were subsequently converted to
$\alpha_{s}^{(5)}(m_{Z})$ employing the 4-loop beta function as well as the
3-loop matching relations at the charm and bottom thresholds. The respective
values for the strong coupling at the Z-scale were quoted in Eq. (4.7) of that
reference as $\alpha_{s}^{(5)}(m_{Z})_{\rm FO}=0.1158$ and
$\alpha_{s}^{(5)}(m_{Z})_{\rm CI}=0.1166$, respectively.

For the conversion, the values $\overline{m}_{c}=1.28$ GeV and
$\overline{m}_{b}=4.2$ GeV were used for the charm and bottom standard running
masses, respectively, and the charm and bottom threshold matching scales were
set to $\mu_{c}=2.0$ GeV and $\mu_{b}=4.0$ GeV, respectively. To reproduce the
conversion with REvolver we define these input values with

      In[]:= mTau = 1.77686;
             {mcmc, mbmb} = {1.28, 4.2};
             matchList = {2.0/mcmc, 4.0/mbmb};
             {amTauFO, amTauCI} = {0.298, 0.304};

and define one Core for each value of $\alpha_{s}^{(3)}(m_{\tau})$

     In[]:= CoreCreate["FO", 5, {3, amTauFO, mTau},
             {{4, mcmc, mcmc}, {5, mbmb, mbmb}}, runAlpha -> 4,
             orderAlpha -> 3, fMatch -> matchList]
            CoreCreate["CI", 5, {3, amTauCI, mTau},
             {{4, mcmc, mcmc}, {5, mbmb, mbmb}}, runAlpha -> 4,
             orderAlpha -> 3, fMatch -> matchList]

where the parameters relevant for running and matching were set with the
optional parameters runAlpha, orderAlpha and fMatch.

The values of $\alpha_{s}^{(5)}(m_{Z})$ can now be easily extracted by

      In[]:= AlphaQCD["FO", mZdef, 5]
             AlphaQCD["CI", mZdef, 5]
      Out[]= 0.11581049250660494
      Out[]= 0.11662117512902259

in full agreement with Eq. (4.7) of Ref. [33].

Similarly, in Ref. [34] $\alpha_{s}^{(3)}(m_{\tau})$ was determined from
$\tau$ decay data and converted to $\alpha_{s}^{(5)}(m_{Z})$ employing 4-loop
evolution and matching. We update the values for the standard running bottom
and charm masses to the ones in the reference paper and define the value of
$\alpha_{s}^{(3)}(m_{\tau})$ as given in Eq. (19) therein as well as the
associated Core

      In[]:= {mcmc, mbmb} = {1.286, 4.164};
             amTau = 0.332;

      In[]:= CoreDeleteAll[]
             CoreCreate["O", 5, {3, amTau, mTau},
              {{4, mcmc, mcmc}, {5, mbmb, mbmb}}, runAlpha -> 4]

where, following Ref. [34], the strong coupling evolution is set to 4-loop
precision. We focus on the central value here for brevity. The value of
$\alpha_{s}^{(5)}(m_{Z})$ can now easily be extracted with

      In[]:= AlphaQCD["O", mZdef, 5]
      Out[]= 0.12019776978833413

agreeing with Eq. (20) of Ref. [34].

#### Flavor number dependence of the QCD scale
$\Lambda_{\mathrm{QCD}}^{(n_{f})}$

To demonstrate the capability of REvolver to extract
$\Lambda_{\mathrm{QCD}}^{(n_{f})}$ in various flavor number schemes we
consider Ref. [35], where $\alpha_{s}^{(5)}(m_{Z})=0.1179$ has been determined
by the ALPHA collaboration and values for $\Lambda_{\mathrm{QCD}}^{(nf)}$ with
$3\leq n_{f}\leq 5$ were determined, see Eqs. (4.14) as well as (5.3) and
(5.4) of that reference.

To reproduce the $\Lambda_{\mathrm{QCD}}^{(n_{f})}$ values with REvolver we
create a Core with massive bottom and charm quarks (with standard running mass
values taken from the PDG [23]) and define the value of
$\alpha_{s}^{(5)}(m_{Z})$ to be the one given in the reference paper

      In[]:= {amZ, mcmc, mbmb} = {0.1179, 1.27, 4.18};

      In[]:= CoreDeleteAll[]
             CoreCreate["O", 5, {5, amZ, mZdef},
              {{4, mcmc, mcmc}, {5, mbmb, mbmb}}]

We now extract $\Lambda_{\mathrm{QCD}}^{(n_{f})}$ in the $n_{f}=5\,,4\,,3$
flavor schemes in the so-called $\overline{\mathrm{MS}}$ scheme [ see Eq. (21)
in App. A.2 ] with

      In[]:= LambdaQCD["O", 5]
             LambdaQCD["O", 4]
             LambdaQCD["O", 3]
      Out[]= 0.20745573124097272
      Out[]= 0.288991643859333
      Out[]= 0.33153068987761314

The numbers are in exact agreement with the values shown in Ref. [35].

####  7.3.2 Top quark running mass at low scales

In Ref. [36] a calibration analysis was provided which suggested that the top
quark running mass $m_{t}(1\,\mathrm{GeV})$ agrees within theoretical
uncertainties with the MC top mass parameter $m_{t}^{\mathrm{MC}}$ of the
Pythia event generator. Specifically, the authors quote a value of
$m_{t}(1\,\mathrm{GeV})=172.82\pm 0.19$ GeV for a calibration with
$m_{t}^{\mathrm{MC}}=173$ GeV, given in Table 1. That analysis was carried out
in the approximation of massless bottom and charm quarks. Here we demonstrate
how REvolver can be used to investigate how finite charm and bottom quark
masses affect the value of the standard running mass $\overline{m}_{t}$
calculated from $m_{t}(1\,\mathrm{GeV})$. As an estimate of the uncertainties
we quote the difference of the values computed using $4$\- and $3$-loop
R-evolution.

To this end we create six Cores: two with massless bottom and charm quarks,
two with a massless charm quark and a massive bottom quark, and two with
massive bottom and charm quarks, employing $\overline{m}_{c}=1.27$ GeV and
$\overline{m}_{b}=4.18$ GeV for the standard running charm and bottom quark
masses. Among the two respective Cores one employs 3-loop and the other 4-loop
precision for the R-evolution equation:

      In[]:= {mcmc, mbmb, mt1} = {1.27, 4.18, 172.82};

      In[]:= CoreDeleteAll[]
             CoreCreate["t4", 6, {5, amZdef, mZdef}, {{5, mt1, 1}},
              runMSR -> 4]
             CoreCreate["t3", 6, {5, amZdef, mZdef}, {{5, mt1, 1}},
              runMSR -> 3]
             CoreCreate["bt4", 6, {5, amZdef, mZdef},
              {{5, mbmb, mbmb}, {4, mt1, 1}}, runMSR -> 4]
             CoreCreate["bt3", 6, {5, amZdef, mZdef},
              {{5, mbmb, mbmb}, {4, mt1, 1}}, runMSR -> 3]
             CoreCreate["cbt4", 6, {5, amZdef, mZdef},
              {{4, mcmc, mcmc}, {5, mbmb, mbmb}, {3, mt1, 1}},
              runMSR -> 4]
             CoreCreate["cbt3", 6, {5, amZdef, mZdef},
             {{4, mcmc, mcmc}, {5, mbmb, mbmb}, {3, mt1, 1}},
             runMSR -> 3]

The Core setups specify that $m_{t}^{\rm MSR}(1\,\mbox{GeV})=172.82$ GeV is
the running top quark mass in the scheme in which all massive flavors above
the scale of $1$ GeV are integrated out. It is easy to extract the top quark
standard running masses and their uncertainties from the Cores above with

      In[]:= MassMS["t4", 6]
             MassMS["t4", 6] - MassMS["t3", 6]
      Out[]= 163.00308510782256
      Out[]= 0.05105270648016358

      In[]:= MassMS["bt4", 6]
             MassMS["bt4", 6] - MassMS["bt3", 6]
      Out[]= 162.94352028117146
      Out[]= 0.050962554040836494

      In[]:= MassMS["cbt4", 6]
             MassMS["cbt4", 6] - MassMS["cbt3", 6]
      Out[]= 162.92268002391415
      Out[]= 0.049529435610025985

where the first numbers refer to the setup with massless bottom and charm
quarks and the last ones to the one where bottom and charm quarks have mass.

We observe that the finite bottom mass lowers the standard top running mass by
around $60$ MeV and the finite charm mass decreases the standard top running
mass by about another $20$ MeV. The perturbative uncertainty of the conversion
amounts to around $50$ MeV in all cases. Together, the finite charm and bottom
masses lower the standard running top mass by about $80$ MeV, which is larger
than the perturbative uncertainty.

####  7.3.3 Charm mass effects for the bottom quark mass

To demonstrate the effect of lighter massive flavors in conversions between
short-distance mass schemes we consider Ref. [37], where effects of the finite
charm quark mass on bottom quark mass determinations from $\Upsilon$ mesons
were (calculated and) examined. In Eq. (99) of that paper the bottom standard
running mass is computed from the 1S mass $m_{b}^{\mathrm{1S}}=4.7$ GeV with
and without the contribution of the charm standard running quark mass
$\overline{m}_{c}=1.5$ GeV. The strong coupling value was set to
$\alpha_{s}^{(4)}(m_{b}^{\mathrm{1S}})=0.216$. We define these parameters with

      In[]:= {mb1S, mcmc, amb1S} = {4.7, 1.5, 0.216};

and create two 4-flavor Cores, one with a massless charm quark and one with
the charm mass given above

      In[]:= CoreDeleteAll[]
             CoreCreate["c0", {4, amb1S, mb1S}]
             CoreCreate["cm", 4, {4, amb1S, mb1S},
              {{4, mcmc, mcmc}}]

We then create two additional Cores by adding the bottom quark 1S mass
$m_{b}^{\mathrm{1S}}=4.7$ GeV to the existing Cores, where, following Ref.
[37], the perturbative order of its relation to the standard running mass is
$\mathcal{O}(\alpha_{s}^{3})$ and the relativistic (“upsilon-expansion”)
counting scheme is adopted.

      In[]:= Add1SMass["c0", "bc0", mb1S, 5, mb1S, "relativistic",
              mb1S, 3]
             Add1SMass["cm", "bcm", mb1S, 5, mb1S, "relativistic",
              mb1S, 3]

We have furthermore set the running mass scale and the renormalization scale
of the strong coupling to $m_{b}^{\mathrm{1S}}$ and the parameter nfConv to
$5$ to specify direct conversion without R-evolution to be in accordance with
the computation carried out in Ref. [37].

Extracting the standard mass from the massless charm Core gives

      In[]:= MassMS["bc0", 5]
      Out[]= 4.190311348710497

in perfect agreement with
$\overline{m}_{b}=(4.7-0.382-0.098-0.030)\,\mathrm{GeV}=4.190$ GeV as quoted
in Ref. [37]. For the standard mass from the massive charm Core we get

      In[]:= MassMS["bcm", 5]
      Out[]= 4.178691644159289

with a sub-MeV difference to the reference value
$\overline{m}_{b}=(4.7-0.382-(0.098+0.0072)-(0.03+0.0049))\,\mathrm{GeV}=4.178$
GeV.

####  7.3.4 Top quark pole masses

In the following we investigate the influence of massive bottom and charm
quarks on the perturbative series relating the top quark standard running mass
and the pole mass.

#### Pole-$\overline{\mathrm{MS}}$ mass series dependence on lighter massive
quarks

In Ref. [2] the dependence of the high-order asymptotic series for the top
quark pole-$\overline{\mathrm{MS}}$ mass relation on finite bottom and charm
quark mass was examined and an algorithm to determine explicit analytic
formulae (accounting for the semianalytical results given in Ref. [38]) were
provided, which are implemented in REvolver. Explicit results for the
asymptotic series coefficients beyond ${\cal O}(\alpha_{s}^{4})$ were provided
in Ref. [2] for the scenarios of massless bottom and charm quarks, massive
bottom and massless charm quarks as well as massive bottom and charm quarks.
To compare the results of REvolver to those of Ref. [2] we set the appropriate
relevant parameters for the standard running charm, bottom and top quark
masses as well as for the strong coupling and create the Cores for the three
scenarios

      In[]:= {mcmc, mbmb, mtmt} = {1.3, 4.2, 163};
             {amZ, mZ} = {0.118, 91.187};

      In[]:= CoreDeleteAll[]
             CoreCreate["tbc", 6, {5, amZ, mZ},
              {{4, mcmc, mcmc}, {5, mbmb, mbmb}, {6, mtmt, mtmt}}]
             CoreCreate["tb", 6, {5, amZ, mZ},
              {{5, mbmb, mbmb}, {6, mtmt, mtmt}}]
             CoreCreate["t", 6, {5, amZ, mZ}, {{6, mtmt, mtmt}}]

The perturbative series relating the standard running mass and the pole mass
can be generated as shown in Sec. 7.2.6 and depicted in the upper left panel
of Fig. 4 for massless bottom and charm quarks. In Fig. 6 we show the result
of repeating this procedure for each of the Cores created above. It is clearly
visible that for each additional accounted massive quark flavor the series
diverges faster and consequently the respective pole mass ambiguity increases.

Figure 6: Top quark pole mass $m^{\mathrm{pole}}_{t}$ as a function of the
perturbative order with bottom and charm quark masses set to zero (left blue),
massless charm quark (middle orange) and massive bottom and charm masses
(right green). The conversion is performed from the $\overline{\mathrm{MS}}$
mass. The central dots are obtained for the default renormalization scales for
the strong coupling $\mu=\overline{m}_{t}$, the error bars represent the scale
variation $\overline{m}_{t}/2\leq\mu\leq 2\,\overline{m}_{t}$.

In accordance with Ref. [2] we read out the value of the asymptotic top quark
pole mass, the associated ambiguity and the order of the minimal correction
term from each Core employing the drange approach, which is based on all
series terms in the superasymptotic region, by executing

      In[]:= MassPoleDetailed["t", 6, mtmt, "drange"]
             MassPoleDetailed["tb", 6, mtmt, "drange"]
             MassPoleDetailed["tbc", 6, mtmt, "drange"]
      Out[]= {173.09649964256326, 0.16178589168021063, 8}
      Out[]= {173.17525229067104, 0.2336983055654116, 7}
      Out[]= {173.17973196293923, 0.26632540751468525, 7}

The three results agree with the numbers given in Table 5 of Ref. [2] within
uncertainties888The $10$ MeV discrepancy in the asymptotic pole mass value is
because in Ref. [2], an additional $\lambda$ variation was applied in the
estimate of the asymptotic coefficients and the quoted coefficients were the
central values of the obtained intervals. REvolver does not support this
functionality..

We now return to Ref. [3]. There, finite charm and bottom quark mass effects
on the value of the top quark asymptotic pole mass were quantified by an
approximate and heuristic method utilizing order-dependent light flavor
decoupling in the strong coupling. A numerical value for the resulting
ambiguity of the asymptotic pole mass was quoted in Eq. (5.1) of that paper
using size of the minimal correction of the asymptotic series.

Fixing the value of the top, bottom and charm quark standard running masses to
the values given in Ref. [3] and creating an associated Core

      In[]:= {mtmt, mbmb, mcmc} = {163.508, 4.2, 1.3};

      In[]:= CoreCreate["tbc2", 6, {5, amZdef, mZdef},
              {{4, mcmc, mcmc}, {5, mbmb, mbmb}, {6, mtmt, mtmt}}]

we obtain

      In[]:= MassPoleDetailed["tbc2", 6, mtmt, mtmt, "min"]
      Out[]= {173.61531203486197, 0.10050070355543994, 7}

where we employ the min method to match the approach used in Ref. [3]. The
asymptotic pole mass value and the associated ambiguity are quoted as
$173.667$ GeV and $108$ MeV, respectively, in that reference and differs by
around $50$ and $8$ MeV, respectively, to the result provided by REvolver. The
difference is mainly related to the exact treatment of light massive flavor
corrections provided by REvolver.

#### Normalization of the pole mass renormalon

Finally, we demonstrate REvolver’s capability to compute the normalization
$N_{1/2}^{(n_{\ell})}$ of the pole mass renormalon for a given number
$n_{\ell}$ of massless quark flavors. Reusing the Cores created in one of the
previous examples, the value of the normalization for $n_{\ell}=3,4,5$ can be
extracted with

      In[]:= N12[#] & /@ {"tbc","tb","t"}
      Out[]= {0.5371120004494861, 0.5055628679671248,
              0.4612967890971816}

in excellent agreement with the corresponding values
$N_{1/2}^{(n_{\ell}=3,4,5)}=\\{0.5370\pm 0.0011,\,0.5056\pm 0.0015,\,0.4616\pm
0.0020\\}$ given for $\mu/\overline{m}=1$ in Table 1 of Ref. [3].

In Eq. (2.21) of Ref. [2] the $\lambda$-parameter was varied in the range
$1/2\leq\lambda\leq 2$, and the associated central value and uncertainty of
the renormalon normalization $N_{1/2}^{(n_{\ell})}$ were determined. The
corresponding values for $n_{\ell}=3,4,5$ can be reproduced by REvolver with

      In[]:= list2 = 2^Range[-1, 1, 2/19];

      In[]:= N3Tab = N12["tbc", #] & /@ list2;
             (Max[N3Tab] + {1, -1} * Min[N3Tab])/2
      Out[]= {0.5256087104548577, 0.011913515642471062}

      In[]:= N4Tab = N12["tb", #] & /@ list2;
             (Max[N4Tab] + {1, -1} * Min[N4Tab])/2
      Out[]= {0.4916415824077933, 0.016342715842279004}

      In[]:= N5Tab = N12["t", #] & /@ list2;
             (Max[N5Tab] + {1, -1} * Min[N5Tab])/2
      Out[]= {0.44606744334202203, 0.023867353429116556}

in perfect agreement with the values $N_{1/2}^{(n_{\ell}=3,4,5)}=\\{0.526\pm
0.012,\,0.492\pm 0.016,\,0.446\pm 0.024\\}$, quoted in Eq. (2.21) of Ref. [2].

##  8 Theoretical input

The perturbative coefficients of the QCD $\beta$-function, that is the
renormalization group equation of the $\overline{\rm MS}$ strong coupling
$\alpha_{s}^{(n_{f})}(\mu)$ in a given $n_{f}$ flavor number scheme, are known
to five loops and can be found in Refs. [39, 40, 41, 42, 43, 44, 45, 46]. If
REvolver is used for scenarios with massless quarks only, it is recommended to
cite those references. The perturbative coefficients for the strong coupling
flavor matching relation at the quark thresholds are known up to four loops
and can be found in Refs. [47, 48, 49, 50]. If REvolver is used for strong
coupling evaluations in different flavor number schemes involving massive
quark threshold matching, these references should be cited.

The perturbative coefficients of the relation between the $\overline{\rm MS}$
mass $\overline{m}_{q}^{(n_{f})}(\mu)$, for the case where $n_{f}-1$ is the
number of flavors lighter than the massive quark $q$, and the pole mass
$m^{\rm pole}_{q}$ for the case that all lighter quarks are massless, are
known to four loops and can be found in Refs. [51, 52, 53, 12]. If REvolver is
used to compute relations between the $\overline{\rm MS}$ mass and the pole
mass when all lighter quarks are massless, these references should be
included. The two- and three-loop corrections coming from the masses of
lighter massive quarks can be found in Refs. [52] and [38, 54], respectively.
It would be appropriate to also cite these references, if REvolver is used to
compute lighter massive flavor effects in the relation between the
$\overline{\rm MS}$ mass and the pole mass of a heavy quark. The perturbative
coefficients of the $\overline{\rm MS}$ mass renormalization group equations
(relevant for renormalization scales above the mass of quark $q$) are known to
five loops and can be found in Refs. [55, 56, 57, 58]. These references should
be cited whenever the code is used to compute the running of the
$\overline{\rm MS}$ mass in any flavor scheme that includes the heavy quark
itself. The perturbative coefficients of the flavor matching relations needed
when the running $\overline{\rm MS}$ mass evolves to scales beyond even
heavier massive quarks are known up to four loops and can be found in Refs.
[47, 59]. These references should be included, if the $\overline{\rm MS}$ mass
is evaluated in different flavor number schemes involving threshold matching
associated to even heavier massive quarks.

The MSR mass $m_{q}^{{\rm MSR},(n_{\ell})}(\mu)$ is derived from the
perturbative relation between the pole mass and the standard running mass
$\overline{m}_{q}$ by integrating out the quark flavor $q$, see Sec. 3. Thus,
all perturbative properties of the MSR mass are derived from expressions known
for the $\overline{\rm MS}$ scheme from the references mentioned in the
previous paragraph. The MSR mass concept was first suggested in Ref. [6] and
allows to consistently consider renormalization group evolution for scales
$\mu$ below the quark mass $m_{q}$, called “R-evolution”. The perturbative
coefficients of (i) the relation between the MSR mass $m_{q}^{{\rm
MSR},(n_{\ell})}(\mu)$, where $n_{\ell}$ is smaller or equal to the number of
flavors lighter than the massive quark $q$, and the pole mass $m^{\rm
pole}_{q}$, (ii) the matching relation of the MSR mass to the $\overline{\rm
MS}$ mass, and (iii) its $n_{\ell}$ flavor-number-dependent renormalization
group equation for the case that all lighter quarks are massless (which are
all known to four loops) can be found in Ref. [1]. The MSR mass scheme
implemented in REvolver is called the “natural” MSR mass scheme in Ref. [1].
The two- and three-loop corrections coming from the masses of lighter massive
quarks including the flavor matching relations, when the MSR mass evolution
crosses a lighter massive quark threshold, can be found in Ref. [2]. The
results for the lighter massive quark effects given in Ref. [2] involve a
parametrization of the results provided in Ref. [38] (in terms of coefficients
$\delta_{2}$ and $\delta_{Q,n}^{(q,q^{\prime},\ldots)}$) that accounts for
theoretical interrelations of the different contributions not considered in
Ref. [38], and allows for a straightforward generalization for the case of
multiple lighter massive quark flavors. The parametrization is also fully
compatible with the recent analytic updates on the three-loop corrections
given Ref. [54] within a few MeV. The expression for $\delta_{2}$ used in
REvolver has been given in Ref. [2]. The concrete expressions for the other
coefficients $\delta_{Q,n}^{(q,q^{\prime},\ldots)}$ have been taken from Ref.
[60]. If REvolver is used to calculate the MSR mass of a massive quark in the
case that all lighter quarks are massless, Refs. [1, 2] should be cited. If
REvolver is used to calculate the MSR mass in different flavor number schemes
involving flavor threshold matching as light massive quark thresholds Ref. [2]
should be cited. The possibility to resum (large) logarithms in the relation
between short-distance quark mass definitions involving scales below the quark
mass scale (e.g. concerning the ratio of the heavy quark mass and a lighter
quark mass) is currently only provided through the MSR mass, which is
therefore the primary short-distance mass scheme for scales below the quark
masses contained in the REvolver core objects. If REvolver is used to resum
such logarithms in the relation of short-distance masses (even not involving
the MSR directly), this MSR mass functionality is used as an intermediate
step, and it is appropriate to cite Ref. [1, 6].

The quark mass dependence of the formulae of the flavor threshold matching
relations for the flavor-number-dependent strong coupling
$\alpha_{s}^{(n_{f})}(\mu)$, the $\overline{\rm MS}$ mass
$\overline{m}_{q}^{(n_{f})}(\mu)$ and the MSR mass $m_{q}^{{\rm
MSR},(n_{\ell})}(\mu)$ implemented in REvolver are expressed in terms of the
standard running mass $\overline{m}_{q}$. Changing the scheme of the quark
mass in these matching relations is not supported in REvolver because the
corresponding numerical impact is tiny and negligible for practical
applications where large logarithmic corrections are properly resummed.

The perturbative coefficients of the relation between the kinetic and pole
mass, assuming lighter massive quarks to be massless, can be found in Refs.
[17, 21]. The 1S mass scheme was first suggested in Refs. [15, 14]. The
perturbative coefficients for the 1S mass scheme have been derived in Refs.
[15, 14, 61, 19, 62, 63, 64]. The perturbative coefficient for the PS mass
scheme, first suggested in Ref. [18], can be obtained from Refs. [18, 19]
using analytic results for the static potential as given in Refs. [65, 66, 67,
68, 69, 70, 71]. The relation connecting the RS and pole mass schemes is
described in Refs. [30, 20] and references therein. Whenever values for any of
these masses are obtained through REvolver, the corresponding references
should be mentioned. The lighter massive flavor corrections in the relation to
the pole mass are known for the kinetic [22], 1S [72, 73] and PS masses [74].
For the 1S mass REvolver currently only accounts for the fixed-order
corrections in the relation to the pole mass. Therefore, if lighter quarks are
considered massive, the appropriate references should be cited.

To compute the RGI mass [13] the standard formula is evaluated as described
below Eq. (28).

##  9 Summary

In this article we have presented REvolver, a C++ library for carrying out
state-of-the-art renormalization group evolution and flavor matching for the
QCD coupling and quark masses, and conversion between the most common quark
mass renormalization schemes. For short-distance quark masses the achievable
precision is at the level of $10$ to $20$ MeV. In addition to similar
libraries that are already available REvolver offers the Core concept, that
allows to define and manage different physical scenarios for coupling and
quark mass values, and to carry out renormalization group summation of
logarithms from scales above and below the quark masses. REvolver supports, in
particular, the summation of logarithms described by the R-evolution equation,
which describes the linear scale evolution characteristic to low-scale short-
distance masses, and it accounts for the flavor threshold corrections that
arise in this linear scale evolution. Furthermore, REvolver provides access to
the asymptotic perturbative relations to the pole mass to in principle any
perturbative order, and it provides quasi-exact solutions to the
renormalization group equations (i.e. exact up to machine-precision) for
complex renormalization scales. The REvolver library can be also accessed
through Mathematica and Python interfaces. We have provided a large number of
examples for the most common applications, which are shipped along with the
library in the form of Mathematica and Jupyter notebooks, as well as a C++
program.

## Acknowledgments

This work was supported in part by FWF Austrian Science Fund under the Project
No. P28535-N27, the Spanish MINECO Ramón y Cajal program (RYC-2014-16022), the
MECD grants FPA2016-78645-P and PID2019-105439GB-C22, the IFT Centro de
Excelencia Severo Ochoa Program under Grant SEV-2012-0249, the EU STRONG-2020
project under the program H2020-INFRAIA-2018-1, grant agreement No. 824093 and
the COST Action CA16201 PARTICLEFACE. CL is supported by the FWF Doctoral
Program “Particles and Interactions” No. W1252-N27. We thank D. Boito and M.
Steinhauser for explanations concerning the content of their respective
articles.

##  Appendix A Algorithms to Solve the Differential Equations

In this appendix some details on the algorithms used in REvolver are provided.

###  A.1 Core Creation

For the creation of a Core, values for the strong coupling and the running
masses of the $n_{q}$ massive quarks at specified scales and in specified
flavor-number schemes, as well as the total number of flavors $n_{T}$, are
provided by the user.999The Core contains $n_{T}-n_{q}$ massless quarks, and
the standard running mass of the $n$-th massive quark is defined with
$n_{T}+n-n_{q}$ active flavors. In the Core the strong coupling and the
running masses can be evaluated in a $n_{f}$ flavor number scheme with
$n_{T}-n_{q}\leq n_{f}\leq n_{T}$, i.e. for $n_{q}+1$ different flavor number
schemes. Let us call the corresponding tuples
$\\{n_{\alpha},\alpha_{s}^{(n_{\alpha})}(\mu_{\alpha})$, $\mu_{\alpha}\\}$ for
the strong coupling and $\\{k_{n},m_{q_{n}}^{(k_{n})}(R_{n}),R_{n}\\}$, with
$\max[k_{n},n_{\alpha}]\leq n_{T}$, $1\leq n\leq n_{q}$, for the $n_{q}$
massive quarks, where $q_{1}$ refers to the lightest and $q_{n_{q}}$ to the
heaviest massive quark. In a first step REvolver determines the values of the
standard running masses $\overline{m}_{q_{n}}$. Subsequently, the following
strong coupling and running mass values at each flavor threshold matching
scale $\mu_{n}$ are obtained: $\alpha_{s}^{(n_{T}+n-1-n_{q})}(\mu_{n})$,
$\alpha_{s}^{(n_{T}+n-n_{q})}(\mu_{n})$ and
$m_{n}^{(n_{T}+n-1-n_{q})}(\mu_{n})$, $m_{n}^{(n_{T}+n-n_{q})}(\mu_{n})$.
These numbers serve as the initial conditions to compute $\alpha_{s}$ and the
running masses $m_{q_{n}}$ at any renormalization scale and in any specified
flavor-number scheme. We have implemented into our code a fast multi-
dimensional recursive algorithm that solves the coupled system of equations to
the specified precision.

The determination of the standard running masses $\overline{m}_{q_{n}}$ in the
algorithm is particularly important since the flavor matching scales $\mu_{n}$
are specified as dimensionless coefficients times $\overline{m}_{q_{n}}$. If
at the $i$-th iteration the standard running masses are
$[\overline{m}_{q_{n}}]_{i}$, the numerical values
$[\overline{m}_{q_{n}}]_{i+1}$ at step $(i+1)$ are computed evolving (and
matching if necessary) the quark masses $m_{q_{n}}^{(k_{n})}(R_{n})$ from the
user-specified $\mu=R_{n}$, for which the numerical value of the mass was
provided, to $\mu=[\overline{m}_{q_{n}}]_{i}$. To do that, before running the
masses, the program computes the values of $\alpha_{s}$ at the various
thresholds $\mu_{n}$ using $[\overline{m}_{q_{n}}]_{i}$ as well. After each
iteration is carried out, the new standard running masses are compared to the
previous ones, and if their largest relative deviation is smaller than the
required precision, the process stops. In the initial step it is assumed that
$[\overline{m}_{q_{n}}]_{i=0}=m_{q_{n}}^{(k_{n})}(R_{n})$. We checked that the
algorithm is very robust and converges quickly even for very extreme starting
conditions such as a charm quark mass defined with $k_{1}=6$ active flavors
taking for $R_{1}$ a very large scale, together with a top quark mass defined
at $R_{3}=1\,$GeV with $k_{3}=3$ active flavors. There is the possibility that
the algorithm does not converge, if a particular sequence of
$[\overline{m}_{q_{n}}]_{i}$-values repeats itself. This is avoided by
adopting, after the $20$-th iteration ($i\geq 21$) the linear combination
$(1-\mathtt{damp})[\overline{m}_{q_{n}}]_{i}+\mathtt{damp}[\overline{m}_{q_{n}}]_{i-1}$,
with damp being a random number between $0$ and $0.2$, as the outcome of the
$i$-th iteration.

To compute the strong coupling flavor matching REvolver takes as exact the
upward relation (at the user-specified loop order)

$\alpha^{(n_{\ell}+1)}_{s}(\mu_{n})=\alpha^{(n_{\ell})}_{s}(\mu_{n})\Biggl{\\{}1+\sum_{i=1}\biggl{[}\frac{\alpha^{(n_{\ell})}_{s}(\mu_{n})}{4\pi}\biggr{]}^{i}\xi_{i}(n_{\ell})\Biggr{\\}},$
(9)

where $\xi_{i}$ may depend on $\log(\overline{m}_{q_{n}}/\mu_{n})$. This means
that the values for $\alpha_{s}^{\\!(n_{\ell})}\\!(\mu_{n})$ and
$\alpha_{s}^{(n_{\ell}+1)}(\mu_{n})$ in REvolver always satisfy exactly Eq.
(9). If $\alpha_{s}^{(n_{\ell})}(\mu_{n})$ is given, the value of
$\alpha_{s}^{(n_{\ell}+1)}(\mu_{n})$ is computed directly from Eq. (9). If, on
the contrary, $\alpha_{s}^{(n_{\ell}+1)}(\mu_{n})$ is given and
$\alpha_{s}^{(n_{\ell})}(\mu_{n})$ shall be obtained, our program numerically
inverts Eq. (9) in an iterative way as described in the following. For
simplicity we omit the argument $\mu_{n}$ and write
$[\alpha^{(n_{\ell})}]_{k}$ for the $k$-th iteration value. The value at the
$(k+1)$-th step the reads

$[\alpha^{(n_{\ell})}]_{k+1}=\frac{\alpha_{s}^{(n_{\ell}+1)}(\mu_{n})}{1+\sum_{i=1}\Bigl{(}\frac{[\alpha^{(n_{\ell})}]_{k}}{4\pi}\Bigr{)}^{\\!i}\,\xi_{i}(n_{\ell})}\,,$
(10)

where for the first iteration
$[\alpha^{(n_{\ell})}]_{k=0}=\alpha_{s}^{(n_{\ell}+1)}(\mu_{n})$ is adopted.
The algorithm converges very quickly, even for
$\alpha_{s}^{(n_{\ell}+1)}(\mu_{n})$ values as large as $0.69$. For smaller
values of the strong coupling, and specially if
$\mu_{n}\simeq\overline{m}_{q_{n}}$, the RHS of Eq. (10) is very weakly
depending on $\alpha^{(n_{\ell})}_{k}$ (since the 1-loop term is very small or
vanishing) and only a few iterations are necessary. If
$\alpha_{s}^{(n_{\ell}+1)}(\mu_{n})$ is larger than $0.69$, there is no
solution to Eq. (9) and the iterative procedure will fail. In those cases the
program will return the leading order solution
$\alpha_{s}^{(n_{\ell})}(\mu_{n})=\alpha_{s}^{(n_{\ell}+1)}(\mu_{n})$. If the
solution exists, the algorithm converges as well, because the slope of the
left-hand-side is smaller than unity in absolute value.

###  A.2 Strong Coupling

For the QCD $\beta$-function at $N$-loop order (where $N$ is the integer value
specified in the variable runAlpha), see Eq. (6), the coefficients $b_{i}$ are
defined by

$\beta^{N}_{\rm{QCD}}(\alpha_{s})\equiv-\frac{\alpha^{2}_{s}}{2\pi}\beta_{0}\biggl{[}1+\sum_{i=1}^{N-1}\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}^{\\!\\!i}\,b_{i}\biggr{]},\qquad
b_{i}\equiv\frac{\beta_{i}}{\beta_{0}}\,.$ (11)

The algorithms implemented in REvolver provide solutions for the strong
coupling differential equation based on the $N$-loop $\beta$-function that are
numerically exact (within the specified precision) and do not involve any
additional approximation related to $N$-loop precision. The associated
coefficients $c^{N}_{n}$, which parametrize the associated inverse QCD
$\beta$-function, are defined by

$\frac{1}{\beta^{N}_{\rm{QCD}}(\alpha_{s})}=-\frac{2\pi}{\beta_{0}}\frac{1}{\alpha^{2}_{s}}\biggl{[}1+\sum_{i=1}c^{N}_{i}\Bigl{(}\frac{\alpha_{s}}{4\pi}\Bigr{)}^{\\!\\!i}\,\biggr{]},$
(12)

where, except for $N=1$, the sum extends to infinity101010The notation where
the upper limit of the sum is omitted signifies that the REvolver algorithms
employ the number of terms mandatory to reach numerical results that are exact
within the specified precision. The analogous notation is used for other
formulae shown below. even though there is a finite limit $N$ in the sum
defining $\beta^{N}_{\rm{QCD}}$ (for $N=1$ one has $c^{N}_{i}=0$ for all
$i\geq 1$). The (infinite set of) $c^{N}_{i}$ coefficients can be computed in
terms of $N-1$ coefficients $b_{1},\ldots b_{N-1}$ using the recursive
relation shown in Eq. (13). The convergence radius of the series is set by the
distance from the origin to the nearest pole in the complex plane, which is
equivalent to the smallest module of the $\beta$-function roots (excluding the
double pole at the origin), which are in general complex numbers. Hence the
convergence radius depends both on the number of flavors and the loop order
$N$. It increases with the number of flavors and decreases with the loop order
$N$, and in all cases is much larger than any physical value of $\alpha_{s}$
that may arise in phenomenological applications.

For convenience we define $b_{0}=c^{N}_{0}=1$, which allows to write

$c^{N}_{n+1}=-\\!\\!\\!\\!\sum_{i=1}^{\min(N-1,n+1)}\\!\\!c^{N}_{n+1-i}\,b_{i}\,.$
(13)

Many analytic formulae implemented in REvolver use the $t$-variable formalism
[6, 1] which is based on the definition

$t=-\frac{2\pi}{\beta_{0}}\frac{1}{\alpha_{s}}\,,$ (14)

which gives the relation

${\rm d}\ln(\mu)=-\hat{b}^{N}\\!(t)\,{\rm d}t\,,$ (15)

for the renormalization scale $\mu$ where

$\hat{b}^{N}\\!(t)=-\frac{2\pi}{\beta_{0}\,t^{2}\beta^{N}_{\rm{QCD}}\bigl{(}-\frac{2\pi}{\beta_{0}}\frac{1}{t}\bigr{)}}\equiv
1+\sum_{i=1}\hat{b}^{N}_{i}\,t^{-i}\,,$ (16)

with $\hat{b}^{N}_{n}=(-1)^{n}c^{N}_{n}/(2\beta_{0})^{n}$, which implies
$\hat{b}^{N}_{0}=1$.111111The $\hat{b}^{N}_{i}$ should not be confused with
the $b_{i}$ coefficients defined in Eq. (11). The $c_{n}^{N}$ and
$\hat{b}_{n}^{N}$ coefficients are obtained numerically whenever needed and
stored in a member vector (whose length is extended if necessary) such that
they are not re-computed again when used later. For the integral of Eq. (15)
it is useful to define the following polynomial in $1/t$ expressed in terms of
$\hat{b}_{i}^{N}$:

$\tilde{G}^{N}\\!(t)\equiv\int_{-\infty}^{t}\text{d}t^{\prime}\biggl{[}\hat{b}^{N}(t^{\prime})-1-\frac{\hat{b}^{N}_{1}}{t^{\prime}}\,\biggr{]}=-\sum_{n=1}\frac{\hat{b}^{N}_{n+1}}{n}\,t^{-n}\,.$
(17)

For $N=1$ one has $\tilde{G}^{N}\\!(t)=0$ while for $N\geq 2$ the sum in
$\tilde{G}^{N}\\!(t)$ always extends to infinity. The exponential of this
expression defines the coefficients $g_{\ell}^{N}$ upon reexpansion in powers
of $1/t$, $e^{\tilde{G}^{N}(t)}\equiv\sum_{\ell=0}g^{N}_{\ell}(-t)^{-\ell}$,
which can be computed with the recursive relation

$g^{N}_{n+1}=\frac{1}{n+1}\sum_{i=0}^{n}(-1)^{i}\,\hat{b}^{N}_{i+2}\,g^{N}_{n-i}\,,$
(18)

where ${g}^{N}_{0}=1$. It is also necessary to define the coefficients that
result from expanding
$e^{-\tilde{G}^{N}(t)}\equiv\sum_{\ell=0}\tilde{g}^{N}_{\ell}(-t)^{-\ell}$ in
an analogous way. They can be computed with a very similar recursive algorithm

$\tilde{g}^{N}_{n+1}=-\frac{1}{n+1}\sum_{i=0}^{n}(-1)^{i}\,\hat{b}^{N}_{i+2}\,\tilde{g}^{N}_{n-i}\,,$
(19)

where $\tilde{g}^{N}_{0}=1$. Once again one has
$g_{n}^{N}=\tilde{g}_{n}^{N}=0$ for $n\geq 1$ when $N=1$. Both $g$ and
$\tilde{g}$ coefficients are computed when first needed and then are
conveniently stored in member vectors for later use.

The QCD scale $\Lambda_{\rm QCD,t}$ in the t-scheme [1] is defined by

$\Lambda^{N}_{\rm{QCD,t}}=\mu\,e^{G^{N}\\!\bigl{(}-\frac{2\pi}{\beta_{0}}\frac{1}{\alpha_{s}(\mu)}\bigr{)}}\,,$
(20)

where $G^{N}\\!(t)=t+\hat{b}^{N}_{1}\log(-t)+\tilde{G}^{N}(t)$ (where we
recall that $\hat{b}_{i>0}^{N}=0$ if $N=1$). The summation for the function
$\tilde{G}^{N}$ in Eq. (17) is terminated after the relative size of the last
term compared to the associated partial sum is smaller than the specified
precision. Since in the evaluation of Eq. (20) the exact solution of the
strong coupling differential equation with $N$-loop $\beta$-function
$\beta^{N}_{\rm QCD}$ (see next section) for $\alpha_{s}(\mu)$ is used, the
numerical value of $\Lambda^{N}_{\rm QCD}$ does (within the specified
precision) not depend on $\mu$. For the evaluation in REvolver the value
$\mu=100$ GeV is always used. The values for the QCD scale are computed for
each number of active flavors and stored as members of the class. For cases in
which $\alpha_{s}$ is very large, the sum $G^{N}$ may not converge within the
specified precision after adding $200$ terms and REvolver will return NaN (see
below why the limit is set to $200$). Again, this only happens for unphysical
$\alpha_{s}$ values.

For the QCD scale in the $\overline{\mathrm{MS}}$ definition, $\Lambda_{\rm
QCD,\overline{\mathrm{MS}}}$ (following [23]), we use the formula

$\Lambda^{N}_{\rm{QCD,\overline{\mathrm{MS}}}}=2^{\hat{b}^{N}_{1}}\,\Lambda^{N}_{\rm{QCD,t}}\,,$
(21)

therefore both schemes coincide for $N=1$.

The algorithm used by REvolver to solve the differential equation for the
strong coupling is based on variable separation. Defining
$\ell_{\mu}\equiv\log(\mu/\mu_{0})$, $\alpha_{\mu}\equiv\alpha_{s}(\mu)$,
$\alpha_{0}\equiv\alpha_{s}(\mu_{0})$ and $a\equiv\alpha_{s}/(4\pi)$ one can
expand $1/\beta_{\rm QCD}^{N}$ as in Eq. (12) and integrate term by term. This
gives

$\ell_{\mu}=\\!\\!\int_{\alpha_{0}}^{\alpha_{\mu}}\\!\\!\\!\frac{\text{d}\alpha}{\beta^{N}_{\rm{QCD}}(\alpha)}=\\!-\frac{1}{2\beta_{0}}\biggl{[}\frac{1}{a_{0}}-\frac{1}{a_{\mu}}+c_{1}^{N}\log\\!\biggl{(}\frac{a_{\mu}}{a_{0}}\biggr{)}+\\!\sum_{i=1}\frac{c^{N}_{i+1}}{i}(a_{\mu}^{i}-a_{0}^{i})\biggr{]}.$
(22)

where the sum on the RHS is carried out until the last term added is smaller
than the specified precision.

We now describe the algorithm to determine $a_{\mu}$ from the algebraic
equation (22). At leading-log (LL), i.e. using the 1-loop $\beta$-function
($N=1$), one has $\beta_{n>0}=0$ and the solution is $a^{\rm
LL}_{\mu}=a_{0}/(1+2\beta_{0}a_{0}\ell_{\mu})$. For $N\geq 2$ the LL solution
is used to rewrite Eq. (22) in the form

$\frac{1}{a_{\mu}}=\frac{1}{a^{\rm{LL}}_{\mu}}+c_{1}^{N}\log\biggl{(}\frac{a_{\mu}}{a_{0}}\biggr{)}+\sum_{i=1}\frac{c^{N}_{i+1}}{i}(a_{\mu}^{i}-a_{0}^{i})\,,$
(23)

which is solved recursively. In this iterative procedure, the $(n+1)$-th
iteration value for the $N$-loop strong coupling is obtained from the $n$-th
iteration one by the relation

$[a_{\mu}]_{n+1}=\frac{1}{\frac{1}{a^{\rm
LL}_{\mu}}+c_{1}^{N}\log\Bigl{(}\frac{[a_{\mu}]_{n}}{a_{0}}\Bigr{)}+\sum_{i=1}\frac{c^{N}_{i+1}}{i}[([a_{\mu}]_{n})^{i}-a_{0}^{i}]}\,.$
(24)

With the initial choice $[a_{\mu}]^{N}_{0}\equiv a^{\rm{LL}}_{\mu}$ it is
ensured that the correct solution is obtained in all cases. The iterative
method converges provided that (i) the solution exists and (ii) the (absolute
value of the) slope of the right-hand-side function is smaller than $1$. The
second condition is always satisfied, and the former will be discussed later
in this paragraph. Since $\beta_{1}>0$ we have that $a_{\mu}^{\rm
LL}<a_{\mu}$, and the solution beyond 1-loop order is approached from below.
The iterative procedure is carried out until the numerical value of $a_{\mu}$
does not change within the specified precision, with a maximum allowed number
of $200$ iterations.

The convergence radius for the infinite sum over $c_{i+1}^{N}$ in Eq. (24) is
the same as for the inverse $\beta$-function in Eq. (12). However, due to
technical limitations related to double-precision floating numbers, REvolver
cannot carry out the sum for $\alpha_{s}$ values arbitrarily close to the
convergence radius within the specified precision. If after adding $200$ terms
this precision is not met, the recursive procedure to obtain $\alpha_{s}$ will
stop and NaN will be returned. This happens for $\alpha_{s}\gtrsim 0.7$, hence
outside the range of any physical application.

###  A.3 $\overline{\rm MS}$ and RGI masses

The $K$-loop renormalization group equation (where $K$ is the integer value
specified in the variable runMSbar) for the $\overline{\rm MS}$ quark mass
$\overline{m}_{q}(\mu)$ has the form

$\frac{\mathrm{\text{d}}\,\overline{m}_{q}(\mu)}{\mathrm{\text{d}}\ln(\mu)}=2\overline{m}_{q}(\mu)\gamma_{m}^{K}(\alpha_{s})=2\overline{m}_{q}(\mu)\sum_{n=0}^{K-1}\gamma_{n}\biggr{[}\frac{\alpha_{s}(\mu)}{4\pi}\biggl{]}^{\\!n+1}\,,$
(25)

which implies logarithmic scale evolution. The REvolver algorithm provides
(within the user specified precision) the exact solution of Eq. (25) employing
the $N$-loop strong coupling described in Sec. A.2. Changing the integration
variable from $\ln(\mu)$ to $\alpha_{s}$ followed by separation of variables,
and using the mass $\overline{m}_{q}(\mu_{0})$ at the scale $\mu_{0}$ as the
initial value, the analytic solution reads

$\\!\\!\\!\\!\overline{m}_{q}(\mu)\\!=\\!\overline{m}_{q}(\mu_{0})\exp\\!\biggr{[}2\\!\\!\int^{\alpha_{\mu}}_{\alpha_{0}}\\!\\!\text{d}\alpha\,\frac{\gamma_{m}^{K}\\!(\alpha)}{\beta^{N}_{\rm{QCD}}(\alpha)}\biggl{]}\equiv\overline{m}_{q}(\mu_{0})\exp\\!\big{[}\tilde{\omega}^{(N,K)}(\mu_{0},\mu)\big{]},$
(26)

which defines the running kernel $\tilde{\omega}^{(N,K)}(\mu_{0},\mu)$. The
REvolver algorithm determines $\tilde{\omega}^{(N,K)}(\mu_{0},\mu)$ from the
Taylor series of the integrand in $\alpha$ where each term is integrated
individually. The result can be written in the form [
$a_{i}=\alpha_{s}(\mu_{i})/(4\pi)$ ]

$\\!\\!\\!\tilde{\omega}^{(N,K)}(\mu_{1},\mu_{2})\\!=\\!-\frac{1}{\beta_{0}}\\!\biggl{[}\gamma_{0}\log\biggl{(}\frac{a_{2}}{a_{1}}\biggr{)}\\!+(1-\delta_{N,1}\delta_{K,1})\\!\sum_{n=1}(a_{2}^{n}-a_{1}^{n})\,d^{(N,K)}_{n}\biggr{]},$
(27)

where the coefficients $d^{(N,K)}_{n}$, defined for $n>0$, are non-zero if
either $N$ or $K$ are larger than $1$. We also have
$d^{(1,K)}_{n}=\gamma_{n}/\beta_{0}$ if $n<K$ and zero otherwise. They are
given by the following expression

$d^{(N,K)}_{n}\equiv\frac{1}{n}\\!\sum^{\min(n,K-1)}_{i=0}c_{n-i}^{N}\,\gamma_{i}\,,$
(28)

where the terms $c^{N}_{i}$ arise in the expansion of the inverse
$\beta$-function (12). REvolver carries out the sum in Eq. (27) until the
relative size of the last computed term with respect to the associated partial
sum is below the user specified precision. The convergence radius of the
infinite sum is the same as for the inverse $\beta$-function in Eq. (12). In
practice, if convergence is found when computing $\alpha_{s}$, the sum for the
computation of $\tilde{\omega}^{(N,K)}$ is going to be convergent as well. If
no convergence is found, NaN is returned.

The RGI-mass $\hat{m}_{q}$ is obtained from the $\overline{\mathrm{MS}}$ mass
$\overline{m}_{q}(\mu)$ by removing its scale-dependence writing
$\tilde{\omega}^{(N,K)}(\mu_{0},\mu)=\tilde{\omega}^{(N,K)}(\mu)-\tilde{\omega}^{(N,K)}(\mu_{0})$
with

$\tilde{\omega}^{(N,K)}(\mu)=-\frac{1}{\beta_{0}}\biggl{[}\gamma_{0}\log\\!\biggl{(}\frac{a_{\mu}}{\pi}\biggr{)}+(1-\delta_{N,1}\delta_{K,1})\sum_{n=1}d^{(N,K)}_{n}a_{\mu}^{n}\biggr{]}.$
(29)

This leads to the definition
$\hat{m}_{q}\equiv\overline{m}_{q}(\mu)\exp[-\tilde{\omega}(\mu)]$. REvolver
uses the $N$-loop exact solution for $a_{\mu}$ and the $K$-loop solution for
$m_{q}(\mu)$, which ensures that the value obtained for the RGI-mass is
$\mu$-independent. In REvolver the computation of $\tilde{\omega}(\mu)$ is
carried out in the same way as for the $\overline{\mathrm{MS}}$ mass explained
above. For the numerical calculation $\mu$ is set to the standard running
mass, $\mu=\overline{m}_{q}$.

###  A.4 MSR mass

The $M$-loop renormalization group equation for the MSR mass $m_{q}^{\rm
MSR}(R)$ (where $M$ is the integer value specified in the variable runMSR) has
the form

$\frac{\text{d}m_{q}^{\rm{MSR}}(R)}{\text{d}\ln(R)}=-R\,\gamma_{R}^{M}[\alpha_{s}(R)]=-R\sum_{n=0}^{M-1}\gamma^{R}_{n}\biggl{[}\frac{\alpha_{s}(R)}{4\pi}\biggr{]}^{n+1}\,,$
(30)

and implies logarithmic and linear scale dependence. The coefficients
$\gamma^{R}_{n}$ are obtained from the series relating the pole and MSR masses
using Eq. (3.2) of Ref. [1]. The REvolver algorithm provides (within the user
specified precision) the exact solution of Eq. (30) employing the $N$-loop
strong coupling described in Sec. A.2. In the following we provide details of
the algorithm’s mathematical derivation.

Switching to the $t$-variable formalism [6, 1] one obtains re-scaled
$\gamma_{R}$ coefficients:

$\gamma_{R}^{M}\\!\biggl{(}-\frac{2\pi}{\beta_{0}}\frac{1}{t}\biggr{)}\equiv\sum_{n=0}^{M-1}\tilde{\gamma}^{R}_{n}(-t)^{-n-1},\qquad\tilde{\gamma}^{R}_{n}=\frac{\gamma^{R}_{n}}{(2\beta_{0})^{n+1}}\,.$
(31)

These can be used to define the series coefficients $S_{j}^{(N,M)}$ in the
expansion of the product function
$\gamma_{R}^{M}\\!(-\frac{2\pi}{\beta_{0}}\frac{1}{t})\hat{b}^{N}\\!(t)e^{-\tilde{G}^{N}\\!(t)}\equiv\sum_{j=0}S^{(N,M)}_{j}(-t)^{-j-1}$
in powers of $1/t$, see Sec. A.2. The coefficients have the form

$S^{(N,M)}_{j}=\\!\\!\sum_{k=0}^{\min(M-1,j)}\\!\\!\tilde{\gamma}^{R}_{k}\sum_{i=0}^{j-k}(-1)^{i}\,\tilde{b}^{N}_{i}\,\tilde{g}^{N}_{j-i-k}\,,$
(32)

which implies that $S^{(1,M)}_{k}=\tilde{\gamma}^{R}_{k}$ for $k<M$ and
$S^{(1,M)}_{k}=0$ for $k\geq M$. The difference of MSR masses at the
renormalization scales $R_{1}$ and $R_{2}$,
$\Delta_{\rm{MSR}}(R_{2},R_{1})\equiv
m_{q}^{\rm{MSR}}(R_{2})-m_{q}^{\rm{MSR}}(R_{1})$, can be computed from Eq.
(30) by integrating over $\alpha_{s}$ (or $t$) using the relation
$R=\Lambda^{N}_{\rm{QCD}}e^{-G^{N}\\!\bigl{[}-\frac{2\pi}{\beta_{0}}\frac{1}{\alpha_{s}^{N}\\!(R)}\bigr{]}}$

$\displaystyle\Delta^{(M,N)}_{\rm{MSR}}(R_{2},R_{1})=-\Lambda^{N}_{\rm{QCD}}\\!\int_{\alpha_{1}}^{\alpha_{2}}\\!\\!\frac{\text{d}\alpha}{\beta^{N}(\alpha)}e^{-G^{N}\\!\bigl{[}-\frac{2\pi}{\beta_{0}}\frac{1}{\alpha}\bigr{]}}\gamma^{M}_{R}(\alpha)$
(33)
$\displaystyle=\Lambda^{N}_{\mathrm{QCD}}\int_{t_{1}}^{t_{2}}\text{d}t\,\hat{b}^{N}(t)e^{-\tilde{G}^{N}\\!(t)}\gamma_{R}\biggl{(}-\frac{2\pi}{\beta_{0}}\frac{1}{t}\biggr{)}e^{-t}(-t)^{-\hat{b}_{1}^{N}}$
$\displaystyle=\Lambda^{N}_{\mathrm{QCD}}\sum_{j=0}S^{(N,M)}_{j}\\!\\!\int_{t_{1}}^{t_{2}}\\!\text{d}t\,e^{-t}(-t)^{-1-\hat{b}_{1}^{N}-j}$
$\displaystyle\equiv\Lambda^{N}_{\mathrm{QCD}}\sum_{j=0}S^{(N,M)}_{j}[I(\hat{b}^{N}_{1},j,t_{2})-I(\hat{b}^{N}_{1},j,t_{1})]\,,$

where in the first line we have used
$\beta^{N}(\alpha)\text{d}R=\Lambda^{N}_{\rm{QCD}}e^{-G^{N}\\!\bigl{[}-\frac{2\pi}{\beta_{0}}\frac{1}{\alpha_{s}^{N}\\!(R)}\bigr{]}}\text{d}\alpha$,
in the second we have switched variables to $t$ and in the third we used the
definition of the coefficients in Eq. (32). Using integration by parts
repeatedly one can write $I(\hat{b}^{N}_{1},j,t)$ in terms of
$I(\hat{b}^{N}_{1},0,t)$

$I(\hat{b}^{N}_{1},j,t)=\frac{1}{(1+\hat{b}^{N}_{1})_{j}}\\!\biggl{[}I(\hat{b}^{N}_{1},0,t)+e^{-t}\sum_{i=1}^{j}(1+\hat{b}^{N}_{1})_{i-1}(-t)^{-i-\hat{b}^{N}_{1}}\biggr{]},$
(34)

with $(a)_{n}=\Gamma(a+n)/\Gamma(a)$ being the Pochhammer symbol. The
exponential contained in $I(\hat{b}^{N}_{1},0,t)$ can be expanded in powers of
$t$ (with infinite convergence radius) and each term can be integrated
separately. In practice, the algorithm computes the difference

$I(\hat{b}^{N}_{1},0,t_{2})-I(\hat{b}^{N}_{1},0,t_{1})=\sum_{i=0}\frac{1}{i!}\frac{(-t_{2})^{i-\hat{b}^{N}_{1}}-(-t_{1})^{i-\hat{b}^{N}_{1}}}{\hat{b}^{N}_{1}-i}\,,$
(35)

terminating the infinite sum once the specified precision is achieved, and the
result for this difference is saved, such that it needs to be computed only
once for each computation of $\Delta^{(M,N)}_{\rm{MSR}}(R_{2},R_{1})$. In case
$\hat{b}_{1}^{N}=0$ the $i=0$ term in Eq. (35) is replaced by
$\ln(t_{1}/t_{2})$. Furthermore, the values of the Pochhammer symbols in the
formula Eq. (34) are computed only once in one Core object and stored for
later use. Likewise, since the sum in Eq. (34) depends on $j$ only through its
upper limit, its value is stored and updated as $j$ grows by one unit.

REvolver computes the (infinite) sum over $j$ in Eq. (33) by truncating it
once the desired precision is achieved. The convergence radius of the $j$ sum
coincides with that of the inverse $\beta$-function explained in Sec. A.2. So,
if convergence is found for the evolution of the strong coupling $\alpha_{s}$,
the sum in Eq. (33) converges as well. If no convergence is found for the sum
or for the strong coupling, NaN is returned. REvolver stores in a member
vector the values of $S_{j}^{(M,N)}$ after they are computed the first time,
such that their values can be recycled in later evaluations of R-evolution.

###  A.5 Asymptotic Pole Mass Coefficients

For the asymptotic formulae for the ${\cal O}(\Lambda_{\rm QCD})$ renormalon-
dominated perturbative coefficients of the relation between the MSR and pole
masses (see Eqs. (2.1) and (2.3) in Ref. [1]), always defined in the flavor
number scheme where all massive quarks are integrated out, we use the
expression

$a_{n}=(2\beta_{0})^{n}\\!\sum_{k=0}^{k_{\rm
max}[n-1]}\\!S^{(N,M)}_{k}\\!\\!\sum_{\ell=0}^{\ell_{\rm
max}[n-1-k]}\\!\\!g^{N}_{\ell}(1+\hat{b}_{1}^{N}+k)_{n-1-\ell-k}\,,$ (36)

where

$\displaystyle k_{\rm max}[m]\,=\,$ $\displaystyle{\rm min}[m,M-1]\,,$ (37)
$\displaystyle\ell_{\rm max}[m]\,=\,$ $\displaystyle{\rm max}[{\rm
min}(m,N-2),0]\,.$

The integers $N$ and $M$ are explained in App. A.2 and A.4. The asymptotic
formula (36) depends on the number of massless quarks (entering the
expressions for the QCD $\beta$-function coefficients $\beta_{i}$ and the
coefficients $S^{(N,M)}_{k}$ and $g^{N}_{\ell}$) and was derived in Ref. [1] [
see Eq. (4.19) in that reference ]. It has the property that it reproduces the
exact coefficients $a_{n}$ up to ${\cal O}(\alpha_{s}^{{\rm min}[M,N-1]})$. In
contrast to the expression given in Ref. [1], Eq. (36) displays concrete
truncation prescriptions for the sums as a function of the loop order
variables.

###  A.6 Renormalon Sum Rule

To extract the pole mass renormalon normalization $N_{1/2}$ from a Core, the
following formula is used:

$N_{1/2}=\frac{\beta_{0}}{2\pi}\sum_{k=0}^{M}\frac{S_{k}^{(N,M)}}{(1+\hat{b}_{1}^{N})_{k}}\;,$
(38)

The formula was derived in Ref. [1], and Eq. (38) is adapted to account for
the optional parameters used in REvolver. In analogy to Eq. (36), $N_{1/2}$ is
always evaluated in the flavor number scheme where all massive quarks are
integrated out. The expression for the coefficients $S^{(N,M)}_{k}$ is given
in Eq. (32) while $b_{1}^{N}$ has been defined in Eq. (16).

## References

  * [1] A. H. Hoang, A. Jain, C. Lepenik, V. Mateu, M. Preisser, I. Scimemi, I. W. Stewart, The MSR mass and the $\mathcal{O}\left({\Lambda}_{\mathrm{QCD}}\right)$ renormalon sum rule, JHEP 04 (2018) 003.  arXiv:1704.01580, doi:10.1007/JHEP04(2018)003.
  * [2] A. H. Hoang, C. Lepenik, M. Preisser, On the Light Massive Flavor Dependence of the Large Order Asymptotic Behavior and the Ambiguity of the Pole Mass, JHEP 09 (2017) 099.  arXiv:1706.08526, doi:10.1007/JHEP09(2017)099.
  * [3] M. Beneke, P. Marquard, P. Nason, M. Steinhauser, On the ultimate uncertainty of the top quark pole mass, Phys. Lett. B 775 (2017) 63–70.  arXiv:1605.03609, doi:10.1016/j.physletb.2017.10.054.
  * [4] K. Chetyrkin, J. H. Kuhn, M. Steinhauser, RunDec: A Mathematica package for running and decoupling of the strong coupling and quark masses, Comput. Phys. Commun. 133 (2000) 43–65.  arXiv:hep-ph/0004189, doi:10.1016/S0010-4655(00)00155-7.
  * [5] F. Herren, M. Steinhauser, Version 3 of RunDec and CRunDec, Comput. Phys. Commun. 224 (2018) 333–345.  arXiv:1703.03751, doi:10.1016/j.cpc.2017.11.014.
  * [6] A. H. Hoang, A. Jain, I. Scimemi, I. W. Stewart, Infrared Renormalization Group Flow for Heavy Quark Masses, Phys. Rev. Lett. 101 (2008) 151602.  arXiv:0803.4214, doi:10.1103/PhysRevLett.101.151602.
  * [7] Wolfram Mathematica, https://www.wolfram.com/mathematica/, accessed: 2020-05.
  * [8] G. Rossum, Python reference manual, Tech. rep., Amsterdam, The Netherlands, The Netherlands (1995).
  * [9] SWIG, http://swig.org/, accessed: 2020-05.
  * [10] T. Kluyver, B. Ragan-Kelley, F. Pérez, B. Granger, M. Bussonnier, J. Frederic, K. Kelley, J. Hamrick, J. Grout, S. Corlay, P. Ivanov, D. Avila, S. Abdalla, C. Willing, J. development team, Jupyter notebooks - a publishing format for reproducible computational workflows, in: F. Loizides, B. Scmidt (Eds.), Positioning and Power in Academic Publishing: Players, Agents and Agendas, IOS Press, 2016, pp. 87–90.
URL https://eprints.soton.ac.uk/403913/

  * [11] Quadpack++, https://github.com/drjerry/quadpackpp, accessed: 2021-01.
  * [12] P. Marquard, A. V. Smirnov, V. A. Smirnov, M. Steinhauser, D. Wellmann, $\overline{\rm MS}$-on-shell quark mass relation up to four loops in QCD and a general SU$(N)$ gauge group, Phys. Rev. D94 (7) (2016) 074025.  arXiv:1606.06754, doi:10.1103/PhysRevD.94.074025.
  * [13] E. G. Floratos, S. Narison, E. de Rafael, Spectral Function Sum Rules in Quantum Chromodynamics. 1. Charged Currents Sector, Nucl. Phys. B 155 (1979) 115–149.  doi:10.1016/0550-3213(79)90359-6.
  * [14] A. H. Hoang, Z. Ligeti, A. V. Manohar, B decay and the Upsilon mass, Phys. Rev. Lett. 82 (1999) 277–280.  arXiv:hep-ph/9809423, doi:10.1103/PhysRevLett.82.277.
  * [15] A. H. Hoang, Z. Ligeti, A. V. Manohar, B decays in the upsilon expansion, Phys. Rev. D59 (1999) 074017.  arXiv:hep-ph/9811239, doi:10.1103/PhysRevD.59.074017.
  * [16] A. H. Hoang, 1S and MS-bar bottom quark masses from Upsilon sum rules, Phys. Rev. D61 (2000) 034005.  arXiv:hep-ph/9905550, doi:10.1103/PhysRevD.61.034005.
  * [17] A. Czarnecki, K. Melnikov, N. Uraltsev, NonAbelian dipole radiation and the heavy quark expansion, Phys. Rev. Lett. 80 (1998) 3189–3192.  arXiv:hep-ph/9708372, doi:10.1103/PhysRevLett.80.3189.
  * [18] M. Beneke, A quark mass definition adequate for threshold problems, Phys. Lett. B434 (1998) 115–125.  arXiv:hep-ph/9804241, doi:10.1016/S0370-2693(98)00741-2.
  * [19] M. Beneke, Y. Kiyo, K. Schuller, Third-order Coulomb corrections to the S-wave Green function, energy levels and wave functions at the origin, Nucl. Phys. B714 (2005) 67–90.  arXiv:hep-ph/0501289, doi:10.1016/j.nuclphysb.2005.02.028.
  * [20] A. Pineda, Determination of the bottom quark mass from the $\Upsilon(1S)$ system, JHEP 06 (2001) 022.  arXiv:hep-ph/0105008, doi:10.1088/1126-6708/2001/06/022.
  * [21] M. Fael, K. Schönwald, M. Steinhauser, Kinetic Heavy Quark Mass to Three Loops, Phys. Rev. Lett. 125 (5) (2020) 052003.  arXiv:2005.06487, doi:10.1103/PhysRevLett.125.052003.
  * [22] M. Fael, K. Schönwald, M. Steinhauser, Relation between the $\overline{\mathrm{MS}}$ and the kinetic mass of heavy quarks, Phys. Rev. D 103 (1) (2021) 014005.  arXiv:2011.11655, doi:10.1103/PhysRevD.103.014005.
  * [23] P. Zyla, et al., Review of Particle Physics, PTEP 2020 (8) (2020) 083C01.  doi:10.1093/ptep/ptaa104.
  * [24] V. Khachatryan, et al., Measurement and QCD analysis of double-differential inclusive jet cross sections in pp collisions at $\sqrt{s}=8$ TeV and cross section ratios to 2.76 and 7 TeV, JHEP 03 (2017) 156.  arXiv:1609.05331, doi:10.1007/JHEP03(2017)156.
  * [25] M. Beneke, M. Jamin, $\alpha_{s}$ and the $\tau$ hadronic width: fixed-order, contour-improved and higher-order perturbation theory, JHEP 09 (2008) 044.  arXiv:0806.3156, doi:10.1088/1126-6708/2008/09/044.
  * [26] D. Boito, M. Jamin, R. Miravitllas, Scheme Variations of the QCD Coupling and Hadronic $\tau$ Decays, Phys. Rev. Lett. 117 (15) (2016) 152001.  arXiv:1606.06175, doi:10.1103/PhysRevLett.117.152001.
  * [27] M. Aaboud, et al., Determination of the strong coupling constant $\alpha_{\mathrm{s}}$ from transverse energy–energy correlations in multijet events at $\sqrt{s}=8~{}\hbox{TeV}$ using the ATLAS detector, Eur. Phys. J. C 77 (12) (2017) 872.  arXiv:1707.02562, doi:10.1140/epjc/s10052-017-5442-0.
  * [28] C. Duhr, F. Dulat, V. Hirschi, B. Mistlberger, Higgs production in bottom quark fusion: matching the 4- and 5-flavour schemes to third order in the strong coupling, JHEP 08 (08) (2020) 017.  arXiv:2004.04752, doi:10.1007/JHEP08(2020)017.
  * [29] P. Marquard, A. V. Smirnov, V. A. Smirnov, M. Steinhauser, Quark Mass Relations to Four-Loop Order in Perturbative QCD, Phys. Rev. Lett. 114 (14) (2015) 142002.  arXiv:1502.01030, doi:10.1103/PhysRevLett.114.142002.
  * [30] C. Peset, A. Pineda, J. Segovia, The charm/bottom quark mass from heavy quarkonium at N3LO, JHEP 09 (2018) 167.  arXiv:1806.05197, doi:10.1007/JHEP09(2018)167.
  * [31] I. I. Y. Bigi, M. A. Shifman, N. G. Uraltsev, A. I. Vainshtein, The Pole mass of the heavy quark. Perturbation theory and beyond, Phys. Rev. D50 (1994) 2234–2246.  arXiv:hep-ph/9402360, doi:10.1103/PhysRevD.50.2234.
  * [32] M. Beneke, V. M. Braun, Heavy quark effective theory beyond perturbation theory: Renormalons, the pole mass and the residual mass term, Nucl. Phys. B426 (1994) 301–343.  arXiv:hep-ph/9402364, doi:10.1016/0550-3213(94)90314-X.
  * [33] D. Boito, M. Golterman, A. Keshavarzi, K. Maltman, D. Nomura, S. Peris, T. Teubner, Strong coupling from $e^{+}e^{-}\to$ hadrons below charm, Phys. Rev. D 98 (7) (2018) 074030.  arXiv:1805.08176, doi:10.1103/PhysRevD.98.074030.
  * [34] P. A. Baikov, K. G. Chetyrkin, J. H. Kuhn, Order $\alpha_{s}^{4}$ QCD Corrections to Z and tau Decays, Phys. Rev. Lett. 101 (2008) 012002.  arXiv:0801.1821, doi:10.1103/PhysRevLett.101.012002.
  * [35] M. Bruno, M. Dalla Brida, P. Fritzsch, T. Korzec, A. Ramos, S. Schaefer, H. Simma, S. Sint, R. Sommer, The $\Lambda$-parameter in 3-flavour QCD and $\alpha_{s}(m_{Z})$ by the ALPHA collaboration, PoS LATTICE2016 (2016) 197.  arXiv:1701.03075, doi:10.22323/1.256.0197.
  * [36] M. Butenschoen, B. Dehnadi, A. H. Hoang, V. Mateu, M. Preisser, I. W. Stewart, Top Quark Mass Calibration for Monte Carlo Event Generators, Phys. Rev. Lett. 117 (23) (2016) 232001.  arXiv:1608.01318, doi:10.1103/PhysRevLett.117.232001.
  * [37] A. H. Hoang, A. V. Manohar, Charm effects in the MS-bar bottom quark mass from Upsilon mesons, Phys. Lett. B483 (2000) 94–98.  arXiv:hep-ph/9911461, doi:10.1016/S0370-2693(00)00595-5.
  * [38] S. Bekavac, A. Grozin, D. Seidel, M. Steinhauser, Light quark mass effects in the on-shell renormalization constants, JHEP 10 (2007) 006.  arXiv:0708.1729, doi:10.1088/1126-6708/2007/10/006.
  * [39] O. V. Tarasov, A. A. Vladimirov, A. Y. Zharkov, The Gell-Mann-Low Function of QCD in the Three Loop Approximation, Phys. Lett. B93 (1980) 429–432.
  * [40] S. A. Larin, J. A. M. Vermaseren, The three-loop QCD $\beta$ function and anomalous dimensions, Phys. Lett. B303 (1993) 334–336.  arXiv:hep-ph/9302208.
  * [41] T. van Ritbergen, J. A. M. Vermaseren, S. A. Larin, The four-loop beta function in quantum chromodynamics, Phys. Lett. B400 (1997) 379–384.  arXiv:hep-ph/9701390, doi:10.1016/S0370-2693(97)00370-5.
  * [42] M. Czakon, The four-loop QCD beta-function and anomalous dimensions, Nucl. Phys. B710 (2005) 485–498.  arXiv:hep-ph/0411261, doi:10.1016/j.nuclphysb.2005.01.012.
  * [43] P. A. Baikov, K. G. Chetyrkin, J. H. Kühn, Five-Loop Running of the QCD coupling constant, Phys. Rev. Lett. 118 (8) (2017) 082002.  arXiv:1606.08659, doi:10.1103/PhysRevLett.118.082002.
  * [44] K. G. Chetyrkin, G. Falcioni, F. Herzog, J. A. M. Vermaseren, Five-loop renormalisation of QCD in covariant gauges, JHEP 10 (2017) 179, [Addendum: JHEP 12, 006 (2017)].  arXiv:1709.08541, doi:10.1007/JHEP10(2017)179.
  * [45] T. Luthe, A. Maier, P. Marquard, Y. Schröder, The five-loop Beta function for a general gauge group and anomalous dimensions beyond Feynman gauge, JHEP 10 (2017) 166.  arXiv:1709.07718, doi:10.1007/JHEP10(2017)166.
  * [46] F. Herzog, B. Ruijl, T. Ueda, J. A. M. Vermaseren, A. Vogt, The five-loop beta function of Yang-Mills theory with fermions, JHEP 02 (2017) 090.  arXiv:1701.01404, doi:10.1007/JHEP02(2017)090.
  * [47] K. G. Chetyrkin, B. A. Kniehl, M. Steinhauser, Decoupling relations to $\mathcal{O}(\alpha_{s}^{3})$ and their connection to low-energy theorems, Nucl. Phys. B510 (1998) 61–87.  arXiv:hep-ph/9708255, doi:10.1016/S0550-3213(98)81004-3,10.1016/S0550-3213(97)00649-4.
  * [48] K. G. Chetyrkin, J. H. Kühn, C. Sturm, QCD decoupling at four loops, Nucl. Phys. B744 (2006) 121–135.  arXiv:hep-ph/0512060, doi:10.1016/j.nuclphysb.2006.03.020.
  * [49] Y. Schröder, M. Steinhauser, Four-loop decoupling relations for the strong coupling, JHEP 01 (2006) 051.  arXiv:hep-ph/0512058, doi:10.1088/1126-6708/2006/01/051.
  * [50] B. A. Kniehl, A. V. Kotikov, A. I. Onishchenko, O. L. Veretin, Strong-coupling constant with flavor thresholds at five loops in the anti-MS scheme, Phys. Rev. Lett. 97 (2006) 042001.  arXiv:hep-ph/0607202, doi:10.1103/PhysRevLett.97.042001.
  * [51] R. Tarrach, The Pole Mass in Perturbative QCD, Nucl. Phys. B183 (1981) 384–396.  doi:10.1016/0550-3213(81)90140-1.
  * [52] N. Gray, D. J. Broadhurst, W. Grafe, K. Schilcher, Three Loop Relation of Quark $\overline{\rm MS}$ and Pole Masses, Z. Phys. C48 (1990) 673–680.  doi:10.1007/BF01614703.
  * [53] K. Melnikov, T. v. Ritbergen, The Three loop relation between the ${\overline{\rm MS}}$ and the pole quark masses, Phys. Lett. B482 (2000) 99–108.  arXiv:hep-ph/9912391, doi:10.1016/S0370-2693(00)00507-4.
  * [54] M. Fael, K. Schönwald, M. Steinhauser, Exact results for ${Z}_{m}^{\mathrm{OS}}$ and ${Z}_{2}^{\mathrm{OS}}$ with two mass scales and up to three loops, JHEP 10 (2020) 087.  arXiv:2008.01102, doi:10.1007/JHEP10(2020)087.
  * [55] J. A. M. Vermaseren, S. A. Larin, T. van Ritbergen, The 4-loop quark mass anomalous dimension and the invariant quark mass, Phys. Lett. B405 (1997) 327–333.  arXiv:hep-ph/9703284, doi:10.1016/S0370-2693(97)00660-6.
  * [56] K. G. Chetyrkin, Quark mass anomalous dimension to $O(\alpha_{s}^{4})$, Phys. Lett. B404 (1997) 161–165.  arXiv:hep-ph/9703278, doi:10.1016/S0370-2693(97)00535-2.
  * [57] P. A. Baikov, K. G. Chetyrkin, J. H. Kühn, Quark Mass and Field Anomalous Dimensions to ${\cal O}(\alpha_{s}^{5})$, JHEP 10 (2014) 076.  arXiv:1402.6611, doi:10.1007/JHEP10(2014)076.
  * [58] T. Luthe, A. Maier, P. Marquard, Y. Schröder, Five-loop quark mass and field anomalous dimensions for a general gauge group, JHEP 01 (2017) 081.  arXiv:1612.05512, doi:10.1007/JHEP01(2017)081.
  * [59] T. Liu, M. Steinhauser, Decoupling of heavy quarks at four loops and effective Higgs-fermion coupling, Phys. Lett. B 746 (2015) 330–334.  arXiv:1502.04719, doi:10.1016/j.physletb.2015.05.023.
  * [60] V. Mateu, P. G. Ortega, Bottom and Charm Mass determinations from global fits to $Q\overline{Q}$ bound states at N3LO, JHEP 01 (2018) 122.  arXiv:1711.05755, doi:10.1007/JHEP01(2018)122.
  * [61] Y. Kiyo, Y. Sumino, Full Formula for Heavy Quarkonium Energy Levels at Next-to-next-to-next-to-leading Order, Nucl. Phys. B889 (2014) 156–191.  arXiv:1408.5590, doi:10.1016/j.nuclphysb.2014.10.010.
  * [62] N. Brambilla, Y. Sumino, A. Vairo, Quarkonium spectroscopy and perturbative QCD: Massive quark loop effects, Phys. Rev. D65 (2002) 034001.  arXiv:hep-ph/0108084, doi:10.1103/PhysRevD.65.034001.
  * [63] A. A. Penin, M. Steinhauser, Heavy quarkonium spectrum at $\mathcal{O}(\alpha_{s}^{5}m_{q})$ and bottom / top quark mass determination, Phys. Lett. B538 (2002) 335–345.  arXiv:hep-ph/0204290, doi:10.1016/S0370-2693(02)02040-3.
  * [64] Y. Kiyo, Y. Sumino, Perturbative heavy quarkonium spectrum at next-to-next-to-next-to-leading order, Phys. Lett. B730 (2014) 76–80.  arXiv:1309.6571, doi:10.1016/j.physletb.2014.01.030.
  * [65] W. Fischler, Quark - anti-Quark Potential in QCD, Nucl. Phys. B129 (1977) 157–174.  doi:10.1016/0550-3213(77)90026-8.
  * [66] A. Billoire, How Heavy Must Be Quarks in Order to Build Coulombic $q\bar{q}$ Bound States, Phys. Lett. B92 (1980) 343–347.  doi:10.1016/0370-2693(80)90279-8.
  * [67] Y. Schroder, The Static potential in QCD to two loops, Phys. Lett. B447 (1999) 321–326.  arXiv:hep-ph/9812205, doi:10.1016/S0370-2693(99)00010-6.
  * [68] A. V. Smirnov, V. A. Smirnov, M. Steinhauser, Fermionic contributions to the three-loop static potential, Phys. Lett. B668 (2008) 293–298.  arXiv:0809.1927, doi:10.1016/j.physletb.2008.08.070.
  * [69] A. V. Smirnov, V. A. Smirnov, M. Steinhauser, Three-loop static potential, Phys. Rev. Lett. 104 (2010) 112002.  arXiv:0911.4742, doi:10.1103/PhysRevLett.104.112002.
  * [70] C. Anzai, Y. Kiyo, Y. Sumino, Static QCD potential at three-loop order, Phys. Rev. Lett. 104 (2010) 112003.  arXiv:0911.4335, doi:10.1103/PhysRevLett.104.112003.
  * [71] R. N. Lee, A. V. Smirnov, V. A. Smirnov, M. Steinhauser, Analytic three-loop static potential, Phys. Rev. D94 (5) (2016) 054029.  arXiv:1608.02603, doi:10.1103/PhysRevD.94.054029.
  * [72] D. Eiras, J. Soto, Light fermion finite mass effects in non-relativistic bound states, Phys. Lett. B491 (2000) 101–110.  arXiv:hep-ph/0005066, doi:10.1016/S0370-2693(00)01004-2.
  * [73] A. H. Hoang, Bottom quark mass from Upsilon mesons: Charm mass effects (2000).  arXiv:hep-ph/0008102.
  * [74] M. Beneke, A. Maier, J. Piclum, T. Rauh, The bottom-quark mass from non-relativistic sum rules at NNNLO, Nucl. Phys. B891 (2015) 42–72.  arXiv:1411.3132, doi:10.1016/j.nuclphysb.2014.12.001.