#  REvolver
Automated running and matching of couplings and masses in QCD

André H. Hoang<EMAIL_ADDRESS>Christopher Lepenik
<EMAIL_ADDRESS>Vicent Mateu<EMAIL_ADDRESS>University of
Vienna, Faculty of Physics, Boltzmanngasse 5, A-1090 Wien, Austria  Erwin
Schrödinger International Institute for Mathematical Physics,
University of Vienna, Boltzmanngasse 9, A-1090 Wien, Austria  Departamento de
Física Fundamental e IUFFyM,
Universidad de Salamanca, E-37008 Salamanca, Spain  Instituto de Física
Teórica UAM-CSIC, E-28049 Madrid, Spain

###### Abstract

In this article we present REvolver, a C++ library for renormalization group
evolution and automatic flavor matching of the QCD coupling and quark masses,
as well as precise conversion between various quark mass renormalization
schemes. The library systematically accounts for the renormalization group
evolution of low-scale short-distance masses which depend linearly on the
renormalization scale and sums logarithmic terms of high and low scales that
are missed by the common logarithmic renormalization scale evolution. The
library can also be accessed through Mathematica and Python interfaces and
provides renormalization group evolution for complex renormalization scales as
well.

###### keywords:

QCD, renormalization group, heavy quarks

††journal: Computer Physics Communications

UWThPh-2020-18, IFT-UAM/CSIC-21-5

PROGRAM SUMMARY

Program Title: REvolver
Licensing provisions: GPLv3 or later
Programming language: C++, Python, Wolfram Language
Program obtainable from: https://gitlab.com/REvolver-hep/REvolver
Operating system: Linux, MacOS, partially Windows
Required RAM: insignificant for a limited number of instances of the Core
class
Number of processors used: one
Running time: fractions of seconds for single commands
Supplementary material: this article, demo programs, doxygen documentation

Nature of problem:
The strong coupling and the quark masses are fundamental parameters of QCD
that are scheme and renormalization-scale dependent. The choice of scheme
depends on the active number of flavors and the range of scales, and is
dictated by the requirements to minimize the size of corrections and to sum
large logarithmic corrections to all orders. For the strong coupling and the
quark masses at high scales, the $\overline{\mbox{MS}}$ scheme with
logarithmic scale dependence is used. For quark masses at low scales, short-
distance mass schemes with linear scale-dependence are used. The REvolver
library provides conversions for the strong coupling and the most common quark
mass schemes, with renormalization scale evolution implemented such that all
types of large logarithmic terms are summed to all orders, accounting for
flavor threshold effects and state-or-the-art correction terms. The pole mass,
which is not a short-distance mass and contains a sizable renormalon
ambiguity, is treated as a derived quantity.

Solution method:
Renormalization group equations are solved for complex-valued scales to
machine precision based on fast-converging iterative algorithms and analytic
all-order expressions. Matching relations for the strong coupling at flavor
thresholds are computed in a way that gives equal results for upward and
downward evolution. Core objects allow to define an arbitrary number of
physical scenarios for strong coupling values and quark mass spectra, where
options for precision and matching scales can be set freely, and values for
quark masses in all common schemes including the pole mass can be extracted.
All REvolver routines are implemented entirely in C++ and can be accessed
through Mathematica and Python interfaces.

###### Contents

  1. 1 Introduction
  2. 2 Terminology
  3. 3 The MSR mass and R-evolution
  4. 4 Setup
    1. 4.1 Installation
    2. 4.2 General Usage
  5. 5 Core Structure
  6. 6 Implemented Functions
    1. 6.1 Constructing a Core and Accessing Scenario Parameters
      1. 6.1.1 C++ / Python only: RunPar and RunParV
      2. 6.1.2 Constructing Core objects with masses
      3. 6.1.3 Constructing Core objects without massive quarks
      4. 6.1.4 Mathematica only: listing and deleting Cores
      5. 6.1.5 Accessing Core parameters
    2. 6.2 Extraction of Masses and Couplings
      1. 6.2.1 Running masses and strong coupling
      2. 6.2.2 Other short-distance masses
      3. 6.2.3 Pole mass
      4. 6.2.4 Norm of the pole mass renormalon ambiguity
      5. 6.2.5 Extracting $\Lambda_{\mathrm{QCD}}$
    3. 6.3 Adding Masses to an Existing Core
  7. 7 Applications and Pedagogical Examples
    1. 7.1 Cores without Massive Quarks
      1. 7.1.1 Strong coupling evolution
      2. 7.1.2 Complex renormalization scales
    2. 7.2 Cores with one Massive Quark
      1. 7.2.1 Strong coupling with a flavor threshold
      2. 7.2.2 Asymptotic pole and low-scale MSR mass
      3. 7.2.3 Bottom and charm quark short-distance masses
      4. 7.2.4 Top quark mass at low scales
      5. 7.2.5 Running masses for complex renormalization scales
      6. 7.2.6 Top quark pole masses
    3. 7.3 Cores with Multiple Massive Quarks
      1. 7.3.1 Strong coupling and $\Lambda_{\rm QCD}$
      2. 7.3.2 Top quark running mass at low scales
      3. 7.3.3 Charm mass effects for the bottom quark mass
      4. 7.3.4 Top quark pole masses
  8. 8 Theoretical input
  9. 9 Summary
  10. A Algorithms to Solve the Differential Equations
    1. A.1 Core Creation
    2. A.2 Strong Coupling
    3. A.3 $\overline{\rm MS}$ and RGI masses
    4. A.4 MSR mass
    5. A.5 Asymptotic Pole Mass Coefficients
    6. A.6 Renormalon Sum Rule

##  1 Introduction

Quark masses are fundamental parameters of quantum chromodynamics(QCD) and
their precise determination in adequate schemes and at appropriate
renormalization scales is of high interest for theoretical as well as
experimental studies of many processes. These can be governed by energy scales
ranging from a few GeV (e.g. for hadronic states) up to several hundred GeV
and even TeV scales (e.g. for particle collisions that take place at the Large
Hadron Collider). One needs to employ the renormalization group evolution
equations to reliably relate the values of quark masses defined at such widely
different energy scales. An interesting situation arises if the dynamical
scale governing the quark mass dependence of an observable is much smaller
than the quark mass itself. In this case, the common running ${\overline{\rm
MS}}$ mass scheme, which obeys a renormalization group equation with
logarithmic scale dependence, cannot be employed for high-precision
applications, because it is only meaningful for scales of the order or larger
than the mass. Rather, so-called low-scale short-distance masses must be used,
which obey renormalization group equations with linear scale dependence. The
numerical impact of the renormalization group evolution is particularly
important for the top quark mass where, due to its large value, significant
scale hierarchies can arise.

Here we present REvolver, a C++ library with routines that provide
renormalization-group resummed conversions between quark mass schemes defined
at different renormalization scales, including scales much lower than the
mass, where low-scale short-distance masses are employed, as well as above the
mass value, where the ${\overline{\rm MS}}$ mass is used. The routines are
based on the creation of so-called Core objects, each of which representing a
certain physical scenario for the heavy quark masses (charm, bottom and top
quarks, as well as hypothetical heavier flavors), the number of massless
quarks and the strong coupling $\alpha_{s}$. In a single session, an (in
principle) arbitrary number of Core objects can be created and managed. Each
Core object can then provide values for the quark masses in the most popular
low-scale short-distance schemes as well as for the ${\overline{\rm MS}}$ mass
and the strong coupling at any (real or complex-valued) renormalization scale
and in any flavor number scheme, consistently accounting for flavor threshold
effects and the resummation of large logarithms of all kinds. The basis of the
quark mass evolution equations for scales below the respective mass is the
renormalization group equation of the natural MSR mass (here simply called the
MSR mass), which was provided in Refs. [1] together with a full treatment of
flavor matching corrections when the evolution crosses the thresholds related
to lighter massive quarks [2]. Furthermore, for each Core object, options can
be set to specify the perturbative precision in the flavor matching and the
renormalization group evolution. Using all available theoretical input in the
literature, see Sec. 8 for a detailed listing, it is possible to relate the
${\overline{\rm MS}}$, MSR and most other low-scale short-distance quark
masses with a theoretical precision of $10$ to $20$ MeV (neglecting any
parametric uncertainties). Quark mass values in the pole scheme cannot be
defined at the same level of precision due to the pole mass renormalon
ambiguity which decreases the accuracy by an order of magnitude [3, 2].
REvolver offers the possibility to set up Core objects using pole masses as an
input or to extract pole mass values from a Core object, but it treats the
pole mass as a derived quantity where the user has to specify the way in which
the pole mass value is defined. Furthermore, REvolver provides various options
to account for the asymptotic higher order corrections of the pole mass and
the pole mass renormalon ambiguity. All REvolver routines are implemented
entirely in C++ and can be accessed through Mathematica and Python interfaces.

There is an existing C++ library accompanied by a Mathematica package called
CRunDec and RunDec [4, 5], respectively, which already provide many
functionalities included in the REvolver library. We have cross checked in
detail that any theoretical (perturbative) input implemented in CRunDec and
RunDec agrees with the corresponding one employed for REvolver. We have
furthermore checked that the numerical output of the routines provided in
CRunDec/RunDec is in agreement with the equivalent routines of REvolver. The
REvolver library, however, exceeds CRunDec/RunDec

  * (i)

by providing the Core concept that allows to automatically create, extend and
manage an arbitrary number of scenarios for strong coupling values, mass
spectra and theory settings, and to extract quark masses and the QCD coupling
in all flavor number schemes and at all scales,

  * (ii)

by accounting for the renormalization group resummation of large logarithms
and lighter massive quark flavor thresholds when dealing with quark masses at
renormalization scales smaller than the quark mass, as well as for low-scale
short-distance masses,

  * (iii)

by giving access to machine-precision numerical routines that provide quasi-
exact solutions of the renormalization group equations for the running masses
and the strong coupling at complex scales, and

  * (iv)

by providing routines to determine the asymptotic series for the pole mass to
an arbitrary order that allow to extract different pole mass definitions and
to quantify the pole mass renormalon ambiguity with various methods.

The emphasis of all REvolver functionalities is to provide integrated and
easy-to-use routines, while maintaining the possibility to deviate from
default settings and specify all available options, useful for high-precision
phenomenological and conceptual QCD studies aiming for uncertainties at the
level of $10$ to $20$ MeV for short-distance masses.

This article is organized as follows: In Sec. 2 essential terminology used for
the description of the REvolver package is explained. In Sec. 3 we succinctly
review the MSR mass and the R-evolution concepts [1, 6] which are essential
for the resummation of the logarithms mentioned above in bullet point (ii).
Section 4 provides general information concerning the REvolver installation
and setting up the C++, Mathematica and Python interfaces. The philosophy of
the Core concept is explained in Sec. 5, and Sec. 6 provides a structured
introduction to all available REvolver routines. In Sec. 7 a sizable number of
pedagogical examples for applications of REvolver routines are provided,
partly using the routine’s default settings, partly using alternative optional
parameter setting, to demonstrate the versatility of REvolver for important
phenomenological applications in the literature. It is recommended that the
user consults the examples shown in this section, which are also collected in
a Mathematica notebook, a Jupyter Notebook using the Python interface, and a
C++ source file provided with the REvolver package. In Sec. 8 important
references are provided which were used as the source for the higher order
corrections implemented for the strong coupling and the mass schemes supported
by REvolver. Here, also a detailed citation recommendation for these higher
order corrections is provided. Section 9 contains a summary. Finally, some
details concerning the algorithms used for Core creation and the quasi-exact
solution of renormalization group equations are given in Appendix A. In
addition, a number of essential formulae for implementation-dependent
quantities are provided which cannot be found in the literature in the form
used in REvolver.

##  2 Terminology

This article employs a particular terminology when referring to
renormalization-scale dependent mass schemes and the pole mass:

  * 1.

Running quark mass in the $n_{f}$-flavor scheme $m^{(n_{f})}_{q}(\mu)$: Refers
to the $\overline{\mathrm{MS}}$ mass if the flavor number $n_{f}$ includes
this massive quark, and the MSR mass otherwise. For example, the running top
quark mass at the scale $\mu$ in the 6-flavor scheme $m_{t}^{(6)}(\mu)$ refers
to the $\overline{\mathrm{MS}}$ mass $\overline{m}_{t}^{(6)}(\mu)$, and the
running top quark mass at the scale $\mu$ in the 4-flavor scheme
$m_{t}^{(4)}(\mu)$ refers to the MSR mass $m_{t}^{{\rm MSR},(4)}(\mu)$.

  * 2.

Standard running mass $\overline{m}_{q}$: Refers to the
$\overline{\mathrm{MS}}$ mass in the flavor number scheme where all lighter
quarks along with this quark are treated dynamically, evaluated at the scale
of this mass. For example, the standard running top mass is the $6$-flavor
$\overline{\mathrm{MS}}$ quark mass evaluated at the scale of this top mass:
$\overline{m}_{t}\equiv\overline{m}_{t}^{(6)}\\!(\overline{m}_{t}^{(6)})$. The
mass dependence of all flavor-threshold corrections is expressed in terms of
the standard running mass.

  * 3.

Asymptotic pole mass: Refers to the pole mass value obtained from the running
mass defined by summing the perturbative series to the order of the minimal
correction.

  * 4.

Order-dependent pole mass: Refers to the pole mass value obtained from the
running mass by truncating the perturbative series at a specified order.

##  3 The MSR mass and R-evolution

The natural MSR mass of a massive quark defined in Ref. [1] (and called just
the MSR mass here) plays a central role in REvolver and is a renormalization
scale and flavor-number-dependent low-scale short-distance mass. It is derived
from the $\overline{\mathrm{MS}}$ mass and treated as the natural extension of
the $\overline{\mathrm{MS}}$ mass for renormalization scales below the mass of
the quark. This combination of the scale-dependent $\overline{\mathrm{MS}}$
and MSR masses extends the well known concept of flavor-number dependent
renormalization group evolution and flavor threshold matching for scales above
the quark mass (where the $\overline{\mathrm{MS}}$ mass scheme is appropriate)
to lower scales. In contrast to usual logarithmic renormalization scale
evolution (as known from the $\overline{\mathrm{MS}}$ masses or the strong
coupling), the MSR mass renormalization group evolution is linear. This is
consequence of the linear dynamical scaling that arises when the off-shell
massive quark quantum fluctuations are integrated out in the nonrelativistic
limit. Together with the flavor number dependent strong
$\overline{\mathrm{MS}}$ coupling, the $\overline{\mathrm{MS}}$ and MSR masses
form the basis of the core concept of REvolver and allow to also resum large
logarithms involving low-scale short-distance mass schemes other than the MSR
mass. This functionality is used by default in the REvolver routines (but can
also be switched off by the user on demand). The MSR mass [1, 2, 6] has
already been used in a number of applications, but we still find it warranted
to briefly review its main concepts in this section. For simplicity, we
consider the case in which all $n_{\ell}$ quarks lighter than $q$ are
massless. The reader is referred to Ref. [2] for the case with massive lighter
quarks.

To define the MSR mass, one starts with the relation between the standard
running mass and pole mass,

$m_{q}^{\rm
pole}-\overline{m}_{q}=\overline{m}_{q}\,\sum_{n=1}^{\infty}\,a_{n}^{\overline{\rm
MS}}(n_{\ell},n_{h})\Biggl{[}\frac{\alpha_{s}^{(n_{\ell}+n_{h})}(\overline{m}_{q})}{4\pi}\Biggr{]}^{\\!n}\,,$
(1)

