aainstitutetext: Theoretical Physics Department, CERN, CH-1211 Geneva 23,
Switzerlandbbinstitutetext: Department of Physics, Brookhaven National
Laboratory, Upton, N.Y., 11973, U.S.A

# An analysis of Bayesian estimates for missing higher orders in perturbative
calculations

Claude Duhr a Alexander Huss a Aleksas Mazeliauskas b Robert Szafron
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

With current high precision collider data, the reliable estimation of
theoretical uncertainties due to missing higher orders (MHOs) in perturbation
theory has become a pressing issue for collider phenomenology. Traditionally,
the size of the MHOs is estimated through scale variation, a simple but ad hoc
method without probabilistic interpretation. Bayesian approaches provide a
compelling alternative to estimate the size of the MHOs, but it is not clear
how to interpret the perturbative scales, like the factorisation and
renormalisation scales, in a Bayesian framework. Recently, it was proposed
that the scales can be incorporated as hidden parameters into a Bayesian
model. In this paper, we thoroughly scrutinise Bayesian approaches to MHO
estimation and systematically study the performance of different models on an
extensive set of high-order calculations. We extend the framework in two
significant ways. First, we define a new model that allows for asymmetric
probability distributions. Second, we introduce a prescription to incorporate
information on perturbative scales without interpreting them as hidden model
parameters. We clarify how the two scale prescriptions bias the result towards
specific scale choice, and we discuss and compare different Bayesian MHO
estimates among themselves and to the traditional scale variation approach.
Finally, we provide a practical prescription of how existing perturbative
results at the standard scale variation points can be converted to 68%/95%
credibility intervals in the Bayesian approach using the new public code MiHO.

††preprint: CERN-TH-2021-058

##  1 Introduction

Perturbative calculations are currently the only way to obtain Standard Model
(SM) predictions for collider experiments. The backbone of almost all
calculations for hadron collider observables is the QCD factorisation theorem,
which relates an observable $\Sigma$ in proton-proton collisions (e.g., a
total cross-section or differential distribution) to the corresponding
partonic quantity,

$\Sigma=\sum_{i,j}\int_{0}^{1}dx_{1}\,dx_{2}\,f_{i}(x_{1},\mu_{F})\,f_{j}(x_{2},\mu_{F})\,\hat{\Sigma}_{ij}(x_{1},x_{2},\mu_{F})+\mathcal{O}\Big{(}\frac{\Lambda_{\text{QCD}}}{Q}\Big{)}\,,$
(1)

where $f_{i}(x,\mu_{F})$ denotes the parton distribution function (PDF) for
parton $i$ and $\hat{\Sigma}_{ij}$ denotes the observable computed on a
collision of the partons $i$ and $j$. Here $\mu_{F}$ is the factorisation
scale, $Q$ denotes a characteristic hard scale of the process, and
$\Lambda_{\text{QCD}}\lesssim 1\text{ GeV}$ is the scale at which QCD becomes
strongly coupled.

The partonic coefficients $\hat{\Sigma}_{ij}$ are computable in perturbation
theory,

$\hat{\Sigma}_{ij}(x_{1},x_{2},\mu_{F})\simeq\sum_{k=0}^{\infty}\alpha_{s}(\mu_{R})^{k_{0}+k}\,\hat{\Sigma}_{ij}^{(k)}(x_{1},x_{2},\mu_{F},\mu_{R})\,,$
(2)

where $k_{0}\geq 0$ is a non-negative integer, $\alpha_{s}(\mu_{R})$ denotes
the renormalised strong coupling constant evaluated at the renormalisation
scale $\mu_{R}$, and $\hat{\Sigma}_{ij}^{(k)}$ is the (next-to-)k-leading
order (NkLO) contribution to the partonic coefficient.111The ‘$\simeq$’-sign
in eq. (2) indicates that the perturbative expansion on the right-hand side is
in general an asymptotic series with zero radius of convergence. For most
phenomenological predictions only very few terms in the perturbative expansion
are known (typically not more than the first 3 or 4 terms), and one does not
have to worry about the non-convergent nature of the series. One can combine
the QCD factorisation theorem in eq. (1) with the perturbative expansion of
the partonic coefficients in eq. (2) to obtain the NnLO approximation for the
observable $\Sigma$ within (fixed-order) perturbation theory

$\Sigma\simeq\Sigma_{n}(\mu_{F},\mu_{R}):=\sum_{k=0}^{n}\Sigma^{(k)}(\mu_{F},\mu_{R})\,.$
(3)

Computations are typically performed in the $\overline{\rm MS}$
renormalisation scheme. Although $\Sigma$ is formally independent of the
choice of the factorisation and renormalisation scales, the truncation of the
perturbative series introduces a residual scale dependence into $\Sigma_{n}$.

Since a perturbative result is always an approximation of the exact value
$\Sigma$, it is important to have robust methods to assess the quality of the
approximation and to quantify the uncertainty attached to it. Three of the
main sources of uncertainty are coming from the values for the strong coupling
constant and other input parameters (we collectively refer to them as the
parametric uncertainty), the PDFs222Note that the PDF uncertainty can itself
be seen as a parametric uncertainty, just like the coupling constants or the
masses. Indeed, there is no strong conceptual difference between the PDF and,
e.g., $\alpha_{s}$ uncertainties, and often they are treated on the same
footing. For our discussion in later sections, we find it convenient to think
of these two sources separately. This does, however, not impact any of our
results, as our primary focus are MHO uncertainties. and the truncation of the
perturbative series. The parametric and PDFs uncertainties are related to
their extraction from experimental data and the fitting methodology used.
Consequently, the PDF and parametric uncertainties are not directly related to
the observable $\Sigma$. They are computed using observable-independent
recipes, cf., e.g., refs. Butterworth:2015oua ; Alekhin:2011sk ; Botje:2011sn
(though care is needed as there may be correlations between $\Sigma$ and the
data used to extract the PDFs and $\alpha_{s}$). The uncertainty resulting
from the truncation of the perturbative series quantifies the effect of the
missing higher orders (MHOs) in perturbation theory. Its estimation is one of
the main topics of this paper. Under the assumption that the PDF, parametric,
and MHO uncertainties are uncorrelated (see Harland-Lang:2018bxd ;
Ball:2018twp ; AbdulKhalek:2019ihb ; Ball:2021icz  for discussion of
correlation between PDF and MHO uncertainties), they can be added in
quadrature to obtain the uncertainty on the NnLO prediction.333For a given
observable there may be other uncertainties one needs to take into account,
e.g., the choice of the renormalisation schemes used for the input parameters,
the truncation of the heavy mass expansion etc. How to estimate these
uncertainties and combine them is often a matter of debate for specific
processes and observables. We do not discuss these issues in this paper as we
intend to keep the discussion general.

Currently, there is no reliable way to estimate the size of the MHOs (other
than computing the next order in perturbation theory). Since the effect of
varying the scales at NnLO is an Nn+1LO effect, one conventionally estimates
the size of the MHO terms by studying the variation of
$\Sigma_{n}(\mu_{F},\mu_{R})$ as the factorisation and renormalisation scales
are varied in a certain range around a central scale $\mu_{F/R,\text{cent}}$
corresponding to a characteristic hard scale of the process. A popular choice
is the so-called _7-point variation_ , where the MHO uncertainty at NnLO is
estimated by the interval

$\Big{[}\min_{(k_{F},k_{R})\in
S_{7}}\\!\\!\\!\\!\\!\Sigma_{n}(k_{F}\,\mu_{F,\text{cent}},k_{R}\,\mu_{R,\text{cent}}),\max_{(k_{F},k_{R})\in
S_{7}}\\!\\!\\!\\!\\!\Sigma_{n}(k_{F}\,\mu_{F,\text{cent}},k_{R}\,\mu_{R,\text{cent}})\Big{]}\,.$
(4)

$(k_{F},k_{R})$ take values in the following set:

$S_{7}=\\{(1/2,1/2),(1/2,1),(1,1/2),(1,1),(1,2),(2,1),(2,2)\\}\,.$ (5)

This differs from the _9-point variation_ where $k_{F},k_{R}\in\\{1/2,1,2\\}$
by the exclusion of large ratios of the scales, i.e., $k_{F}/k_{R}=4$ or
$1/4$. We stress that the procedure of estimating the MHO uncertainty from
scale variation is completely ad hoc. There is no reason other than experience
from NLO and NNLO calculations why the obtained uncertainty interval at NnLO
should capture the actual size of the Nn+1LO corrections $\Sigma^{(n+1)}$ or
the true value of $\Sigma$. Indeed, the higher-order terms probed by varying
the scales are universal and predictable from renormalisation group methods.
They do not take into account the size of the genuine higher-order
corrections, which are process-dependent and independent of the perturbative
scales. With advances in precision calculations over the last 5 years, NNLO
results are becoming routinely available for many LHC processes, and even N3LO
corrections to colourless final states have been recently computed. It is,
therefore, time to critically assess the MHO uncertainties and to look for
alternative prescriptions. A number of attempts have been made in the past to
estimate or predict the truncation errors and MHOs, see, e.g., refs
Ellis:1995jv ; Gardi:1996iq ; Ellis:1996zn ; Weniger:1997zz ; Ellis:1997sb ;
Karliner:1998ge ; Caprini:1998wg ; Cvetic:1998zm ; Elias:2000iw ;
Caprini:2000js ; Jentschura:2000iw ; Beneke:2008ad ; David:2013gaa ;
Furnstahl:2015rha ; Boito:2016pwf ; Jamin:2016ihy ; Melendez:2017phj ;
Caprini:2019kwp ; Melendez:2019izc ; Costin:2020hwg ; Boito:2020hvu ;
Drischler:2020hwi . Various prescriptions for reducing or eliminating scale
ambiguity have also been discussed broadly in the literature, e.g., refs.
Politzer:1981vc ; Pennington:1981sq ; Stevenson:1982wn ; Grunberg:1982fw ;
Brodsky:1982gc ; Cvetic:1997ca ; Maxwell:2000mm ; Brodsky:2011ig ;
Kataev:2014jba ; Kataev:2016aib ; Czakon:2016dgf ; Garkusha:2018mua ;
Wu:2019mky ; Boito:2018rwt ; Brodsky:2019tmf ; Chawdhry:2019uuv ;
Chishtie:2020cen ; DiGiustino:2020fbk .

Here we focus on a promising, widely applicable approach to estimate the MHO
contributions using Bayesian inference on the known fixed-order computations.
This idea was pioneered by Cacciari and Houdeau (CH) in ref. Cacciari:2011ze ,
and developed further in refs. Forte:2013mda ; Bagnaschi:2014wea ;
Bonvini:2020xeo . Bayesian inference is a powerful method to construct
probability distributions in which Bayes’ theorem is used to update the
probability as new information becomes available iteratively. The procedure
can be summarised as follows: one starts from some model and the prior
distribution of the model parameters. These inputs are purely subjective. The
next step is a construction of a posterior distribution for model parameters
based on available data. This is the inference step, where one updates prior
belief as new data become available. Finally, one constructs the posterior
distribution for the next data point by marginalising over the posterior
distribution for the model parameters. The main idea of the CH-method it to
use Bayesian inference on the known perturbative orders to model the size of
the MHOs. The original CH-model can be applied only after one specifies the
choices for the perturbative scales. In this paper, we follow up on the recent
work in ref. Bonvini:2020xeo , where the CH-approach was extended by the
postulate that unphysical scales should be treated on par with other hidden
model parameters. According to the lore of Bayesian inference, one expects
that the posterior distribution for the perturbative scales peaks at the
‘optimal value’ for the scale choice, and the dependence on the scales is
eliminated by marginalisation. However, this is only one example of possible
prescriptions, as there is no reason, physical or mathematical, to interpret
the perturbative scales as hidden model parameters. One of the main objectives
of this paper is to thoroughly scrutinise the consequences of this postulate
and to explore alternative approaches.

This paper is organised as follows: In section 2 we discuss two prescriptions
of how to incorporate scale dependence into a model for Bayesian inference of
fixed-order calculations, and we investigate implicit biases inherent to them.
In section 3 we recap the mathematical formalism behind Bayesian inference at
fixed scale and introduce a model to account for asymmetric probability
distributions. In section 4 we study extensively the model performance for
representative perturbative calculations without scale dependence or at fixed
scale choice. Finally, in section 5 we study the performance of the two
prescriptions to incorporate scale information for an extensive set of
collider processes known at NNLO and/or N3LO. We conclude with a discussion in
section 6.

##  2 Bayesian methods for observables in quantum field theory

###  2.1 Incorporating scale information into Bayesian inference

In this section we present a way to estimate the MHO uncertainty for an
observable $\Sigma$ via Bayesian methods. To be concrete, let us assume that
we have computed the first $n+1$ orders $\Sigma^{(k)}$, $0\leq k\leq n$, in
the perturbative expansion of $\Sigma$. Our goal is to construct a probability
distribution ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ for $\Sigma$ given the
sequence of perturbative coefficients
${\bm{\Sigma}}_{n}:=(\Sigma^{(0)},\Sigma^{(1)},\ldots,\Sigma^{(n)})$. This
distribution can be used to construct credibility intervals (CIs) for
$\Sigma$, which may serve as the uncertainty estimate. In particular, we
define the $x$% CI as
$\text{CI}_{x}=[\Sigma_{x}^{\text{low}},\Sigma_{x}^{\text{up}}]$ where

$\int_{-\infty}^{\Sigma_{x}^{\text{low}}}\\!d\Sigma\,{\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})=\int^{\infty}_{\Sigma_{x}^{\text{up}}}\\!d\Sigma\,{\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})=\frac{1-\frac{x}{100}}{2}\,.$
(6)

We will mainly consider $x\%\approx 68\%$ and $x\%\approx 95\%$ CIs.

The original approach of ref. Cacciari:2011ze , known as the _Cacciari-
Houdeau_ (CH) method, allows one to construct a probability
${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ in the case where the
${\bm{\Sigma}}_{n}$ are _numbers_. In perturbation theory, however, the
perturbative coefficients ${\bm{\Sigma}}_{n}$ are _functions_ of the
perturbative scales, which are arbitrary quantities introduced by perturbation
theory without any physical interpretation. Clearly,
${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ should be independent of the scale
choice, but the value of the NnLO prediction may depend strongly on this
choice, and there is no reason to prefer one scale over another (as long as
the scale does not differ too much from the hard scale of the process).

Reference Bonvini:2020xeo  introduced a possible way of incorporating
perturbative scales into a Bayesian approach to estimating MHOs (see also ref.
Bagnaschi:2014wea  for an alternative approach). The approach of ref.
Bonvini:2020xeo 444Reference Bonvini:2020xeo  excludes the factorisation scale
from this discussion. We will come back to this point in section 2.1.3. to
include the information on the perturbative scales into a Bayesian model is
not the only one. Broadly speaking, we identify two different prescriptions of
how one can incorporate information on perturbative scales into the
construction of a probability distribution:

Prescription 1: Scale-marginalisation (sm). Despite the perturbative scales
not being physical, one assigns a Bayesian degree of belief to them. In this
way the scales become (hidden) parameters of a Bayesian model. It is then
possible to construct ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ through
Bayesian inference and marginalisation over the perturbative scales. This is
the approach advocated in ref. Bonvini:2020xeo .

Prescription 2: Scale-averaging (sa). Since the perturbative scales are devoid
of any physical meaning, they cannot be treated as hidden parameters of a
Bayesian model, i.e., they cannot be treated as random variables with a
statistical interpretation. Instead, the result of a perturbative calculation
through NnLO should be interpreted as a family of independent predictions for
the true value $\Sigma$ parametrised by the scales. The probability
distribution ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ is then defined by
adding coherently the probabilities computed for each family member.

As we will show later, the two prescriptions of obtaining scale-independent
probabilities differ not only in their philosophy, but in practice they can
produce different results due to implicit biases towards some scales over
others. In the next section we present the mathematical framework to construct
the probability distribution ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$. In
order to keep the notations concise, we focus first on the simplified
situation where the perturbative coefficients ${\bm{\Sigma}}_{n}$ depend on a
single perturbative scale $\mu$ (which we may think of as the renormalisation
scale for now).

####  2.1.1 Prescription 1: scale-marginalisation

We start by reviewing the _scale-marginalisation (sm)_ prescription of ref.
Bonvini:2020xeo . Since $\mu$ is a parameter of the model, the probability
distribution ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ is obtained by
marginalising over $\mu$;

$\begin{split}{\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})\,&=\int
d\mu\,{\mathcal{P}}(\Sigma,\mu|{\bm{\Sigma}}_{n})=\int
d\mu\,{\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n},\mu)\,{\mathcal{P}}(\mu|{\bm{\Sigma}}_{n})\,.\end{split}$
(7)

Let us interpret the probabilities in the integrand. For the first factor, we
apply Bayes’ formula to write

$\begin{split}{\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n},\mu)&\,=\frac{{\mathcal{P}}(\Sigma,{\bm{\Sigma}}_{n}|\mu)}{{\mathcal{P}}({\bm{\Sigma}}_{n}|\mu)}=\frac{P(\Sigma,{\bm{\Sigma}}_{n}(\mu))}{P({\bm{\Sigma}}_{n}(\mu))}=P(\Sigma|{\bm{\Sigma}}_{n}(\mu))\,.\end{split}$
(8)

The interpretation of this equation is as follows: Since $\mu$ is a parameter
of the model, whenever we evaluate a conditional probability _for a fixed
value of $\mu$_, we can replace the perturbative coefficients
${\bm{\Sigma}}_{n}$ (which are functions) by their values
${\bm{\Sigma}}_{n}(\mu)$ for fixed $\mu$ (which are numbers). In the following
we always indicate by ${\mathcal{P}}$ a probability distribution which takes
are arguments the perturbative coefficients ${\bm{\Sigma}}_{n}$ (which are
functions), while $P$ denotes a probability distribution on perturbative
coefficients computed at a fixed scale $\mu$. For now we assume that $P$ is
given. In section 3 we discuss a framework how a distribution $P$ can be
constructed.

The second factor in the integrand of eq. (7) is the posterior distribution
for $\mu$, and we have

$\begin{split}{\mathcal{P}}(\mu|{\bm{\Sigma}}_{n})&\,=\frac{{\mathcal{P}}(\mu,{\bm{\Sigma}}_{n})}{{\mathcal{P}}({\bm{\Sigma}}_{n})}=\frac{{\mathcal{P}}({\bm{\Sigma}}_{n}|\mu)\,P_{0}(\mu)}{\int
d\mu^{\prime}\,{\mathcal{P}}({\bm{\Sigma}}_{n}|\mu^{\prime})\,P_{0}(\mu^{\prime})}=\frac{P({\bm{\Sigma}}_{n}(\mu))\,P_{0}(\mu)}{\int
d\mu^{\prime}\,P({\bm{\Sigma}}_{n}(\mu^{\prime}))\,P_{0}(\mu^{\prime})}\,,\end{split}$
(9)

where $P_{0}(\mu)$ is an (arbitrary) prior distribution on the perturbative
scale $\mu$. In ref. Bonvini:2020xeo  it was argued that a physically well-
motivated choice is a flat distribution in the logarithm of the scale,

$P_{0}(\mu):=\frac{1}{2\mu\,\log F}\,\Theta\left(\log
F-\left|\log\frac{\mu}{\mu_{0}}\right|\right)\,,$ (10)

where $F>1$ is a constant and $\Theta$ denotes the Heaviside step function;

$\Theta(x)=\left\\{\begin{array}[]{ll}1\,,\,\quad\text{ if }x>0\,,\\\
0\,,\,\quad\text{ otherwise }\,.\end{array}\right.$ (11)

The central scale $\mu_{0}$ is in principle arbitrary, but assumed to be of
the order of the hard scale of the process. Putting everything together, we
see that ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ can be expressed in terms
of the probability distribution $P$ and the prior $P_{0}$ on the scale:

$\displaystyle{\mathcal{P}}$
$\displaystyle(\Sigma|{\bm{\Sigma}}_{n})\equiv{\mathcal{P}}_{\text{sm}}(\Sigma|{\bm{\Sigma}}_{n}):=\frac{\int
d\mu\,P(\Sigma|{\bm{\Sigma}}_{n}(\mu))\,P({\bm{\Sigma}}_{n}(\mu))\,P_{0}(\mu)}{\int
d\mu^{\prime}\,P({\bm{\Sigma}}_{n}(\mu^{\prime}))\,P_{0}(\mu^{\prime})}\,.$
(12)

Equation (12) allows one to relate
${\mathcal{P}}_{\text{sm}}(\Sigma|{\bm{\Sigma}}_{n})$ to two other probability
distributions (in addition to the prior on $\mu$). The first one is the
probability distribution $P(\Sigma|{\bm{\Sigma}}_{n}(\mu))$ for the (true)
value of the observable $\Sigma$ given the $n+1$ first perturbative orders
calculated at the scale $\mu$. The second one is the probability distribution
$P({\bm{\Sigma}}_{n}(\mu))$, which is related through Bayes’ theorem to
$P(\Sigma^{(n)}(\mu)|{\bm{\Sigma}}_{n-1}(\mu))$. The latter computes the
probability distribution for the NnLO corrections $\Sigma^{(n)}(\mu)$
evaluated at the scale $\mu$, given the $n$ first perturbative orders
evaluated at that scale. The two distributions
$P(\Sigma|{\bm{\Sigma}}_{n}(\mu))$ and
$P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}_{n}(\mu))$ are the answers to two
different mathematical questions, and there is a priori no reason why they
should be related. However, if perturbation theory is applicable, then the
true value of $\Sigma$ should be well approximated by $\Sigma_{n+1}(\mu)$ (at
least for reasonable choices of $\mu$). We therefore expect that also the
distribution $P(\Sigma|{\bm{\Sigma}}_{n}(\mu))$ is well approximated by
$P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}_{n}(\mu))$. Substituting
$\Sigma^{(n+1)}(\mu)\approx\Sigma-\Sigma_{n}(\mu)$ in the latter, we expect
that, if perturbation theory is applicable, then

$P(\Sigma|{\bm{\Sigma}}_{n}(\mu))\approx
P(\Sigma-\Sigma_{n}(\mu)|{\bm{\Sigma}}_{n}(\mu))\,.$ (13)

For a more detailed derivation of this approximation in a Bayesian framework
we refer to ref. Bonvini:2020xeo . Using this approximation, we can cast eq.
(12) into the form:

$\displaystyle{\mathcal{P}}_{\text{sm}}$
$\displaystyle(\Sigma|{\bm{\Sigma}}_{n})\approx\frac{\int
d\mu\,P(\Sigma-\Sigma_{n}(\mu)|{\bm{\Sigma}}_{n}(\mu))\,P({\bm{\Sigma}}_{n}(\mu))\,P_{0}(\mu)}{\int
d\mu^{\prime}\,P({\bm{\Sigma}}_{n}(\mu^{\prime}))\,P_{0}(\mu^{\prime})}\,.$
(14)

We stress that, unlike indicated by the notation in eq. (13), the probability
distributions in the left and right-hand sides of eqs. (12) and (13) are
actually _distinct_ : $P(\Sigma|{\bm{\Sigma}}_{n}(\mu))$ is the conditional
distribution for the true value $\Sigma$ of the observable, while
$P(\Sigma-\Sigma_{n}(\mu)|{\bm{\Sigma}}_{n}(\mu))$ is the distribution for the
value of the Nn+1LO correction
$\Sigma^{(n+1)}(\mu)\approx\Sigma-\Sigma_{n}(\mu)$ computed at the scale
$\mu$. We will always assume that available NnLO results sufficiently well
approximate $\Sigma$ and that we are in the perturbative regime where eq. (13)
holds, so we do not distinguish the notation for the two distributions.

####  2.1.2 Prescription 2: scale-averaging

In this section we present an alternative way to include scale information
into the model which does not require the perturbative scale $\mu$ to be
promoted to a model parameter with a statistical or physical meaning. We call
this method the _scale-averaging (sa)_ prescription. The starting point is to
interpret $\mu$ as a parameter defining a family of perturbative predictions
which approximate $\Sigma$ at NnLO in perturbation theory. For each member of
this family, i.e., for each fixed value $\mu$ of the scale, we can compute the
probability $P(\Sigma|{\bm{\Sigma}}_{n}(\mu))$ (as before we assume again that
$P$ is given). We then combine this family of predictions into a single
distribution ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ as a weighted average
of probabilities

$\begin{split}{\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})\equiv{\mathcal{P}}_{\text{sa}}(\Sigma|{\bm{\Sigma}}_{n})&\,:=\int
d\mu\,w(\mu)\,P(\Sigma|{\bm{\Sigma}}_{n}(\mu))\,.\end{split}$ (15)

The form of the weighting function $w(\mu)$ is arbitrary, but it must be
chosen such that ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ satisfies all the
axioms of a normalised probability distribution. In particular, $w(\mu)$ must
itself be normalised:

$\int d\mu\,w(\mu)=1\,.$ (16)

Since one expects that the perturbative scale $\mu$ should be chosen close to
the characteristic hard scale of the process in order to ensure good
perturbative convergence, we postulate a window function around a fixed scale
$\mu_{0}\sim Q$:

$w(\mu):=\frac{1}{2\mu\,\log F}\,\Theta\left(\log
F-\left|\log\frac{\mu}{\mu_{0}}\right|\right)\,,$ (17)

where $F>1$ is a constant. We note that our choice for $w(\mu)$ is identical
to the prior on the scale in the scale marginalisation in eq. (10). The
interpretation, however, is different: in the previous section the scale $\mu$
is a model parameter with a statistical interpretation that leads to eq. (12).
In eq. (15), the scale has no statistical interpretation, and is simply a
parameter defining a family of distributions, which are then combined into a
single distribution with weighting function $w(\mu)$ according to eq. (15).
With the weight given in eq. (17), this corresponds to a uniform average in
the logarithm of the scale over the interval $\mu_{0}/F<\mu<F\mu_{0}$. Using
the approximation in eq. (13) we obtain the final expression

${\mathcal{P}}_{\text{sa}}(\Sigma|{\bm{\Sigma}}_{n})\approx\int
d\mu\,w(\mu)\,P(\Sigma-\Sigma_{n}(\mu)|{\bm{\Sigma}}_{n}(\mu))\,.$ (18)

####  2.1.3 Extension to multiple perturbative scales

In the case where the perturbative coefficients ${\bm{\Sigma}}_{n}$ are
independent of the unphysical perturbative scale $\mu$, i.e.,
${\bm{\Sigma}}_{n}(\mu)={\bm{\Sigma}}_{n}(\mu_{0})$ for some $\mu_{0}$, the
distributions ${\mathcal{P}}$ reduce to $P$:

${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})=P(\Sigma|{\bm{\Sigma}}_{n}(\mu_{0}))\,,\quad\text{if}\quad{\bm{\Sigma}}_{n}(\mu)={\bm{\Sigma}}_{n}(\mu_{0})\,.$
(19)

To extend the discussion to more scales $\mu_{i}$, $1\leq i\leq s$ we
introduce a multivariate prior $P_{0}(\mu_{1},\ldots,\mu_{s})$ or weighting
function $w(\mu_{1},\ldots,\mu_{s})$, and we integrate over each scale
$\mu_{i}$ separately.

We need to comment on hadronic observables, which always depend on the
factorisation scale $\mu_{F}$. In ref. Bonvini:2020xeo  the factorisation
scale was explicitly excluded from the marginalisation procedure described in
section 2.1.1. Indeed, the factorisation scale also appears in the PDFs. It is
not clear what is the correlation and how the PDF uncertainty impacts the
probability distribution ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$, which is
supposed to represent the probability distribution for the true value of the
observable and should supersede and incorporate all other uncertainties. Here
we aim for something much simpler. We want to use the distribution
${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ to construct credibility intervals
that serve as estimators for the MHO uncertainty. In this sense, the role of
${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ is similar to scale variation as a
tool to estimate the size of the MHOs. Like in the case of scale variation, we
will assume that the MHO uncertainty is independent of the PDF uncertainty,
and the two can be added in quadrature.

###  2.2 Connection to prescriptions for choosing the perturbative scales

In the previous section we have presented two methods to construct a scale-
independent probability distribution
${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$. The two methods differ not only by
a different philosophy of how the unphysical perturbative scales are included,
but can also lead to different results for the inferred probability
distributions. Here we show that under rather general conditions the position
of the peak of ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ as a function of
$\Sigma$ will be biased towards two different scales $\mu_{\text{FAC}}$ and
$\mu_{\text{PMS}}$:

  * •

The point of _fastest apparent convergence_ (FAC) Grunberg:1980ja ;
Stevenson:1981vj  of an observable computed at NnLO is the scale
$\mu_{\text{FAC}}$ such that

$\Sigma_{n}(\mu_{\text{FAC}})=\Sigma_{n-1}(\mu_{\text{FAC}})\,,$ (20)

where we used the notation $\Sigma_{n}(\mu):=\Sigma_{n}(\mu,\mu)$. Note that
the value of $\mu_{\text{FAC}}$ depends on the order at which it is computed,
and there may be more than one solution to eq. (20) at a given order in
perturbation theory.

  * •

The _principle of minimal sensitivity_ (PMS) Stevenson:1981vj  aims at
selecting a scale $\mu_{\text{PMS}}$ at NnLO that maximises the stability of
the observable under changes of the scale

$\frac{\partial}{\partial\mu}\Sigma_{n}(\mu)\big{|}_{\mu=\mu_{\text{PMS}}}=0\,.$
(21)

Just like the FAC point, the PMS point depends on the order on which it is
computed and may not be unique at a given order in perturbation theory.

Consider now a broad class of Bayesian inference models satisfying the
following conditions:

  1. 1.

We work with the prior in eq. (10) or the weighting function in eq. (17).

  2. 2.

$P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}_{n}(\mu))$ is _symmetric_ , i.e., it is
insensitive to the signs of the correction:

$P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}_{n}(\mu))=P(|\Sigma^{(n+1)}(\mu)|\,\big{|}\,|\Sigma^{(0)}(\mu)|,\ldots,|\Sigma^{(n)}(\mu)|)\,.$
(22)

  3. 3.

The distributions $P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}_{n}(\mu))$ have a
single peak in $\Sigma^{(n+1)}(\mu)$ for each fixed value of $\mu$, and the
peak is reached for some $\mu\in[\mu_{0}/F,\mu_{0}F]$.

  4. 4.

The peak is more and more pronounced as $n$ increases, i.e., as more
information on the progression of the perturbative series becomes known.