where $n_{f}=n_{\ell}+n_{h}$ is the number of active flavors, with $n_{\ell}$
being the number of massless quarks and $n_{h}=1$ referring to the quark $q$.
Since the MSR mass is employed for renormalization scales below $m_{q}$, one
integrates out the virtual heavy quark loops by setting $n_{h}=0$. This allows
to define a renormalization scale $R$ smaller than $m_{q}$ and the corrections
to the pole mass having a linear dependence on $R$ to implement a consistent
nonrelativistic scaling behavior. The MSR mass is thus defined by furthermore
setting $\overline{m}_{q}\to R$:

$\displaystyle m_{q}^{\rm pole}-m_{q}^{{\rm
MSR},(n_{\ell})}(R)=R\sum_{n=1}^{\infty}a_{n}^{\overline{\mathrm{MS}}}(n_{\ell},0)\\!\Biggl{[}\frac{\alpha_{s}^{(n_{\ell})}(R)}{4\pi}\Biggr{]}^{\\!n}\,.$
(2)

In contrast to the $\overline{\mathrm{MS}}$ mass, which has only logarithmic
dependence on the scale $\mu$, the MSR mass has an additional linear
dependence on $R$. The $\overline{\mathrm{MS}}$ and MSR masses can be related
perturbatively and unambiguously through Eqs. (1) and (2) because the pole
mass in both equalities is identical.111This is consistent since the two
series on the RHS of Eqs. (1) and (2) have the same leading linear and mass-
independent renormalon ambiguity. The resulting perturbative series for the
difference of two MSR masses at different renormalization scales $R_{1}$ and
$R_{2}$ is renormalon free, as long as it is expressed in powers of the strong
coupling at the same renormalization scale. As a result, for disparate values
of $R_{1}$ and $R_{2}$ large logarithms will appear. These logarithms can be
consistently summed up with the renormalization group equation

$-\\!\frac{\rm d}{{\rm d}R}m_{q}^{{\rm
MSR},(n_{\ell})}(R)=\gamma_{R}^{M}[\alpha_{s}^{(n_{\ell})}(R)]=\sum_{n=0}^{\infty}\gamma_{n}^{R}(n_{\ell})\\!\Biggl{[}\frac{\alpha_{s}^{(n_{\ell})}(R)}{4\pi}\Biggr{]}^{\\!n+1},$
(3)

In contrast to the logarithmic renormalization group equations for the
$\overline{\mathrm{MS}}$ mass and the strong coupling, it shows a linear power
scaling and has therefore been dubbed as the R-evolution equation. The
anomalous dimension coefficients $\gamma_{n}^{R}(n_{\ell})$ can be calculated
from the relation [ see Eq. (6) for the definition of the QCD $\beta$-function
coefficients $\beta_{i}$ ]

$\displaystyle\gamma_{n}^{R}(n_{\ell})=a^{\overline{\mathrm{MS}}}_{n+1}(n_{\ell},0)-2\sum_{j=0}^{n-1}(n-j)\,\beta_{j}\,a^{\overline{\mathrm{MS}}}_{n-j}(n_{\ell},0)\,.$
(4)

The coefficients of the R-evolution equation have the following explicit form:

$\displaystyle\gamma_{0}^{R}(n_{\ell})$
$\displaystyle={\textstyle\frac{16}{3}}\,,$ (5)
$\displaystyle\gamma_{1}^{R}(n_{\ell})$
$\displaystyle=96.1039-9.55076\,n_{\ell}\,,$
$\displaystyle\gamma_{2}^{R}(n_{\ell})$
$\displaystyle=1595.75-269.953\,n_{\ell}-2.65945\,n_{\ell}^{2}\,,$
$\displaystyle\gamma_{3}^{R}(n_{\ell})$ $\displaystyle=(12319.\pm
417.)-(9103.\pm 10.)\,n_{\ell}+610.264\,n_{\ell}^{2}-6.515\,n_{\ell}^{3}\,.$

The uncertainties in the ${\cal O}(\alpha_{s}^{4})$ coefficient arise from the
numerical uncertainties in the relation between the $\overline{\mathrm{MS}}$
and the pole masses at this order. In App. A.4 we present an efficient
algorithm to exactly integrate Eq. (3).

Adopting appropriate values for $R$ and $n_{\ell}$ the MSR mass $m_{q}^{{\rm
MSR},(n_{\ell})}(R)$ can be related in a renormalon-free manner to any other
low-scale short-distance mass without the appearance of large logarithms and
can thus be used to also resum potentially large R-evolution logarithms in the
relation of other low-scale short-distance mass schemes. In the presence of
massive quarks with masses lighter than $m_{q}$, the MSR mass has an
$n_{\ell}$-dependent renormalization group evolution and flavor threshold
corrections in close analogy to the renormalization group evolution of the
strong coupling and the $\overline{\mathrm{MS}}$ mass. This allows for the
resummation of large logarithms involving the masses of the lighter massive
quarks. For details we refer to Ref. [2].

##  4 Setup

There are three ways to access the functionalities of the REvolver library:

  * 1.

via the C++ library directly, which might be most suitable for extensive
automated tasks and to interface with other libraries and codes,

  * 2.

via the Wolfram Mathematica [7] interface (using WSTP / MathLink), which is
suitable for interactive tasks and for using in parallel with other
Mathematica features,

  * 3.

via the Python [8] interface (generated using SWIG [9]) for usage in scripts
and interactive execution in Jupyter notebooks [10].

###  4.1 Installation

Note that slightly more detailed instructions for installing the code,
including Windows-specific commands, are given in the README.md file provided
with the source code. Here we only describe the installation procedure for
Linux and MacOS, and only for REvolver itself (not for CMake and other
auxiliaries).

For the compilation of REvolver a C++11 compatible compiler is needed. The
recommended (and tested) choices are gcc on Linux, Apple Clang on MacOS and
MinGW on Windows. It is expected that REvolver compiles on other platforms and
with different compilers as well, although this has not been tested and we do
not provide any specific instructions.

We provide a CMake script with various options controlling which interfaces
and demonstration codes are built. To use the script, at least version 3.1 of
CMake is required.

If REvolver is to be used via Mathematica, Wolfram Mathematica is required in
a version which supports WSTP or MathLink.

For the Python interface, at least version 3 of Python has to be installed,
including the development packages. We note that the Python interface is
currently supported only on Linux and MacOS.

After downloading the code to the local hard drive, open a command line
interface, navigate to the directory code/ and run the commands

      $ mkdir build
      $ cd build/

In the next step, the CMake script will be executed. Depending on which
interfaces are to be prepared, various flags can be set:

  * 1.

wolfr determines if the Mathematica interface is prepared (default: OFF),

  * 2.

py determines if the Python interface is prepared (default: OFF),

  * 3.

`cpp_demo` determines if the C++ demo executable is built (default: ON).

The static REvolver C++ library is always compiled. Note that to compile the
C++ demo executable the library Quadpack++ [11] is used. However, no
additional steps are required by the user since the library is provided with
REvolver. Note that the `Quadpack++` library is only compiled if the flag
`cpp_demo` is set to `ON`.

To execute the CMake script with default flags and compile the code, one has
to run the terminal commands

      $ cmake ..
      $ make install

or in general

      $ cmake [(-D <flag>={ON|OFF})...] ..
      $ make install

where <flag> is a placeholder for one of the flags listed above. For example,
to prepare the Mathematica interface, but not the compilation of the C++ demo
code one would use

      $ cmake -D wolfr=ON -D cpp_demo=OFF ..
      $ make install

The directory code/build/ can be safely removed after the compiling is done.
The resulting libraries and executables can be found at the following
locations:

  * 1.

the static C++ library file:
code/lib/libREvolver.a

  * 2.

the MathLink / WSTP executable, ready to be loaded in a Mathematica notebook:
code/bin/REvolver

  * 3.

the Python module file and dynamic library, ready to be imported in a Python
script:
code/pyREvolver/lib/pyREvolver.py
`code/pyREvolver/lib/_pyREvolver.so`

  * 4.

the C++ demo executable:
code/bin/examples

###  4.2 General Usage

####  C++ Interface

To use the REvolver C++ static library, the respective header file has to be
included which is done via

      #include REvolver.h

and the library has to be properly linked when compiling the code. After
including the header file, the implemented classes and routines are available
in the namespace revo and accessible with the scope resolution prefix revo::
unless the instruction using namespace revo; has been invoked, such that the
resolution prefix is not necessary.

For a demonstration, we refer to the source file code/examples/examples.cpp
and the related executable code/bin/examples (if `cpp_demo=ON` was set).

####  Mathematica Interface

To load the WSTP executable in a Mathematica notebook, execute

      Install["<path to executable>/REvolver"]

with `<path to executable>` referring to the directory path where the
executable is located. For future convenience it might be useful to execute

      CopyFile["<path to executable>/REvolver",
               $UserBaseDirectory <> "/Applications/REvolver"]

which copies the executable to the user base directory of Mathematica, making
it possible to load the executable with the command

      Install["REvolver"]

in the future.

For a demonstration of how to load REvolver in a Mathematica notebook and a
general overview of the available functions, see the demonstration notebook
code/examples/examples.nb provided with the package.

####  Python Interface

The module can be loaded in a Python script or Jupyter notebook with the usual
syntax

      import pyREvolver

assuming that pyREvolver.py and the shared library file (*.so) are located in
the same folder as the script or notebook, or have been added to the module
search path with the following command

      import sys
      sys.path.append(’<path to pyREvolver.py>’)

For a demonstration of how to load and use the module, see the Jupyter
notebook code/examples/examples.ipynb provided with the package.

##  5 Core Structure

As described in the introduction, all functionalities of the library are
centered around instances of the class revo::Core (simply called “Core
objects” or “Cores” in the following) each representing a certain physical
scenario for the quark mass spectrum and the strong coupling and from which
numerical values for quark masses and the strong coupling in specified schemes
and at specified scales can be extracted. In principle, the number of Core
objects defined at the same time is only limited by the available memory,
regardless of the interface used.

Figure 1: Schematic structure of the class revo::Core in C++: the class
contains an instance of the class revo::Alpha as a member providing
functionalities related to the strong coupling such as running and matching,
and an instance of the class revo::Mass providing functionalities related to
the quark masses such as running, matching and conversion. revo::Mass uses
revo::Alpha for coupling evolution, and revo is the namespace where all
relevant REvolver classes are defined.

The schematic structure of the class revo::Core in C++ is depicted in Fig. 1:
the class has objects of the classes revo::Alpha and revo::Mass as members.
The class revo::Alpha has various member functions related to the strong
coupling like running, matching and for obtaining the QCD scale
$\Lambda_{\mathrm{QCD}}$. The class revo::Mass has member functions related to
the evolution and matching of the running masses as well as the extraction of
quark mass values in specified schemes. It uses an instance of the class
revo::Alpha to obtain the necessary coupling values. The respective member
objects of the revo::Core class can be accessed through the member functions
revo::Core::alpha() and revo::Core::masses(), respectively. The class
revo::Core itself represents the frame to access these member functions and
provides additional functionalities related to setting up a physical coupling
and quark mass spectrum scenario, and extending an existing scenario by adding
additional heavier massive quarks.

Although in principle possible, helper classes such as revo::Mass and
revo::Alpha are not meant to be used outside Core objects. All available
functionalities can (and should) be accessed through Core objects.

When using the Mathematica interface, the specific structure of the classes
are not relevant since the wrapper hides most details to fit into the Wolfram
language syntax. To preserve the possibility to have multiple Cores defined at
the same time in Mathematica, a unique name has to be specified for each Core
instance, which is referred to when extracting mass and coupling values or
when extending scenarios.

For concrete usage and examples we refer to Secs. 6 and 7.

##  6 Implemented Functions

In the following descriptions and examples we will assume that the namespace
revo has been introduced in the C++ code with the instruction

      using namespace revo;

such that the scope resolution prefix revo:: can be omitted, and that in
Python the module was loaded with

      from pyREvolver import *

to keep code snippets uncluttered.

The syntax and interface structure in Python is the same as in C++, with a few
exceptions:

  * 1.

If a constant of an enum class type has to be provided as an input, the scope
resolution operator :: has to be exchanged with `_`, e.g. the enumerator MSbar
of type MScheme has to be provided using `MScheme_MSbar` instead of
MScheme::MSbar.

  * 2.

The C++ function Mass::mPole allows for two optional pointer-type inputs to
provide the possibility of accessing several output values (see Sec. 6.2.3).
In Python, instead, within the class Mass the additional member function
mPoleDetailed is provided which returns a tuple of values. In Mathematica the
same functionality is provided by the function MassPoleDetailed.

  * 3.

C++ specific syntax cannot be used, e.g. initializing an std::vector with an
initializer list.

We will treat the C++ and Mathematica interfaces on an equal footing, always
stating the C++ function prototypes and definitions first with the Mathematica
ones following. We will then briefly describe the inputs and outputs, and in
most cases give short examples. If not stated otherwise, the related Python
syntax is the same as in C++. Also, to focus on the essential functionalities
first, we will present all commands without optional parameters at the
beginning and describe additional options in a second step. Note that the
given function prototypes do not always correspond exactly to the ones present
in the source codes to make the descriptions more transparent, e.g. for
template functions in C++ or type restrictions in Mathematica.