Note that these assumption imply that the peak of
$P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}_{n}(\mu))$ is located at
$\Sigma^{(n+1)}(\mu)=0$, for all $\mu$. The aforementioned assumptions are
reasonable, and they are satisfied for some of the models considered in ref.
Bonvini:2020xeo , in particular by the geometric model and the scale-variation
model of ref. Bonvini:2020xeo  (see also section 3.2). In appendix A we show
that any model satisfying the assumptions above will also satisfy the
following property for the peak of the distribution
${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$:

  * •

If there is a single $\mu_{\text{FAC}}\in[\mu_{0}/F,\mu_{0}F]$, then
${\mathcal{P}}_{\text{sm}}(\Sigma|{\bm{\Sigma}}_{n})$ peaks at
$\Sigma=\Sigma_{n}(\mu_{\text{FAC}})$, i.e., at the NnLO prediction for the
observable evaluated at the point of fasted apparent convergence.

  * •

If there is a single $\mu_{\text{PMS}}\in[\mu_{0}/F,\mu_{0}F]$, then
${\mathcal{P}}_{\text{sa}}(\Sigma|{\bm{\Sigma}}_{n})$ peaks at
$\Sigma=\Sigma_{n}(\mu_{\text{PMS}})$, i.e., at the NnLO prediction for the
observable evaluated at the point of minimal sensitivity.

As a consequence, the sensitivity of the probability distributions and CIs on
the choice of the scale variation interval (given by $F$) depends on the
presence (or absence) of FAC and PMS points in the interval for the sm- and
sa-prescriptions respectively. We will illustrate this on concrete examples in
section 5.

Let us conclude this section with an important comment about the sm-
prescription. In the sm-prescription, the perturbative scales are treated as
model parameters with a statistical interpretation. Bayesian inference will
then select the ‘best’ value for the perturbative scales according to the
model, see also the discussion in ref. Bonvini:2020xeo . At high enough
perturbative orders, the dependence on the priors is expected to be mild, so
that the conclusions derived from the inference are robust and (to a certain
extent) prior-independent. The results of this section indicate that such a
Bayesian interpretation of the perturbative scales should be taken with a
grain of salt: For symmetric models with the prior in eq. (10) (which applies
in particular to some of the models of ref. Bonvini:2020xeo ), as soon as the
range $[\mu_{0}/F,\mu_{0}F]$ contains an FAC point, the scale selected by the
inference is biased: it will _automatically_ select the value of the
observable at the FAC point as the most likely value. In other words, such a
model will _always_ select the FAC point as its ‘best’ estimate for the
scales, independently of the actual progression of the perturbative series.
Note that the previous discussion does not apply to the sa-prescription,
because the perturbative scales are not assigned a probabilistic
interpretation (in particular, the weight $w(\mu)$ in eq. (17) is not a
probability distribution). In fact we will see examples in section 5 where the
CIs for the sa-prescription increase monotonically with $F$.

###  2.3 Connection to MHO estimates from scale variation

One may wonder about the practicality of the prescriptions in sections 2.1.1
and 2.1.2, because evaluating the integrals of the scale could be rather slow,
as it may require to evaluate the perturbative approximation of $\Sigma$ for
many values of the scales $\mu_{i}$. In many examples of higher-order
computations for the LHC, the dependence on the factorisation and
renormalisation scales is rather smooth. We can then approximate the integrals
over scale in eqs. (14) and (18) by a low-order quadrature rule.

We focus for now on observables depending on a single scale $\mu$, and we
comment on multiple scales later. It is particularly convenient to choose the
quadrature points to coincide with the ones used in the conventional scale
variation, namely $\mu/\mu_{0}\in\\{1/2,1,2\\}$.555We note that the functional
dependence on $\mu_{R}$ is determined by the running of the coupling constant
and the 9-point $\mu_{F},\mu_{R}$ grid can also be reconstructed from the
7-points given in eq. (5) in some cases. For a dynamical scale choice, this
reconstruction is typically only feasible for distributions that allow for the
determination of the numerical value of the scale that contributes to a
specific bin within the uncertainties of the bin size.  For the three-point
Gauss-Legendre quadrature rule this can be achieved by using the value of
$F=2^{\sqrt{5/3}}\approx 2.45$ for $w(\mu)$ in eq. (17). Then the integral in
eq. (18) can be approximated as

$\int d\mu\,w(\mu)\,f(\mu)\approx
w_{-1}\,f(\mu_{0}/2)+w_{0}\,f(\mu_{0})+w_{+1}\,f(2\mu_{0})\,,$ (23)

where the weights $w_{i}$ for the Gauss-Legendre quadrature are

$w_{0}=\frac{8}{18}\text{~{}~{}~{}and~{}~{}~{}}w_{\pm 1}=\frac{5}{18}\,.$ (24)

This quadrature achieves high accuracy for smooth functions in the interval
$-\log F<\log\mu/\mu_{0}<\log F$ (it is in fact exact if $f(\mu)$ is
polynomial in $\log\mu$ up to degree 5; note that hadronic observables are
never of this form, because some logarithms have been resummed to all orders
into the PDFs and $\alpha_{s}(\mu)$). Similar formulas can of course be
derived for quadrature rules of higher order and multi-dimensional $\mu$
integrals. The advantage of the three-point rule in eq. (23) is that it allows
one to easily convert the available data from the scale-variation approach
into an approximation for the Bayesian probability distribution, which is
scale-marginalised or scale-averaged over the interval
$\mu_{0}/F<\mu<\mu_{0}F$ with $F\approx 2.45$. This approximation is valid as
long as the three-point quadrature rule gives a good approximation of the
integral. In applications, the dependence on the scale is usually smooth, so
that low-order quadrature rules are expected to give reasonable
approximations. We will illustrate this on examples in section 5. One could
also restrict the integration region to the interval $[\mu_{0}/2,2\mu_{0}]$,
i.e. $F=2$, by using a quadrature rule where the nodes coincide with the end-
points of the integration, e.g., the Gauss-Lobatto quadrature rule (accurate
at cubic order). We compare the $\log\mu/\mu_{0}$ ranges of different
quadrature rules using the same scale variation points in the left panel of
figure 1.

Figure 1: Left: The integration regions for different quadrature rules using
scale variation points. Right: sketch of probability distributions at scale
variation points when the width is smaller than scale variation interval
(shaded) and when it is larger (dashed lines).

The quadrature rule in eq. (23) allows us to obtain an approximate bound on
the probability that the true value lies in the range obtained from
conventional scale variation. We will first discuss the case of an observable
depending on a single perturbative scale $\mu$, and generalise the result
later. Define the lower and upper values of the scale variation around some
central scale $\mu_{0}$ as

$\begin{split}\Sigma_{\text{min}}&\,:=\Sigma_{n}(k_{\textrm{min}}\mu_{0})=\min_{k\in\\{1/2,1,2\\}}\Sigma_{n}(k\mu_{0})\,,\\\
\Sigma_{\text{max}}&\,:=\Sigma_{n}(k_{\textrm{max}}\mu_{0})=\max_{k\in\\{1/2,1,2\\}}\Sigma_{n}(k\mu_{0})\,.\end{split}$
(25)

Consider the sa-prescription and a model such that

  1. 1.

the model is symmetric,

  2. 2.

$P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}(\mu))$ peaks only at
$\Sigma^{(n+1)}(\mu)=0$, for all values of $\mu$,

  3. 3.

the three-point quadrature in eq. (23) gives a good approximation of the
integral.

If these assumptions are satisfied, then we can show that the probability that
$\Sigma\in[\Sigma_{\text{min}},\Sigma_{\text{max}}]$ is bounded from above:

${\mathcal{P}}_{\text{sa}}(\Sigma\in[\Sigma_{\text{min}},\Sigma_{\text{max}}]|{\bm{\Sigma}}_{n})\lesssim
1-\frac{1}{2}(w_{\textrm{min}}+w_{\textrm{max}})\,,$ (26)

where $w_{\textrm{min}}$ is the quadrature weight associated to the scale
$\mu=k_{\textrm{min}}\mu_{0}$ through eq. (23), and similarly for
$w_{\textrm{max}}$. The $\lesssim$ sign reminds us that the bound is only as
accurate as the assumptions above, in particular the validity of eq. (23). To
understand the origin of this approximate bound, we first note that

${\mathcal{P}}_{\text{sa}}(\Sigma\in[\Sigma_{\text{min}},\Sigma_{\text{max}}]|{\bm{\Sigma}}_{n})=\Phi(\Sigma_{\text{max}})-\Phi(\Sigma_{\text{min}})\,,$
(27)

where $\Phi(\Sigma)$ is the cumulative distribution
$\Phi(\Sigma):=\int_{-\infty}^{\Sigma}d\Sigma^{\prime}\,{\mathcal{P}}_{\text{sa}}(\Sigma^{\prime}|{\bm{\Sigma}}_{n})$.
Using eqs. (18) and (23), we express $\Phi(\Sigma)$ as a sum of cumulative
distributions
$\varphi(\Sigma,\mu):=\int_{-\infty}^{\Sigma}d\Sigma^{\prime}\,P(\Sigma^{\prime}-\Sigma_{n}(\mu)|{\bm{\Sigma}}_{n}(\mu))$
for a fixed choice of the scale:

$\begin{split}\Phi(\Sigma_{\text{min}})&\,\approx
w_{-1}\,\varphi(\Sigma_{\text{min}},\mu_{0}/2)+w_{0}\,\varphi(\Sigma_{\text{min}},\mu_{0})+w_{+1}\,\varphi(\Sigma_{\text{min}},2\mu_{0})\\\
&\,\gtrsim
w_{\textrm{min}}\,\varphi(\Sigma_{\text{min}},k_{\textrm{min}}\mu_{0})\,,\end{split}$
(28)

where the last step follows from the fact that all three terms are positive.
An illustration of the three functions
$P(\Sigma-\Sigma_{n}(k\mu_{0})|{\bm{\Sigma}}_{n}(k\mu_{0}))$,
$k\in\\{1/2,1,2\\}$, is shown in figure 1 for two scenarios: when the
distributions are narrow compared to the scale variation interval
$[\Sigma_{\text{min}},\Sigma_{\text{max}}]$ (shaded areas) and when the
distributions are wider than $[\Sigma_{\text{min}},\Sigma_{\text{max}}]$
(dashed lines). Note that we do not know which value of $k$ corresponds to
which distribution in figure 1. However, from the assumptions above it follows
that only one of them peaks at $\Sigma_{\text{min}}$, namely the one
corresponding to the value of $k=k_{\textrm{min}}$. Since the distribution is
by assumption symmetric around $\Sigma=\Sigma_{\text{min}}$, we must have
$\varphi(\Sigma_{\text{min}},k_{\textrm{min}}\mu_{0})=\frac{1}{2}$. Hence, we
obtain

$\Phi(\Sigma_{\text{min}})\gtrsim\frac{1}{2}\,w_{\textrm{min}}\,.$ (29)

Using a similar argument, we can conclude

$\Phi(\Sigma_{\text{max}})\lesssim 1-\frac{1}{2}\,w_{\textrm{max}}\,.$ (30)

Equation (26) then follows upon inserting the bounds for
$\Phi(\Sigma_{\text{min}})$ and $\Phi(\Sigma_{\text{max}})$ into eq. (27).

We can ask when the approximate bound in eq. (26) is saturated. It is easy to
see that the bound for the cumulative distribution in eq. (29) is saturated if
the three distributions
$P(\Sigma-\Sigma_{n}(k\mu_{0})|{\bm{\Sigma}}_{n}(k\mu_{0}))$,
$k\in\\{1/2,1,2\\}$ are narrow compared to their separation, because in that
scenario we have
$\Phi(\Sigma_{\textrm{min}})\approx\frac{1}{2}\,w_{\textrm{min}}$. (see figure
1). A similar argument holds for the bound in eq. (30),
$\Phi(\Sigma_{\textrm{max}})\approx 1-\frac{1}{2}\,w_{\textrm{max}}$. Hence,
we conclude that the approximate bounds on the probability for
$[\Sigma_{\text{min}},\Sigma_{\text{max}}]$ are saturated for observables and
models for which the distributions
$P(\Sigma-\Sigma_{n}(k\mu_{0})|{\bm{\Sigma}}_{n}(k\mu_{0}))$,
$k\in\\{1/2,1,2\\}$ are well separated.

Quadrature rule | $\beta_{1}$ | $\beta_{2}$
---|---|---
Gauss-Legendre | $13/18\approx 72\%$ | $\approx 92\%$
Gauss-Lobatto | $\phantom{11}5/6\approx 83\%$ | $\approx 97\%$
Equidistant | $\phantom{11}2/3\approx 66\%$ | $\approx 89\%$
Table 1: Approximate bounds on the probability for interval
$[\Sigma_{\text{min}},\Sigma_{\text{max}}]$.

We can repeat exactly the same argument for observables depending on multiple
scales $(\mu_{1},\ldots,\mu_{s})$ by using the same quadrature rule $s$ times
for scale-averaging for each scale $\mu_{i}$ in eq. (28). Each cumulative
distribution $\varphi\left(\Sigma,\mu_{1},\ldots,\mu_{s}\right)$ in the sum is
multiplied by a product of $s$ quadrature weights. It is then straightforward
to show that the bound becomes

${\mathcal{P}}_{\text{sa}}(\Sigma\in[\Sigma_{\text{min}},\Sigma_{\text{max}}]|{\bm{\Sigma}}_{n})\lesssim\beta_{s}\,,\textrm{
with
}\beta_{s}:=1-\frac{1}{2^{s}}\,(w_{\textrm{min}}+w_{\textrm{max}})^{s}\,.$
(31)

Let us discuss some implications of the approximate bound in eq. (31). Under
the assumptions spelled out above, eq. (31) establishes a connection between
the interval obtained by varying the scale up and down by a factor of $2$ and
the Bayesian probability distribution constructed from the sa-prescription.
The value of the approximate bound depends on the quadrature weights used to
approximate the integral. For the three quadrature rules illustrated in figure
1, the values of the bound are shown in table 1 for $s=1,2$. Note that
different quadratures correspond to different intervals over which the scale-
averaging is done. In addition, the accuracy of the approximation in eq. (23)
also depends on the order of the quadrature rule. Therefore the approximate
bounds in table 1 depend on how well the assumptions stated at the beginning
of this section are satisfied. In particular, it _does not_ allow one to
attach a definite probability that the true value lies in the range given by
the scale variation. However, we can show that the $(100\beta_{s})\%$
credibility interval
$\text{CI}_{100\beta_{s}}=[\Sigma_{100\beta_{s}}^{\text{low}},\Sigma_{100\beta_{s}}^{\text{up}}]$
will contain $[\Sigma_{\text{min}},\Sigma_{\text{max}}]$ and the two intervals
(approximately) agree whenever the approximate bound is saturated. Indeed,
from the definition of the credibility interval in eq. (6) it follows that

$\Phi(\Sigma_{100\beta_{s}}^{\text{low}})=\frac{1}{2}(1-\beta_{s})\textrm{ and
}\Phi(\Sigma_{100\beta_{s}}^{\text{up}})=\frac{1}{2}(1+\beta_{s})\,.$ (32)

Equations (29) and (30) and their generalisations for $s>1$ then imply:

$\Phi(\Sigma_{100\beta_{s}}^{\text{low}})\lesssim\Phi(\Sigma_{\text{min}})\leq\Phi(\Sigma_{\text{max}})\lesssim\Phi(\Sigma_{100\beta_{s}}^{\text{up}})\,,$
(33)

and so the monotonicity of the cumulative distribution implies

$\Sigma_{100\beta_{s}}^{\text{low}}\lesssim\Sigma_{\text{min}}\leq\Sigma_{\text{max}}\lesssim\Sigma_{100\beta_{s}}^{\text{up}}\,,$
(34)

and the claim follows.

The previous discussion has an important consequence: Assume that we are in a
situation with $s=2$ scales. Then we conclude from table 1 that the 9-point
scale variation interval satisfies

$[\Sigma_{\text{min}},\Sigma_{\text{max}}]\subset\textrm{CI}_{100\beta_{2}}\approx\textrm{CI}_{89-97}\,.$
(35)

Clearly, the 7-point interval is also contained in the same interval. In the
case where the approximate bound is saturated, we expect that for two scales
$(\mu_{F},\mu_{R})$ the intervals obtained from scale variation are comparable
to $\approx 95\%$ credibility intervals. In practise, we often find that
observables are much less sensitive to one of the scales, namely, $\mu_{F}$,
and we are effectively dealing with $s=1$ case, where $\beta_{1}\approx
66$–$83\%$. We will illustrate in section 5 that the conventional scale
variation intervals are often similar in size to CI68 using Bayesian methods.
The arguments presented in this section show that this is not a coincidence.
Nevertheless we reiterate that the derivation of these bounds relies on
specific assumptions. In particular, the sm-prescription and models with
asymmetric probability distributions do not need comply with these
requirements and the bounds can be broken. However, as we will see in the
examples in section 5, in many cases the size of the CIs is similar for
symmetric and asymmetric models.

##  3 Modelling probability distributions at fixed scales

In the previous section we have discussed two ways to define a probability
distribution ${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ in the context of
perturbation theory. In both approaches the computation of the distribution
${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ is reduced to an integral over the
perturbative scale $\mu$. The integrand involves, on the one hand, the prior
$P_{0}(\mu)$ in eq. (10) or the weighting function $w(\mu)$ in eq. (17), and,
on the other hand, the probability $P(\Sigma|{\bm{\Sigma}}_{n}(\mu))$
evaluated at a fixed value $\mu$ of the perturbative scale. The latter is
related through eq. (13) to the probability distribution
$P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}_{n}(\mu))$ for the size of the next
perturbative order at the scale $\mu$, given the $n+1$ first orders computed
at the same scale. Various models for
$P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}_{n}(\mu))$ have been presented in the
literature, in particular the CH-method of ref. Cacciari:2011ze  and its
generalisations of ref. Bonvini:2020xeo . All theses examples define symmetric
distributions in the sense of eq. (22). In section 2.2 we have argued that
symmetric distributions always lead to very special distributions
${\mathcal{P}}(\Sigma|{\bm{\Sigma}}_{n})$ after integration over $\mu$, with
peaks that are correlated to the FAC and PMS points. This may introduce
unwanted bias into the model.

The goal of this section is to introduce a model that is not symmetric, which
is a necessary condition for distributions not necessarily peaked at the FAC
or PMS points. We start by presenting a framework to define different models
for $P(\Sigma^{(n+1)}(\mu)|{\bm{\Sigma}}_{n}(\mu))$. This framework contains,
in particular, the CH-method of ref. Cacciari:2011ze  and its generalisations
of ref. Bonvini:2020xeo  as special cases. We then use this framework to
define a new model that takes into account the signs of the corrections and
leads to asymmetric distributions.

###  3.1 The mathematical setup

Consider the set consisting of elements $({\bm{\delta}},{\mathbf{p}})$, where
${\mathbf{p}}\in\mathbb{R}^{p}$ is a vector and
${\bm{\delta}}=(\delta_{0},\delta_{1},\ldots)$ is an infinite sequence of real
numbers representing the infinite series $\sum_{k=0}^{\infty}\delta_{k}$. We
prefer to think of the sequence ${\bm{\delta}}$ rather than the infinite
series, because the order of the $\delta_{k}$’s matters.

We make the following assumptions:

  1. 1.

The parameters ${\mathbf{p}}$ are hidden, and we do not have access to them.
Instead, they are distributed according to a known prior distribution
$P_{0}({\mathbf{p}})$.

  2. 2.

The sequence ${\bm{\delta}}$ and the hidden parameters ${\mathbf{p}}$ are not
independent. For a given value of ${\mathbf{p}}$, we know the probability that
the first $n+1$ elements of ${\bm{\delta}}$ take a certain value, i.e., for
each $n\geq 0$, if
${\bm{\delta}}_{n}:=(\delta_{0},\delta_{1},\ldots,\delta_{n})$, we know the
conditional probability $P^{(n)}({\bm{\delta}}_{n}|{\mathbf{p}})$.

  3. 3.

The conditional probabilities $P^{(n)}({\bm{\delta}}_{n}|{\mathbf{p}})$ for
different values of $n$ are related by marginalisation:

$P^{(n-1)}({\bm{\delta}}_{n-1}|{\mathbf{p}})=\int
d\delta_{n}\,P^{(n)}({\bm{\delta}}_{n}|{\mathbf{p}})\,,\text{ for all }n>0\,.$
(36)

We refer to a choice of $(P_{0},P^{(0)},P^{(1)},P^{(2)},\ldots)$ as a _model_
for sequences $({\bm{\delta}},{\mathbf{p}})$. We always assume that all
probabilities are normalised, and we will often keep the superscript on
$P^{(n)}$ implicit and use the simplified notation
$P({\bm{\delta}}_{n}|{\mathbf{p}}):=P^{(n)}({\bm{\delta}}_{n}|{\mathbf{p}})$.

Assume that we are given the first $n+1$ terms ${\bm{\delta}}_{n}$ of
${\bm{\delta}}$. The main idea is that, within a given model, we can use
Bayesian inference to obtain a conditional probability distribution
$P(\delta_{n+1}|{\bm{\delta}}_{n})$ for the next term $\delta_{n+1}$ in the
sequence ${\bm{\delta}}$, given the first $n+1$ terms ${\bm{\delta}}_{n}$. The
probability distribution $P(\delta_{n+1}|{\bm{\delta}}_{n})$ is given by
Bayes’ formula:

$P(\delta_{n+1}|{\bm{\delta}}_{n})=\frac{P({\bm{\delta}}_{n+1})}{P({\bm{\delta}}_{n})}\,.$
(37)

The probability distributions for $P({\bm{\delta}}_{n})$ and
$P({\bm{\delta}}_{n+1})$ are themselves obtained from Bayes’ formula:

$P({\bm{\delta}}_{n})=\int
d^{p}{\mathbf{p}}\,P({\bm{\delta}}_{n}|{\mathbf{p}})\,P_{0}({\mathbf{p}})\,.$
(38)

We refer to a model as _symmetric_ if it satisfies (cf. eq. (22))

$P(\delta_{0},\ldots,\delta_{n})=P(|\delta_{0}|,\ldots,|\delta_{n}|)\,.$ (39)

A symmetric model always predicts a vanishing median and expectation value for
$\delta_{n+1}$ given ${\bm{\delta}}_{n}$,

$\langle\delta_{n+1}\rangle_{{\bm{\delta}}_{n}}=\delta_{n+1}^{\text{median}}=0\,,\quad\text{
if the model is symmetric}\,,$ (40)

and all CIs computed from a symmetric model are symmetric around
$\delta_{n+1}=0$.

###  3.2 The geometric model

The CH model assumes a universal bound $c$ for all coefficients
$\delta_{k}(\mu)=\Sigma^{(k)}(\mu)/\alpha_{s}^{k}(\mu)$ for a fixed value of
$\mu$. In ref. Bonvini:2020xeo  a variation of the CH model, called the
_geometric model_ , was introduced to take into account the fact that the
expansion parameter does not need to be equal to $\alpha_{s}$.666In ref.
Bagnaschi:2014wea  another modification of the CH model was introduced to take
this effect into account. We do not explore the method of ref.
Bagnaschi:2014wea  in this paper. The geometric model considers a sequence of
coefficients $\delta_{k}(\mu)=\Sigma^{(k)}(\mu)/\Sigma^{(0)}(\mu)$ (so that
$\delta_{0}(\mu)=1$), whose magnitude is bounded from above by a geometric
series. Consequently, the model has two hidden parameters
${\mathbf{p}}=(a,c)$, and the coefficients of the sequence
${\bm{\delta}}(\mu)$ must satisfy the constraint (for simplicity we suppress
the dependence on the scale):

$\frac{|\delta_{k}|}{a^{k}}\leq c\,,\quad\text{ for all }k\geq 0\,.$ (41)

The constraints for $k=0,1$ imply $c>1$ and $a>0$. Since we do not make any
other assumption about the coefficients, we choose a flat distribution that
takes into account eq. (41) Bonvini:2020xeo :

$\begin{split}P_{\text{geom}}({\bm{\delta}}_{n}|a,c)&\,=\prod_{k=0}^{m}\frac{1}{2a^{k}c}\,\Theta\left(c-\frac{|\delta_{k}|}{a^{k}}\right)\\\
&\,=\frac{1}{(2c)^{n+1}a^{n(n+1)/2}}\,\Theta\left(c-\max\left(\frac{|\delta_{0}|}{a^{0}},\ldots,\frac{|\delta_{n}|}{a^{n}}\right)\right)\,.\end{split}$
(42)

The model is clearly symmetric. For the prior on the hidden parameters, $a$
and $c$ are assumed independent, $P_{0}(a,c)=P_{0}(a)P_{0}(c)$. The prior on
$c$ is given by

$P_{0}(c)=\frac{\varepsilon}{c^{1+\varepsilon}}\,\Theta(c-1)\,.$ (43)

Here $\varepsilon>0$ is a constant introduced to obtain a normalisable prior.
The prior on $a$ is Bonvini:2020xeo :

$P_{0}(a)=(1+\omega)\,(1-a)^{\omega}\,\Theta(a)\,\Theta(1-a)\,,$ (44)

where $\omega>0$ suppresses the prior at $a=1$. Note that $\varepsilon$ and
$\omega$ are constants, and they are not considered hidden parameters of the
model. For comparisons with the asymmetric model defined in the next section,
we will use $\epsilon=0.1$ and $\omega=1$. The form of eq. (42) and the priors
are simple enough that the integral in eq. (38) can be performed in closed
form, see ref. Bonvini:2020xeo .

Let us conclude this section with some comments. First, we note that the
geometric model contains the original CH model as a special case, where
$a=\alpha_{s}(\mu)$, i.e., the prior for $a$ is

$P_{0,\textrm{CH}}(a)=\delta(a-\alpha_{s}(\mu))\,,$ (45)

so that

$P_{\textrm{CH}}({\bm{\delta}}_{n}|c)=\int
da\,P_{\text{geom}}({\bm{\delta}}_{n}|a,c)\,P_{0,\textrm{CH}}(a)=P_{\text{geom}}({\bm{\delta}}_{n}|\alpha_{s}(\mu),c)\,.$
(46)

Second, so far we have chosen the coefficients $\delta_{k}(\mu)$ to be real
numbers. We can assume that this sequence was obtained by evaluating some
sequence of functions $\Sigma^{(k)}$ at some fixed value $\mu$, i.e.,
$\delta_{k}(\mu)=\Sigma^{(k)}(\mu)/\Sigma^{(0)}(\mu)$. Reference
Bonvini:2020xeo  proposes a variant where the probability distribution takes
input from other values of the functions $\Sigma^{(k)}$. The _scale variation
model_ of ref. Bonvini:2020xeo  considers sequences ${\bm{\delta}}(\mu)$ whose
coefficients satisfy a bound of the form:

$\frac{|\delta_{k}(\mu)|}{r_{k}(\mu)}\leq c\,,\quad\text{ for all }k\geq 0\,,$
(47)

where $r_{k}(\mu)$ encodes information on the derivative of $\Sigma_{k}(\mu)$
with respect to $\mu$. The probability distribution for _fixed_ $\mu$ is
related to $P_{\textrm{CH}}$ via the relation (with $\alpha_{s}$ replaced by
unity):

$P_{\text{sc.-var.}}({\bm{\delta}}_{n}(\mu)|c)=P_{\text{CH}}\left(\frac{\delta_{0}(\mu)}{r_{0}(\mu)},\ldots,\frac{\delta_{n}(\mu)}{r_{n}(\mu)}\Big{|}c\right)\,\prod_{k=0}^{n}\frac{1}{r_{k}(\mu)}\,.$
(48)

In this case the prior for $c$ is chosen in ref. Bonvini:2020xeo  to be a
Poisson distribution.

###  3.3 An asymmetric geometric model: the $abc$-model

The symmetry in the geometric model can be traced back to the fact that the
constraint in eq. (41) defining the model only depends on the magnitude of the
coefficients, and it is insensitive to their signs. Alternative models that
lead to asymmetric distributions and include information about coefficient’s
signs are discussed in appendix B of ref. Bonvini:2020xeo  (see in particular
B.2 and B.5). Here, we consider a new model where the upper and lower bounds
are different in magnitude. More precisely, we have three hidden parameters
${\mathbf{p}}=(a,b,c)$, and we consider a sequence of coefficients
$\delta_{k}=\Sigma^{(k)}/\Sigma^{(0)}$ such that (again, we suppress the
dependence on the scale $\mu$)

$b-c\,\leq\frac{\delta_{k}}{a^{k}}\leq b+c\,,\quad\text{ for all }k\geq 0\,.$
(49)

We refer to this model as the _asymmetric geometric model_ , or the _$abc$
-model_. The conditional probability distribution takes a form very similar to
the one in the geometric model in eq. (42):

$\displaystyle P_{abc}({\bm{\delta}}_{n}|a,b,c)$
$\displaystyle\,=\prod_{k=0}^{m}\frac{1}{2|a|^{k}c}\,\Theta\left(c-\left|\frac{\delta_{k}}{a^{k}}-b\right|\right)$
(50)
$\displaystyle\,=\frac{1}{(2c)^{n+1}|a|^{n(n+1)/2}}\,\Theta\left(c-\max\left(\left|\frac{\delta_{0}}{a^{0}}-b\right|,\ldots,\left|\frac{\delta_{n}}{a^{n}}-b\right|\right)\right)\,.$

Like in the geometric model, our choice of priors assumes that $a$ is
independent of the parameters $b$ and $c$. However, since $b$ and $c$ jointly
define the upper and lower bounds, we do not see any reason to consider them
independent. Instead, we use a joint probability for $(b,c)$:

$P_{0}(a,b,c)=P_{0}(a)\,P_{0}(b,c)\,,$ (51)

where the factors in the right-hand side are given by:

$\begin{split}P_{0}(a)&\,=\frac{1}{2}(1+\omega)\,(1-|a|)^{\omega}\,\Theta(1-|a|)\,,\\\
P_{0}(b,c)&\,=\frac{\varepsilon\,\eta^{\epsilon}}{2\xi
c^{2+\varepsilon}}\,\Theta(c-\eta)\,\Theta(\xi c-|b|)\,.\end{split}$ (52)

Unlike in the geometric model, we do not constrain $a$ to be positive and only
require that $|a|<1$. The prior $P_{0}(b,c)$ is characterised by three
constants $\epsilon,\eta$ and $\xi$. As in eq. (43), $\epsilon>0$ is used to
make the prior normalisable. The second parameter $\eta>0$ determines the
minimal difference between the upper and lower bounds of the series, see eq.
(49). Finally, $\xi>0$ restricts the allowed range of $b$, with $\xi=0$
forcing $b$ to vanish and we recover a symmetric geometric model.

Figure 2: The shaded area indicates the ranges for $b$ and $c$ allowed by the
prior $P_{0}(b,c)$, in eq. (52) and the inference from the universal input
$\delta_{0}=1$ in eq (50) for $\xi>1$ (left) and $\xi<1$ (right). The position
of $c=\eta>0$ line is arbitrary.

We note that the constraints in eq. (49) for $k=0$ imply $c>|b-1|$. Combining
this with the prior in eq. (52), we obtain $c>\max\\{\eta,\frac{1}{1+\xi}\\}$.
In figure 2 we show the ranges for $b$ and $c$ allowed by the prior in eq.
(52) and $c>|b-1|$ for the two cases $\xi>1$ and $\xi<1$. We can clearly see
from figure 2 that if $\eta<\frac{1}{1+\xi}$, then $\eta$ no longer affects
the model inference. Similarly if $\xi>\frac{\eta+1}{\eta}>1$, i.e.,
$\eta>\frac{1}{\xi-1}$, then the inference becomes independent of $\xi$. Every
subsequent term $\ \delta_{k}$ in the series adds the constraint
$c>\left|\frac{\delta_{k}}{a^{k}}-b\right|$, whose position depends on the
values of $\delta_{k}$ and $a$.

We conclude this section by stressing that the quantities $\omega$, $\xi$,
$\epsilon$ and $\eta$ are not hidden parameters of the model. Rather, they are
(positive) constants defining a family of $abc$-models. In the next section we
will study the model sensitivity to the choice of these constants. However,
one has to keep in mind that this only explores the robustness for the
specific model defined by eqs. (49), (50) and (52). There is an infinite-
dimensional functional space of possible models and priors that one might
consider. We believe that the $abc$-model offers a good compromise between the
flexibility of the model (it gives asymmetric probability distribution and can
be equally well applied to monotonic and simple alternating series) and the
transparency and simplicity of the model.

####  3.3.1 Sensitivity to the choice of priors

Figure 3:  Top left panel: The probability distribution from the $abc$-model
for $S_{n+1}^{\text{est}}$ for different values of $n$ for the geometric
series with $x=0.7$. Top right panel: The same as the left panel, but we show
the probability for the scaled deviation from the known correction
$(S_{n+1}^{\text{est}}-S_{n+1})/|S_{n+1}-S_{n}|$. Bottom left panel: The
median (plus), 68% CI (errorbox) and 95% CI (errorbar) for the posterior of
$S_{n+1}^{\text{est}}$ , computed from the $abc$ (blue) and geometric (red)
models using information on the previous orders. The exact values of $S_{n}$
are shown as black circles. Bottom right panel: The same as the left panel,
but the exact $S_{n}$ value is subtracted from $S_{n+1}^{\text{est}}$ and the
difference is normalised by $|S_{n+1}-S_{n}|$.

Let us illustrate the $abc$-model and its sensitivity to the choice of priors
on the example of the geometric series $\delta_{k}=x^{k}$, for $|x|<1$. Then
there exists $a=x$ for which the region of values of $(b,c)$ excluded by eq.
(50) collapses to $c>|1-b|$. We see from figure 2 that $\eta>0$ prevents $c$
from reaching zero. Setting $\eta\to 0$ and $\xi\to\infty$ turns off any
constraints by the prior $P_{0}(b,c)$, and the $abc$-model can describe such a
geometric series arbitrary well with $a\to x,b\to 1,c\to 0$ in eq. (50).

We compute the conditional probability for the estimated next term in the
series $\delta_{n+1}^{\text{est}}$:

$P_{abc}(\delta_{n+1}^{\text{est}}|{\bm{\delta}}_{n})=\frac{P_{abc}({\bm{\delta}}_{n+1})}{P_{abc}({\bm{\delta}}_{n})}\,,$
(53)

where

$\displaystyle P_{abc}({\bm{\delta}}_{n})$ $\displaystyle=\int
da\,db\,dc\,P_{abc}({\bm{\delta}}_{n}|a,b,c)P_{0}(a)P_{0}(b,c)$ (54)
$\displaystyle=\int da\frac{P_{0}(a)}{|a|^{n(n+1)/2}}\int
db\,\frac{\epsilon\,\eta^{\epsilon}}{2^{n+2}\xi(n+2+\epsilon)}\max\left(\eta,\max_{0\leq
k\leq
n}\left|\frac{\delta_{k}}{a^{k}}-b\right|,\frac{|b|}{\xi}\right)^{-(n+2+\epsilon)}$

For fixed $a$, the integral over $b$ can be done analytically and expressed in
terms of maximum and minimum functions of $\frac{\delta_{k}}{a^{k}}$, $\eta$
and $\xi$. In appendix B we show how to evaluate the distribution
$P_{abc}({\bm{\delta}}_{n})$ analytically in terms of Gauss’ hypergeometric
function, allowing for a fast and efficient numerical implementation of the
model.