For a more detailed and technical documentation of the full functionality and
interface structure of the C++ library, please consider reading the online
doxygen documentation (see https://revolver-hep.gitlab.io/REvolver).

A detailed documentation of the functions accessible via the Mathematica
interface is available through the Mathematica internal documentation.

###  6.1 Constructing a Core and Accessing Scenario Parameters

In the following we describe how to construct Cores in the various interfaces
and how to read out their scenario parameters. The scenario parameters of a
Core uniquely reflect its physical scenario. They include the total flavor
number, the flavor number scheme, value as well as scale of the strong
coupling specified at Core creation, the running masses at reference scales,
the flavor matching scales, and the parameters that specify the precision of
the theoretical input and scheme choices. The latter include the perturbative
orders of renormalization group equations and threshold matching relations,
the lambda parameters setting variations in renormalization group equations,
the variation related to the uncertainty of the perturbative 4-loop
pole-$\overline{\mathrm{MS}}$ mass coefficient, and the coefficients of the
QCD $\beta$-function. All scenario parameters, except for coupling and quark
mass values, acquire default values if not specified at Core creation.

####  6.1.1 C++ / Python only: RunPar and RunParV

In the C++ and Python interfaces, the RunPar struct

      struct RunPar {
        int nf;
        double value;
        double scale;
      };

is used to collect the parameters of the running coupling and masses. RunPar
structs contain the active number of flavors nf specifying the flavor number
scheme, the parameter (coupling or mass) value value and the respective
renormalization scale scale. All numbers referring to quantities with
dimensions of energy handled by REvolver (e.g. masses, renormalization scales
or $\Lambda_{\rm QCD}$) are understood in GeV units.

The related type RunParV is an alias for std::vector<RunPar>, i.e. a
collection of RunPars.

In the following C++ example we define the RunPar structs alphaPar and
alphaPar2 specifying flavor number schemes, values and renormalization scales
for the strong coupling, and the RunParVs mPar and mPar2 containing three
RunPars each, specifying values for running masses of charm, bottom and top
quarks. alphaPar sets a realistic value for the strong coupling
$\alpha_{s}^{(5)}(m_{Z})=0.1181$, while alphaPar2 contains parameters to
specify the strong coupling $\alpha_{s}^{(4)}(4.2\,\mathrm{GeV})$. mPar
defines standard running masses with realistic values of charm, bottom and top
quarks, namely $\overline{m}_{c}=1.3$ GeV, $\overline{m}_{b}\\!=\\!4.2$ GeV,
and $\overline{m}_{t}\\!=\\!163$ GeV, while mPar2 defines values for
different flavor number schemes and scales, specifically
$m^{(6)}_{c}\\!(163.0\,\mathrm{GeV})=\overline{m}_{c}^{(6)}\\!(163.0\,\mathrm{GeV})$,
$m_{b}^{(4)}(4.2\,\mathrm{GeV})=m_{b}^{{\rm MSR},(4)}(4.2\,\mathrm{GeV})$, and
$m_{t}^{(5)}(4.2\,\mathrm{GeV})=m_{t}^{{\rm MSR},(5)}(4.2\,\mathrm{GeV})$. In
the Python example we only define alphaPar and mPar for brevity.

####  C++ example

      RunPar alphaPar = {5, 0.1181, 91.187};
      RunPar alphaPar2 = {4, 0.22491680889566054, 4.2};

      RunParV mPar;
      mPar.push_back({4, 1.3, 1.3});
      mPar.push_back({5, 4.2, 4.2});
      mPar.push_back({6, 163.0, 163.0});

      RunParV mPar2;
      mPar2.push_back({6, 0.6173718176865822, 163.0});
      mPar2.push_back({4, 4.20502733598667, 4.2});
      mPar2.push_back({5, 172.37293079716443, 4.2});

####  Python example

      alphaPar = RunPar(5, 0.1181, 91.187)

      mPar = RunParV(3)
      mPar[0] = RunPar(4, 1.3, 1.3)
      mPar[1] = RunPar(5, 4.2, 4.2)
      mPar[2] = RunPar(6, 163.0, 163.0)

      mPar2 = RunParV(3)
      mPar2[0] = RunPar(6, 0.6173718176865822, 163.0)
      mPar2[1] = RunPar(4, 4.20502733598667, 4.2)
      mPar2[2] = RunPar(5, 172.37293079716443, 4.2)

####  6.1.2 Constructing Core objects with masses

The prototypes for the functions constructing Cores in C++ and Mathematica,
respectively, are

      Core::Core(int nTot, const RunPar& alphaPar,
                 const RunParV& mPar);

      CoreCreate[CoreName_String, nTot_Integer, alphaPar_List,
                 mPar_List]

with the mandatory input nTot, specifying the total number of quark flavors in
the scenario, as well as the input parameters for the strong coupling and the
quark masses. In C++, the coupling and mass parameters are given by RunPar
structs and std::vectors RunParV, respectively, which are described in Sec.
6.1.1, while in Mathematica, the parameter collections are given by lists and
lists of lists, respectively. The argument CoreName in Mathematica specifies
the user-defined unique name of the created Core instance. The given masses
must be sorted in increasing order with respect to their standard running mass
values starting with the lightest. The number of massless quarks in a Core is
equal to nTot minus the number of elements in mPar. To construct a Core
without massive quarks, see Sec. 6.1.3.

####  C++ example

The instructions

      Core core1(6, alphaPar, mPar);
      Core core2(6, alphaPar2, mPar2);

construct two Core objects named core1 and core2, respectively, with a total
flavor number of $6$, and the parameters determining the strong coupling and
masses contained in alphaPar, alphaPar2, mPar and mPar2 as defined in the
example of Sec. 6.1.1. These are the minimal set of parameters that have to be
specified to create Core objects.

####  Mathematica example

To construct the same Cores in Mathematica one can use

      alphaPar = {5, amZdef, mZdef};
      mPar = {{4, 1.3, 1.3}, {5, 4.2, 4.2}, {6, 163.0, 163.0}};
      CoreCreate["core1", 6, alphaPar, mPar]

      alphaPar2 = {4, 0.22491680889566054, 4.2};
      mPar2 = {{6, 0.6173718176865822, 163.0},
               {4, 4.20502733598667, 4.2},
               {5, 172.37293079716443, 4.2}};
      CoreCreate["core2", 6, alphaPar2, mPar2]

using the predefined parameters amZdef = 0.1181 for the strong coupling and
mZdef = 91.187 for the Z-boson mass.

#### Optional parameters

The Core constructor allows to set a number of optional parameters to control
the flavor matching scales, the perturbative order of matching relations and
renormalization group equations, to perform scale variation of the
renormalization group equations, and to vary the 4-loop
pole-$\overline{\mathrm{MS}}$ mass coefficient within its error band. The
values of these optional parameters are a defining property of the physical
scenario represented by a Core object and respected by all functionalities
related to numerical values of the strong coupling and the running masses.

The full C++ constructor prototype is

      Core::Core(int nTot, const RunPar& alphaPar,
                 const RunParV& massPar,
                 const doubleV& fMatch = doubleV(),
                 int runAlpha = kMaxRunAlpha,
                 double lambdaAlpha = 1.0,
                 int orderAlpha = kMaxOrderAlpha,
                 int runMSbar = kMaxRunMSbar,
                 double lambdaMSbar = 1.0,
                 int orderMSbar = kMaxOrderMSbar,
                 int runMSR = kMaxRunMSR,
                 double lambdaMSR = 1.0,
                 int orderMSR = kMaxOrderMSR,
                 double msBarDeltaError = 0.0);

with doubleV being an alias for std::vector<double> set by REvolver. If one of
the optional parameters shown in the constructor above is explicitly
specified, all parameters appearing prior in the argument list must be
specified as well. The values kMaxRunAlpha, kMaxOrderAlpha, kMaxRunMSR,
kMaxOrderMSR, kMaxRunMSbar and kMaxOrderMSbar are predefined constants
representing the respective defaults. In the Mathematica interface, the
optional parameters of the same name can be set individually via the options
parameter syntax, i.e. by adding opt->val after the last regular function
input, as shown in the examples below.

The meaning of the optional parameters is as follows:

  * 1.

fMatch: a vector / list containing elements {f1, f2, ...}, where fn specifies
that the flavor matching scale $\mu_{n}$ for the n-th lightest massive quark
threshold is fn times the standard running mass:
$\mu_{n}=\texttt{fn}\times\overline{m}_{n}$. Default: all fn are set to $1.0$.
(Note that the mass dependence of the flavor threshold corrections is
expressed in terms of the standard running masses as well. The specification
to use a different mass scheme to parameterize the flavor threshold
corrections is not supported.)

  * 2.

runAlpha: the loop order used for the running of the strong coupling. Default:
highest available order which is $5$.

  * 3.

lambdaAlpha: a parameter probing the renormalization scale dependence of the
QCD $\beta$-function. With respect to the perturbative series of the
$\beta$-function truncated at the order set by runAlpha (used for the value
$1.0$), a (runAlpha$+1$) order term is estimated from renormalization scale
variation (when a value different from $1.0$ is specified): the estimate is
obtained by expanding the original perturbative series
$\beta[\alpha_{s}(\mu)]$ (truncated at order runAlpha) in terms of
$\alpha_{s\\!}(\mathtt{lambdaAlpha}\\!\times\\!\\!\mu)$, truncating at order
runAlpha. The result is expanded in $\alpha_{s}(\mu)$ truncating again at
order (runAlpha$+1$). A variation around $1.0$ of ${\cal O}(\pm 10\%)$ leads
to an adequate uncertainty estimation for the known lower orders of the
$\beta$-function, so that variations exceeding this range should be avoided.
Default: $1.0$.

  * 4.

orderAlpha: the loop order used for the strong coupling flavor threshold
matching relations. The renormalization scale dependence of these matching
relations is precise to loop order orderAlpha and independent of the value
specified for runAlpha. Default: highest available order which is $4$.

  * 5.

runMSbar: loop order used for the $\overline{\mathrm{MS}}$ mass running.
Default: highest available order which is $5$.

  * 6.

lambdaMSbar: a parameter probing the renormalization scale dependence of the
anomalous dimension $\gamma$ of the $\overline{\mathrm{MS}}$ mass in analogy
to the parameter lambdaAlpha. The variation is performed by expanding the
original series for $\gamma[\alpha_{s}(\mu)]$ (truncated at order runMSbar) in
terms of $\alpha_{s}(\mathtt{lambdaMSbar}\times\mu)$, truncating at order
runMSbar. A variation around $1.0$ of ${\cal O}(\pm 10\%)$ leads to an
adequate uncertainty estimation for the known lower orders of $\gamma$, so
that variations exceeding this range should be avoided. Default: $1.0$.

  * 7.

orderMSbar: loop order used for the flavor threshold matching relations of the
$\overline{\mathrm{MS}}$ masses. The renormalization scale dependence of these
matching relations is precise to loop order orderMSbar and independent of the
values specified for runAlpha and orderAlpha. Default: highest available order
which is $4$.

  * 8.

runMSR: the loop order used for the MSR mass running. Default: highest
available order which is $4$.

  * 9.

lambdaMSR: a parameter probing the renormalization scale dependence of the
anomalous dimension $\gamma^{R}$ of the MSR mass in analogy to the parameters
lambdaAlpha and lambdaMSbar. The variation is performed by expanding the
original $\gamma^{R}[\alpha_{s}(R)]$ series (truncated at order runMSR) in
terms of $\alpha_{s}(\mathtt{lambdaMSR}\times R)$, truncating at order runMSR.
A variation around $1.0$ by factors of around $0.5$ and $2$ leads to an
adequate uncertainty estimation for the known lower orders of the
$\gamma^{R}$, so that variations exceeding this range should be avoided.
Default: $1.0$.

  * 10.

orderMSR: the loop order used for the flavor threshold matching relations of
the MSR masses associated to the massive quark itself and all lighter massive
quarks. The renormalization scale dependence of these matching relations is
precise to loop order orderMSR and independent of the value specified for
runMSR, runAlpha and orderAlpha. Default: highest available order which is
$4$.

  * 11.

msBarDeltaError: controls the error of the 4-loop coefficient in the
pole-$\overline{\mathrm{MS}}$ mass relation. Should be varied between $-1$ and
$1$ to scan the standard deviation as quoted in Ref. [12]. Default: $0.0$.

  * 12.

precisionGoal: the parameter setting the relative precision of all convergent
infinite sums and iterative algorithms. The input value is clipped to the
range $[10^{-6},10^{-15}]$. Default: $10^{-15}$ which we refer to as machine
precision. The default should be adequate for most applications, but a lower
precision goal may be specified for improving speed.

The quark mass dependence of all flavor matching relations (for the strong
coupling and the running masses) is expressed in terms of the corresponding
standard running masses $\overline{m}_{q}$. Changing this to an arbitrary mass
scheme is not supported in REvolver. This concerns flavor threshold matching
as well as perturbative reexpansions of the strong coupling in other flavor
number schemes. The resulting numerical differences are, however, tiny and
smaller than the corresponding perturbative uncertainties.

####  C++ example

The instruction

      Core core3(6, alphaPar, mPar, {2.0, 1.0, 1.0});

constructs a Core object named core3 with a total flavor number of $6$ and the
parameters specifying the strong coupling and masses contained in alphaPar and
mPar, respectively, as defined in the example of Sec. 6.1.1. The matching
scale of the flavor threshold related to the lightest massive particle is
$2\times\overline{m}_{c}=2.0\times 1.3$ GeV.

####  Mathematica example

The command

      CoreCreate["core3", 6, alphaPar, mPar,
                 fMatch->{2.0, 1.0, 1.0}]

has the same effect as the analogous C++ example, using the lists alphaPar and
mPar defined in the previous example of this section.

####  6.1.3 Constructing Core objects without massive quarks

The functions with the prototypes

      Core::Core(const RunPar& alphaPar);

      CoreCreate[CoreName_String, alphaPar_List]

are used to construct a Core object with massless quarks only and without
specifying any optional parameters. The parameters are analogous to the
massive case described in Sec. 6.1.2.

#### Optional parameters

The full C++ constructor prototype for a Core with only massless quarks is

      Core::Core(const RunPar& alphaPar,
                 int runAlpha = kMaxRunAlpha,
                 double lambdaAlpha = 1.0,
                 const doubleV& beta = doubleV());

where the optional variables runAlpha and lambdaAlpha are analogous to the
massive case described in Sec. 6.1.2 and can be set in Mathematica using
option parameters. With the optional input beta (which is a constant reference
to a std::vector<double> of arbitrary length) one can specify an arbitrary
number of custom $\beta$-function coefficients. Their default values are the
common $\overline{\mathrm{MS}}$ QCD $\beta$-function coefficients up to $5$
loops with all higher order coefficients set to zero. These are defined based
on the $\beta$-function form

$\mu\frac{{\rm d}\alpha_{s}(\mu)}{{\rm d}\mu}=\frac{{\rm
d}\alpha_{s}(\mu)}{{\rm d}\ln(\mu)}=\beta_{\rm
QCD}(\alpha_{s}(\mu))=-2\alpha_{s}(\mu)\\!\sum_{n=0}^{\infty}\beta_{n}\biggl{[}\frac{\alpha_{s}(\mu)}{4\pi}\biggr{]}^{\\!n+1}\,,$
(6)

where the elements of the C++ container `beta` correspond to the ordered list
of coefficients $\beta_{0},\cdots,\beta_{n}$. Note that adding masses to Core
objects with custom $\beta$-function coefficients is not supported and that,
depending on the choice of runAlpha not all coefficients specified by the user
may be used.

In Mathematica the functionality of custom QCD $\beta$-function coefficients
can be used with

      CoreCreate[CoreName_String, alphaPar_List, beta_List]

with `beta` being the list of $\beta$-function coefficients and the additional
option parameters already explained before.

####  6.1.4 Mathematica only: listing and deleting Cores

      CoreList[]
      CoreDelete[CoreName_String]
      CoreDelete[CoreNames_List]
      CoreDeleteAll[]

These commands list the names of the Cores currently defined, and delete
specific or all Cores, respectively. The argument of CoreDelete is a string
referring to a Core name or a list containing several Core names.

#### Example

      In[]:= CoreList[]
      Out[]= {core1, core2, core3}

      In[]:= CoreDelete["core3"]
             CoreList[]
      Out[]= {core1, core2}

where the Core named core3 has been deleted from memory. We explicitly show
In[] and Out[] to separate in- from out-put and assumed that definitions from
previous examples are still valid.

####  6.1.5 Accessing Core parameters

A Core represents a certain physical scenario for the strong coupling and the
quark mass spectrum that also depends on the theoretical approximations and
conventions implemented (with optional/default parameters) specified at the
time the Core was created. The scenario parameters of a Core (coupling, quark
masses, theoretical approximations and conventions) unambiguously specify a
given scenario and can be accessed by dedicated routines. Note that the
scenario parameters of a Core for the strong coupling depend on the way how
the strong coupling was specified when the Core has been created. Therefore it
is possible to create two physically equivalent Cores with differing scenario
parameters for the strong coupling.

While in C++ the scenario parameters of Core objects are returned from
separate functions, some Mathematica commands print collections of them. The
respective function prototypes in C++ are

      int            Core::nTot() const;
      const RunPar&  Alpha::defParams() const;
      const doubleV& Core::standardMasses() const;
      const doubleV& Core::fMatch() const;
      int            Core::getOrder(OrderPar para) const;
      double         Core::getLambda(LambdaPar para) const;
      double         Core::msBarDeltaError() const;
      const doubleV& Core::betaCoefs(int nf) const;

returning the total number of flavors, the defining RunPar related to the
coupling, an std::vector of the standard running masses in increasing order,
an std::vector with the fn factors specifying the flavor matching scales, the
perturbative orders used for coupling and mass evolutions, the lambda scaling
parameters set for coupling and mass evolution, the variation parameter of the
4-loop coefficient in the pole-$\overline{\mathrm{MS}}$ evolution, and the
$\beta$-function coefficients in the nf-flavor number scheme.

The function inputs of the types defined as

      enum class OrderPar {
        runAlpha,
        orderAlpha,
        runMSbar,
        runMSR,
        orderMSbar,
        orderMSR
      };

and

      enum class LambdaPar { lambdaAlpha, lambdaMSbar, lambdaMSR };

i.e. the enumerators of type OrderPar and LambdaPar, respectively, govern to
which evolution (coupling, $\overline{\mathrm{MS}}$, MSR) or matching
procedure the output of the functions getOrder and getLambda refers to.

In Mathematica, the parameters discussed above can be extracted using the
functions

      CoreParams[CoreName_String]
      CoreParamsDetail[CoreName_String]
      BetaCoefs[CoreName_String]

CoreParams returns from the specified Core a list containing the total number
of flavors, the flavor number scheme, value and renormalization scale of the
strong coupling specified at Core creation, and the standard running masses in
increasing order. CoreParamsDetail returns the coupling and mass values at all
flavor matching scales (in the flavor schemes above as well as below the
corresponding threshold), the flavor number scheme, value as well as scale of
the strong coupling specified at Core creation, and all optional parameters
set by CoreCreate and described in the C++ description above. BetaCoefs
returns the $\beta$-function coefficients in all relevant flavor number
schemes with the normalization as given in Eq. (6).

In addition to just printing the parameters, CoreParamsDetail allows for an
optional parameter to which the parameters are saved as a nested list. The
corresponding function prototype is

      CoreParamsDetail[CoreName_String, output_Symbol]

where output is the symbol in which the list is stored.

####  C++ example

The instructions

      core1.nTot();
      core1.getOrder(OrderPar::runAlpha);

return the total number of flavors and the perturbative order used in the
running of the strong coupling, respectively, for the Core named core1. They
correspond to the values (int)6 and (int)5, respectively, given the definition
of core1 from the C++ example of Sec. 6.1.2.

####  Mathematica example

In the following we show how the output of the functions CoreParams and
BetaCoefs looks like for the Core named core1 given in the Mathematica example
of Sec. 6.1.2:

      In[]:= CoreParams["core1"]
      Out[]= {6, {5, 0.1181, 91.187}, 1.3, 4.2, 163.}

showing the list containing the total number of flavors, the flavor number
scheme, value and renormalization scale of the strong coupling given at Core
creation, and the standard running masses in increasing order; and

      In[]:= BetaCoefs["core1"]
      Out[]= nl = 3: {9.,64.,643.833,12090.4,130378.}
             nl = 4: {8.33333,51.3333,406.352,8035.19,58310.6}
             nl = 5: {7.66667,38.6667,180.907,4826.16,15470.6}
             nl = 6: {7.,26.,-32.5,2472.28,271.428}

showing the $\beta$-function coefficients in the relevant flavor number
schemes.

The output of the function CoreParamsDetail is more extensive than that of
CoreParams, giving

      In[]:= CoreParamsDetail["core1"]
      Out[]=

$\hskip 25.0pt\begin{pmatrix}\verb|[core1,
6]|&\texttt{m1-threshold}&\texttt{m2-threshold}&\texttt{m3-threshold}\\\
\texttt{aS-values}&\verb|{0.385234, 0.383676}|&\verb|{0.224917,
0.224684}|&\verb|{0.108577, 0.108555}|\\\ \texttt{m1-values}&\verb|{1.30636,
1.3}|&\verb|{0.947058, 0.945337}|&\verb|{0.617575, 0.617372}|\\\
\texttt{m2-values}&\verb|{4.64086, 4.62912}|&\verb|{4.20503,
4.2}|&\verb|{2.7438, 2.7429}|\\\ \texttt{m3-values}&\verb|{172.819,
172.807}|&\verb|{172.383, 172.373}|&\verb|{163.032, 163.}|\\\ \end{pmatrix}$

             aS input: {5,0.1181,91.187}
             Matching f-factors: {1.,1.,1.}
             runAlpha: 5
             lambdaAlpha: 1.
             orderAlpha: 4
             runMSbar: 5
             lambdaMSbar: 1.
             orderMSbar: 4
             runMSR: 4
             lambdaMSR: 1.
             orderMSR: 4
             msBarDeltaError: 0.
             precisionGoal: 1.*10^-15

for the Core named core1. The set of printed values represents the full
physical content of a Core. Creating a new Core from this output in an
arbitrary way leads to a physically equivalent core. So Cores are created in a
self-consistent way. This is made possible because REvolver’s algorithm to
solve the renormalization group evolution for the strong coupling provides
(machine precision) exact solutions and because its algorithm for the flavor
matching is self-consistent, see the routines described in Sec. 6.2.1.
Therefore, all massive quarks consistently affect each others flavor number
scheme dependent running values depending on the parameters specified at Core
creation.

For a demonstration, consider the output of CoreParamsDetail of the Core named
core2

      In[]:= CoreParamsDetail["core2"]
      Out[]=

$\hskip 25.0pt\begin{pmatrix}\verb|[core2,
6]|&\texttt{m1-threshold}&\texttt{m2-threshold}&\texttt{m3-threshold}\\\
\texttt{aS-values}&\verb|{0.385234, 0.383676}|&\verb|{0.224917,
0.224684}|&\verb|{0.108577, 0.108555}|\\\ \texttt{m1-values}&\verb|{1.30636,
1.3}|&\verb|{0.947058, 0.945337}|&\verb|{0.617575, 0.617372}|\\\
\texttt{m2-values}&\verb|{4.64086, 4.62912}|&\verb|{4.20503,
4.2}|&\verb|{2.7438, 2.7429}|\\\ \texttt{m3-values}&\verb|{172.819,
172.807}|&\verb|{172.383, 172.373}|&\verb|{163.032, 163.}|\\\ \end{pmatrix}$

             aS input: {4,0.224917,4.2}
             Matching f-factors: {1.,1.,1.}
             runAlpha: 5
             lambdaAlpha: 1.
             orderAlpha: 4
             runMSbar: 5
             lambdaMSbar: 1.
             orderMSbar: 4
             runMSR: 4
             lambdaMSR: 1.
             orderMSR: 4
             msBarDeltaError: 0.
             precisionGoal: 1.*10^-15

which is (apart from the Core name and the strong coupling specifications at
Core creation) exactly the same. In fact, the numbers used to create core2 in
the examples of Sec. 6.1.2 have been taken from the output of
CoreParamDetails["core1"].

###  6.2 Extraction of Masses and Couplings

The following section presents the routines to extract the values for the
strong coupling and the running quark masses in any flavor number scheme at
any scale, as well as values in other quark mass schemes from a given Core.

####  6.2.1 Running masses and strong coupling

The functions with the prototypes

      double Mass::mMS(int nfIn, double scale) const;
      double Alpha::operator()(double scale) const;

      MassMS[CoreName_String, nfIn_Integer, scale_Real]
      AlphaQCD[CoreName_String, scale_Real]

are used to extract a running mass (i.e. the $\overline{\mathrm{MS}}$ mass if
the flavor number scheme includes this massive quark, the MSR mass otherwise)
and strong coupling at a specific renormalization scale scale from a Core
where the optional parameters valid in the creation of the Core are respected.
The functions furthermore use an automatic matching convention, which means
that the flavor number scheme of the output is deduced automatically from
scale, with flavor threshold matching at
$\mu_{n}=\mathtt{fn}\times\overline{m}_{n}$, see Sec. 6.1.2. nfIn specifies
for which quark the running mass is returned and refers to the number of
dynamical flavors of the associated standard running mass.

In C++, as indicated by the scope resolution prefixes of Mass::mMS and
Alpha::operator(), these functions are not by themselves members of the class
Core, but members of the classes Mass and Alpha respectively. Consequently, in
practice they are accessed through the Core member functions Core::alpha and
Core::masses (see Sec. 5 the following examples).

Note that matching at flavor thresholds is always done from below to above the
thresholds, meaning that matching between $n_{\ell}$ and $n_{\ell}+1$ flavor
schemes always uses the perturbative expansion in $\alpha_{s}^{(n_{\ell})}$.
If this is not possible directly, e.g. when $\alpha_{s}^{(n_{\ell}+1)}$ is
given and $\alpha_{s}^{(n_{\ell})}$ still needs to be determined, the solution
is computed iteratively. This convention is strictly applied everywhere, and
in particular for the determination of the Core parameters described in Sec.
6.1.5. For the computation of the renormalization group evolution equations
(at any specified order) we use algorithms which are exact, i.e. they provide
results with machine precision. This, in combination with the flavor matching
convention, has the advantage that Cores are created in a self-consistent way.
This means that a Core that is created from the coupling and the masses at any
renormalization scales extracted from an existing core Core will lead to a
Core that is physically equivalent within machine precision if all optional
parameters are set to equivalent values.

Both functions above are also available in a version permitting complex-valued
input for scale, resulting in a complex-valued output. The function prototypes
in this case are

      std::complex<double>
        Mass::mMS(int nfIn, std::complex<double> scale);
      std::complex<double>
        Alpha::operator()(std::complex<double> scale);

      MassMS[CoreName_String, nfIn_Integer, scale_Complex]
      AlphaQCD[CoreName_String, scale_Complex]

and the flavor number scheme of the output is determined from the automatic
matching conventions based on the absolute value of scale.

Note that for applications where negative real values of scale are expected,
the user of the C++ or Python interfaces should declare scale as complex-
valued explicitly from the start so that the intrinsic C++ function
evaluations can be performed. Otherwise, in such a case NaN will be returned.
Using the Mathematica interface, the complex-valued declaration is
automatically employed if scale is negative real (or complex). In the way C++
treats the branch cuts of complex-valued functions involved in the
calculations, this corresponds to adding an infinitesimally small positive
imaginary part to scale.

In case the input parameter scale is not provided when executing Mass::mMS and
MassMS, the routines return the respective standard running mass.

####  C++ example

The instruction

      core1.masses().mMS(6, 20.0);

returns the mass value of the heaviest of the six quarks defined in the Core
named core1, referred to by the flavor number $6$ of the corresponding
standard running mass flavor scheme, at the scale $20.0$ GeV with automatic
flavor matching (performed at the standard running mass of the heaviest quark
in core1). Referring to the corresponding quark mass as $m_{t}$, the returned
value (double)171.046 corresponds to
$m_{t}^{\mathrm{MSR},(5)}(20.0\,\mathrm{GeV})$.

The instruction

      core1.alpha()(10.0);

returns the value of the strong coupling defined in the Core named core1 at
the scale $10.0$ GeV. Due to automatic matching the returned value
(double)0.178468 refers to $\alpha_{s}^{(5)}(10.0\,\rm{GeV})$. Note that the
Core member functions Core::alpha and Core::masses have been used to access
member functions of the classes Alpha and Mass, see Sec. 5.

####  Mathematica example

The Mathematica command

      In[]:= MassMS["core1", 5, 30.0]
      Out[]= 3.199542552851507

returns the mass value of the next-to-heaviest of the six quarks defined in
the Core named core1, referred to by the flavor number $5$, at the scale
$30.0$ GeV. Referring to the corresponding quark mass as $m_{b}$, due to
automatic matching, the shown output corresponds to the
$\overline{\mathrm{MS}}$ mass $\overline{m}_{b}^{(5)}(30.0\,\mathrm{GeV})$.
The command

      In[]:= AlphaQCD["core1", -10.0 + 0.1 I]
      Out[]= 0.11619771241600231 - 0.08307724827828394 I

returns the value of the strong coupling defined in the Core named core1 at
the complex scale $\mu=(-10.0+0.1\,i)$ GeV. Due to the automatic matching
convention, the flavor number scheme is automatically chosen to be $5$ since
$|\mu|=10.0005$ GeV exceeds the standard running bottom mass in core1.
Consequently the output refers to
$\alpha_{s}^{(5)}((-10.0+0.1\,i)\,\mathrm{GeV})$.

#### Optional parameters

Automatic matching at the flavor thresholds can be overruled using the
additional input nfOut, which specifies the flavor number scheme. The function
prototypes are

      double Mass::mMS(int nfIn, double scale,
                       int nfOut = kDefault) const;
      double Alpha::operator()(double scale,
                               int nfOut = kDefault) const;

      MassMS[CoreName_String, nfIn_Integer, scale_Real,
             nfOut_Integer:kDefault]
      AlphaQCD[CoreName_String, scale_Real, nfOut_Integer:kDefault]

and analogously for complex-valued input. In C++ as well as in Mathematica,
the value kDefault is a predefined constant, internally set to $-1$,
specifying that the respective default values will be used.

####  6.2.2 Other short-distance masses

From a given Core, the extraction of values of quark masses in a number of
other short-distance quark mass schemes is supported. This includes the
renormalization group invariant (RGI) scheme [13] and the following low-scale
short-distance mass schemes: 1S (1S) [14, 15, 16], kinetic (Kin) [17],
potential subtracted (PS) [18, 19], and renormalon subtracted222We implement
only the “unprimed” version of the RS mass, which has a finite ${\cal
O}(\alpha_{s})$ term in its relation to the pole mass. (RS) [20]. Note that
the optional parameters setting the loop orders of the conversion formulae
used for the extraction of these masses are independent of the loop order
parameters specified during Core creation (runAlpha, orderAlpha, runMSbar,
orderMSbar, runMSR, orderMSR), see Sec. 6.1.2.

The functions with the prototypes

      double Mass::mRGI(int nfIn) const;
      double Mass::m1S(int nfIn) const;
      double Mass::mKin(int nfIn, double scaleKin) const;
      double Mass::mPS(int nfIn, double muF) const;
      double Mass::mRS(int nfIn, double scaleRS) const;

      MassRGI[CoreName_String, nfIn_Integer]
      Mass1S[CoreName_String, nfIn_Integer]
      MassKin[CoreName_String, nfIn_Integer, scaleKin_Real]
      MassPS[CoreName_String, nfIn_Integer, muf_Real]
      MassRS[CoreName_String, nfIn_Integer, scaleRS_Real]

extract from a Core a quark mass value in the specified short-distance scheme.
The variable nfIn is again the specifier for the massive quark, referring to
the number of dynamical flavors of the associated standard running mass. The
variables scaleKin,333Following Refs. [21, 22], the default value for the
kinetic mass intrinsic scale (where the default log resummed conversion
between the running and the kinetic masses is carried out, see below) is set
to be twice its renormalization scale, $2\times$scaleKin. muF and scaleRS
specify the renormalization scales of the kinetic, potential subtracted and
renormalon subtracted masses, respectively. Note that the 1S, Kin, PS and RS
masses are defined in (nfIn$-1$)-flavor schemes, while the RGI mass is defined
in the nfIn flavor scheme.

In C++, all routines outlined above are member functions of the class Mass, as
indicated by the scope resolution prefix Mass::, and are accessed through the
Core member function Core::masses (see Sec. 5 and the following examples).

Without specifying any additional optional parameters, the extraction of the
low-scale short-distance masses is done in the following default way: first,
the running mass in the ($\mathtt{nfIn}-1$)-flavor scheme is determined at the
intrinsic scale of the low-scale short-distance mass using R-evolution. The
intrinsic scales of the Kin, PS and RS schemes are the respective values of
$2\times$scaleKin, muf and scaleRS, while for the 1S scheme it is the inverse
Bohr radius $M_{q,B}$ determined with the routine Mass::mBohr described below
with default setting. Subsequently, the running mass is converted to the low-
scale short-distance mass using both masses’ perturbative relation to the pole
mass employing the strong coupling at the intrinsic scale, where the pole mass
is then consistently eliminated to the corresponding order. The resulting
relation is renormalon-free avoiding large logarithms. Whenever known, finite
mass effects stemming from lighter massive quarks are taken into account in
this relation, which is up to three loops for the pole mass relations of the
1S, PS and Kin444For the kinetic mass REvolver adopts the definition of
lighter massive quark corrections given in Ref. [22] where these quark mass
corrections are absent in the flavor number scheme below the corresponding
threshold and exclusively come from the flavor number decoupling relations of
the strong coupling above the threshold. schemes. For the RS mass, finite
lighter quark mass effects have not been explicitly specified in the
literature. Note that the lighter quark mass effects in the perturbative
relation between the pole and running ($\overline{\mathrm{MS}}$ or MSR) masses
are known to three loops. In the conversion between the running and the RS
masses these lighter quark mass corrections are set to zero coherently
everywhere to avoid upsetting the renormalon cancellation. Effective methods
to simulate finite quark mass effects (e.g. by enforcing a change in the
flavor number scheme of the strong coupling) are not implemented as regular
REvolver functionalities, but can still be realized, as we show in the
examples given in Sec. 7.

REvolver provides the routine Mass::mBohr to calculate the heavy quarkonium
inverse Bohr radius $M_{q,B}$ for a massive quark $q$. In the default setting,
$M_{q,B}$ is the root of the function
$f(x)=C_{F}\alpha_{s}^{(n_{\ell})}(x)m_{q}^{\mathrm{MSR},(n_{\ell})}(x)-x$,
and determined using an algorithm based on a modified version of Dekker’s
method. The corresponding routine is a global function and always uses machine
precision. The inverse Bohr radius for a massive quark $q$ can be extracted
from a given Core through the functions with the following prototypes

      double Mass::mBohr(int nfIn, double nb = 1.0) const;

      MBohr[CoreName_String, nfIn_Integer, nb_Real:1.0]

where nfIn is the specifier for the massive quark and nb refers to an optional
rescaling factor such that the inverse Bohr radius is determined from the
equality
$M_{q,B}^{\mathtt{nb}}=\mathtt{nb}\,C_{F}\alpha_{s}^{(n_{\ell})}(M_{q,B}^{\mathtt{nb}})m_{q}^{\mathrm{MSR},(n_{\ell})}(M_{q,B}^{\mathtt{nb}})$.
This option is useful for scale variation e.g. in the context of higher
excited heavy quarkonium states. For the computation of the intrinsic scale
used for the 1S mass scheme conversions, $\mathtt{nb}$ is used with the
default setting $\mathtt{nb}\,$=1.

####  C++ example

The instruction

      core1.masses().m1S(6);

returns the mass value in GeV units of the heaviest of the six quarks defined
in the Core named core1, referred to by the flavor number $6$, in the 1S
scheme, which is (double)171.517.

####  Mathematica example

The Mathematica command

      In[]:= MassPS["core1", 5, 2.0]
      Out[]= 4.521091787631138

returns the mass value of the next-to-heaviest of the six quarks defined in
the Core named core1, referred to by the flavor number $5$, in the potential
subtracted scheme at the scale $2.0\,$GeV, i.e.
$m_{b}^{\mathrm{PS}}(2\,\mathrm{GeV})$, referring to that quark mass as
$m_{b}$.

#### Optional parameters

The functions responsible for extracting quark mass values in short-distance
mass schemes other than the running mass allow for several optional parameters
whose nature is tied to the respective schemes. The full prototypes are given
by

      double Mass::mRGI(int nfIn,
                        int order = kMaxRunMSbar) const;
      double Mass::m1S(int nfIn, int nfConv = kDefault,
                       double scale = kDefault,
                       Count1S counting = Count1S::Default,
                       double muA = kDefault,
                       int order = kMaxOrder1s) const;
      double Mass::mKin(int nfIn, double scaleKin,
                        int nfConv = kDefault,
                        double scale = kDefault,
                        double muA = kDefault,
                        int order = kMaxOrderKinetic) const;
      double Mass::mPS(int nfIn, double muF,
                       int nfConv = kDefault,
                       double scale = kDefault,
                       double muA = kDefault,
                       double rIR = 1,
                       int order = kMaxOrderPs) const;
      double Mass::mRS(int nfIn, double scaleRS,
                       int nfConv = kDefault,
                       double scale = kDefault,
                       double muA = kDefault,
                       int order = kMaxRunMSR,
                       int nRS = kMaxRunAlpha - 1,
                       double N12 = kDefault) const;

      MassRGI[CoreName_String, nfIn_Integer,
              order_Integer:kMaxRunMSbar]
      Mass1S[CoreName_String, nfIn_Integer,
             nfConv_Integer:kDefault,
             scale_Real:kDefault,
             counting_String:"default",
             muA_Real:kDefault,
             order_Integer:kMaxOrder1s]
      MassKin[CoreName_String, nfIn_Integer, scaleKin_Real,
              nfConv_Integer:kDefault,
              scale_Real:kDefault,
              muA_Real:kDefault,
              order_Integer:kMaxOrderKinetic]
      MassPS[CoreName_String, nfIn_Integer, muF_Real,
             nfConv_Integer:kDefault,
             scale_Real:kDefault,
             muA_Real:kDefault,
             rIR_Real:1.0,
             order_Integer:kMaxOrderPs]
      MassRS[CoreName_String, nfIn_Integer, scaleRS_Real,
             nfConv_Integer:kDefault,
             scale_Real:kDefault,
             muA_Real:kDefault,
             order_Integer:kMaxRunMSR,
             nRS_Integer:kMaxRunAlpha - 1,
             N12_Real:kDefault]

where the variables kDefault, kMaxRunMSbar, kMaxRunMSR, kMaxOrder1s,
kMaxOrderKinetic and kMaxOrderPs as well as the values "default" and
Count1S::Default in Mathematica and C++, respectively, are predefined
constants representing the respective default values. In C++ and Mathematica,
if one of the optional parameters shown above is explicitly specified, all
parameters appearing prior in the argument list must be specified as well.

In the following the meaning of the optional parameters is explained.

The extraction of the RGI mass has only one optional parameter, which is
order, specifying the loop order of the $\beta$-function and
$\overline{\mathrm{MS}}$ mass anomalous dimensions (which are taken equal)
entering the conversion formula from the standard running mass, see Eq. (27).
The default value is the highest available order, which is $5$.

Considering the functions extracting low-scale short-distance masses, there
are some optional parameters affecting the formulae used for the conversion
computations:

  * 1.

nfConv: specifies the flavor number scheme of the running mass that is used in
the conversion formula. For nfConv = nfIn - 1 the MSR scheme is used; for
nfConv = nfIn the $\overline{\mathrm{MS}}$ scheme is used. The strong coupling
is always employed in the nfIn - 1 flavor scheme. Default: nfIn - 1.

  * 2.

scale: specifies the scale of the running mass from which the conversion is
determined. Default: intrinsic scale of the low-scale mass.

  * 3.

muA: specifies the scale of the strong coupling used for the conversion.
Default: intrinsic scale of the low-scale mass.

  * 4.

order: specifies how many perturbative orders are used in the conversion.
Default: Highest available order for the low-scale mass.

Note that for very heavy quarks (such as the top quark) specifications of the
parameters nfConv, scale and muA that differ from their defaults can lead to
large logarithmic corrections in the conversions involving low-scale short-
distance masses. On the other hand, for the case of lighter massive quarks
(such as the bottom and especially the charm quark) the default settings may
lead to unphysically low scales causing perturbative instabilities. The
parameters nfConv, scale and muA should therefore be used with some care.

Furthermore, for some low-scale masses there are additional optional
parameters related to their definition and specific properties:

  * 1.

counting: specifies the order counting used for the conversion computation to
extract the 1S mass. The available options are Count1S::Nonrelativistic and
Count1S::Relativistic (enumerations of type Count1S) in C++, and
"nonrelativistic" and "relativistic" in Mathematica. For non-relativistic
counting it is assumed that scale is of the same order as the inverse Bohr
radius, which is much smaller than the quark mass value and both are counted
as $\mathcal{O}(m_{n}\alpha_{s})$; for the relativistic counting it is assumed
that scale is of the same order as the quark mass, and the inverse Bohr radius
is counted as $\mathcal{O}(m_{n})$, see Sec. 5.2 of Ref. [1]. If one chooses
to convert from the $\overline{\mathrm{MS}}$ mass, only the relativistic
counting is supported. Note that here the optional parameter order always
specifies how many non-zero terms are taken into account in the conversion
series, e.g. in the case of non-relativistic counting $\mathtt{order}=1$
refers to $\mathcal{O}(m_{n}\alpha_{s}^{2})$, while it refers to
$\mathcal{O}(m_{n}\alpha_{s})$ in the relativistic counting. (So for
$\mathtt{order}=1$ the conversion formula is the same in both counting
schemes.) To simplify language, we refer to a perturbative term in the 1S-pole
mass relation as $n$-loop, if it is combined with the $n$-loop coefficient in
the relation between the running mass and the pole mass, independent of the
employed counting. Default: Count1S::Default and "default" for C++ and
Mathematica, respectively, which corresponds to "nonrelativistic" for an MSR
mass input and to "relativistic" for $\overline{\mathrm{MS}}$.

  * 2.

rIR: specifies the ratio of the IR subtraction scale $\mu_{\mathrm{IR}}$ with
respect to the scale muF employed in the 4-loop term of the pole-PS mass
relation, $\mathtt{rIR}=\mu_{\mathrm{IR}}/\mu_{F}$, see Secs. 4.5.3 and 5.1 of
Ref. [1]. Default: $1.0$, which corresponds to the definition given in Ref.
[19].

  * 3.

nRS: specifies the number of terms (used in the perturbative construction of
the Borel function and the normalization N12) for the calculation of the
coefficients of the pole-RS mass perturbation series, see Ref. [20]. Default:
highest available order which is $4$.

  * 4.

N12: specifies the pole mass renormalon normalization constant employed in the
pole-RS mass relation. Default: value computed by the sum rule formula
employed by the routine Mass::N12 (C++) or N12 (Mathematica), see Sec 6.2.4,
summing up nRS terms in the sum rule series using all available information on
the QCD $\beta$-function and anomalous dimension of the MSR mass.

####  C++ example

The instruction

      core1.masses().mKin(5, 2.0, 4, 5.0, 3.0, 2);

returns the kinetic mass value of the next-to-heaviest of the six quarks
defined in the Core named core1, referred to by the flavor number 5 of the
corresponding standard running mass flavor scheme. The intrinsic kinetic mass
scale is set to $2\,$GeV, the conversion formula is applied to the MSR mass in
the $4$ flavor scheme at $5\,$GeV and with the renormalization scale of the
strong coupling set to $3\,$GeV. Terms up to $\mathcal{O}(\alpha_{s}^{2})$ are
included in the conversion formula. The returned value is (double)4.3105.

####  Mathematica example

The Mathematica command

      In[]:= MassRS["core1", 6, 20.0, 6, 163.0, 60.0, 3]
      Out[]= 170.48477730509723

returns the RS mass value of the heaviest of the six quarks defined in the
Core named core1, referred to by the flavor number 6. The intrinsic scale of
the RS mass is set to $20\,$GeV and the conversion formula is applied to the
$\overline{\mathrm{MS}}$ mass in the $6$ flavor scheme. The
$\overline{\mathrm{MS}}$ mass renormalization scale is specified to be
$163\,$GeV, while the renormalization scale of the strong coupling is set to
$60\,$GeV. The last parameter shown above specifies that three perturbative
orders are used in the conversion. For nRS and N12 no inputs are specified,
consequently the default setting is applied.

####  6.2.3 Pole mass

Quark mass values in the pole mass scheme can be extracted from a given Core.
The pole mass scheme suffers from a renormalon ambiguity so that there are
several options to quote a value. REvolver supports two ways of extracting
pole quark mass values, accessible through the functions corresponding to the
prototypes

      double Mass::mPoleFO(int nfIn, int nfConv, double scale,
                           double muA, int order) const;
      double Mass::mPole(int nfIn, double scale) const;

      MassPoleFO[CoreName_String, nfIn_Integer, nfConv_Integer,
                 scale_Real, muA_Real, order_Integer]
      MassPole[CoreName_String, nfIn_Integer, scale_Real]

All parameters shown have to be specified by the user. For both routines the
variable nfIn is the specifier for the quark whose pole mass value is
extracted, referring to the number of dynamical flavors of the associated
standard running mass, while scale specifies the scale of the running mass for
which the conversion formula to the pole mass is employed. Furthermore, the
loop orders of the various components entering the conversion formula used in
the computation of the pole mass value are provided by the settings of the
loop order parameters specified during Core creation (runAlpha, orderAlpha,
runMSbar, orderMSbar, runMSR, orderMSR), see Sec. 6.1.2. Here the parameters
runMSR and runAlpha are particularly important because they also set the loop
order of the exact coefficients accounted for in the perturbative series
between the MSR and the pole masses, see Eq. (36) in Sec. A.5. Because in
REvolver matching is always applied at flavor thresholds, the latter series
are the only ones being affected by the renormalon.

The functions Mass::PoleFO and MassPoleFO, respectively, return the pole mass
using the conversion formula up to the order set by order from the running
mass in the flavor number scheme specified by nfConf. The input parameter muA
specifies the renormalization scale of the strong coupling used in the
conversion formula. The input order allows, in principle, an arbitrarily large
integer. For $\mathtt{order}\leq\min(\texttt{runMSR},\,\texttt{runAlpha}-1)$
the exact perturbative coefficients up to this order are utilized in the
relation between MSR and pole mass. For
$\mathtt{order}>\min(\texttt{runMSR},\,\texttt{runAlpha}-1)$ a renormalon-
based asymptotic formula for the coefficients of the asymptotic series is
used, see Eq. (36) in Sec. A.5 for details.

The functions Mass::Pole and MassPole, respectively, return the value of the
asymptotic pole mass, i.e. the conversion is done by summing the perturbative
series for the conversion formula relating running and pole mass to the order
of minimal correction, again potentially employing the already mentioned
asymptotic formula if the order of minimal correction is larger than
$\min(\texttt{runMSR},\,\texttt{runAlpha}-1)$. The asymptotic pole mass is
always determined from the MSR mass in the flavor number scheme where all
massive quarks are integrated out regardless of the value of scale. The method
used for the calculation of the asymptotic value is the minimal correction
method ("min") described below.

####  C++ example

The instruction

      core1.masses().mPole(5, 2.0);

returns the asymptotic pole mass of the next-to-heaviest of the six quarks
defined in the Core named core1 in GeV units, referred to by the flavor number
$5$, using the minimal correction method, which is (double)4.91150. The scale
employed for the running mass at conversion is $2\,$GeV.

####  Mathematica example

The Mathematica command

      In[]:= MassPoleFO["core1", 6, 5, 20.0, 20.0, 16]
      Out[]= 1190.2576448418606

returns the order-dependent pole mass value of the heaviest of the six quarks
defined in the Core named core1, referred to by the flavor number 6. The
conversion formula is applied to the respective MSR mass with $5$ active
flavors at the scale $20\,$GeV, choosing the same scale for the
renormalization scale of the strong coupling. The perturbative series terms
are summed up to ${\cal O}(\alpha_{s}^{16})$.555The value obtained with this
command is unphysically large due to the asymptotic (non-convergent) nature of
the series that relate the pole mass and short-distance masses. This behavior
is the origin of the ambiguity of the pole mass.

#### Optional parameters

The functions Mass::mPoleFO and MassPoleFO in C++ and Mathematica,
respectively, are completely general and do not support any additional
optional parameters.

The full prototypes of the functions responsible for accessing the asymptotic
pole mass are

      double Mass::mPole(int nfIn, double scale,
                         double muA = kDefault,
                         PoleMethod method = PoleMethod::Default,
                         double f = 1.25,
                         double* ambiguity = nullptr,
                         int* nMin = nullptr) const;

      MassPole[CoreName_String, nfIn_Integer, scale_Real,
               muA_Real:kDefault,
               method_String:"min",
               f_Real:1.25]
      MassPoleDetailed[CoreName_String, nfIn_Integer, scale_Real,
                       muA_Real:kDefault,
                       method_String,
                       f_Real:1.25]

where MassPoleDetailed in Mathematica is a function returning a list,
containing the asymptotic pole mass, the associated renormalon ambiguity and
the order of the minimal correction term, providing the functionality
corresponding to the C++ routine Mass::mPole (with specified optional pointer
parameters ambiguity and nMin), see the descriptions below. In Python there is
the analogous function Mass.MassPoleDetailed.

The meaning of the optional input parameters is as follows:

  * 1.

muA: specifies the renormalization scale of the strong coupling in the
conversion formula. The default is the value of scale.

  * 2.

method: specifies the method used to obtain the asymptotic pole mass and the
ambiguity. In C++ the available options are PoleMethod::Min for the minimal
correction term method, PoleMethod::Range for the range method and
PoleMethod::DRange for the corresponding discrete version, see the explanation
below. In Mathematica these options correspond to the string inputs "min",
"range" and "drange", respectively. The input method is optional if only the
asymptotic pole mass is returned; the defaults in that case are
PoleMethod::Min (C++) and "min" (Mathematica). The input method is not
optional if the pole mass ambiguity is also returned, i.e. for
MassPoleDetailed in Mathematica and when the optional pointer to ambiguity is
given in C++.

  * 3.

f: specifies a constant larger than unity multiplying the minimal correction
for the methods "range" and "drange". Default: $1.25$

  * 4.

ambiguity: a pointer to double. If specified, the value of the pole mass
ambiguity is saved in the variable pointed to (C++ only).

  * 5.

nMin: a pointer to int. If specified, the order of the minimal correction term
is saved in the variable pointed to (C++ only).

Note that for the Mathematica routine MassPoleDetailed the input variable
method does not have a default and must always be specified (because it always
returns a value for the pole mass). Calling the function with $4$ arguments
means that values are specified for the four input parameters CoreName, nfIn,
scale, and method, while the variables muA and f are set to their default
values.

The minimal correction method ("min") to obtain the asymptotic pole mass value
and its ambiguity refers to the method suggested in Ref. [3] where the
ambiguity is determined from the size of the minimal correction based on a
quadratic function fitted to the smallest correction and the two neighboring
corrections. However, in contrast to the procedure described in Ref. [3],
REvolver accounts for the mass effects of lighter massive quarks by the exact
expressions given in Ref. [2] instead of including them in an approximate way
by flavor number scheme modifications of the strong coupling. Also, for
coefficients of order higher than
$\min(\texttt{runMSR},\,\texttt{runAlpha}-1)$ the asymptotic formula of Eq.
(36) in Sec. A.5 is employed.

The drange (“discrete range”, "drange") choice refers to the method suggested
in Ref. [2], where the pole mass value and its ambiguity are computed from the
range in orders around the minimal term where the corrections are smaller than
f times the minimal correction. The method range ("range") refers to a
continuous generalization which is analogous but provides smoother results.
Here the order-dependent discrete-valued individual perturbative coefficients
of the relation between the pole and running masses, as well as the related
cumulant, are made continuous by a cubic interpolation. From these functions,
the asymptotic pole mass value and the ambiguity are determined in analogy to
the "drange" method. For all methods (including "range"), the returned order
of the minimal correction term is an integer and refers to the original series
without any interpolation.

####  C++ example

With the instructions

      double ambiguity;
      int nMin;
      core1.masses().mPole(6, 10.0, 10.0, PoleMethod::Range, 1.25,
                           &ambiguity, &nMin);

first the variables nMin and ambiguity are initialized and pointers to them
are passed in the call of the function Mass::mPole. The instruction returns
the asymptotic pole mass value (double)173.107 for the heaviest of the $6$
quarks of the Core named core1, obtained from the conversion formula for the
running mass at the scale $10.0\,$GeV for $3$ active flavors [ $m_{t}^{{\rm
MSR},(3)}(10,\mbox{GeV})$ ] using the range method. The values of the pole
mass ambiguity (double)0.179811, and the order of the minimal correction term
(int)4 are stored in the variables ambiguity and nMin, respectively.

####  Mathematica example

The Mathematica command

      In[]:= MassPoleDetailed["core1", 6, 10.0, "min"]
      Out[]= {173.09681693689498, 0.1307735468951421, 4}

corresponds to the previous commands given in C++, however, here the minimal
correction method is used. muA is automatically set to default value which is
10.0 in this case. The output list entries correspond to the asymptotic pole
mass value, the ambiguity and the order of the smallest correction term,
respectively.

####  6.2.4 Norm of the pole mass renormalon ambiguity

The pole mass renormalon normalization constant can be accessed using the
functions with the following prototypes

      double Mass::N12(double lambda = 1.0) const;
      double Mass::P12(double lambda = 1.0) const;

      N12[CoreName_String, lambda_Real:1.0]
      P12[CoreName_String, lambda_Real:1.0]

where the two functions correspond to the normalization conventions $P_{1/2}$
and $N_{1/2}=\Gamma(1+\hat{b}_{1})\beta_{0}P_{1/2}/(2\pi)$ as described in
Ref. [1]. The routines employ the renormalon sum rule formula shown in Eq.
(38) of App. A.6. The number of massless quarks for which the normalization is
determined is tied to the Core used to extract the normalization. The optional
input parameter lambda is a scaling parameter to estimate the uncertainty of
the output ($1.0$ by default). The number of terms summed up in the sum rule
formula is set by runMSR specified at Core construction. The QCD
$\beta$-function coefficients entering the sum rule formula are used up to
runAlpha loop order, all higher order coefficients are set to zero, so that
$\hat{b}_{1}=0$ when runAlpha $=1$.

####  6.2.5 Extracting $\Lambda_{\mathrm{QCD}}$

The functions corresponding to the prototypes

      double Alpha::lambdaQCD(int nf) const;

      LambdaQCD[CoreName_String, nf_Integer]

return from a given Core the QCD scale
$\Lambda_{\mathrm{QCD}}^{(\mathtt{nf})}$ in the $\overline{\mathrm{MS}}$
definition (see Ref. [23]) for the QCD coupling in the nf flavor scheme,
utilizing the exponential formula given in Eq. (20) of App. A.2. The
parameters runAlpha and lambdaAlpha specified at Core creation set the number
of coefficients of the $\beta$-function used. The possible values for nf range
from the number of massless quarks to the total number of quarks in the
specified Core.

####  C++ example

The instruction

      core1.alpha().lambdaQCD(3);

returns $\Lambda_{\mathrm{QCD}}^{(3)}$ in the $\overline{\mathrm{MS}}$
definition in GeV units from the scenario encoded in core1, corresponding to
(double)0.335547

#### Optional parameters

The definition in which $\Lambda_{\mathrm{QCD}}$ is extracted can be chosen
with an optional parameter. The full function prototypes are

      double Alpha::lambdaQCD(
          int nf, LambdaConvention convention =
                      LambdaConvention::MSbar) const;

      LambdaQCD[CoreName_String, nf_Integer,
                convention_String:"MSbar"]

In addition to the conventional $\overline{\mathrm{MS}}$ definition according
to the PDG 2020 [23], which is the default, the “t-scheme” is supported
referring to the definition based on the $t$-variable notations of Ref. [1].
To choose the t-scheme the optional parameter convention has to be set to
LambdaConvention::tScheme and "tScheme" in C++ and Mathematica respectively.
See App. A.2 for more details on the different conventions.

###  6.3 Adding Masses to an Existing Core

The following section presents routines to extend existing Core scenarios by
adding one heavier massive quark. The mass value of the additional quark can
be given as a running mass, in the pole mass scheme (order-dependent as well
as asymptotic), or in any of the short-distance schemes described in Sec.
6.2.2, with the same obligatory and optional parameters. Unless the mass of
the additional quark is already given in the running mass scheme, the quark
mass value is first converted to the running mass (specified by the
parameters). The value for the running mass is then added to the scenario. The
optional parameter fnQ specifies the matching scale for the new threshold,
analogous to the individual entries of fMatch at Core creation, see Sec.
6.1.2.

The routines employed for converting to the running mass are the exact inverse
of those used for extracting mass values in the respective scheme. This is
achieved by numerically inverting the respective relations with iterative
algorithms.

The related function prototypes, including all relevant optional parameters
are

      void Core::addMsMass(int nf, double mass, double scale,
                           double fnQ = 1.0);
      void Core::addPoleMass(double mPole, double scale,
                             double muA = kDefault,
                             PoleMethod method =
                               PoleMethod::Default,
                             double f = 1.25,
                             double fnQ = 1.0);
      void Core::addPoleMassFO(double mPole, int nfConv,
                               double scale, double muA, int order,
                               double fnQ = 1.0);
      void Core::addPSMass(double mPS, double muF,
                           int nfConv = kDefault,
                           double scale = kDefault,
                           double muA = kDefault,
                           double rIR = 1,
                           int order = kMaxOrderPs,
                           double fnQ = 1.0);
      void Core::add1SMass(double m1S,
                           int nfConv = kDefault,
                           double scale = kDefault,
                           Count1S counting = Count1S::Default,
                           double muA = kDefault,
                           int order = kMaxOrder1s,
                           double fnQ = 1.0);
      void Core::addRGIMass(double mRGI,
                            int order = kMaxRunMSbar,
                            double fnQ = 1.0);
      void Core::addRSMass(double mRS, double scaleRS,
                           int nfConv = kDefault,
                           double scale = kDefault,
                           double muA = kDefault,
                           int order = kMaxRunMSR,
                           int nRS = kMaxRunAlpha - 1,
                           double N12 = kDefault,
                           double fnQ = 1.0);
      void Core::addKinMass(double mKin, double scaleKin,
                            int nfConv = kDefault,
                            double scale = kDefault,
                            double muA = kDefault,
                            int order = kMaxOrderKinetic,
                            double fnQ = 1.0);

in C++ and

      AddMSMass[CoreName_String, NewCoreName_String, nf_Integer,
                mass_Real, scale_Real]
      AddPoleMass[CoreName_String, NewCoreName_String, mPole_Real,
                  scale_Real,
                  muA_Real:kDefault,
                  method_String:"min",
                  f_Real:1.25]
      AddPoleMassFO[CoreName_String, NewCoreName_String,
                    mPole_Real, nfConv_Integer, scale_Real,
                    muA_Real, order_Integer]
      AddPSMass[CoreName_String, NewCoreName_String, mPS_Real,
                muF_Real,
                nfConv_Integer:kDefault,
                scale_Real:kDefault,
                muA_Real:kDefault,
                rIR_Real:1,
                order_Integer:kMaxOrderPs]
      Add1SMass[CoreName_String, NewCoreName_String, m1S_Real,
                nfConv_Integer:kDefault,
                scale_Real:kDefault,
                counting_String:"default",
                muA_Real:kDefault,
                order_Integer:kMaxOrder1s]
      AddRGIMass[CoreName_String, NewCoreName_String, mRGI_Real,
                 order_Integer:kMaxRunMSbar]
      AddRSMass[CoreName_String, NewCoreName_String, mRS_Real,
                scaleRS_Real,
                nfConv_Integer:kDefault,
                scale_Real:kDefault,
                muA_Real:kDefault,
                order_Integer:4,
                nRS_Integer:kMaxRunAlpha - 1,
                N12_Real:kDefault]
      AddKinMass[CoreName_String, NewCoreName_String, mKin_Real,
                scaleKin_Real,
                nfConv_Integer:kDefault,
                scale_Real:kDefault,
                muA_Real:kDefault,
                order_Integer:kMaxOrderKinetic]

in Mathematica, where fnQ can be set by the option parameter syntax (fnQ ->
<value>).

In Mathematica, adding a new heavier quark results in creating a new Core with
a (not yet assigned) name specified in the obligatory input argument
NewCoreName. This new Core is based on the Core with the name specified in
CoreName containing only lighter massive quarks. The scenario parameters of
the Core CoreName are passed on to the new Core named NewCoreName. If it turns
out that the standard running mass associated to the mass to be added to a
Core is not the heaviest one in the new configuration, an error is returned
and the new Core is not created.

In C++ the new heavier quark is added directly to the Core from which the
member function is called. This saves resources and copying the old Core can
be easily done applying the assignment operator =, see the C++ example below.

####  C++ example

With the instructions

      RunParV mPar4;
      mPar4.push_back({4, 1.3, 1.3});
      mPar4.push_back({5, 4.2, 4.2});
      Core core4(5, alphaPar, mPar4);
      Core core5 = core4;
      core5.add1SMass(171.51726494075493);

in the first line, the std::vector<RunPar> mPar4 is declared, and subsequently
filled in the next two lines according to Sec. 6.1.1 with information on two
massive quarks. Together with alphaPar taken from the example in Sec. 6.1.1,
it is the input for creating the Core named core4 in the third line,
containing a total of $5$ quark flavors. In the fifth line a copy of core4,
named core5, is produced. In the last line a heavier quark with a mass value
in the 1S scheme is added to core5 which contains a total number of $6$
flavors. The input value for the 1S mass is taken from the C++ example of Sec.
6.2.2 where the 1S mass value was extracted from core1. Therefore core5 and
core1 contain physically equivalent scenarios. core4 has not been modified by
adding the additional heavier quark and consequently its scenario still
contains $5$ quark flavors.

####  Mathematica example

The command

      In[]:= mPar4 = {{4, 1.3, 1.3}, {5, 4.2, 4.2})};
             CoreCreate["core4", 5, alphaPar, mPar4]
             Add1SMass["core4", "core5", 171.51726494075493]

has the same effect as the C++ example. First a list mPar4 is defined,
containing the standard running masses of two quarks, which are used as an
input for creating core4 with a total number of $5$ flavors. In the next line
the new Core named core5 is created, which is a copy of core4, to which a
heavier quark is added with a mass value specified in the 1S scheme.

##  7 Applications and Pedagogical Examples

In this section we demonstrate some of the features provided by REvolver in a
number of concrete examples as they may arise in practical applications. We
give the examples in terms of Mathematica code since, due to its interactive
nature, it is especially suitable for that purpose. We provide all examples
given here as a Mathematica notebook, as Python code in form of a Jupyter
Notebook as well as C++ code, together with the library. For sake of clarity,
all digits of the Mathematica output are displayed.

###  7.1 Cores without Massive Quarks

####  7.1.1 Strong coupling evolution

#### Strong coupling from inclusive jet cross sections

In this first application we demonstrate one possible way to make use of the
strong coupling evolution in REvolver by reproducing and analyzing some
results of Ref. [24], where the CMS collaboration carried out a strong
coupling measurement from inclusive jet cross sections in different $p_{T}$
bins based on $8$ TeV LHC data. In that publication, the evolution of the
strong coupling from different values of $Q$ (the average $p_{T}$ in the bins)
to $m_{Z}$ was carried out with $n_{f}=5$ active flavors and 2-loop accuracy.
In the following, we focus on the result
$\alpha^{(5)}_{s}(Q)=0.0822^{+0.0034}_{-0.0031}$ for $Q=1508.04$ GeV shown in
Table 5 of that article.