In the top left panel of figure 3 we show the probability distributions for
the partial sums $S_{n}:=\sum_{k=0}^{n}\delta_{k}=\sum_{k=0}^{n}x^{k}$ for
$x=0.7$ and $n\leq 7$ using the $abc$-model with parameter values
$(\epsilon,\omega,\xi,\eta)=(0.1,1,2,0.1)$ (see the discussion below for the
choice of these values). For $n=0$, the probability distribution is symmetric
and centred around $S_{0}=1$. For $n>0$ the distributions are clearly not
symmetric and become more and more peaked as $n$ increases. In the top right
panel we show the probability distributions for the scaled deviation from the
known correction
$(S_{n+1}^{\text{est}}-S_{n+1})/|S_{n+1}-S_{n}|=(\delta_{n+1}^{\text{est}}-\delta_{n+1})/|\delta_{n+1}|$.
This allows us to compare different orders without the suppression of the
expansion parameter. In this plot the $S_{n+1}$ value corresponds to zero on
the $x$-axis, while $S_{n}$ corresponds to $\pm 1$ (depending on the sign of
$\delta_{n+1}$). Again, for $n=0$ the distribution is centred around the
initial value $S_{0}$, but for each subsequent order, the distribution shifts
towards the true value. We note that the shape of the distribution does not
change significantly beyond $n=3$, so that the narrowness of the distributions
in the left panel is solely due to the higher power of the expansion
parameter. In the simple case of a geometric series, the scaled distributions
are controlled by the prior $P_{0}(b,c)$. For $\eta=0.1$ and $\xi=2$, the
minimum of $c$ is at
$(b,c)=(\frac{\xi}{1+\xi},\frac{1}{1+\xi})\approx(0.6,0.3)$ (see figure 2),
which gives a good approximation for the centre position and the (half) width
of distributions seen in figure 3. As explained earlier, increasing $\xi$ and
reducing $\eta$ can narrow down the posterior distribution around the exact
geometric series value $x^{n+1}$.

In the bottom panels of figure 3 we show the 68% and 95% CIs for the $abc$
(blue) and geometric models (red) (we use the same values
$(\omega,\epsilon)=(1,0.1)$ for both models). On the left we show the CIs for
$S_{n+1}^{\text{est}}$, while on the right we again scale to the known size of
the next term in the series. We see that for $n>0$ the 68% CIs computed with
the $abc$-model encompass the next value of $S_{n}$. This is not the case for
the geometric model, even though the 68% CIs are larger in magnitude. This is
due to the fact that the geometric model is symmetric around the last known
order, thereby not accommodating the fact that all the coefficients
$\delta_{k}$ are positive. In contrast, the CIs of the $abc$-model are
systematically shifted towards the exact values for $n>3$, and the 95% CIs
even exclude the previous value $S_{n}$. However, one must remember that these
features of the $abc$-model may not always be desired in practical
applications. and can be controlled by the priors (see below).

For completeness, in appendix C we have also included a similar plot showing
the alternating series with $x=-0.7$. Here we just note that $abc$-model
correctly predicts the negative sign of the expansion parameter and remains
consistently shifted towards the next value in the geometric series
$\delta_{n+1}$. In all other respects the distributions and the CIs are
identical to the ones shown in figure 3.

Figure 4: The dependence of the probability distribution
$1/S_{0}P_{abc}((S_{n+1}^{\text{est}}-S_{n})/S_{0}|{\bm{\delta}}_{n})$ for
geometric series for $n=1$ (dashed) and $n=3$ (solid) on the parameters
$(\epsilon,\omega,\eta,\xi)$ for 4 selected values, with the others held fixed
at $(\epsilon,\omega,\eta,\xi)=(0.1,1,0.1,2)$.

In figure 4 we study the dependence of
$P_{abc}(\delta_{n+1}|{\bm{\delta}}_{n})$ on the choice of the values for
$(\epsilon,\omega,\eta,\xi)$ that define the priors on $(a,b,c)$ in eq. (52).
We vary each parameter independently, with the other held fixed at their
default values $(\epsilon,\omega,\xi,\eta)=(0.1,1,2,0.1)$. We see that, as
expected, the dependence on $\epsilon$, $\omega$, and $\eta$ reduces as $n$
increases. The $\xi$-dependence, however, does not disappear completely with
increasing $n$. This is due to the fact that $\xi$ controls the amount of
asymmetry of $P_{abc}(\delta_{n+1}|{\bm{\delta}}_{n})$. For small values of
$\xi$, the distribution approaches a symmetric distribution around the last
known value $S_{n}$, while for large $\xi$ it becomes much more peaked around
the next value $S_{n+1}$. Empirically, we find that our choice $\xi=2$ allows
distributions to be shifted towards the next term, but keeping the width of
the distribution of the same size as the shift. Next we look at $\omega$,
which controls the suppression of the prior for $a$ around $|a|=1$ ($\omega=0$
corresponds to a flat prior). Because here we consider a rather large value of
the expansion parameter $x=0.7$, increasing $\omega$ slightly reduces the
expectation value of $a$ and shifts the distribution to smaller values of
$S_{n+1}^{\text{est}}$. Our default choice of $\omega=1$ biases the hidden
model parameter inference towards small $|a|$ values. In the bottom left panel
of figure 4 we see that larger values of $\eta$ shift the posterior
distribution towards the previous order. $\eta$ plays a similar role to $\xi$
in controlling the allowed minimum value of $c$. We use a small value of
$\eta=0.1$, which is in fact superseded by the constraint from $\xi$ (see
figure 2). Finally, in contrast to $\eta$, for larger values of $\epsilon$ the
model becomes slightly more peaked and shifted towards the next term in the
series. We choose a small value $\epsilon=0.1$ to maintain fairly flat
distribution in the prior of $\log c$.

##  4 Quantities without explicit scale dependence

The goal of this section and the next is to illustrate the concepts and models
from the previous sections on various examples of observables in QFT. We start
by investigating several quantities whose perturbative expansion does not have
an explicit dependence on the perturbative scales, or for which scale
dependence is not relevant. This allows us to illustrate the application of
the $abc$-model to genuine perturbative expansions from QFT. The sequence
${\bm{\delta}}_{n}$ is the sequence of the $n+1$ first perturbative
coefficients, normalised such that $\delta_{0}=1$

$\delta_{k}:=\Sigma^{(k)}/\Sigma_{0}\,.$ (55)

The rescaling introduces a Jacobian into eq. (13), which now takes the form

$P(\Sigma|{\bm{\Sigma}}_{n})\approx\frac{1}{\Sigma_{0}}P_{abc}\left(\frac{\Sigma}{\Sigma_{0}}-\sum_{k=0}^{n}\delta_{k}\Big{|}{\bm{\delta}}_{n}\right).$
(56)

We then use this probability distribution to compute CIs for $\Sigma$ within
the $abc$-model using perturbative input through NnLO. These intervals
estimate the size of the missing Nn+1LO terms, and so they serve as measures
of the MHO uncertainty at NnLO. To assess the validity of this procedure, we
show in each case how the size of the CIs computed at different perturbative
orders compares to the actual size of the next order whenever it is available.
By abuse of notation, we will use both
$P_{abc}(\Sigma_{n+1}^{\text{est}}|{\bm{\delta}}_{n})$ and
$P_{abc}(\delta_{n+1}^{\text{est}}|{\bm{\delta}}_{n})$, where the relation to
eq. (56) has to be inferred from the argument.

Before discussing explicit examples of perturbative expansions in QFT, let us
make some general comments. First, we stress that the approximation in eq.
(56) is only valid if we are in the perturbative regime. It is in particular
not valid once we enter the regime where the perturbative series starts to
diverge. In that situation the right-hand side of eq. (56) may still represent
a valid probability distribution for the next perturbative coefficient, but
this distribution is no longer a good approximation of the left-hand side,
which is the probability distribution for the true value of $\Sigma$. We note
that for most applications only very few perturbative orders are available,
and this will not be a problem. Second, it is possible to be in a situation
where the model happens to describe the input data extremely well. The
$abc$-model strongly favours a perfect geometric series. If the available
perturbative data closely follow a geometric progression up to a certain order
by some numerical coincidence, then the resulting distributions can be
strongly peaked. In such a scenario, the 68% or 95% CIs may even exclude the
highest available order $\Sigma_{n}$. This is a _design feature_ of the
$abc$-model that can be controlled by the priors (see section 3.3).

We stress that regardless of how well the truncated perturbative series is
described by a chosen model, there is no guarantee that the actual next order
in perturbative QFT or the measured experimental value of $\Sigma$ will fall
into a specific choice of CI. Therefore, the choices of the model, the priors
and the CI used to estimate the MHO uncertainty of an NnLO computation should
be carefully scrutinised for each new observable separately.

###  4.1 The cusp anomalous dimension

Our first example is the light-like cusp anomalous dimension. It is known
through four loops in QCD and $\mathcal{N}=4$ Super Yang-Mills (SYM)
Brandt:1981kf ; Korchemsky:1987wg ; Moch:2004pa ; Grozin:2014hna ;
Grozin:2015kna ; Henn:2016men ; vonManteuffel:2016xki ; Davies:2016jie ;
Lee:2016ixa ; Boels:2017skl ; Boels:2017ftb ; Moch:2017uml ; Moch:2018wjh ;
Grozin:2018vdn ; Lee:2019zop ; Henn:2019rmi ; Bruser:2019auj ;
vonManteuffel:2019wbj ; Henn:2019swt ; vonManteuffel:2020vjv . In the planar
limit of the $\mathcal{N}=4$ SYM theory, it is even known to all orders in
perturbation theory from integrability Beisert:2006ez . The cusp anomalous
dimension is therefore an ideal candidate to study the applicability of
Bayesian techniques to derive MHO uncertainties for perturbative field theory
quantities.

##### Planar $\mathcal{N}=4$ SYM.

The cusp anomalous dimension in planar $\mathcal{N}=4$ SYM admits a
perturbative expansion:

$\gamma_{K}(\lambda)=\sum_{k=1}^{\infty}\gamma_{K}^{(k)}\,\lambda^{k}\,.$ (57)

Here the expansion parameter is the ’t Hooft coupling
$\lambda=\frac{g_{s}^{2}\,N_{c}}{8\pi^{2}}$, where $g_{s}$ is the Yang-Mills
coupling and $N_{c}$ is the number of fundamental SU$(N_{c})$ colours. Note
that $\mathcal{N}=4$ SYM is a conformal field theory (even away from the
planar limit), and so the coupling is not UV renormalised and does not depend
on any renormalisation scale. The first few expansion coefficients are shown
in table 2. It was observed that the series in eq. (57) is convergent if
$|\lambda|\lesssim\frac{1}{8}$.

$k$ | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
$\gamma_{K}^{(k)}$ | $4$ | $-6.57$ | $23.81$ | $-117.17$ | $662.62$ | $-4028.78$ | $25618.7$
$\gamma^{(k)}_{\rm QCD}$ | $1$ | $1.73$ | $2.80$ | $2.89$ | $-$ | $-$ | $-$
Table 2: The first few terms in the perturbative series of the cusp anomalous
dimension for planar $\mathcal{N}=4$ SYM, eq. (57), and QCD, eq. (59).  Figure
5: Top left panel: The probability distribution from the $abc$-model for the
cusp anomalous dimension $(\gamma_{K})_{n+1}^{\text{est}}$ in planar
$\mathcal{N}=4$ SYM for $\lambda=0.06$ and for different values of $n$. Top
right panel: The same distributions normalised to the exact Nn+1LO correction.
Bottom left panel: the median (plus), 68% CI (errorbox) and 95% CI (errorbar)
for the posterior of $(\gamma_{K})_{n+1}^{\text{est}}$ , computed from the
$abc$ (blue) and geometric (red) models using information on the previous
orders. The exact value of $(\gamma_{K})_{n}$ is shown as black circles.
Bottom right panel: CIs scaled to the exact Nn+1LO correction. Figure 6:  The
probability distributions and CIs for the cusp anomalous dimension
$(\gamma_{K})_{n+1}^{\text{est}}$ in planar $\mathcal{N}=4$ SYM for
$\lambda=0.155$ and for different values of $n$. See the caption of figure 6
for details.

In figures 5 and 6 we show the probability distributions obtained from the
$abc$-model for the partial sums

$(\gamma_{K})_{n}:=\sum_{k=1}^{n+1}\gamma_{K}^{(k)}\,\lambda^{k}\,,$ (58)

with $\lambda=0.06$ and $\lambda=0.155$. The series ${\bm{\delta}}_{n}$ is
normalised according to eq. (55).

For $\lambda=0.06<\frac{1}{8}$, the perturbative series (which is alternating)
is expected to converge. In figure 5 we clearly see a convergent pattern in
the left top and bottom panels. The distributions quickly become very peaked
around the limiting value with very small $68\%$ and $95\%$ CIs. In the right
panels, we show the distributions and CIs normalised to the exact Nn+1LO
correction. This effectively cancels the suppression due to increasing power
of the expansion parameter $\lambda$. We see that after such a rescaling, the
distributions do not change significantly for $n>2$ (up to a trivial
reflection around the vertical axis, because of the alternating nature of the
expansion). We note that in contrast to the case of an exact geometric series
(see figure 3), here the posterior distribution for
$\delta_{n+1}^{\text{est}}$ is centred around the half of the exact value
$|\delta_{n+1}|$. As the series progresses, the scaled distribution reaches a
peak for $n=5$. After that it becomes steeper and broader. In the bottom right
panel of figure 5, we see that the 95% CIs shrink more than the 68% CIs. For
comparison we also show the CIs for the geometric model, which are
significantly wider, but are no better in capturing the next order than the
narrower CIs from the $abc$-model. In both cases the 95% CIs appear to give a
reasonable estimate for the size of $\delta_{n+1}$, even though the series is
clearly not a perfect geometric series.

For $\lambda=0.155$ the series is asymptotic and starts to diverge after the
seventh perturbative order, i.e. $n>6$. This behaviour is reflected in the
probability distributions in the left panel of fig. 6: they become peaked
around the limiting value up to the inclusion of the sixth order (the sixth
and seventh order corrections are numerically almost identical). After that a
plateau develops which becomes broader as we increase the perturbative order.
Consequently, the width of the CIs increases steadily after the seventh order
(bottom left panel). However, the rescaled distributions (right panels) look
essentially identical to the convergent $\lambda=0.06$ case. This is not
surprising, because the two series differ only by the value of the expansion
parameter, which is scaled away in the right panels. Therefore, like before
the $95\%$ CIs give a good estimate of the size of the next order corrections,
even in the case of an asymptotic series. We note, however, that for $n>6$ we
enter the regime where the next order is not expected to give a good
approximation of the true value of $\gamma_{K}$. Hence, eq. (56) is not valid,
and so the credibility intervals are not expected to give a good estimate of
the uncertainty attached to the perturbative result.

##### The cusp anomalous dimension in QCD.

In QCD, the cusp anomalous dimension is known up to four loops Brandt:1981kf ;
Korchemsky:1987wg ; Moch:2004pa ; Grozin:2014hna ; Grozin:2015kna ;
Henn:2019swt . It depends on the representation of the Wilson loop and the
matter content of the theory. Here we focus on the fundamental representation,
and we work with 5 massless flavours. In this case we have

$\Gamma^{\rm QCD}_{\rm
cusp}(\mu_{R})=C_{F}\sum_{k=1}^{\infty}\left(\frac{\alpha_{s}(\mu_{R})}{\pi}\right)^{k}\gamma^{(k)}_{\rm
QCD}\,,$ (59)

with $C_{F}=4/3$ and the numerical values of coefficients $\gamma^{(k)}_{\rm
QCD}$ are summarised in table 2. We note that the scale dependence enters
solely through the coupling constant $\alpha_{s}$. Without any context, the
scale is completely arbitrary and there exists no physically motivated choice
for the central scale without specifying the physical process.

Figure 7:  Top left panel: The probability distribution from the $abc$-model
for the QCD cusp anomalous dimension
$(\Gamma^{\text{QCD}}_{\text{cusp}})_{n+1}^{\text{est}}$ evaluated with
$\alpha_{s}(\mu_{R}=5\,\text{GeV})\approx 0.213$ and for different values of
$n$. Top right panel: The same distributions normalized to the exact Nn+1LO
correction. Bottom left panel: the median (plus), 68% CI (errorbox) and 95% CI
(errorbar) for the posterior of
$(\Gamma_{\text{cusp}}^{\text{QCD}})_{n+1}^{\text{est}}$ , computed from the
$abc$ (blue) and geometric (red) models using information on the previous
orders. The exact values of $(\Gamma_{\text{cusp}}^{\text{QCD}})_{n}$ are
shown as black circles. Bottom right panel: CIs scaled to the exact Nn+1LO
correction.

In figure 7 we show (with $\alpha_{s}(\mu_{R}=5\,\text{GeV})\approx 0.213$)
the probability distributions and CIs for the partial sums $(\Gamma^{\rm
QCD}_{\rm
cusp})_{n}=C_{F}\sum_{k=1}^{n+1}\left(\frac{\alpha_{s}}{\pi}\right)^{k}\gamma^{(k)}_{\rm
QCD}$ obtained from the $abc$-model. In the top left panel we see that the
series is rapidly convergent, and the probability distributions become
progressively more peaked. From the rescaled distribution in the top right
panel we see that the probability distribution for $n=2$ peaks at the exact
value of the correction. We do not show scaled results for $n=3$, because the
exact five-loop result is not known yet. In the bottom left panel we show the
68% and 95% CIs for the $abc$-model and the geometric model. We see that,
unlike in the case of planar $\mathcal{N}=4$ SYM, the $68\%$ CIs always
include the next correction, thus giving a reliable estimate of the MHOs (at
least for the first few orders considered here). We emphasise that this is
particularly non-trivial at four loop order, where Casimir scaling is violated
for the first time, and a new quartic Casimir enters. Estimates of the size of
the four-loop corrections based on Padé approximants from the first three
orders were within $5\%$ of the four-loop contributions that admit Casimir
scaling for $n_{f}=5$, but were 57% off from the complete result Moch:2005ba .
The $68\%$ CI obtained from the first three loop orders in the $abc$-model
correctly captures the size of the four-loop result, including the
contribution from the quartic Casimir.

From the bottom right plot we see that the median of the posterior
distribution for the $abc$ model (marked by a plus) is closer to the exact
value than for the geometric model, although the $abc$-model appears to expect
larger higher-order correction, and the distributions are skewed towards
larger corrections. For $n=4$ the CIs are not visible in the bottom left panel
due to the suppression of the expansion parameter. Using the $abc$-model we
estimate that CIs for MHO corrections given the 4-loop order result for
$\mu_{R}=5\,\text{GeV}$, i.e.,
$\Sigma=\Gamma^{\text{QCD}}_{\text{cusp}}-(\Gamma^{\text{QCD}}_{\text{cusp}})_{4}$,
is

$\text{CI}_{68}=C_{F}\left(\frac{\alpha_{s}(\mu_{R})}{\pi}\right)^{5}[2.1,9.5],\quad\text{CI}_{95}=C_{F}\left(\frac{\alpha_{s}(\mu_{R})}{\pi}\right)^{5}[-0.38,21],$
(60)

where we scaled out $C_{F}\left(\frac{\alpha_{s}(\mu_{R})}{\pi}\right)^{5}$
for visibility. We see that the $abc$-model at 68% credibility level would
expect a positive higher-order correction. In contrast, the geometric model
estimates a symmetric interval CI${}_{68}=[-5.4,5.4]$ with the same
normalisation .

###  4.2 On-shell and $\overline{\text{MS}}$-scheme quark masses

Our next example is the relation between the heavy quark mass $m$ in the on-
shell scheme and the mass $\overline{m}$ in the $\overline{\text{MS}}$-scheme.
The bare mass $m_{0}$ is related to the on-shell and $\overline{\text{MS}}$
quark masses via the renormalisation factors $Z_{m}$

$m_{0}=Z_{m}^{\overline{\rm MS}}(\mu_{R})\,\overline{m}(\mu_{R})=Z_{m}^{\rm
OS}\,m\,.$ (61)

For heavy quarks it is possible to perturbatively compute these
renormalisation factors. They can be used to define the scheme-conversion
factor between the two schemes,

$z_{m}(\mu_{R})=\frac{Z_{m}^{\rm OS}}{Z_{m}^{\overline{\rm
MS}}(\mu_{R})}=\frac{\overline{m}(\mu_{R})}{m}=\sum_{k=0}^{\infty}\left(\frac{\alpha_{s}(\mu_{R})}{\pi}\right)^{k}z_{m}^{(k)}(\mu_{R})\;.$
(62)

For heavy quarks $z_{m}(\mu_{R})$ is currently known up to four loops
Melnikov:2000qh ; Marquard:2015qpa ; Marquard:2016dcn . We can express the on-
shell mass as a series calculated in perturbative QFT

$m=z^{-1}_{m}(\mu_{R})\overline{m}(\mu_{R})=\sum_{k=0}^{\infty}m^{(k)},$ (63)

In practice, one does not consider arbitrary $\mu_{R}$ but rather work with a
self-consistent definition of the $\overline{m}$ mass
$\overline{m}(\overline{m})$. Thus for the quark masses we consider the scale
as fixed number rather than as a free parameter.

Figure 8:  Top left panel: The probability distribution from the $abc$-model
for the on-shell top quark mass $(m_{t})^{\text{est}}_{n+1}$ evaluated at
$\mu_{R}=\overline{m}_{t}$ and for different values of $n$. Top right panel:
The same distributions normalised to the exact Nn+1LO correction. Bottom left
panel: the median (plus), 68% CI (errorbox) and 95% CI (errorbar) for the
posterior of $(m_{t})_{n+1}^{\text{est}}$ , computed from the $abc$ (blue) and
geometric (red) models using information on the previous orders. The exact
values of $(m_{t})_{n}$ are shown as black circles. Bottom right panel: CIs
scaled to the exact Nn+1LO correction.

In figure 8 we show the probability distributions and CIs of the on-shell top
quark mass $(m_{t})_{n+1}^{\text{est}}=\sum_{k=0}^{n+1}m^{(k)}$ given the
first $n+1$ perturbative orders of the scheme-conversion factor
$z^{k}_{m_{t}}(\mu_{R})$ evaluated at $\mu_{R}=\overline{m}_{t}$. In the top
left panel we see that the distributions become more and more peaked as
perturbative information becomes available. On the top right we see that the
$abc$-model correctly estimates the positive sign of the corrections, but the
estimated corrections are in general smaller than the exact value. In the
bottom left panel we see that the CIs shrinks rapidly as $n$ increases. In the
bottom right we see that, unlike in the case of the QCD anomalous dimension,
the $68\%$ CIs of both the $abc$ and geometric models do not include the next
perturbative order, although the 95% CIs are sufficiently wide to include it.
We note that the $abc$-model distributions are asymmetric, with a tail towards
positive perturbative corrections showing that, although the $abc$-model does
not peak at the exact value of NnLO contribution, it (correctly) anticipates
that larger positive corrections are more likely. This is a result of the fact
that the perturbative coefficients $z_{m}^{(k)}$ are all negative (at least
through four loops). Using the $abc$-model we estimate the CIs for the MHOs
given the 4-loop result, i.e., $\Sigma=m_{t}-(m_{t})_{4}$, to be (in GeV):

$\text{CI}_{68}=[0.008,0.046]\,,\quad\text{CI}_{95}=[-0.027,0.112]\,.$ (64)

It is interesting to compare the size of the CIs obtained from the $abc$-model
to the range five-loop estimates (in GeV) obtained in ref. Kataev:2018gle :

$\displaystyle(m_{t})_{0}\delta_{5}^{\text{ECH}}$
$\displaystyle=0.073\,,\quad$ $\displaystyle(m_{t})_{0}\delta_{5}^{\text{ECH
direct}}$ $\displaystyle=0.083\,,\quad$
$\displaystyle(m_{t})_{0}\delta_{5}^{\text{FL}}$ $\displaystyle=0.126\,,\quad$
$\displaystyle(m_{t})_{0}\delta_{5}^{\text{FL},M\rightarrow\bar{m}}$
$\displaystyle=0.086\,,\quad,$ $\displaystyle(m_{t})_{0}\delta_{5}^{r-n}$
$\displaystyle=0.112\,.$ (65)

We see that the $68\%$ CI from the $abc$-model are lower than the estimates of
ref. Kataev:2018gle , but the $95\%$ CI includes all values of ref.
Kataev:2018gle , except one.

For completeness, in appendix C, we include plots for the bottom and charm
quarks. The scaled probability distributions and CIs look similar to the ones
in the right panels of figure 8. We note, however, that the MHO estimate using
Bayesian inference is less reliable for lighter quarks, as in this case
corrections become large and the perturbative expansion itself fails. This can
be seen especially in the case of the charm quark mass.

Figure 9:  The sensitivity of the probability distribution
$1/(m_{t})_{0}P_{abc}(((m_{t})_{n+1}^{\text{est}}-(m_{t})_{n})/(m_{t})_{0}|{\bm{\delta}}_{n})$
for on-shell top quark mass for $n=1$ (dashed) and $n=3$ (solid) on the
parameters $(\epsilon,\omega,\eta,\xi)$ for 4 selected values with the others
held fixed at $(\epsilon,\omega,\eta,\xi)=(0.1,1,0.1,2)$.

In figure 9 we show the dependence of the probability distributions for the
$(m_{t})_{n+1}^{\text{est}}$ for $n=1,3$ from the $abc$-model on the
parameters $(\eta,\xi,\epsilon,\omega)$. As already discussed in section 3.3,
we see that the dependence on $(\eta,\epsilon,\omega)$ reduces as more
perturbative information becomes available, while $\xi$ controls the amount of
asymmetry of the distribution. However, after the fourth perturbative order,
the distributions become very narrow and peaked, and the shape is largely
independent of the values of $(\eta,\xi,\epsilon,\omega)$ if $\xi\geq 2$.

###  4.3 Electron anomalous magnetic moment and positronium hyperfine-
splitting

In this section we illustrate how Bayesian techniques work in the context of
two classical results in QED which are known to relatively high order in the
coupling, namely the electron anomalous magnetic moment and the positronium
hyperfine splitting. QED corrections are typically believed to be rapidly
convergent due to the smallness of the QED coupling
$\alpha=1/137.035\,998\,995(85)$ Mohr:2018hvt . Further, these results are
evaluated in the on-shell renormalisation scheme, thus they do not depend on
any unphysical scale.

Let us start by discussing the electron anomalous magnetic moment. The QED
corrections take the form

$\frac{g-2}{2}\equiv
a_{e}=\sum_{k=1}^{\infty}\left(\frac{\alpha}{\pi}\right)^{k}C_{e}^{(2k)}\,.$
(66)

Here we focus only on the mass-independent part of the corrections, which are
known up to 5 loops Schwinger:1948iu ; Petermann:1957hs ; Sommerfield:1957zz ;
Laporta:1996mq ; Aoyama:2012wj ; Aoyama:2014sxa ; Laporta:2017okg ;
Aoyama:2017uqe ; Volkov:2018jhy ; Aoyama:2019ryr . The numerical values for
the coefficients $C_{e}^{2k}$ are collected in table 3.

$k$ | 1 | 2 | 3 | 4 | 5 |
---|---|---|---|---|---|---
$C_{e}^{(2k)}$ | $0.5$ | $-0.328$ | $1.181$ | $-1.912$ | $6.675(192)$ |
$ij$ | $00$ | $10$ | $21$ | $20$ | $32$ | $31$
$C_{\rm HFS}^{(ij)}$ | $0.583$ | $-1.24$ | $-2.06$ | $-3.88$ | $-8.64$ | $14.97$
Table 3:  Known coefficients of mass-independent corrections to electron
anomalous magnetic moment and to positronium hyperfine splitting. Figure 10:
Top left panel: The probability distribution from the $abc$-model for
$(a_{e})_{n+1}^{\text{est}}$ and for different values of $n$. Top right panel:
The same distributions normalised to the exact Nn+1LO correction. Bottom left
panel: the median (plus), 68% CI (errorbox) and 95% CI (errorbar) for the
posterior of $(a_{e})_{n+1}^{\text{est}}$ , computed from the $abc$ (blue) and
geometric (red) models using information on the previous orders. The exact
values of $(a_{e})_{n}$ are shown as black circles. Bottom right panel: CIs
scaled to the exact Nn+1LO correction.

Figure 10 illustrates the $abc$-model estimates for the probability
distributions and the CIs for the partial sums
$(a_{e})_{n+1}^{\text{est}}=\sum_{k=1}^{n+1}\left(\frac{\alpha}{\pi}\right)^{k}C_{e}^{(2k)}$.
Due to the very small values of the expansion parameter, the absolute
distributions and CIs are hardly visible for $n>1$. The distribution for $n=1$
is rather broad, resulting in CIs that overestimate the actual size of the
higher-order corrections. At $n=2$, the 68% CI captures the size of the next
correction, while for $n=3$ only the $95$% CI includes the next order.

Figure 11: The probability distribution and the CIs for the hyperfine-
splitting $\Delta E$ in units of the electron mass $m_{e}$. See the caption of
figure 10.

Our second QED example is the positronium hyperfine-splitting. This is a non-
perturbative quantity and admits a double expansion in $\alpha$ and
$\ln\alpha$:

$\Delta E=m_{e}\alpha^{4}\sum_{i,j=0}^{\infty}C^{(ij)}_{\rm
HFS}\left(\frac{\alpha}{\pi}\right)^{i}\ln^{j}\alpha\,.$ (67)

Only terms with $j<i$ are non-zero with the exception of the coefficient
$C_{00}$. This introduces a natural ordering among the coefficients, which are
ordered first by $i$ then by $j$. For the known terms this ordering is also
equivalent to ordering the contributions by their size, since
$\ln\alpha\sim-5$. The currently known values Karplus:1952wp ; Lepage:1977gd ;
karshenboim1993new ; Czarnecki:1998zv ; Kniehl:2000cx  are shown in table 3.

The presence of the logarithmic corrections in the series clearly violates the
geometric progression of the series assumed in both the geometric and
$abc$-models. Our Bayesian inference for such series is therefore less
reliable as an estimate of the MHOs in comparison to a purely perturbative
quantity such as $a_{e}$. This can be also seen in figure 11: At the lowest
two perturbative orders the CIs overestimate the size of the perturbative
corrections. As the logarithmically enhanced corrections appear, both
geometric and $abc$ models fail to provide reliable estimates of the next-
order correction. This is expected, since the models consider only a single
expansion parameter $a$ and expect the next-order correction to be suppressed
by an additional power of $a$.

In principle, we can consider even more complicated QED bound-state
observables, for example, the bound electron $g$-factor Pachucki:2005px ;
Pachucki:2017xfd ; Czarnecki:2017kva ; Czarnecki:2020kzi . In this case, one
has a triple expansion in powers of $(Z\alpha)$, $\alpha$ and $\ln(Z\alpha)$,
where $Z$ is atomic charge. We expect that such series require more attuned
models that take into account non-perturbative terms and double or triple
expansions. We conclude this section by observing that QED corrections require
dedicated studies within Bayesian framework that we leave for future works.

##  5 Hadronic observables with explicit scale dependence

In this section we discuss the application of Bayesian techniques to estimate
MHOs for a selection of hadronic observables known to NNLO, or even N3LO, in
the strong coupling constant. In all cases, we compute scale-independent
probability distributions and CIs using the scale-marginalisation (sm) and
scale-averaging (sa) prescriptions introduced in sections 2.1.1 and 2.1.2.
Throughout this section we work with the geometric and $abc$-models of section
3. We start by discussing inclusive cross-sections and discuss differential
distributions towards the end of this section. In contrast to the examples
studied in ref. Bonvini:2020xeo , we take into account the variation of the
factorisation scale.

###  5.1 Higgs production in gluon fusion

We start by discussing the inclusive Higgs production cross-section in gluon
fusion:

$(\sigma_{ggH})_{n}(\mu_{F},\mu_{R})=\sum_{k=0}^{n}\sigma^{(k)}_{ggH}(\mu_{F},\mu_{R})\,.$
(68)

The inclusive cross-section is known through NLO in QCD including all finite
top-mass effects Georgi:1977gs ; Dawson:1990zj ; Djouadi:1991tka ;
Graudenz:1992pv ; Spira:1995rr ; Harlander:2005rq ; Aglietti:2006tp ;
Aglietti:2004nj ; Anastasiou:2006hc ; Anastasiou:2009kn , and through N3LO in
an effective theory expansion where the top quark is infinitely heavy
Harlander:2002wh ; Anastasiou:2002yz ; Ravindran:2003um ; Anastasiou:2015ema ;
Anastasiou:2016cez ; Mistlberger:2018etf . All results in this section were
obtained with iHixs 2 Dulat:2018rbf , neglecting finite top-mass effects, with
the `PDF4LHC15_nnlo_100` PDF set and the on-shell top quark mass
$m_{t}=172.5\,{\rm GeV}$. We choose as value for the central scale
$\mu_{0}=m_{H}/2=62.5\,\text{GeV}$ (both for scale-variation results and in
the priors in eqs. (10) and (17)).

First, in figure 12 we show the probability distributions and CIs for
$(\sigma_{ggH})_{n}^{\text{est}}$ for fixed $\mu_{R}=\mu_{F}=\mu_{0}$.

Figure 12:  Top left panel: The probability distributions from the $abc$-model
for the inclusive gluon-fusion Higgs production cross-section
$(\sigma_{ggH})^{\text{est}}_{n+1}$ evaluated at $\mu_{R}=\mu_{F}=m_{H}/2$ and
for different values of $n$. Top right panel: The same distributions
normalised to the exact Nn+1LO correction. Bottom left panel: the median
(plus), 68% CI (errorbox) and 95% CI (errorbar) for the posterior of
$(\sigma_{ggH})_{n+1}^{\text{est}}$, computed from the $abc$ (blue) and
geometric (red) models using information on the previous orders. The exact
values of $(\sigma_{ggH})_{n}$ are shown as black circles. Bottom right panel:
CIs scaled to the exact Nn+1LO correction.