In a first step the relevant values can be defined by

    In[]:= {aQCentral, Q} = {0.0822, 1508.04};
           {aQMin, aQMax} = aQCentral + {-0.0031, 0.0034};

and the respective Cores for the central, upper and lower strong coupling
values, setting the strong-coupling evolution to 2-loop order, can be created
by

      In[]:= CoreCreate["central2", {5, aQCentral, Q},
               runAlpha -> 2];
             CoreCreate["min2", {5, aQMin, Q}, runAlpha -> 2];
             CoreCreate["max2", {5, aQMax, Q}, runAlpha -> 2];

The central value of $\alpha^{(5)}_{s}(m_{Z})$ can then be extracted by

      In[]:= amZCentral2 = AlphaQCD["central2", mZdef]
      Out[]= 0.11616452350227859

coinciding with $\alpha_{s}^{(5)}(m_{Z})=0.1162$, as given in Ref. [24].
Likewise the quoted uncertainties ${}^{+0.0070}_{-0.0062}$ are easily
reproduced by executing

     In[]:= (AlphaQCD[#, mZdef] & /@ {"max2","min2"}) - amZCentral2
     Out[]= {0.007009865309335922, -0.006168760659529524}

Note that in the Mathematica interface the variable mZdef has the predefined
value $91.187$ GeV.

For comparison, we also employ 5-loop evolution for the strong coupling,
resulting in equivalent numbers after rounding:

      In[]:= CoreCreate["central5", {5, aQCentral, Q},
               runAlpha -> 5];
             CoreCreate["min5", {5, aQMin, Q}, runAlpha -> 5];
             CoreCreate["max5", {5, aQMax, Q}, runAlpha -> 5];

      In[]:= amZCentral5 = AlphaQCD["central5", mZdef]
             AlphaQCD[#, mZdef] & /@ {"max5", "min5"} - amZCentral5
      Out[]= 0.11624523920392597
      Out[]= {0.007030395615252941, -0.006184136846029031}

Note that the specification of the running order runAlpha -> 5 at Core
creation is not mandatory since 5-loop running is the default. We show the
specification to be explicit.

We now have a look at the perturbative uncertainty of these results. One
approach to estimate the perturbative uncertainty of the given 2-loop result
is to consider the difference of the values obtained by 2-loop and 1-loop
evolution. This can be easily computed using REvolver by creating a new Core
with specified 1-loop coupling evolution:

      In[]:= CoreCreate["central1", {5, aQCentral, Q},
               runAlpha -> 1];
             amZCentral2 - AlphaQCD["central1", mZdef]
      Out[]= 0.0017745561167654411

The obtained conservative perturbative error estimate is about $25\%$ of the
stated experimental error. For the highest available perturbative order for
running, $5$ loop, this proportion shrinks to $0.005\%$:

      In[]:= CoreCreate["central4", {5, aQCentral, Q},
               runAlpha -> 4];
             amZCentral5 - AlphaQCD["central4", mZdef]
      Out[]= 3.573226484421266*^-7

Figure 2:  Values of $\alpha^{(5)}_{s}(m_{Z})$ over the $\beta$-function
scaling parameter $\lambda$, obtained by evolving the strong coupling down
from $\alpha^{(5)}_{s}(1508.04\,\mathrm{GeV})$ with different loop orders in
the evolution equation. The value of $\alpha^{(5)}_{s}(m_{Z})$ computed by
including one more perturbative order in the evolution with $\lambda=1$ is
subtracted to highlight for which ranges of $\lambda$ both results for
$\alpha^{(5)}_{s}(m_{Z})$ coincide or come close. The results indicate that a
variation for $\lambda$ of about $10\%$ around $1$ is a reasonable range to
estimate the perturbative uncertainty.

A different approach to estimate the perturbative error of
$\alpha_{s}^{(5)}(m_{Z})$ is to vary the $\beta$-function scaling parameter
$\lambda$ (controlled by the optional parameter lambdaAlpha) as described in
Sec. 6.1.2. In Fig. 2 the value of $\alpha^{(5)}_{s}(m_{Z})$ is shown,
obtained by evolving the strong coupling down from
$\alpha^{(5)}_{s}(1508.04\,\mathrm{GeV})$ with different loop orders in the
evolution equation and with varying values of $\lambda$; the value of
$\alpha^{(5)}_{s}(m_{Z})$ obtained by including one more perturbative order in
the evolution with $\lambda=1$ (which corresponds to the standard form of the
$\beta$-function) is subtracted. We observe that, to reproduce a value near to
the one obtained with one more perturbative order included in the running,
$\lambda$ has to be varied by about $10\%$ around $1$, as already described in
Sec. 6.1.2.

Assuming this to be also the appropriate range to estimate the perturbative
uncertainty of the 5-loop result, the error can be estimated by scanning over
$\lambda$ values. Employing REvolver for this task, we first define a list of
$20$ logarithmically distributed values of $\lambda$ in the appropriate range

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

and create one Core for each value in $\lambda$. For convenience, the Core
names are set to the associated $\lambda$ values

      In[]:= CoreCreate[ToString[#], {5, aQCentral, Q},
                 lambdaAlpha -> #] & /@ lamList;

The span of values of $\alpha^{(5)}_{s}(m_{Z})$ corresponding to the range in
$\lambda$ can then be obtained with

      In[]:= aLamList = AlphaQCD[ToString[#], mZdef] & /@ lamList;

The central value of this range and the associated error are consequently
given by

      In[]:= (Max[aLamList] + {1, -1} * Min[aLamList])/2
      Out[]= {0.11624517272452983, 1.0224197660724244*^-7}

giving the same order of magnitude as the conservative approach. The variation
in $\lambda$ of $\alpha^{(5)}_{s}(m_{Z})$ when employing 5-loop running is
also depicted in Fig. 3.

Figure 3:  Variation of $\alpha^{(5)}_{s}(m_{Z})$ with $\lambda$, obtained by
evolving the strong coupling down from
$\alpha^{(5)}_{s}(1508.04\,\mathrm{GeV})$ with 5-loop accuracy. The value of
$\alpha_{s}(m_{Z})$ for $\lambda=1$ is subtracted. Varying $\lambda$ in a
range of about $10\%$ around $1$ the value of $\alpha_{s}(m_{Z})$ changes by
$10^{-7}$.

####  7.1.2 Complex renormalization scales

In this application we illustrate the REvolver functionality to determine the
strong coupling $\alpha_{s}(\mu)$ at complex scales $\mu$. To this end we
consider the analyses in Refs. [25, 26], where the perturbative QCD
corrections $\delta^{(0)}$ to the inclusive $\tau$ hadronic width were
considered in fixed-order perturbation theory (FOPT) as well as contour-
improved perturbation theory (CIPT). In the second reference a scheme
different from the usual $\overline{\mathrm{MS}}$ definition is used for
$\alpha_{s}$. We exploit this fact to demonstrate the ability of REvolver to
deal with user-defined $\beta$-function coefficients.

#### Hadronic $\tau$ decay

In CIPT, determining $\delta^{(0)}$ involves integrating over powers of the
strong coupling $\alpha_{s}^{(3)}(\mu)$ multiplied with a kinematic weight
function in the complex $\mu^{2}$ plane along a circle with radius
$m_{\tau}^{2}$

$\displaystyle\delta^{(0)}_{\mathrm{CI}}$
$\displaystyle=\sum_{n=1}^{\infty}c_{n,1}J_{n}^{a}(m_{\tau}^{2})\,,$
$\displaystyle J_{n}^{a}(m_{\tau}^{2})$ $\displaystyle=\frac{1}{2\pi
i}\oint_{|x|=1}\\!\frac{\mathrm{d}x}{x}(1-x)^{3}(1+x)\biggl{[}\frac{\alpha_{s}^{(3)}(\sqrt{-m_{\tau}\,x})}{\pi}\biggr{]}^{n}\,,$
(7)

where the coefficients $c_{n,1}$ are given in Eqs. (2.13), (2.15) and (3.10)
of Ref. [25]. In that article, 4-loop running for the strong coupling was
used, as well as the input value $\alpha_{s}^{(3)}(m_{\tau})=0.34$. To
reproduce the numbers for $\delta^{(0)}_{\mathrm{CI}}$ given in Eq. (3.9) of
Ref. [25] we first define variables accounting for the employed input values

     In[]:= {nfa, aTau, mTau} = {3, 0.34, 1.77686};
            {c11, c21, c31, c41, c51} =
              {1, 1.640, 6.371, 49.076, 283};

and a function to compute $J_{n}^{a}(m_{\tau}^{2})$

      In[]:= Ja[n_, core_] := 1/(2 Pi)
                NIntegrate[(1 - E^(I phi))^3 (1 + E^(I phi))
                  (AlphaQCD[core, Sqrt[-mTau^2 E^(I phi)]]/Pi)^n,
                 {phi, 0, 2 Pi}];

where we change variables via $x=\exp(i\phi)$. Furthermore, we define a Core
with 4-loop running for the coupling

      In[]:= CoreCreate["CI4", {nfa, aTau, mTau}, runAlpha -> 4]

All requirements are now set to reproduce the numbers given in Eq. (3.9) of
Ref. [25]

      In[]:= Ja[1, "CI4"]*c11
             Ja[2, "CI4"]*c21
             Ja[3, "CI4"]*c31
             Ja[4, "CI4"]*c41
             Ja[5, "CI4"]*c51
      Out[]= 0.14789839179248082 + 0. I
      Out[]= 0.02968556861749222 - 1.3583618952083769*^-18 I
      Out[]= 0.0121854520515486 + 4.947098114161149*^-18 I
      Out[]= 0.008592183988059407 - 8.468364884304664*^-19 I
      Out[]= 0.0037863076541974016‘ - 9.919281576456942*^-19 I

with full agreement.

We can now easily inspect the corrections induced by the 5-loop running of the
strong coupling by creating an appropriate Core

      In[]:= CoreCreate["CI5", {nfa, aTau, mTau}, runAlpha -> 5];

and evaluating

      In[]:= Ja[1, "CI5"]*c11
             Ja[2, "CI5"]*c21
             Ja[3, "CI5"]*c31
             Ja[4, "CI5"]*c41
             Ja[5, "CI5"]*c51
      Out[]= 0.14775354880331787 - 8.834874115176436*^-18 I
      Out[]= 0.02960219692750483 - 1.811149193611169*^-18 I
      Out[]= 0.012122455789576415 + 3.4080009230887917*^-18 I
      Out[]= 0.008521503950040327 - 1.6936729768609328*^-18 I
      Out[]= 0.0037376807938461127 - 1.5260433194549141*^-18 I

We observe a small negative $\mathcal{O}(1\%)$ shift in the individual
coefficients. Note that the explicit specification of the running order
runAlpha -> 5 at Core creation is not mandatory since 5-loop running is the
default. We included it to be explicit.

#### Hadronic $\tau$ decay with $C$-scheme strong coupling

We now turn to Ref. [26] for which we reproduce the value quoted for
$\delta^{(0)}_{\mathrm{CI}}$ in Eq. (22). To obtain this number the 3-flavor
strong coupling, the related $\beta$-function coefficients and the
coefficients $c_{n,1}$ were converted to a class of schemes for the strong
coupling where the $\beta$-function adopts the following exact all order form

$\hat{\beta}(\hat{\alpha}_{s})=-2\hat{\alpha}_{s}\frac{\frac{\hat{\alpha}_{s}}{4\pi}\beta_{0}}{1-\frac{\hat{\alpha}_{s}}{4\pi}\frac{\beta_{1}}{\beta_{0}}}\,=-2\hat{\alpha}_{s}\\!\sum_{i=0}\hat{\beta}_{i}\biggl{(}\frac{\hat{\alpha}_{s}}{4\pi}\biggr{)}^{\\!\\!i+1}\,,$
(8)

where $\hat{\beta}_{i}=\beta_{0}(\beta_{1}/\beta_{0})^{i}$. Within that class
of schemes one needs to fix a parameter $C$ to uniquely specify the strong
coupling. This parameter was set to $C=-1.246$ in Eq. (22) of Ref. [26] with
the argument that the unknown 5-loop coefficient can then be neglected
resulting in the value
$\delta^{(0)}_{\mathrm{CI}}(\hat{\alpha}(m_{\tau}),C=-1.246)=0.1840\pm
0.0062$. The quoted uncertainty refers to the size of the 4-loop correction
term. Solving Eq. (6) of the reference paper to convert the quoted value of
the $\overline{\mathrm{MS}}$ strong coupling
$\alpha_{s}^{(3)}(m_{\tau})=0.316$ to the scheme described above leads to
$\hat{\alpha}_{s}(m_{\tau},C=-1.246)=0.477$. The transformed coefficients
$c_{n,1}$ can be extracted from Eq. (12) therein. REvolver allows the order-
by-order specification of user-defined $\beta$-functions, so we expand
$\hat{\beta}(\hat{\alpha}_{s})$ up to $\mathcal{O}(\hat{\alpha}_{s}^{11})$.

We define the relevant input values related to $\hat{\alpha}_{s}$, the
coefficients $c_{n,1}$ and $\hat{\beta}_{i}$ for $n_{f}=3$ with

      In[]:= {nfa, aTau, mTau} = {3, 0.477, 1.77686};
             {c11, c21, c31, c41} =
               {1, 1.640 + 2.25 c, 7.682 + 11.38 c + 5.063 c^2,
                 61.06 + 72.08 c + 47.4 c^2 + 11.39 c^3} /.
                c -> -1.246;
             betaHat = 9. * (64/9)^Range[0, 9];

and the related Core

      In[]:= CoreCreate["Hat", {nfa, aTau, mTau}, betaHat]

Reusing the function `Ja[n_, core_]` as defined in the previous example we
obtain

      In[]:= Ja[1, "Hat"]*c11 + Ja[2, "Hat"]*c21 +
              Ja[3, "Hat"]*c31 + Ja[4, "Hat"]*c41
      Out[]= 0.18403340012158337 - 1.0177988675705949*^-17 I

for $\delta^{(0)}_{\mathrm{CI}}(\hat{\alpha}(m_{\tau}),C=-1.246)$ and

      In[]:= Ja[4, "Hat"]*c41
      Out[]= 0.006222359703784199 - 2.951274747776552*^-19 I

for the last correction term. Both the central value and the size of the last
correction term are in perfect agreement with Ref. [26].

#### Cauchy integral theorem for the strong coupling

It is worth mentioning that Cauchy’s integral formula can be utilized to check
the numerical quality of the coupling evolution routine implemented in
REvolver. For example, valuating $\alpha_{s}^{(3)}(2\,\mathrm{GeV})$ directly
as well as by employing Cauchy’s integral formula with a radius of $1$ GeV
gives equivalent results up to machine precision

      In[]:= mu0 = 2;
             aDirect = AlphaQCD["CI5", mu0]
             aResidue =
              1/(2 Pi) NIntegrate[AlphaQCD["CI5", mu0 + E^(I phi)],
                {phi, 0, 2 Pi}, PrecisionGoal -> 10]
             aResidue - aDirect
      Out[]= 0.3169005366613899
      Out[]= 0.31690053666139 + 1.1043592643970545*^-17 I
      Out[]= 1.1102230246251565*^-16 + 1.1043592643970545*^-17 I

For all practical purposes the solution for the strong coupling evolution
provided by REvolver based on a given QCD $\beta$-function can be considered
as exact.

###  7.2 Cores with one Massive Quark

In the previous sample applications the impact of flavor thresholds was not
considered. In the following we discuss examples where threshold effects
associated to one massive quark are accounted for.

####  7.2.1 Strong coupling with a flavor threshold

#### Strong coupling from multijet events

We consider Ref. [27], where the strong coupling value was determined by the
ATLAS collaboration from transverse energy-energy correlations in multijet
events based on $8$ TeV LHC data. The measurements were made for different
values of the transverse momentum sum $H_{\mathrm{T2}}$ of the two leading
jets. All values of $H_{\mathrm{T2}}$ considered in that analysis were much
larger than the top quark mass, consequently the associated respective strong
coupling values $\alpha_{s}^{(6)}(Q)$ are defined in the 6-flavor scheme. In
Tabs. 2 and 3 of Ref. [27] values for $\alpha_{s}^{(6)}(Q)$ and the associated
results for $\alpha_{s}^{(5)}(m_{Z})$ are quoted. For the strong coupling
running an approximate analytic 2-loop solution of the evolution equation was
used accounting for continuous matching at the top quark mass (which is
correct for 1-loop matching when the matching scale is at the top quark mass).
In the following, we focus on the associated results
$\alpha_{s}^{(5)}(m_{Z})=0.1186^{+0.0090}_{-0.0047}$ and
$\alpha_{s}^{(6)}(810\,\mathrm{GeV})=0.0907^{+0.0052}_{-0.0026}$, where for
simplicity, we added the respective upper and lower uncertainties in
quadrature. We reproduce the asymmetric uncertainties for
$\alpha_{s}^{(6)}(810\,\mathrm{GeV})$ for the given $\alpha_{s}^{(5)}(m_{Z})$
range and start by defining the relevant parameters for the 5-flavor coupling
$\alpha_{s}^{(5)}(m_{Z})$, where we set the top quark standard running mass to
$\overline{m}_{t}=163$ GeV

      In[]:= {amZCentral, Q, mtmt} = {0.1186, 810, 163};
             {amZMax, amZMin} = amZCentral + {0.0090, -0.0047};

Next we create the Cores for our evaluation (including a `CoreDeleteAll[]` to
remove older cores)

      In[]:= CoreDeleteAll[]
             CoreCreate["central2", 6, {5, amZCentral, mZdef},
              {{6, mtmt, mtmt}}, runAlpha -> 2]
             CoreCreate["max2", 6, {5, amZMax, mZdef},
              {{6, mtmt, mtmt}}, runAlpha -> 2]
             CoreCreate["min2", 6, {5, amZMin, mZdef},
              {{6, mtmt, mtmt}}, runAlpha -> 2]

To obtain the $6$-flavor strong coupling values at $Q=810$ GeV including the
error range we simply evaluate

      In[]:= aQCentral2 = AlphaQCD["central2", Q, 6]
             AlphaQCD[#, Q, 6] & /@ {"min2", "max2"} - aQCentral2
      Out[]= 0.09079931707610696
      Out[]= {-0.002757013433150532, 0.005129916277931454}

agreeing very well with the result given in Ref. [27].

We investigate this setup with two-loop evolution further by determining the
uncertainty related to varying the matching scale accounting for the top
threshold corrections at one loop. This is easily done using REvolver by
creating a set of Cores with a range of f-parameters. To this end we first
define a table containing $20$ logarithmically scaled f-parameters in the
range $[1/2,2]$, which corresponds to matching scales between one half and
twice the standard running top mass. Subsequently the table is used to create
the corresponding set of Cores

      In[]:= list2 = 2^Range[-1, 1, 2/19];
      In[]:= CoreCreate[ToString[#] <> "run2", 6,
                 {5, amZCentral, mZdef}, {{6, mtmt, mtmt}},
                 runAlpha -> 2, orderAlpha -> 1,
                 fMatch -> {#}] & /@ list2;

Finally, we can determine a list of the corresponding strong coupling values
$\alpha_{s}^{(6)}(810\,\mathrm{GeV})$ and compute the central value as well as
the error range

      In[]:= aQList2 = AlphaQCD[ToString[#] <> "run2", Q, 6] & /@
                list2;
             (Max[aQList2] + {1, -1} * Min[aQList2])/2
      Out[]= {0.09081555469961156, 0.00009936201139525841}

The perturbative uncertainties associated to the threshold corrections are
about $3\%$ of the quoted experimental error.

Employing 5-loop instead of 2-loop running for the strong coupling and 4-loop
matching corrections we obtain

      In[]:= CoreCreate["central5", 6, {5, amZCentral, mZdef},
              {{6, mtmt, mtmt}}, runAlpha -> 5]
             CoreCreate["max5", 6, {5, amZMax, mZdef},
              {{6, mtmt, mtmt}}, runAlpha -> 5]
             CoreCreate["min5", 6, {5, amZMin, mZdef},
              {{6, mtmt, mtmt}}, runAlpha -> 5]
      In[]:= aQCentral5 = AlphaQCD["central5", Q, 6]
             AlphaQCD[#, Q, 6] & /@ {"min5", "max5"} - aQCentral5
      Out[]= 0.09078701609518454
      Out[]= {-0.00275523573924355, 0.0051259314534172346}

for $\alpha_{s}^{(6)}(810\,\mathrm{GeV})$ and its upper and lower uncertainty
and

      In[]:= CoreCreate[ToString[#] <> "run5", 6,
                 {5, amZCentral, mZdef}, {{6, mtmt, mtmt}},
                 runAlpha -> 5, fMatch -> {#}] & /@ list2;
      In[]:= aQList5 = AlphaQCD[ToString[#] <> "run5", Q, 6] & /@
                list2;
             (Max[aQList5] + {1, -1} * Min[aQList5])/2
      Out[]= {0.09078682756270738, 1.9445996186917558*^-7}

for the central value and error estimate derived from varying the matching
scale. Within the experimental uncertainties as quoted in Ref. [27] using
2-loop evolution and continuous matching is perfectly adequate.

#### Strong coupling from inclusive jet cross sections

We return to the analysis of Ref. [24] and investigate the impact of a top
quark threshold on the values given there, staying with 2-loop running as
employed in that reference. In the following we create Cores with the given
range of values of $\alpha_{s}^{(5)}(m_{Z})=0.1162^{+0.0070}_{-0.0062}$ and
determine the strong coupling $\alpha_{s}^{(6)}(1508.04\,\mathrm{GeV})$ in the
6-flavor scheme instead of the 5-flavor coupling
$\alpha_{s}^{(5)}(1508.04\,\mathrm{GeV})$ determined in Ref. [24]. After
defining the parameters

      In[]:= {amZCentral, Q, mtmt} = {0.1162, 1508.04, 163};
             {amZMax, amZMin} = amZCentral + {0.007, -0.0062};

and creating the respective Cores, accounting for the standard running top
quark mass $\overline{m}_{t}=163$ GeV and default 4-loop matching

      In[]:= CoreCreate["central2", 6, {5, amZCentral, mZdef},
               {{6, mtmt, mtmt}}, runAlpha -> 2];
             CoreCreate["min2", 6, {5, amZMin, mZdef},
               {{6, mtmt, mtmt}}, runAlpha -> 2];
             CoreCreate["max2", 6, {5, amZMax, mZdef},
               {{6, mtmt, mtmt}}, runAlpha -> 2];

we obtain

      In[]:= aQCentral2 = AlphaQCD["central2", Q, 6]
             (AlphaQCD[#, Q, 6] & /@ {"min2","max2"}) - aQCentral2
      Out[]= 0.08406136347372324
      Out[]= {-0.003262229612787812, 0.0035626775539831235}

for the central value and the upper and lower uncertainties, respectively.
Comparing these values to
$\alpha_{s}^{(5)}(1508.04\,\mathrm{GeV})=0.0822^{+0.0034}_{-0.0031}$, quoted
in Ref. [24], we observe a positive shift of about $0.002$ in the central
value. This already amounts to about $60\%$ of the given experimental error,
illustrating that flavor-thresholds effects can lead to significant changes.

####  7.2.2 Asymptotic pole and low-scale MSR mass

To demonstrate the REvolver functionalities related to mass conversions
accounting for flavor threshold effects, we investigate how much the
asymptotic top quark pole mass $m_{t}^{\rm pole}$ as well as the top MSR mass
at $2$ GeV $m_{t}^{\rm MSR}(2\,\mbox{GeV})$ would change if the strong
coupling value quoted in Ref. [24]
$\alpha_{s}(1508.04\,\mbox{GeV})=0.0822^{+0.0034}_{-0.0031}$ would be
interpreted as a $6$-flavor compared to a $5$-flavor result. After defining
the relevant constants and Cores (we keep using the values of Q and mtmt
already defined in the previous section)

      In[]:= aQ = 0.0822;
             CoreCreate["5", 6, {5, aQ, Q}, {{6, mtmt, mtmt}}];
             CoreCreate["6", 6, {6, aQ, Q}, {{6, mtmt, mtmt}}];

where we have employed default highest-order precision for coupling and mass
evolution and flavor threshold matching, we first extract
$m_{t}^{\mathrm{MSR}}(2\,\mathrm{GeV})$ from both Cores and then determine
their difference:

      In[]:= mMSR25 = MassMS["5", 6, 2]
             mMSR26 = MassMS["6", 6, 2]
             mMSR25 - mMSR26
      Out[]= 172.45858180514585
      Out[]= 172.10074770514723
      Out[]= 0.3578340999986267

For the asymptotic pole masses we obtain

      In[]:= mPole5 = MassPole["5", 6, mtmt]
             mPole6 = MassPole["6", 6, mtmt]
             mPole5 - mPole6
      Out[]= 172.8732843267365
      Out[]= 172.4868549739421
      Out[]= 0.3864293527944085

where we adopted $\overline{m}_{t}=163$ GeV as the conversion scale for the
asymptotic pole mass determination. The differences we obtain amount to $358$
MeV and $386$ MeV, both of which exceed the uncertainty of the current world
average for direct top mass measurements [23].

####  7.2.3 Bottom and charm quark short-distance masses

#### Bottom $\overline{\mathrm{MS}}$ mass at high scales

In Ref. [28] the partonic cross section for Higgs production via bottom quark
fusion was presented. This involves the bottom running mass at the Higgs scale
$m_{H}=125.09$ GeV, quoted to be
$m_{b}^{(5)}(m_{H})\equiv\overline{m}_{b}^{(5)}(m_{H})=2.79$ GeV using 4-loop
$\overline{\mathrm{MS}}$-mass running, $\alpha_{s}^{(5)}(m_{Z})=0.118$ and the
bottom quark standard running mass $\overline{m}_{b}=4.18$ GeV, see Table 2
therein. We can reproduce this relation in REvolver by defining the relevant
parameters

      In[]:= {mH, amZ, mbmb} = {125.09, 0.118, 4.18};

and creating a corresponding Core

      In[]:= CoreDeleteAll[]
             CoreCreate["b", 5, {5, amZ, mZdef}, {5, mbmb, mbmb},
                runMSbar -> 4]

The value of $m_{b}^{(5)}(m_{H})$ can now be extracted with

      In[]:= MassMS["b", 5, mH]
      Out[]= 2.78854676339097

showing perfect agreement with Ref. [28].

#### Bottom quark PS and 1S masses

In Ref. [29] the conversion between the $\overline{\mathrm{MS}}$ and various
low-scale short-distance masses was carried out with 4-loop fixed-order
formulae. To illustrate the functionalities of REvolver for mass conversions
we now consider the values of the PS and 1S bottom quark masses given in Table
II of that article. We start defining the input parameters
$\alpha_{s}^{(5)}(m_{Z})=0.1185$, $\overline{m}_{b}=4.163$ GeV and
$\mu_{f}=2.0$ GeV for the strong coupling, the standard running bottom mass
and the bottom PS mass renormalization scale, respectively, as specified in
Ref. [29],

      In[]:= {amZ, mbmb, mufB} = {0.1185, 4.163, 2.0};

and the related Core

      In[]:= CoreCreate["b", 5, {5, amZ, mZdef},
              {{5, mbmb, mbmb}}, fMatch -> {2.0}]

Following Ref. [29] we set the bottom quark matching scale to twice the
standard running bottom mass. Comparing the 4-loop value for the bottom PS
mass $m_{b}^{\rm PS}(2\,\mbox{GeV})=4.483$ GeV quoted in Ref. [29] and
obtained by REvolver we obtain

      In[]:= mbPSRef = 4.483;
             mbPSREvo = MassPS["b", 5, mufB, 5, mbmb, mbmb]
             mbPSREvo - mbPSRef
      Out[]= 4.484037803646092
      Out[]= 0.0010378036460920725

where, following Ref. [29], the conversion is carried out using the fixed-
order relation between the standard running and PS masses using the standard
running mass $\overline{m}_{b}$ as the renormalization scale for the strong
coupling.

The small difference of about $1$ MeV results from the fact that in REvolver
all conversions between mass schemes are based on formulae starting from the
value of the running mass, while in Ref. [29] the conversion was obtained the
other way around, i.e. the running mass was computed starting from a value in
the PS scheme. The difference is naturally covered by the perturbative
uncertainty and not relevant for practical purposes.

In the previous example we have used REvolver to determine the PS from the
standard running mass $\overline{m}_{b}$. The conversion in the opposite
direction can be achieved using the REvolver functionality to add a heavier
mass to an existing Core, see Sec. 6.3. To determine the standard running mass
we first create a Core containing $4$ massless flavors

      In[]:= anf4mu3 = AlphaQCD["b", 3.0, 4];
             CoreCreate["O4", {4, anf4mu3, 3.0}]

which requires the 4-flavor strong coupling
$\alpha_{s}^{(4)}(3\,\mathrm{GeV})$ as an input, which here we obtain from the
previously created Core named b.