We see that for $n=1,2$, the CIs for both the $abc$ and geometric models
include the next order at 68% credibility level. We note, however, that for
$\mu_{0}=m_{H}/2$ the N3LO corrections are particularly small, and the CIs of
the $abc$-model overshoot the computed value. We will see below that for other
scale choices, the CIs of the $abc$-model are centred around the N3LO values.

Figure 13: The inclusive Higgs production cross-section as a function of the
renormalisation scale (with $\mu_{F}=m_{H}/2$ held fixed) at different orders
$n=0,1,2,3$. The errorbars and erroboxes correspond to the NNLO 95% and 68%
CIs for the geometric model (light green) and the $abc$-model (light blue),
given the $n<3$ results.

Next, in figure 13 we plot the renormalisation scale dependence of
$(\sigma_{ggH})_{n}$ at different orders (for this plot only we keep
$\mu_{F}=\mu_{0}$ fixed). We see that in the conventional scale variation
window, $F=2$, the cross-section depends quite strongly on the scale, except
for the highest order. The N3LO corrections vanish for $\mu_{R}\approx
0.8\mu_{0}$, and the lines representing the NNLO and N3LO results as a
function of the scale cross, i.e., we reach a FAC point. The FAC point at NNLO
is at a smaller value, $\mu_{R}\approx 0.3\mu_{0}$. Also, the crossing points
visible in figure 13 approximately coincide with the local maxima for $n=2,3$,
i.e., the FAC and PMS points are surprisingly close to each other.

In figure 13 we further show the CIs in the geometric and $abc$-models with
LO, NLO and NNLO cross-sections (at fixed scale) used as an input. We observe
that the CIs for the geometric model are centred around the $n=2$ line due to
the symmetry of the geometric model, but the $abc$-model assigns a higher
probability to a positive correction. Indeed, for $\mu_{R}>2\mu_{0}$ the CIs
are almost centred around the N3LO line. For $\mu_{R}<0.8\mu_{0}$, the N3LO
corrections become negative, while the $abc$-model still anticipates a
positive correction. We note that the CIs for the geometric model shrink
dramatically as one approaches the NNLO FAC point at $\mu_{R}=0.3\mu_{0}$,
i.e., large corrections are deemed unlikely by the geometric model close to
that scale. Once the NNLO corrections become negative, the $abc$-model also
becomes centred around the NNLO result, i.e., the $abc$-model no longer
expects a positive correction. The CIs in this region do not shrink as
dramatically as for the geometric model, although generally the $abc$-model
leads to narrower CIs than the geometric model.

####  5.1.1 Sensitivity to the scale interval and scale prescription

We begin by analysing the dependence of our sm- and sa-prescriptions on the
range $F$ that is chosen in eqs. (14) and (18), i.e., $\mu_{0}/F\leq\mu\leq
F\mu_{0}$.

Figure 14: Left: The dependence of the CIs on the size $F$ of the scale
interval ($\mu_{0}/F\leq\mu\leq F\mu_{0}$) for gluon-fusion Higgs cross-
section for the geometric and $abc$-models using the sm-prescription Right:
The same plot using the sa-prescription.

Figure 14 shows the CIs for the values $F=2,\,4,\,10$ using the geometric and
$abc$-models. The left and right panels summarise the results from the sm- and
sa-prescriptions respectively. For $n=2$ the sm-prescription exhibits a
somewhat smaller dependence on the choice of $F$, while the CIs for the sa-
prescription grow with $F$. This can be understood from the discussion in
section 2.2: for symmetric models like the geometric model, the sm-
prescription tries to adapt to the point where the higher-order corrections
are minimised, i.e., the FAC point. Once this point is covered by the range in
the marginalisation, a further increase does not have a substantial impact on
the uncertainties. From figure 13 we see that this is precisely what happens
for the geometric model for $n=2$, when $F$ is increased from 2 to 4, and even
for $F=10$ the CIs change only very little. The discussion of section 2.2 does
not apply to the $abc$-model, which is asymmetric. Recall from figure 13 that
the $abc$-model does not have such dramatic reduction in the size of the CIs
when reaching an FAC point. Therefore it is also less biased to the inclusion
of FAC point into the integration range for $\mu$ in the sm-prescription. For
$n=3$ the second FAC point at very small scales becomes relevant only for the
largest $F$ values.

The sa-prescription, on the other hand, is biased to the regions where the
scale dependence is flat, i.e., it is biased towards the PMS point. However,
for this particular process we observe the peculiar situation that FAC and PMS
points are very close to each other (at least for the orders we considered).
Due to the coherent addition of probabilities, for the sa-prescription
increasing the range of the $\mu$-integration inevitably will also increase
the CIs.

Let us conclude this discussion with an important point. From a Bayesian
perspective, one expects that the predictions of the model become independent
of the priors once enough data have become available. In particular, in the
sm-prescription, the perturbative scales are treated as model parameters, and
so one expects that at high enough orders the model predictions should only
mildly depend on the prior in eq. (10), and the choice of $F$, in agreement
with our findings in figure 14. However, it would be premature to conclude
that the probability distributions are independent of the prior and $F$: we
have shown in section 2.2 that for the geometric model in the sm-prescription,
the preferred scale is the FAC point. The fact that the CIs are insensitive to
the choice of $F$ is likely to be related to the fact that the distributions
will be highly peaked at the values of the cross-section at the FAC point,
which is not necessarily related to the prior independence. Note that the
argument of section 2.2 does not apply to the $abc$-model in the sm-
prescription (because the $abc$-model is not symmetric).

####  5.1.2 Accuracy of Gauss-Legendre quadrature rule

As already discussed in section 2.3, performing the numerical integral over
the scale can become computationally expensive. Quadrature rules, on the other
hand, allow one to approximate the full integral using only a small number of
input points. In particular, a set of points for a 9-point scale variation may
be readily turned into an integral over $(\mu_{R},\mu_{F})$.

Figure 15: Left: The coloured lines show the probability distributions for the
$abc$-model at different orders using the sa-prescription and the 3-point
Gauss-Legendre quadrature on
$\mu_{F}/\mu_{0},\mu_{R}/\mu_{0}\in\\{1/2,1,2\\}$. The black dashed line
indicates the “exact” numerical integration in the interval with $F=2.45$.
Right: The CIs for the geometric and $abc$ models at different orders using
the sm- and sa-prescriptions using the Gauss-Legendre quadrature. The result
of the “exact” numerical integral are overlaid using black dashed lines. For
comparison we also include the 9-point and 7-point variation intervals.

The left panel of figure 15 illustrates the accuracy of this approximation for
the probability distributions in the $abc$-model using the sa-prescription.
The distributions for different values of $n$ (in colour) are obtained using
the 3-point Gauss-Legendre quadrature for both $\mu_{R}$ and $\mu_{F}$. The
results obtained by a finely-spaced equidistant numerical integration over the
same scale interval with $F\approx 2.45$ is overlaid on top (black dashed
lines). We see that, except for some features around the peak, the quadrature
rules reproduce the “exact” numerical integration very well. In particular,
this approximation has very little effect on the CIs, which are shown in the
right panel of figure 15, where we show the 68% and 95% CIs for the geometric
and $abc$ models using both the sm- and sa-prescriptions. The dashed black
lines indicate the results of the finely-spaced numerical integration. We note
that after scale integration the probability distributions are no longer
symmetric, even if the distributions at fixed scale are. Overall, we observe
that the Gauss–Legendre quadrature rule is able to reproduce both the 68% and
95% CIs from the “exact” numerical integration to a very good degree.
Therefore, in the subsequent sections we will only show the results obtained
using the Gauss-Legendre quadrature.

For Higgs production in gluon fusion both the sm- and sa-prescriptions produce
very similar CIs shown in the right panel of figure 15. The CIs are also very
similar in size for both the geometric and the $abc$-models. The results for
the latter are, however, systematically shifted upwards due to the sequence of
positive corrections at lower orders. For comparison we also display the scale
variation intervals using 9- and 7-point rules, which have no probabilistic
interpretation. We observe that the scale variation intervals are slightly
smaller, but rather close, to the 68% CIs obtained from the Bayesian methods.
This illustrates the arguments of section 2.3 in a practical example. In
particular, thanks to the slow $\mu_{F}$ variation, the scale variation
intervals are close to saturating the $s=1$ bound in eq. (34).

Finally, we conclude with the summary of 68% and 95% CIs and scale-variation
intervals for the highest known order $n=3$ of $(\sigma_{ggH})_{n}$:

model | prescription |  $\text{CI}_{68}$ (pb) |  $\text{CI}_{95}$ (pb) | 7-point (pb) | 9-point (pb)
---|---|---|---|---|---
$abc$ | sa | [47.8, 51.2] | [44.4, 56.5] | [46.3, 48.3] | [46.0, 48.3]
$abc$ | sm | [48.3, 51.1] | [46.1, 55.8]
geo | sa | [46.0, 49.2] | [41.6, 53.0]
geo | sm | [47.0, 49.2] | [42.7, 53.0]

We note that as the N3LO cross-section reaches a local maximum near the
central scale $\mu_{0}$, traditional scale variation produces a one-sided
interval. The geometric model has approximately symmetric intervals around the
central value, while the $abc$-model anticipates positive MHO terms. Although
all 68% CIs and the scale-variation intervals are numerically close in size,
we stress that they have completely different statistical meanings. Scale-
variation intervals have no probabilistic interpretation and cannot be
associated with a particular level of credibility. In contrast, the Bayesian
CIs express a mathematically rigorous degree of belief for the MHOs from the
posterior distribution for a given model and scale prescription.

###  5.2 Vector-Boson Fusion Higgs and di-Higgs production

Our next two examples are Higgs production in vector-boson fusion (VBF), and a
similar process for double Higgs production. The inclusive VBF cross-section
is known up to N3LO. We denote VBF cross-sections for single Higgs
Bolzoni:2010xr ; Bolzoni:2011cu ; Dreyer:2016oyx ; Dreyer:2020urf  and double-
Higgs production Liu-Sheng:2014gxa ; Dreyer:2018qbw ; Dreyer:2020urf  by
($0\leq n\leq 3$)

$\begin{split}(\sigma_{\rm
VBF-H})_{n}(\mu_{F},\mu_{R})&\,=\sum_{k=0}^{n}\sigma^{(k)}_{\rm
VBF-H}(\mu_{F},\mu_{R})\,,\\\ (\sigma_{\rm VBF-
HH})_{n}(\mu_{F},\mu_{R})&\,=\sum_{k=0}^{n}\sigma^{(k)}_{\rm VBF-
HH}(\mu_{F},\mu_{R})\,.\end{split}$ (69)

As before the centre-of-mass energy is $\sqrt{s}=13\,\rm TeV$. The central
scale is given by the vector boson momentum Han:1992hr  and we take into
account the dependence on both factorisation and renormalisation scales.
Computations were performed with the proVBFH code proVBFH .

Figure 16: The 68% and 95% CIs for the VBF cross-sections for Higgs and di-
Higgs production for the geometric and $abc$-models using the sa- and sm-
prescriptions. The scale variation intervals using 7 and 9 points are shown
for comparison.

In the left panel of figure 16 we display the CIs for different models and
prescriptions for single Higgs VBF production. For $n<2$ the Bayesian approach
gives a larger uncertainty (68% CIs) than the traditional scale variation.
Because the NLO correction is negative, the $abc$-model anticipates an
alternating series, and consequently the CIs for $n=1$ for the $abc$-model are
positively shifted compared to the NLO result. However, the NNLO corrections
are again negative, and for $n=2$ all studied models and prescriptions give
very similar 68% CIs, although the $abc$-model has much larger 95% CIs than
the geometric model. For $n=3$ the 68% CIs shrink even further and become
somewhat smaller than the scale variation intervals. For the single Higgs VBF
cross-section $(\sigma_{\text{VBF-H}})_{n}$ at $n=3$ these CIs are:

model | prescription |  $\text{CI}_{68}$ (fb) |  $\text{CI}_{95}$ (fb) | 7 point (fb) | 9 point (fb)
---|---|---|---|---|---
$abc$ | sa | [3.9306, 3.9357] | [3.9287, 3.9478] | [3.9304, 3.9367] | [3.9304, 3.9367]
$abc$ | sm | [3.9304, 3.9337] | [3.9290, 3.9430]
geo | sa | [3.9305, 3.9343] | [3.9287, 3.9385]
geo | sm | [3.9304, 3.9324] | [3.9293, 3.9355]

We note that the sm-prescription gives much smaller CIs than the sa-
prescription. In fact, the 95% CIs of the scale-marginalised geometric model
is smaller and does not contain the scale-variation interval, demonstrating
that the bounds discussed in section 2.3 do not necessarily apply to the sm-
prescription. In contrast the 95% CIs for the geometric model in sa-
prescription contain the scale variation intervals, as expected.

In the right panel of figure 16 we display the CIs for different models and
prescriptions for di-Higgs VBF production. We observe very good convergence of
the cross-section, and correspondingly the CIs from Bayesian inference shrink
rapidly. We observe that the sa-prescription gives larger CIs than the sm-
prescription, which is due to the presence of an FAC point in the scale
integration interval. We also see a significant difference between the 7- and
9-point scale variation intervals. The 9-point interval for $n=2$ is larger
than the 68% CI for the geometric model with the sa-prescription, but still
within the 95% CI. This is a sign that both renormalisation and factorisation
scale dependencies are comparable and the one dimensional ($s=1$) bound
estimated in section 2.3 is not applicable in this case. For $n=2$ and $n=3$
the 95% CI for scale-marginalised geometric model is even smaller than 9-point
scale variation interval, while sa-prescription 95% still encompasses it, as
expected. The CIs for di-Higgs VBF cross-section $(\sigma_{\text{VBF-
HH}})_{n}$ at $n=3$ are:

model | prescription |  $\text{CI}_{68}$ (fb) |  $\text{CI}_{95}$ (fb) | 7 point (fb) | 9 point (fb)
---|---|---|---|---|---
$abc$ | sa | [1.7237, 1.7248] | [1.7229, 1.7280] | [1.7232, 1.7247] | [1.7232, 1.7258]
$abc$ | sm | [1.7241, 1.7248] | [1.7238, 1.7252]
geo | sa | [1.7237, 1.7247] | [1.7229, 1.7263]
geo | sm | [1.7239, 1.7247] | [1.7232, 1.7249]

We see that sa-prescription gives generally more conservative estimates for
the MHOs than either the sm-prescription or scale variation.

###  5.3 Drell-Yan processes

Drell-Yan-type processes play a crucial role in hadron collider phenomenology.
The inclusive cross-section for off-shell photon production with virtuality
$Q^{2}$, $pp\to\gamma^{(*)}$, is known though N3LO Altarelli:1978id ;
Hamberg:1990np ; Duhr:2020seh :

$\left(\sigma_{\rm DY-
NC}\right)_{n}(Q^{2},\mu_{F},\mu_{R})=\sum_{k=0}^{n}\sigma^{(k)}_{\rm DY-
NC}(Q^{2},\mu_{F},\mu_{R})\,.$ (70)

N3LO corrections are also known for the charged-current Drell-Yan process
$pp\to W^{\pm}\to\ell^{\pm}\nu_{\ell}$ Altarelli:1979ub ; Aurenche:1980tp ;
KubarAndre:1978uy ; Hamberg:1990np ; Anastasiou:2003ds ; Melnikov:2006di ;
Duhr:2020sdp . We focus here on the lepton charge asymmetry defined as

$A_{W}(Q^{2})=\frac{\frac{d\sigma_{\rm DY-W^{+}}}{dQ^{2}}-\frac{d\sigma_{\rm
DY-W^{-}}}{dQ^{2}}}{\frac{d\sigma_{\rm DY-W^{+}}}{dQ^{2}}+\frac{d\sigma_{\rm
DY-W^{-}}}{dQ^{2}}}\,,$ (71)

and the $n$-th order approximation is

$\left(A_{W}\right)_{n}(Q^{2},\mu_{F},\mu_{R})=\sum_{k=0}^{n}A_{W}^{(k)}(Q^{2},\mu_{F},\mu_{R})\,.$
(72)

For both processes we only consider here $Q=m_{W}$ with the central scale set
to $\mu_{0}=Q=m_{W}$.

In the left panel of figure 17 we show results for the neutral current Drell-
Yan cross-section. We see that for $n=2,3$ the sa-prescription produces 68%
CIs which are similar in size to the regular scale variation intervals. The
CIs from the sm-prescription for $n=2$ are much smaller, again due to the
presence of an FAC point in the integration region for $\mu$. The N3LO
corrections are known to be sizeable, and they are not covered by the
conventional 7-point scale variation at $n=2$ Duhr:2020seh . The CIs for the
neutral-current Drell Yan cross-section $(\sigma_{\text{DY-NC}})_{n}$ at $n=3$
are:

model | prescription |  $\text{CI}_{68}$ (nb) |  $\text{CI}_{95}$ (nb) | 7 point(nb) | 9 point (nb)
---|---|---|---|---|---
$abc$ | sa | [45.6, 46.6] | [44.8, 49.0] | [45.6, 46.4] | [45.5, 46.4]
$abc$ | sm | [45.9, 46.5] | [45.1, 48.3]
geo | sa | [45.5, 46.4] | [44.6, 47.2]
geo | sm | [45.8, 46.3] | [45.0, 46.9]

We observe that the 68$\%$ CIs are similar in size among themselves, and to
the scale-variation intervals, but the CIs from the $abc$-model are slightly
shifted upwards in the anticipation of a positive MHO correction.

Figure 17: The 68% and 95% CIs for the neutral-current Drell-Yan cross-section
and charged-current lepton-charge asymmetry for the geometric and $abc$-models
using the sa- and sm-prescriptions. The scale variation intervals using 7 and
9 points are shown for comparison.

In the right panel of figure 17 we show results for the lepton charge
asymmetry for $\mu_{0}=Q=m_{W}$. The perturbative expansion for
$A_{W}(m_{W}^{2})$ is quickly convergent with only a mild scale dependence,
because some corrections cancel in the ratio. The perturbative coefficients
feature a monotonic increase with the perturbative order, and the $abc$-model
correctly anticipates positive contributions from MHOs. The CIs from the
$abc$-model are slightly smaller than for the geometric model. We do not
observe significant differences between the sm- and sa-prescriptions, except
for $n=3$, where scale-marginalisation gives more aggressive CIs. We note that
the traditional 7-point scale variation intervals for $n=0,1$ fail to include
the next correction, but for $n=2$ they are similar to the 68% CIs obtained
from Bayesian inference. The results for the CIs for $A_{W}(m_{W}^{2})$ at
N3LO are:

model | prescription | $\text{CI}_{68}$ | $\text{CI}_{95}$ | 7 point | 9 point
---|---|---|---|---|---
$abc$ | sa | [0.1489, 0.1494] | [0.1486, 0.1497] | [0.1488, 0.1493] | [0.1485, 0.1494]
$abc$ | sm | [0.1491, 0.1494] | [0.1488, 0.1497]
geo | sa | [0.1487, 0.1493] | [0.1484, 0.1495]
geo | sm | [0.1490, 0.1494] | [0.1485, 0.1495]

We observe that the sm-prescription produces smaller CIs than the sa-
prescription, indicating the presence of FAC point, which also pulls the CIs
upwards with respect to the N3LO result. Scale-averaged CIs are more
conservative and more centred around the N3LO results for both the geometric
and $abc$-models.

###  5.4 Deep inelastic scattering

The Bayesian procedure of estimating MHOs can be readily applied to
differential spectra by considering individual points/bins of the associated
observable.777Note that here no correlations are considered, i.e., cross-
section entries in each bin of the distribution are assumed to be independent.
The incorporation of correlations is left for future studies. In the
following, we will investigate how the geometric and $abc$-models perform for
various differential predictions. Results will be further contrasted between
the sm- and sa-prescriptions, and compared to the 9-point scale variation
intervals. Like for inclusive observables, we will perform scale integrals
using the Gauss–Legendre method from section 2.3. This setup allows us to
accommodate already existing results based on 7- or 9-point scale variation,
without the need for a full re-computation of these predictions. In this
section we begin our discussion of differential spectra with the DIS process.
We will present results for several benchmark processes at the LHC in the
subsequent sections.

Figure 18: Deep inelastic scattering distribution $d\sigma/dQ^{2}$ up to
N${}^{\text{3}}$LO with $\mu_{0}=Q$ as the central scale choice. Left:
absolute distributions and 9-point scale variation envelopes normalised to the
NLO prediction at $\mu_{0}$. Right: 68% CIs for the $abc$ and geometric models
using the scale-marginalisation (top) and scale-averaging (bottom)
prescriptions (also normalised to NLO).

Predictions up to N${}^{\text{3}}$LO for the DIS process were obtained from
the calculation in ref. Currie:2018fgr , where we consider the distribution
$d\sigma/dQ^{2}$ with respect to the virtuality $Q^{2}$ of the intermediate
photon. The results are presented in figure 18. The left panels show the
absolute distributions (top) and the standard scale variation envelopes
normalised by the NLO prediction evaluated at the central scale (bottom). The
four right panels contrast the MHO estimates (68% CIs) obtained using the
$abc$ and geometric models with the sa- (top two panels) and sm-prescriptions
(bottom two panels). All CIs are normalised to the NLO prediction evaluated at
the central scale.888For clarity, we do not show 68% CIs at LO, nor any of the
95% CIs, which can be very large. Bands of different colours correspond to
estimates at different perturbative orders: grey (LO), blue (NLO), red (NNLO)
and dark grey (N3LO). Overall, we can observe that all models and
prescriptions perform well in capturing the progression of the perturbative
series with either overlapping or touching uncertainty bands. At NLO the
Bayesian approaches give rise to 68% CIs which are slightly larger than the
9-point scale variation interval, but are comparable in size at NNLO and
N${}^{\text{3}}$LO. In the case of the sm-prescription shown in the top two
panels on the right in figure 18, the CI bands narrow dramatically at NLO for
$Q^{2}\sim 400~{}\mathrm{GeV}$, as well as at NNLO and N${}^{\text{3}}$LO at
the lower end of the spectrum. This can be traced back to the presence of a
FAC points that are strongly emphasised by the sm-prescription. In the case of
the geometric model, we see that both scale prescriptions move the median of
the NLO probability distribution towards the direction of the higher-order
corrections, which slightly improves the overlap between the CIs at NNLO and
N3LO. For the $abc$-model, on the other hand, the central value at NLO moves
in the opposite direction, with strongly asymmetric uncertainty bands that
extend further towards the LO prediction. This is because after the negative
NLO correction, the $abc$-model anticipates a positive NNLO correction, as
would be the case for an alternating series. However, for this process the
NNLO correction is again negative. The geometric model, which ignores the
signs of the corrections, appears to perform better. We note that at higher
orders the $abc$-model captures the perturbative progression well with only
slightly larger CIs than the geometric model.

###  5.5 Di-photon production

We now consider predictions for di-photon production at the LHC at
$\sqrt{s}=8~{}\mathrm{TeV}$ using the results of ref. Gehrmann:2020oec . The
production of a pair of photons, ${p}\,{p}\to\gamma\gamma$, provides a clean
and well-measured final state that is an important background in Higgs
measurements and can be used in the search for New Physics resonances. One
striking feature of this process is the apparent slow convergence of the
perturbative QCD series, as partially induced by the staggered transverse-
momentum cuts on the photons. As a consequence, lower-order predictions are
found to be inadequate to describe the measurement, and only starting from
NNLO a sensible comparison to the data can be performed.

Figure 19: Di-photon production up to NNLO with $\mu_{0}=m_{\gamma\gamma}$ as
the central scale choice. Left: absolute $d\sigma/dm_{\gamma\gamma}$
distributions and 9-point scale variation envelopes normalised to the NLO
prediction at $\mu_{0}$. Right: 68% CIs for the $abc$ and geometric models
using the scale-marginalisation (top) and scale-averaging (bottom)
prescriptions (also normalised to NLO).

Results for the di-photon invariant-mass distribution
$d\sigma/dm_{\gamma\gamma}$ are shown in figure 19 following the same pattern
as in figure 18. We note that the sizeable higher-order corrections put the
respective predictions far outside from the traditional scale variation
intervals of the previous order. Moreover, the scale variation envelope at
NNLO is found to be of similar size, or even larger, than at NLO, with no
convincing sign of convergence. The Bayesian models, on the other hand,
produce 68% CIs that are generally larger, and the CIs do not increase when
going to higher orders. With the large positive corrections prohibiting any
FAC points, the two scale prescriptions produce CIs that are almost identical.
While the CI bands between NLO and NNLO touch in the case of the geometric
model, for the $abc$-model the NNLO CIs are outside of the NLO band for
$m_{\gamma\gamma}\gtrsim 200~{}\mathrm{GeV}$. This is likely due to large
positive corrections encountered at each perturbative order. The $abc$-model
anticipates this pattern to continue at higher orders and shifts the
probability distribution towards larger cross-section values. Finally, in view
of the noticeable differences for the size and position of the bands at the
highest order (NNLO) using different approaches, we can conclude that
N${}^{\text{3}}$LO corrections will be highly relevant for this process to
obtain more robust MHO estimates, which the traditional scale variation is
likely substantially underestimating.

###  5.6 Gauge boson production in association with a jet

The production of electroweak gauge bosons that recoil against hard QCD
emissions is among the most important standard candle processes at the LHC
with a wide range of applications that span from detector calibration,
precision QCD studies, to New Physics searches. This process is now known up
to NNLO for all gauge boson types. Here we focus on $W^{+}$+jet production
using the predictions from ref. Gehrmann-DeRidder:2017mvr .

Figure 20: $W^{+}$+jet production up to NNLO with $\mu_{0}=E_{T,W}$ as the
central scale choice. Left: absolute $d\sigma/dp_{T,W}$ distributions and
9-point scale variation envelopes normalised to the NLO prediction at
$\mu_{0}$. Right: 68% CIs for the $abc$ and geometric models using the scale-
marginalisation (top) and scale-averaging (bottom) prescriptions (also
normalised to NLO).

Figure 20 presents the results for the inclusive transverse momentum
distribution $d\sigma/dp_{T,W}$ of the $W$ boson following the same layout as
in figure 18. We can observe a good perturbative behaviour for this process
with NNLO corrections that are well captured by the traditional scale
variation at one order lower. The Bayesian approaches give rise to 68% CIs
that are wider at NLO, while the size of the bands at NNLO are comparable
between all approaches. The sm- and sa-prescriptions give rise to almost
identical CIs, similar to the observations made for the di-gamma process. In
the case of the $abc$-model, we can appreciate the flexibility of the
asymmetric probability distribution. The positive offset at NLO centres the
CIs at this order around the NNLO band. Finally, we note that once NNLO
corrections are included, all the different prescriptions of estimating MHOs
considered here give almost identical bands indicating that the uncertainties
assigned at this order are very robust.

###  5.7 Higgs boson production in association with a jet

The study of the Higgs boson and its properties is among the highest
priorities of the LHC programme and the transverse momentum distribution, in
particular, plays an important role in this endeavour not only for boosted
analyses but also to resolve the details of the production dynamics for this
process. We focus here on Higgs boson production in association with a jet at
$\sqrt{s}=13~{}\mathrm{TeV}$ in the four-lepton decay channel, $H\to 4\ell$,
and use the predictions for the “ATLAS II” setup computed in ref. Chen:2019wxf
.

Figure 21: Higgs+jet production up to NNLO with $\mu_{0}=\tfrac{1}{2}E_{T,H}$
as the central scale choice. Left: absolute $d\sigma/dp_{T,H}$ distributions
and 9-point scale variation envelopes normalised to the NLO prediction at
$\mu_{0}$. Right: 68% CIs for the $abc$ and geometric models using the scale-
marginalisation (top) and scale-averaging (bottom) prescriptions (also
normalised to NLO).

The results up to NNLO for the transverse momentum distribution
$d\sigma/dp_{T,H}$ of the Higgs boson are shown in figure 21. We observe large
positive corrections at each perturbative order with the exception of the
lowest $p_{T}$ bins. In contrast to the other differential distributions
considered so far, we notice that the Bayesian models are not giving rise to
substantially larger uncertainty bands (68% CI) at NLO compared to the
traditional scale variation envelope, but they are similar in size. Going to
NNLO, the uncertainty estimates are found to be comparable between the
different approaches, with the sm-prescription giving somewhat narrower bands.
This compatibility between the different approaches again supports the
robustness of the uncertainties assigned to these predictions. With the
monotonic pattern of positive corrections at each order, the $abc$-model is
again able to adapt to this feature in the intermediate–high $p_{T}$ range,
capturing to a large extent the NNLO corrections in the NLO shift. This is
more pronounced in the sm-prescription suggesting that the central scale
choice of $\mu_{0}=\tfrac{1}{2}E_{T,H}$ is probing regimes in the scale
variation that are close a FAC point.

###  5.8 Inclusive jet production

Jet production is the most prominent process at hadron colliders and an
essential tool for the study of QCD. The inclusive transverse momentum
distribution $d\sigma/dp_{T,j}$ is particularly interesting in the context of
theory uncertainties, due to its sensitivity on the specific choice of the
central scale that further persists at NNLO due to instabilities in the sub-
leading jet distribution that impacts the inclusive spectrum Currie:2018xkj .
In order to gain insights into this sensitivity and how the different
prescriptions for assigning MHO uncertainties compare in these scenarios, we
consider two different central scale choices: the jet transverse momentum
$p_{T,j}$ and the scalar sum of the transverse momenta of the partons
$\hat{H}_{T}=\sum_{i\in\text{partons}}p_{T,i}$, which are shown in figures 22
and 23, respectively.

Figure 22: Inclusive jet production up to NNLO with $\mu_{0}=p_{T,j}$ as the
central scale choice. Left: absolute $d\sigma/dp_{T,j}$ distributions and
9-point scale variation envelopes normalised to the NLO prediction at
$\mu_{0}$. Right: 68% CIs for the $abc$ and geometric models using the scale-
marginalisation (top) and scale-averaging (bottom) prescriptions (also
normalised to NLO).  Figure 23: Inclusive jet production up to NNLO with
$\mu_{0}=\hat{H}_{T}$ as the central scale choice. Left: absolute
$d\sigma/dp_{T,j}$ distributions and 9-point scale variation envelopes
normalised to the NLO prediction at $\mu_{0}$. Right: 68% CIs for the $abc$
and geometric models using the scale-marginalisation (top) and scale-averaging
(bottom) prescriptions (also normalised to NLO).

We begin our discussion with the scale choice $\mu_{0}=p_{T,j}$, which makes
the NLO corrections artificially small. From the results presented in figure
22, we see that in the low-$p_{T}$ region all approaches give rise to sizeable
NNLO CIs. This is most clearly visible for the $abc$-model, where no reduction
in the size of the CI band is observed as we move from NLO to NNLO. Neither of
the bands overlap, because the $abc$-model expects a positive NNLO correction,
which is not true for this central scale choice. The geometric model, instead,
is able to assign more conservative MHO uncertainties at NLO that allow for
some overlap between the last two orders. Overall, we can conclude that the
Bayesian models considered here are not able to “un-do” an inappropriate scale
choice through either the sm- or sa-prescriptions.

More specifically, no prescription based on a simple re-scaling from a central
scale choice will be able to adapt between different _dynamical_ scales that
encapsulate a different kinematic dependence. This becomes most apparent in
observables that are sensitive to infrared physics, because it is the
cancellation between real and virtual corrections that strongly impacts such a
sensitivity. For two dynamical scale choices that have a different dependence
on the kinematics, the difference in the scales assigned to the real and
virtual corrections can vary substantially. Such a feature can never be
captured by a global scaling, as done in the sa- and sm-prescriptions
considered here, that varies the scales for the real and virtual contributions
in the same manner.

Using $\hat{H}_{T}$ as the central scale choice as advocated in ref.
Currie:2018xkj , we observe in figure 23 a stable progression of the series
throughout with the next order, typically lying within the uncertainty
estimate of the previous order. The 68% CIs are slightly larger for the
Bayesian approaches at NLO, while at NNLO all uncertainty estimates are
mutually compatible. This further supports this scale as an appropriate choice
given its robustness with respect to the different prescriptions in assigning
theory uncertainties. Finally, it is interesting to note that the $abc$-model
is able to adapt to the sizeable positive corrections in the tail of the
distribution, thus giving a noticeable positive shift at NLO that reduces the
overall impact of NNLO corrections.

##  6 Discussion

The estimation of MHOs is a vital aspect of collider phenomenology.
Traditional approaches are based on scale variation, an ad hoc prescription to
obtain uncertainty intervals for MHOs without a clear probabilistic
underpinning. In this paper, we have analysed Bayesian approaches to estimate
the MHO uncertainties from the progression of the perturbative expansion of a
physical observable. Building on the recent work in ref. Bonvini:2020xeo , we
have focused in particular on the question of how the dependence on the