# Effects of Renormalon Scheme and Perturbative Scale Choices
on Determinations of the Strong Coupling from $e^{+}e^{-}$ Event Shapes

Guido Bell<EMAIL_ADDRESS>Theoretische Physik 1, Center for
Particle Physics Siegen, Universität Siegen,
Walter-Flex-Strasse 3, 57068 Siegen, Germany  Christopher Lee<EMAIL_ADDRESS>Theoretical Division, Los Alamos National Laboratory, P.O. Box 1663, MS B283,
Los Alamos, NM 87545, USA  Yiannis Makris  Theoretical Division, Los Alamos
National Laboratory, P.O. Box 1663, MS B283,
Los Alamos, NM 87545, USA  INFN Sezione di Pavia, via Bassi 6, I-27100 Pavia,
Italy  Jim Talbert<EMAIL_ADDRESS>DAMTP, University of Cambridge,
Wilberforce Road, Cambridge, CB3 0WA, United Kingdom  Theoretical Division,
Los Alamos National Laboratory, P.O. Box 1663, MS B283,
Los Alamos, NM 87545, USA  Bin Yan<EMAIL_ADDRESS>Theoretical Division,
Los Alamos National Laboratory, P.O. Box 1663, MS B283,
Los Alamos, NM 87545, USA  Institute of High Energy Physics, Chinese Academy
of Sciences, Beijing 100049, China

###### Abstract

We study the role of renormalon cancellation schemes and perturbative scale
choices in extractions of the strong coupling constant $\alpha_{s}(m_{Z})$ and
the leading non-perturbative shift parameter $\Omega_{1}$ from resummed
predictions of the $e^{+}e^{-}$ event shape thrust. We calculate the thrust
distribution to N3LL′ resummed accuracy in Soft-Collinear Effective Theory
(SCET) matched to the fixed-order $\mathcal{O}(\alpha_{s}^{2})$ prediction,
and perform a new high-statistics computation of the
$\mathcal{O}(\alpha_{s}^{3})$ matching in EERAD3, although we do not include
the latter in our final $\alpha_{s}$ fits due to some observed systematics
that require further investigation. We are primarily interested in testing the
phenomenological impact sourced from varying amongst three renormalon
cancellation schemes and two sets of perturbative scale profile choices. We
then perform a global fit to available data spanning center-of-mass energies
between 35-207 GeV in each scenario. Relevant subsets of our results are
consistent with prior SCET-based extractions of $\alpha_{s}(m_{Z})$, but we
are also led to a number of novel observations. Notably, we find that the
combined effect of altering the renormalon cancellation scheme and profile
parameters can lead to few-percent-level impacts on the extracted values in
the $\alpha_{s}-\Omega_{1}$ plane, indicating a potentially important
systematic theory uncertainty that should be accounted for. We also observe
that fits performed over windows dominated by dijet events are typically of a
higher quality than those that extend into the far tails of the distributions,
possibly motivating future fits focused more heavily in this region. Finally,
we discuss how different estimates of the three-loop soft matching coefficient
$c_{\tilde{S}}^{3}$ can also lead to measurable changes in the fitted
$\\{\alpha_{s},\Omega_{1}\\}$ values.

††preprint: LA-UR-23-26138††preprint: SI-HEP-2023-24††preprint: P3H-23-079

##  I Introduction

Electron-positron event shapes are amongst the oldest and most established
observables for testing the predictions of perturbative QCD Dasgupta and Salam
(2004). They typically measure geometric configurations of hadronic final-
state momentum distributions in detectors which, thanks to the absence of
initial-state color-charged particles, are free of many uncertainties
associated to hadron-hadron collisions. They are therefore generally assumed
to be theoretically clean observables, permitting high-order perturbative
calculations as well as precision extractions of the fundamental parameter of
massless QCD—the strong coupling constant $\alpha_{s}$.

It is to this end that multiple experimental and theoretical collaborations
have pursued $e^{+}e^{-}$ event shapes, utilizing available data from the
Large Electron-Positron collider (LEP) and other experiments. Some of these
results Bethke  _et al._ (2009); Dissertori  _et al._ (2009, 2010); Abbate
_et al._ (2011); Abbiendi  _et al._ (2011); Schieck  _et al._ (2013); Gehrmann
_et al._ (2013); Hoang  _et al._ (2015a); Kardos  _et al._ (2018); Verbytskyi
_et al._ (2019) contribute along with analogous studies utilizing lattice
simulations, $\tau$-decays, PDF fits, heavy quarkonia decays, etc. to the 2022
Particle Data Group (PDG) world average, which is found to be
$\alpha_{s}(m_{Z})=0.1179\pm 0.0009$ Workman and Others (2022). Thanks to the
desirable properties of $e^{+}e^{-}$ event shapes and the associated precision
of theoretical predictions now available—which exist up to the N3LL′ resummed
perturbative accuracy using Soft-Collinear Effective Theory (SCET) Bauer  _et
al._ (2000, 2001, 2002); Beneke  _et al._ (2002) and
$\mathcal{O}(\alpha_{s}^{3})$ accuracy in fixed-order QCD Gehrmann-De Ridder
_et al._ (2007); Weinzierl (2009)—determinations based on observables such as
thrust Abbate  _et al._ (2011) and $C$-parameter Hoang  _et al._ (2015a) are
among the most precise quoted in the PDG.

One notices, however, that these extractions tend to yield significantly lower
values of $\alpha_{s}$ as compared to the PDG world average, which motivates
further investigation. It is generally observed in the literature that lower
values of $\alpha_{s}$ emerge from event-shape predictions that account both
for the resummation of logarithms in perturbation theory and a semi-analytic
description of non-perturbative power corrections from hadronization.
Recently, this discrepancy has motivated studies of models of power
corrections in the three-jet region of event-shape distributions Luisoni  _et
al._ (2021); Caola  _et al._ (2022a, b) and associated impacts on the
determination of $\alpha_{s}$ Luisoni  _et al._ (2021); Nason and Zanderighi
(2023), indicating additional systematic uncertainties from such effects. In
the current work, we follow a somewhat orthogonal approach that addresses the
question whether there exist theoretical uncertainties in the _dijet_
factorized prediction for the thrust distribution that have a similarly non-
negligible impact on the $\alpha_{s}$ determination. We stress that the
purpose of this work is to identify such systematic effects that remain to be
better understood or controlled in future fits, and that we do not aim at a
competitive $\alpha_{s}$ extraction here.

While multiple explanations for this tension could exist, a definitive answer
will at the very least require a better understanding of the impact of (A)
non-perturbative (NP) physics effects and (B) perturbative scale uncertainties
on the error estimation of the extracted results. For example, regarding
${\bf{(A)}}$, the SCET thrust and $C$-parameter fit results mentioned above
were performed in a two-dimensional plane in $\alpha_{s}$ and $\Omega_{1}$,
where the latter is a universal parameter encoding the leading NP shift in the
tail regions of the differential distributions Webber (1994); Dokshitzer and
Webber (1995); Korchemsky (1998); Korchemsky and Sterman (1999); Lee and
Sterman (2007):

$\frac{d\sigma}{de}(e)\underset{\text{NP}}{\longrightarrow}\frac{d\sigma}{de}\Big{(}e-c_{e}\frac{\Omega_{1}}{Q}\Big{)}\,,$
(1)

where $e$ denotes a generic event shape whose dominant power corrections come
from the soft sector, and where $c_{e}$ is a calculable, observable-dependent
scaling coefficient given by $c_{\tau}=2$ for thrust and $c_{C}=3\pi$ for
$C$-parameter. While (1) is only valid in a region of $e$ large enough to
justify an operator product expansion of the soft shape function that
describes the NP physics, its implication for fitting analyses is evident:
shifts in theory distributions due to differing values of $\alpha_{s}$ can
partly be compensated by shifts from differing values of $\Omega_{1}$. As a
result, extractions are often given in terms of ellipses in the combined
$\alpha_{s}-\Omega_{1}$ plane, whose sizes correspond to a given statistical
confidence level and whose slopes encode the correlation between the two fit
parameters $\alpha_{s}$ and $\Omega_{1}$. Furthermore, often embedded in the
definition of $\Omega_{1}$ is an implicit scheme for simultaneously cancelling
the leading soft infrared renormalon appearing in the perturbative dijet soft
function and certain hadronization model parameter(s) Hoang and Stewart
(2008). This renormalon cancellation scheme amounts to a choice, and as
recently argued in Bachu  _et al._ (2021), many reasonable choices exist. We
stress that when we make the same choices as in prior SCET-based thrust
analyses Abbate  _et al._ (2011); Hoang  _et al._ (2015a), we obtain results
entirely compatible with these works. But here we will also show that multiple
well-defined scheme choices may well lead to a significant spread of the fit
values, implying a potential systematic theory uncertainty that should be
accounted for, while also motivating further strategies to reduce or avoid
such systematic effects.

As part of our investigations we will also find that $\alpha_{s}$ extractions
that are constrained to a more restricted domain of the observable—in
particular the dominantly dijet region of differential distributions—typically
result in better-quality fits than those performed up to or even beyond the
multi-jet threshold. This represents a complementary perspective to the one
that was recently presented in Nason and Zanderighi (2023), which advocated
$\alpha_{s}$ fits in the far-tail regions, where three-jet and multi-jet
events dominate, and where recent advances in understanding three-jet NP power
corrections have focused. While the latter developments are important to
improve the theoretical understanding of event-shape distributions, we note
that these anyalses are currently model-dependent, and therefore not on the
same footing as the corresponding dijet studies. In particular, the
universality of the three-jet NP power corrections is currently not well
understood.

In this work we study the thrust distribution $T=1-\tau$ Brandt  _et al._
(1964); Farhi (1977), a canonical $e^{+}e^{-}$ event shape defined as

$\tau=\frac{1}{Q}\,\underset{i}{\sum}\;|{\bf
p}_{\perp}^{i}|\;e^{-|\eta_{i}|}\,,$ (2)

with $Q$ the center-of-mass (c.o.m.) energy of the collision, and ${\bf
p}_{\perp}^{i}$ and $\eta_{i}$ respectively denoting the transverse momentum
and rapidity of the $i$th final-state particle measured with respect to the
thrust axis. Thrust is amongst the best studied $e^{+}e^{-}$ event shapes, a
fact in part due to the ability to calculate its spectra at
N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{3})$ resummed and matched accuracy
Becher and Schwartz (2008); Abbate  _et al._ (2011). As already mentioned,
thrust has been utilized for prior effective-field-theory-based extractions of
the strong coupling, and this provides a controlled environment for us to
study the impacts of certain assumptions embedded in the theoretical
framework, and to compare with prior analyses. Indeed, in addition to
revealing the novel physics points mentioned above, our study provides the
first independent cross-check of the results in Abbate  _et al._ (2011); Hoang
_et al._ (2015a), despite a number of minor systematic differences in our
analysis that will be described in detail in Sec. V.4 below.

The remainder of the paper develops as follows: in Sec. II we will review the
dijet SCET factorization theorem employed to predict the thrust distribution,
the various perturbative ingredients required therein, and the matching to the
fixed-order QCD prediction up to $\mathcal{O}(\alpha_{s}^{2})$ accuracy. We
also describe a new high-statistics calculation of the three-loop
$\mathcal{O}(\alpha_{s}^{3})$ remainder function using the public EERAD3 code
Gehrmann-De Ridder  _et al._ (2014), and highlight some systematic concerns
about it that lead us to not include the matching to this order in our final
extraction code. In Sec. III we review the further factorization of the soft
function into a non-perturbative shape function and a fully perturbative
component, as well as the formalism for achieving the leading infrared
renormalon cancellation between the two. We will then introduce the three
different renormalon cancellation schemes we intend to study, and in Sec. IV
we review two alternative sets of profile functions as well as our method of
varying them to estimate the overall perturbative uncertainty. Finally, in
Sec. V we present our method of extracting $\\{\alpha_{s},\Omega_{1}\\}$, show
our results for these quantities in different schemes, discuss their
implications, and also touch on the uncertainties associated with the three-
loop soft constant $c_{\tilde{S}}^{3}$, which is the only N3LL′ ingredient
that is currently not known exactly. We conclude in Sec. VI. App. A and B
respectively collect details on the fixed-order expansion of the resummed
cross section and the order-by-order expansion of the renormalon cancellation
formulae, while App. C provides more details of an analysis of one of the
renormalon schemes introduced in Sec. III that does not play a central role in
the comparison of schemes in the main body of this paper nor in its
conclusions.

##  II Perturbative Treatment

The bulk of the theory we implement is thoroughly presented in prior
literature, and so for brevity we only outline the core features of our
calculation in what follows, leaving many background details to those
references (see e.g. Abbate  _et al._ (2011); Almeida  _et al._ (2014); Bell
_et al._ (2019)).

We predict the singular part of the integrated (or cumulative) thrust
distribution

$\sigma_{c}(\tau)=\int_{0}^{\tau}d\tau^{\prime}\frac{d\sigma}{d\tau^{\prime}}$
(3)

with methods from SCET, for which factorization theorems for dijet
$e^{+}e^{-}$ event shapes are well established, and match to fixed-order QCD
to account for non-singular contributions. Hence the overall perturbative
cross section consists of two parts,

$\sigma_{c}^{PT}(\tau)=\sigma_{c,\text{sing}}(\tau)+\sigma_{c,\text{ns}}(\tau)\,,$
(4)

with the non-singular matching contribution $\sigma_{c,\text{ns}}(\tau)$
implicitly defined in (11) below. We will now address both of these
contributions in turn.

###  II.1 Singular Contribution

The singular and resummed contribution $\sigma_{c,\text{sing}}(\tau)$,
normalized by the Born cross section $\sigma_{0}$, reads111We use the form
derived in Bell  _et al._ (2019), which depends at a given order of
perturbation theory on the factorization scales $\mu_{H,J,S}$, but is
otherwise explicitly independent of the renormalization scale $\mu$ at _every_
order of perturbation theory.

$\displaystyle\frac{\sigma_{c,\text{sing}}(\tau)}{\sigma_{0}}=e^{\widetilde{K}(\mu_{H},\mu_{J},\mu_{S};Q)+K_{\gamma}(\mu_{H},\mu_{J},\mu_{S})}\biggl{(}\frac{1}{\tau}\biggr{)}^{\Omega(\mu_{J},\mu_{S})}$
$\displaystyle\times
H(Q^{2},\mu_{H})\widetilde{J}\Bigl{(}\partial_{\Omega}+\ln\frac{\mu_{J}^{2}}{Q^{2}\tau},\mu_{J}\Bigr{)}^{2}\;\widetilde{S}\Bigl{(}\partial_{\Omega}+\ln\frac{\mu_{S}}{Q\tau},\mu_{S}\Bigr{)}$
$\displaystyle\times\frac{e^{\gamma_{E}\Omega}}{\Gamma(1-\Omega)}\,,$ (5)

where the evolution kernels $\tilde{K}$, $K_{\gamma}$, and $\Omega$ are given
by

$\displaystyle\widetilde{K}(\mu_{H},\mu_{J},\mu_{S};Q)$
$\displaystyle\equiv-\kappa_{H}\widetilde{K}_{\Gamma}(\mu,\mu_{H};Q)$ (6)
$\displaystyle\\!\\!\\!\\!-2j_{J}\kappa_{J}\widetilde{K}_{\Gamma}(\mu,\mu_{J};Q)-\kappa_{S}\widetilde{K}_{\Gamma}(\mu,\mu_{S};Q)\,,$
$\displaystyle K_{\gamma}(\mu_{H},\mu_{J},\mu_{S})$ $\displaystyle\equiv
K_{\gamma_{H}}(\mu,\mu_{H})+2K_{\gamma_{J}}(\mu,\mu_{J})$
$\displaystyle\quad+K_{\gamma_{S}}(\mu,\mu_{S})\,,$
$\displaystyle\Omega(\mu_{J},\mu_{S})$
$\displaystyle\equiv-2\kappa_{J}\,\eta_{\Gamma}(\mu,\mu_{J})-\kappa_{S}\,\eta_{\Gamma}(\mu,\mu_{S})\,.$

The hard function $H$ arises from the matching of the effective theory to QCD
and it describes virtual corrections to $e^{+}e^{-}\rightarrow q\bar{q}$
scattering, while collinear radiation along the jet directions is encoded in
the jet functions $\widetilde{J}$ (with the tilde denoting an evaluation in
Laplace space). Finally, background low-energetic radiation that communicates
between the two jets is described by the (Laplace-space) soft function
$\widetilde{S}$. All of these functions are evaluated in (5) at an associated
‘natural’ scale $\mu_{H,J,S}$ at which their logarithmic corrections are
minimized and their perturbative series well-behaved. In SCET the resummation
is achieved via renormalization group (RG) evolution, which generates the
kernels $\widetilde{K}$, $K_{\gamma}$, $\Omega$ in (6) defined in terms of
cusp ($\Gamma_{\rm cusp}$) and non-cusp ($\gamma_{F}$) anomalous dimensions,

$\displaystyle\widetilde{K}_{\Gamma}(\mu,\mu_{F};Q)$
$\displaystyle\equiv\int_{\mu_{F}}^{\mu}\frac{d\mu^{\prime}}{\mu^{\prime}}\Gamma_{\rm
cusp}[\alpha_{s}(\mu^{\prime})]\ln\frac{\mu^{\prime}}{Q}\,,$ $\displaystyle
K_{\gamma_{F}}(\mu,\mu_{F})$
$\displaystyle\equiv\int_{\mu_{F}}^{\mu}\frac{d\mu^{\prime}}{\mu^{\prime}}\gamma_{F}[\alpha_{s}(\mu^{\prime})],$
$\displaystyle\eta_{\Gamma}(\mu,\mu_{F})$
$\displaystyle\equiv\int_{\mu_{F}}^{\mu}\frac{d\mu^{\prime}}{\mu^{\prime}}\Gamma_{\rm
cusp}[\alpha_{s}(\mu^{\prime})].$ (7)

The corresponding values for the parameters $\kappa_{F}$ and $j_{F}$ are given
for the thrust distribution by

$\displaystyle j_{H}$ $\displaystyle=1,\,\,\,\,\,$ $\displaystyle
j_{J}=2,\,\,\,\,\,$ $\displaystyle j_{S}=1,$ $\displaystyle\kappa_{H}$
$\displaystyle=4,\,\,\,\,\,\,$ $\displaystyle\kappa_{J}=-2,\,\,\,\,\,$
$\displaystyle\kappa_{S}=4.$ (8)

Given the perturbative expansions of the cusp and non-cusp anomalous
dimensions, one can solve (7) order by order, achieving approximate analytic
expressions for the evolution kernels that resum the logarithmic corrections
to the thrust distribution to a given accuracy.222The systematic uncertainties
of this approximation have been studied e.g. in Billis  _et al._ (2020); Ebert
(2022). The order to which each individual component of (5) must be calculated
can be read off e.g. from Table 6 of Bell  _et al._ (2019). As our goal
consists in resumming the thrust distribution to N3LL′ accuracy, we need the
non-cusp anomalous dimensions $\gamma_{F}$ and the fixed-order hard, jet, and
soft functions $H,\widetilde{J},\widetilde{S}$ up to
$\mathcal{O}(\alpha_{s}^{3})$ corrections, while $\mathcal{O}(\alpha_{s}^{4})$
ingredients are required for the cusp anomalous dimension $\Gamma_{\rm cusp}$
and the QCD $\beta$-function.

Of these, only the three-loop matching correction $c_{\tilde{S}}^{3}$ to the
Laplace-space soft function $\tilde{S}$ is currently unknown, although an
extraction of this quantity with the EERAD3 event-generator results from Monni
_et al._ (2011) was performed in Brüser  _et al._ (2018), which found

$c_{\tilde{S}}^{3}\big{|}_{\tt EERAD3}\leavevmode\nobreak\ =-19988\pm 5440\,,$
(9)

where we have added the uncertainties quoted in Brüser  _et al._ (2018)
linearly. However, in Sec. V.5 below we will also examine scenarios with a
Padé approximated value for $c_{\tilde{S}}^{3}$ that has been used in earlier
studies of the thrust distribution (in particular in Abbate  _et al._ (2011))
with

$c_{\tilde{S}}^{3}\big{|}_{\rm Pad\acute{e}}\leavevmode\nobreak\
=\leavevmode\nobreak\ 691\pm 1000\,.$ (10)

The goal of that exercise consists in understanding how sensitive $\alpha_{s}$
determinations are to this unknown N3LL′ coefficient in different perturbative
schemes that we will introduce below. In any case, despite the somewhat large
systematic difference between these two estimates of $c_{\tilde{S}}^{3}$, and
the remaining large uncertainty in (9), we will still label results that
depend on these estimates as N3LL′ in what follows, in accord with previous
analyses Abbate  _et al._ (2011); Hoang  _et al._ (2015a).

###  II.2 Non-singular Contribution

In the far-tail region of the distribution, where resummation effects are less
important, it is necessary to match the singular SCET cross section from (5)
to the full QCD prediction. Obtaining the additional non-singular component in
(4) is then a matter of determining remainder coefficients $r_{c}^{i}$ by
subtracting off the singular SCET prediction from the QCD result, at a given
order in perturbation theory,

$\displaystyle\frac{\sigma^{PT}_{c}(\tau)}{\sigma_{0}}-\frac{\sigma_{\text{c,sing}}(\tau)}{\sigma_{0}}=r_{c}(\tau)$
$\displaystyle\simeq\theta(\tau)\,\biggl{\\{}\bar{r}_{c}^{1}(Q,\tau)+\bar{r}_{c}^{2}(Q,\tau)+\bar{r}_{c}^{3}(Q,\tau)\biggr{\\}}\,,$
(11)

where
$\bar{r}_{c}^{i}(Q,\tau)\equiv\big{(}\alpha_{s}(Q)/(2\pi)\big{)}^{i}\,r_{c}^{i}(\tau)$.
The singular cross section is obtained from (5) by expanding
$F=H,\widetilde{J},\widetilde{S}$ at $\mu_{H}=\mu_{J}=\mu_{S}=Q$ to fixed
order in $\alpha_{s}$, multiplying out these expressions, and then inverse
Laplace-transforming the result. This process is straightforward given the
well-known expressions for the fixed-order coefficients $F_{n}$ defined by

$F(L_{F},\mu_{F})=\sum_{n=0}^{\infty}\left(\frac{\alpha_{s}(\mu_{F})}{4\pi}\right)^{n}F_{n}(L_{F})\,,$
(12)

which depend on (logarithmically weighted) anomalous dimensions
$\gamma_{F}^{n}$, $\Gamma_{F}^{n}$, as well as unweighted singular constants
$c_{F}^{n}$, briefly discussed above for $F=\widetilde{S}$. In App. A we have
collected explicit expressions for these functions up to
$\mathcal{O}(\alpha_{s}^{3})$, in both momentum and Laplace space. In terms of
these quantities, the singular cross section $\sigma_{\text{c,sing}}(\tau)$ is
fully determined, and one only needs to obtain the QCD result
$\sigma^{PT}_{c}(\tau)$ to determine the remainder coefficients. At
$\mathcal{O}(\alpha_{s})$ this can be done analytically (from, e.g., Ellis
_et al._ (2011)), while we have implemented the relevant
$\mathcal{O}(\alpha_{s}^{2})$ matching from Bell  _et al._ (2019) in this
work, which used EVENT2 Catani and Seymour (1997) to obtain $r_{c}^{2}$. This
then leaves the $\mathcal{O}(\alpha_{s}^{3})$ remainder coefficient
$r_{c}^{3}$.

Figure 1: Differential remainder function $r^{3}(\tau)/(2\pi)^{3}$ as
extracted from EERAD3. The data points are the EERAD3 data, and the red line
is a combined interpolation for large $\tau$ and a fit to a basis of
logarithmic functions for small $\tau$ as explained in the text.

While the differential analogue to $r_{c}^{3}$, $r^{3}=dr_{c}^{3}/d\tau$, is
available from the corresponding SCET matching performed in Abbate  _et al._
(2011), we have opted to perform our own extraction using the EERAD3
generator, which yields a binned approximation to the differential thrust
distribution away from $\tau=0$,

$\frac{1}{\sigma_{0}}\frac{d\sigma_{c}}{d\tau}\Big{|}_{\tt
EERAD3}\simeq\frac{1}{\sigma_{0}}\frac{d\sigma_{c}}{d\tau}\Big{|}_{\tau>0}=B(\tau)+r(\tau)\,,$
(13)

where $B(\tau)$ is a singular function that is reproduced by the SCET
calculation. Critically, the unknown constant $c_{\tilde{S}}^{3}$ mentioned
above does not appear in this expression, as it only contributes to the
coefficient of the $\delta(\tau)$ term. This can be seen explicitly in (58) in
App. A, which can also be used to obtain $r^{3}$ using the differential
analogue of (11). The resulting object can then be integrated to obtain the
cumulative remainder function $r_{c}^{3}$.

Following this procedure, we have computed the remainder function $r^{3}$ in
EERAD3 with $1.5\times 10^{10}$ events with the internal infrared cutoff
parameter set to $10^{-7}$. The result is shown in Fig. 1. We note that this
number of events is considerably greater than in prior determinations
appearing in the literature. With these statistics, we are able to learn about
the behavior of the fixed-order prediction of EERAD3 in a bit more detail than
previously possible, especially in the small $\tau$ region. The remainder
$r^{3}$ is obtained by subtracting all the singular contributions predicted by
(64) down to the single-logarithmic coefficient, i.e. the logarithms with
coefficients $\sigma_{3k}$ with $k\geq 1$. However, we have noticed that the
coefficient $\sigma_{31}$ of the single logarithm at three loops is itself not
predicted accurately by EERAD3—it is off by nearly a factor of two. Thus, the
result in Fig. 1 likely contains an uncancelled singular contribution. We have
gone ahead with the exercise of fitting the data in Fig. 1 to a basis of sub-
leading (power-suppressed) logarithmic functions as was done in Abbate  _et
al._ (2011), which is represented by the red line for $\tau<0.2$ (the red line
for $\tau>0.2$ is simply a direct interpolation of the EERAD3 data). In this
way we could integrate $r^{3}$ to obtain a cumulative remainder function
$r_{c}^{3}$. However, because of the uncancelled singular contribution in the
EERAD3 data, we do not believe we can obtain a reliable prediction for
$r_{c}^{3}$ at this time, nor can we use these data to obtain a reliable
extraction of the unknown three-loop soft constant $c_{\tilde{S}}^{3}$ using
the procedure described in Hoang and Kluth (2008). These observations should
not, however, cast doubt on previous extractions of $\alpha_{s}$ using EERAD3
calculations of $r^{3}$, as the ‘missing logarithmic effect’ is less
pronounced at larger values of $\tau$ where the $\alpha_{s}$ fits are
performed, and the statistical uncertainties for the lower number of events
that were used in those works likely encompass the error, anyway.

Nevertheless, as a result of these observations, we choose to match only to
the $\mathcal{O}(\alpha_{s}^{2})$ remainder function $r_{c}^{2}$ in this work,
effectively setting $r_{c}^{3}$ to zero. Thus we label our results as
N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$. We stress, however, that our
study of renormalon cancellation scheme and perturbative scale profile
choices, and the conclusions drawn therefrom, remain unaffected from this
restriction. We have verified this by performing extractions with $r^{3}$ in
Fig. 1 turned on or off. We defer further study of $r^{3}$ and its integral
$r_{c}^{3}$ to a future publication. We also note that the more recent
$\alpha_{s}$ determinations in Kardos  _et al._ (2018); Verbytskyi  _et al._
(2019) used CoLoRFulNNLO Del Duca  _et al._ (2016) to determine the
$\mathcal{O}(\alpha_{s}^{3})$ remainder function.

With this set of remainder functions $r^{i}_{c}$ at hand, the final
perturbative prediction for the cumulative thrust distribution is given by

$\displaystyle\sigma^{PT}_{c}(\tau)$
$\displaystyle=\frac{\sigma_{c,\text{sing}}(\tau;\mu_{H},\mu_{J},\mu_{S})}{\sigma_{0}}+\frac{\alpha_{s}(\mu_{ns})}{2\pi}r_{c}^{1}(\tau)$
$\displaystyle+\left(\frac{\alpha_{s}(\mu_{ns})}{2\pi}\right)^{2}\left[r_{c}^{2}(\tau)+\beta_{0}r_{c}^{1}(\tau)\ln\frac{\mu_{ns}}{Q}\right]\,,$
(14)

where we have made the dependence of the singular cross section on the hard,
jet and soft scales $\mu_{H,J,S}$ explicit, and we have allowed for an
independent scale $\mu_{ns}$ that can be varied to probe the perturbative
uncertainty of the non-singular contribution. In Sec. IV we will discuss how
we vary these scales to estimate unknown higher-order corrections in all
sectors of our calculation.

##  III Non-Perturbative Treatment

###  III.1 Gapped Shape Function

A complete theoretical treatment of the thrust distribution must also take
into account non-perturbative effects due to hadronization. These effects are
encoded in the dijet soft function in (5) describing low-energy wide-angle
radiation between the two jets. We use a model for the soft function
Korchemsky (1998); Korchemsky and Sterman (1999); Hoang and Stewart (2008),

$S(k,\mu_{S})=\int
dk^{\prime}\,S_{PT}(k-k^{\prime},\mu_{S})\,f_{\text{mod}}(k^{\prime}-2\overline{\Delta})\,,$
(15)

where $f_{\text{mod}}$ is a non-perturbative shape function that modifies the
perturbative prediction for the soft function, and $\overline{\Delta}$ is a
gap parameter modeling the minimum soft momentum in a final state due to
hadronization. This implies the following formula for the cross section
itself:

$\displaystyle\frac{1}{\sigma_{0}}\,\sigma_{c}(\tau)$ $\displaystyle=\int
dk\,\sigma_{c}^{PT}\Bigl{(}\tau-\frac{k}{Q};\mu_{H},\mu_{J},\mu_{S},\mu_{ns}\Bigr{)}$
$\displaystyle\times f_{\text{mod}}\bigl{(}k-2\overline{\Delta}\bigr{)}\,.$
(16)

In this prediction, the shape function is convolved with both the (resummed)
singular and the non-singular parts of the cross section in (II.2)—as in
previous treatments Abbate  _et al._ (2011); Hoang  _et al._ (2015a); Bell
_et al._ (2019)—in order to smooth the transition from the resummation region
to the fixed-order region. This is an important point to keep in mind in the
following discussion, as any modifications we make to the non-perturbative gap
parameter will affect both regions.

For the gapped shape function in (15), we adopt the form

$f_{\text{mod}}(k)=\frac{1}{\lambda}\left[\sum_{n=0}^{\infty}\,b_{n}\,f_{n}\left(\frac{k}{\lambda}\right)\right]^{2}\,,$
(17)

where the functions $f_{n}(x)$ form a complete orthonormal basis composed of
Legendre polynomials Ligeti  _et al._ (2008). For the present analysis, we set
$b_{0}=1$ and $b_{n}=0$ for $n>0$, since we are only interested in the tail
region of the distribution, where the first moment $\lambda$ of the shape
function is the only relevant parameter. In the absence of the non-
perturbative gap parameter $\overline{\Delta}$, one then recovers the leading
shift from (1) when the shape function is convolved with the perturbative
cross section. With a non-zero $\overline{\Delta}$, on the other hand, the
relationship between $\overline{\Omega}_{1}$ and the shape function is
modified to

$2\overline{\Omega}_{1}=2\overline{\Delta}+\int dk\,k\,f_{\text{mod}}(k)\,,$
(18)

with $\overline{\Delta}\sim\Lambda_{\rm{QCD}}$, and where it is understood
that the $\lambda$ parameter in (17) is given by
$\lambda=2\left(\overline{\Omega}_{1}-\overline{\Delta}\right)$. Note,
however, that our notation in (18) has slightly changed from (1), as the
presence of the barred notation ($\overline{\Omega}_{1},\overline{\Delta}$)
indicates objects defined in an $\overline{\text{MS}}$-like renormalization
scheme, where the soft function $S_{PT}$ has been calculated. However, both
the perturbative soft function and the gap parameter $\overline{\Delta}$
exhibit renormalon ambiguities Hoang and Stewart (2008), which must be
cancelled to obtain stable predictions.

To do so, one redefines the gap parameter in (15) as

$\overline{\Delta}=\Delta(\mu_{\delta},\mu_{R})+\delta(\mu_{\delta},\mu_{R}),$
(19)

where $\Delta$ is renormalon free and $\delta$ is chosen to cancel the
ambiguity in the perturbative soft function. We will review the procedure for
calculating $\delta$ in various schemes in Sec. III.2 below. In these schemes,
the perturbative series for $\delta$ is calculated starting from the
perturbative soft function that is renormalized at a ‘reference’ scale
$\mu_{\delta}$. Moreover, the condition imposed on the Laplace-space soft
function is evaluated at an argument proportional to the inverse of a second
‘subtraction’ scale $\mu_{R}$. The dependence on these two scales, as well as
the details on the chosen subtraction scheme to compute $\delta$, enters the
final prediction for the cross section.

Without repeating the details of the derivation, we quote here the final
result for our prediction of the renormalon-free cross section Abbate  _et
al._ (2011); Hoang  _et al._ (2015b); Bell  _et al._ (2019):

$\displaystyle\frac{1}{\sigma_{0}}\,\sigma_{c}(\tau)$ $\displaystyle=\int
dk\,\sigma_{c}^{PT}\Bigl{(}\tau-\frac{k}{Q};\mu_{H},\mu_{J},\mu_{S},\mu_{ns}\Bigr{)}$
(20)
$\displaystyle\times\Bigl{[}e^{-2\delta(\mu_{\delta},\mu_{R})\frac{d}{dk}}f_{\text{mod}}\bigl{(}k-2\Delta(\mu_{\delta},\mu_{R})\bigr{)}\Bigr{]}\,.$

The scheme chosen to compute the renormalon subtraction series, and thus the
definition of the gap parameter $\Delta$, must be specified in computing the
cross section. The scheme includes a definition of the reference scale
$\mu_{\delta}$, and in case it is not equal to $\mu_{S}$, the terms in
brackets should be re-expanded in powers of $\alpha_{s}(\mu_{S})$ Bachu  _et
al._ (2021).

Independent of the chosen renormalon cancellation scheme is the practical
approach for evaluating (20). We do so by collecting terms with the same
explicit powers of $\alpha_{s}$, such that the cross section necessary for
N3LL′ accuracy can be expanded as

$\frac{\sigma_{c}(\tau)}{\sigma_{0}}=\sigma_{c}^{[0]}(\tau)+\sigma_{c}^{[1]}(\tau)+\sigma_{c}^{[2]}(\tau)+\sigma_{c}^{[3]}(\tau)\,,$
(21)

whose components are given order-by-order by

$\sigma_{c}^{[i]}(\tau)=\int
dk\,\left[\sum_{n=0}^{i}\,\sigma_{c}^{\text{N${}^{n}$LO}}\Bigl{(}\tau-\frac{k}{Q}\Bigr{)}\,f_{\text{mod}}^{(i-n)}(k-2\Delta)\right]\,,$
(22)

where we have suppressed the scale dependence in both
$\sigma_{c}^{\text{N${}^{n}$LO}}$ (the purely perturbative resummed and
matched cross section from (II.2), with any fixed-order prefactors truncated
to NnLO accuracy, following the notation used in Eq. (4.38) of Bell  _et al._
(2019)) and the shape function $f_{\text{mod}}$ for brevity. The component
notation $f_{\text{mod}}^{(i)}$ represents the fact that in practice we expand
out the renormalon-corrected shape function as

$\displaystyle e^{-2\delta\frac{d}{dk}}f_{\text{mod}}(k-2\Delta)$
$\displaystyle=\sum_{i}f_{\text{mod}}^{(i)}(k-2\Delta)\,,$ (23)

where the coefficients $f_{\text{mod}}^{(i)}$ refer to an expansion in
$\alpha_{s}(\mu_{S})$, given in (67). Explicit expressions for these
coefficients can be found order-by-order in App. B, where it is clear that
these terms _depend_ on the particular renormalon cancellation scheme. It also
follows that in the presence of a gapped and renormalon-corrected shape
function, the actual shift of the differential distribution will no longer be
a constant, as in (1). Indeed, the scale dependence of both $\Delta$ and
$\delta$, which will be discussed in upcoming sections, leads to a
$\tau$-dependent ‘effective’ shift, which can be calculated as

$\zeta_{\text{eff}}(\tau)\equiv\int
dk\,k\left[\sum_{i}f_{\text{mod}}^{(i)}(k-2\Delta)\right]\,.$ (24)

The behavior of this effective shift will guide in large part our
considerations below on a set of renormalon cancellation schemes for the
thrust distribution.

###  III.2 Renormalon Cancellation Schemes

As mentioned above, both the perturbative soft function $S_{PT}$ and the
subtraction term $\delta$ suffer from renormalon ambiguities associated to
infrared poles in their (all-order) Borel-series representations. As in prior
studies, we adopt a formalism to cancel them against one another, thereby
rendering the overall cross section free of the leading soft renormalon. A
generalized set of schemes achieving this cancellation was presented in Bachu
_et al._ (2021), each defined by imposing a condition on the soft function in
Laplace space to render it free of the leading renormalon:

$\frac{d^{n}}{d(\ln\nu)^{n}}\ln\left[\widetilde{S}_{PT}(\nu,\mu_{\delta})\,e^{-2\nu\delta(\mu_{\delta},\mu_{R})}\right]_{\nu=\xi/\mu_{R}}=0\,,$
(25)

where $\nu$ is the Laplace-space variable, $\mu_{\delta}$ the reference
renormalization scale where $\widetilde{S}_{PT}$ is evaluated, and $\mu_{R}$
the renormalon subtraction scale, with the condition on the $n$-th derivative
of $\widetilde{S}_{PT}$ imposed at the Laplace-space argument
$\nu=\xi/\mu_{R}$. The parameters $n$, $\xi$, along with the choice of
$\mu_{\delta}$, define the renormalon cancellation scheme being used to define
the subtraction term $\delta$ and thus the gap parameter $\Delta$. From (25)
one immediately obtains an expression for the subtraction term $\delta$ in
terms of the perturbative soft function,

$\delta(\mu_{\delta},\mu_{R})=\frac{\mu_{R}}{2\xi}\frac{d^{n}}{d(\ln\nu)^{n}}\ln\widetilde{S}_{PT}(\nu,\mu_{\delta})\big{|}_{\nu=\xi/\mu_{R}}\,.$
(26)

Here it is clear that we need to control the RG evolution of the subtraction
term, both in terms of the reference scale $\mu_{\delta}$ and the subtraction
scale $\mu_{R}$, which is governed by the anomalous dimensions

$\displaystyle\gamma_{\Delta}\left[\alpha_{s}(\mu_{\delta})\right]$
$\displaystyle=\frac{d}{d\ln\mu_{\delta}}\Delta(\mu_{\delta},\mu_{R})=-\frac{d}{d\ln\mu_{\delta}}\delta(\mu_{\delta},\mu_{R})\,,$
$\displaystyle\gamma_{R}\left[\alpha_{s}(\mu_{R})\right]$
$\displaystyle=-\frac{d}{d\mu_{R}}\Delta(\mu_{R},\mu_{R})=\frac{d}{d\mu_{R}}\delta(\mu_{R},\mu_{R})\,.$
(27)

The $\mu_{\delta}$-evolution is well-known Hoang and Kluth (2008), and is
given in terms of the cusp evolution function $\eta_{\Gamma}$ defined in (7).
The second equation is the so-called ‘$R$-evolution’ equation Hoang and Kluth
(2008); Hoang  _et al._ (2008), and its anomalous dimension $\gamma_{R}$
explicitly depends on the subtraction scheme defined in (26). Note that
$\gamma_{R}$ is computed at $\mu_{\delta}=\mu_{R}$. We give the solution for
the $\mu_{\delta}$\- and $R$-evolved gap parameter to three-loop order in
(69).

We now address the freedom in defining renormalon cancellation schemes which,
as is evident in (26), depend on the choice of the derivative rank (the
parameter $n$), the overall normalization (the parameter $\xi$), and the
reference scale $\mu_{\delta}$. As long as an appropriate hierarchy of scales
is maintained, one may choose these parameters freely. In Abbate  _et al._
(2011); Hoang  _et al._ (2015a), for example, the so-called _R(-Gap) Scheme_
was defined by

R Scheme:
$\displaystyle\\{n,\xi,\mu_{\delta},\mu_{R}\\}=\\{1,e^{-\gamma_{E}},\mu_{S},R\\}\,,$
(28)

where the functional form of the profile $R$ will be given in Sec. IV.1 below.
Exact expressions for $n=1$ for both the subtraction term $\delta$ and the
anomalous dimension $\gamma_{R}$ are given in App. B up to
$\mathcal{O}(\alpha_{s}^{3})$.

In Fig. 2 we plot the effective shift (24) of the differential distribution
that results from various scheme choices across the $\tau$-domain relevant for
the $\alpha_{s}$ fits. Specifically, the red and green curves refer to the R
scheme for two different perturbative scale choices (‘2010’ and ‘2018’) that
we will introduce in Sec. IV.1 below. A notable feature of Fig. 2 is that the
effective shift of the cross section _grows_ as $\tau$ increases, even into
the region of multi-jet configurations, where the applied formalism for the
description of NP corrections, which is based on a dijet factorization
theorem, loses its validity. We stress that the growth of $\zeta_{\text{eff}}$
in (24) is a direct consequence of the $R$ evolution of the gap parameter
$\Delta$, predicted by (III.2). Its lasting effect in the large $\tau$ region
of the cross section, in particular, is a result of the choice made in (16) to
convolve the shape function with both the singular and non-singular parts of
the cross section. As argued above, this choice was made in order to smooth
the transition between the dijet and multi-jet regions. Interestingly, we note
that similar effects of a growing non-perturbative shift are also mimicked by
models for power corrections to the three-jet region, recently studied in
Caola  _et al._ (2022b), although the quantitative effect found in that work
is more pronounced. In the context of renormalon cancellation-induced effects
in a _dijet_ factorization formula, however, the growth of the effective shift
$\zeta_{\text{eff}}$ well into the multi-jet region can appear surprising.

These observations led us to explore the space of renormalon cancellation
schemes that may yield an effective shift whose growth is mitigated for large
$\tau$ values, while still achieving the cancellation of the leading soft
renormalon. As the growth of (24) is related to the RG evolution of the gap
parameter $\Delta(\mu_{\delta},\mu_{R})$, and therefore the scale profiles of
Sec. IV.1, this motivated us to look for other choices of the reference scale
$\mu_{\delta}$ and/or the subtraction scale $\mu_{R}$ to achieve this goal.

To that end, we also define the following _R ⋆ Scheme_:

R⋆ Scheme:
$\displaystyle\\{n,\xi,\mu_{\delta},\mu_{R}\\}=\\{1,e^{-\gamma_{E}},R^{\star},R^{\star}\\}\,,$
(29)

where we have chosen to identify the reference scale $\mu_{\delta}$ not with
the soft scale $\mu_{S}$, but rather with the subtraction scale $\mu_{R}$
itself. This eliminates right away any logarithms of $\mu_{\delta}/\mu_{R}$ in
the renormalon subtraction term $\delta$, as evident from (76) and (81),
though there will still be logarithms of $\mu_{S}/\mu_{R}$ when re-expanding
the subtractions in $\alpha_{s}(\mu_{S})$. In this particular scheme, however,
these logarithms do not show up until $\mathcal{O}(\alpha_{s}^{3})$, since the
one-loop subtraction vanishes exactly, see (81a).

The scheme choice (29) gives us some freedom to tweak the behavior of the
subtraction scale $R^{*}$ without inducing too large logarithms with respect
to the soft scale $\mu_{S}$. For instance, we may _freeze_ the growth of
$R^{*}$ as a function of $\tau$ with an appropriate choice of the profile
functions in Sec. IV.1. In particular, we will choose a functional form for
this new scale, $\mu_{\delta}=\mu_{R}=R^{\star}$, to be a piecewise function:

$R^{\star}(\tau)=\begin{cases}R(\tau)&R(\tau)<R_{\rm{max}}\\\
R_{\rm{max}}&R(\tau)\geq R_{\rm{max}}\,.\end{cases}$ (30)

where $R(\tau)$ is the profile function one might have chosen for the R
scheme, and the value $R_{\rm{max}}$ is the value of $R(\tau)$ at the point
$\tau=t_{1}$, which marks the transition between the non-perturbative and
resummation regions, as will be discussed in more detail in Sec. IV.1.
Importantly, due to the choice $\mu_{\delta}=R^{\star}$ all logarithms in
$\mu_{\delta}/R$ appearing in the subtraction term $\delta$ are turned off in
the (dominantly non-perturbative) domain with $\tau<t_{1}$, while for
$\tau>t_{1}$ we allow for small logarithmic contributions in $\mu_{\delta}/R$.
In contrast to this, imposing the modified profile in (30) directly in the R
scheme would have induced large logarithmic corrections in $\mu_{\delta}/R$,
since the reference scale $\mu_{\delta}$ traces the soft scale $\mu_{S}$ in
this case.

The resulting effective shift in the R⋆ scheme is plotted in black and blue in
Fig. 2, respectively, for the two perturbative scale choices that we mentioned
before. In particular, one observes a leveling off in comparison to the
standard R scheme, for both scale choices, as intended. This figure also
illustrates that the range of $n=1$ scheme choices we consider corresponds to
a variation of power corrections in the multi-jet region of order 10%.

Figure 2: The effective shift (24) of the cross section due to different
renormalon cancellation and perturbative scale profile schemes, evaluated
using ‘central’ profiles at $Q=m_{Z}$. The flat, dashed line corresponds to
the constant shift of (1), while the vertical gray lines correspond to the two
fitting windows discussed below. All functions are calculated at
N3LL′$+\mathcal{O}(\alpha_{s}^{2})$ accuracy.

As a third and entirely independent scheme, we also consider an instance with
$n=0$, which we refer to as the _R 0 Scheme_:

R0 Scheme:
$\displaystyle\\{n,\xi,\mu_{\delta},\mu_{R}\\}=\\{0,2\pi,\mu_{S},R\\}\,.$ (31)

We also give explicit expressions for $\delta$ and $\gamma_{R}$ in the $n=0$
case in App. B, where one notices that the cancellation ingredients
generically appear with one higher logarithmic power than in the $n=1$ case.
The overall effective shift of the R0 Scheme, with varying scale choices, is
also shown in Fig. 2 in brown and cyan, where one notices that while both
curves increase with $\tau$ with a much steeper slope, their overall magnitude
is noticeably reduced with respect to the $n=1$ schemes for all $\tau$ (notice
that the curves are multiplied by 1.25 in the figure). These effects are
jointly associated to the additional power of logarithms appearing in the
$n=0$ cancellation terms with respect to $n=1$ counterparts (cf. App. C), and
to the $\mathcal{O}(10)$ difference in the normalization factor $\xi$ between
R(⋆) ($e^{-\gamma_{E}}$) and R0 ($2\pi$) scheme definitions we have
chosen.333Note that the gap parameter $\Delta$ appearing in
$\zeta_{\text{eff}}$ in (24) evaluated at the reference scale $R_{\Delta}$
remains the same between the $R^{0}$ and $R^{(\star)}$ schemes:
$\Delta(R_{\Delta},R_{\Delta})=0.1$ GeV. However, due to less stable
perturbative convergence and harder-to-control uncertainties in this scheme,
we will not consider it as a candidate for our main analysis below.
Regardless, for completeness, we do show some results in this scheme in Sec.
V.5 and App. C.

The important point to emphasize is that there are a number of consistent
schemes that can be used to cancel the leading soft renormalon of the
perturbative soft function $S_{PT}$ and the gap parameter $\Delta$, and in
upcoming sections we will study the impact of this choice on the numerical
extraction of $\\{\alpha_{s},\Omega_{1}\\}$. There are, of course, other
considerations one can use to determine if one scheme is preferable or better
behaved than another, e.g. perturbative stability and convergence, and, if so,
one could thus argue that the theoretical uncertainty coming from varying
schemes can be removed (as we have so chosen for $n=0$). However, in the
absence of such arguments, the variation of the subtraction scheme ought to be
considered to be its own source of systematic theoretical uncertainty.

##  IV Estimating Theory Uncertainties

###  IV.1 Profile Functions

The resummed, matched and NP-corrected cross section in (20) depends on a set
of scales characterizing the physics hierarchies present in our factorization
framework, and these should smoothly transition across the full $\tau$ domain
we study. Those in the perturbative cross section $\sigma_{c}^{PT}$ should of
course be chosen to live at values that minimize the logarithms present in the
hard, jet and soft functions in the $\tau$ domain most sensitive to
resummation effects, i.e. the tail region of the distribution. For smaller
values of $\tau$, towards the peak of the distribution, the full shape
function of (20) becomes necessary to describe the non-perturbative physics at
play, and we force all scales to plateau at some value $\mu_{0}$ just above
$\Lambda_{\rm{QCD}}$. For larger values of $\tau$, on the other hand, in the
far-tail region where matching to fixed-order QCD is required, we merge the
scales onto the hard scale $\mu_{H}$. All of these mergers are achieved with
particular choices of profile functions Ligeti  _et al._ (2008); Abbate  _et
al._ (2011) and, as mentioned above, demonstrating the impact of this choice
on the $\alpha_{s}$ extractions is one of our central messages in this work.
To that end we will present the explicit forms for the profiles designed in
Bell  _et al._ (2019) (‘2018 profiles’) and Abbate  _et al._ (2011) (‘2010
profiles’) in what follows.

### 2018 Profiles

In Bell  _et al._ (2019) we designed a set of profile functions for a
generalized class of event shapes that encompasses thrust—the angularity
distributions. Specifically for thrust the profiles, which were inspired by
those presented in Hoang  _et al._ (2015b, a), were chosen as:

$\displaystyle\mu_{H}$ $\displaystyle=e_{H}Q\,,$ $\displaystyle\mu_{S}(\tau)$
$\displaystyle=\left[1+e_{S}\,\theta(t_{3}-\tau)\left(1-\frac{\tau}{t_{3}}\right)^{2}\right]\mu_{\rm{run}}(\tau)\,,$
$\displaystyle\mu_{J}(\tau)$
$\displaystyle=\left[1+e_{J}\,\theta(t_{3}-\tau)\left(1-\frac{\tau}{t_{3}}\right)^{2}\right]\,\sqrt{\mu_{H}\,\mu_{\rm{run}}(\tau)}\,,$
$\displaystyle\mu_{ns}(\tau)$
$\displaystyle=\begin{cases}\frac{1}{2}\big{(}\mu_{H}+\mu_{J}(\tau)\big{)}&n_{s}=1\\\
\mu_{H}&n_{s}=0\\\ \frac{1}{2}\big{(}3\mu_{H}-\mu_{J}(\tau)\big{)}&n_{s}=-1\\\
\end{cases}\,,$ $\displaystyle\mu_{R}(\tau)$
$\displaystyle=R(\tau)\equiv\mu_{S}(\tau)\,\,\,\,\,\text{with}\,\,\,\,\,\mu_{0}\rightarrow
R_{0},$ (32)

where the function $\mu_{\rm{run}}(\tau)$ ensures that the profiles evolve
smoothly over the full $\tau$ domain. Its specific form is given by

$\mu_{\rm{run}}(\tau)=\begin{cases}\mu_{0}&\tau\leq t_{0}\\\
\zeta\Big{(}\tau;\\{t_{0},\mu_{0},0\\},\\{t_{1},0,2r\mu_{H}\\}\Big{)}&t_{0}\leq\tau\leq
t_{1}\\\ 2r\mu_{H}\tau&t_{1}\leq\tau\leq t_{2}\\\
\zeta\Big{(}\tau;\\{t_{2},0,2r\mu_{H}\\},\\{t_{3},\mu_{H},0\\}\Big{)}&t_{2}\leq\tau\leq
t_{3}\\\ \mu_{H}&\tau\geq t_{3}\end{cases}$ (33)

where $\zeta$ controls the interpolation between different regions, with

$\displaystyle\zeta\left(\tau;\\{t_{0},y_{0},r_{0}\\},\\{t_{1},y_{1},r_{1}\\}\right)$
(34) $\displaystyle\quad=\begin{cases}a+r_{0}(\tau-t_{0})+c(\tau-
t_{0})^{2}&\tau\leq\frac{t_{0}+t_{1}}{2}\\\ A+r_{1}(\tau-t_{1})+C(\tau-
t_{1})^{2}&\tau\geq\frac{t_{0}+t_{1}}{2}\end{cases}$

and the various coefficients therein are determined by the continuity
requirement of this function and its first derivative,

$\displaystyle a$ $\displaystyle=y_{0}+r_{0}t_{0}\,,$ $\displaystyle c$
$\displaystyle=2\,\frac{A-a}{(t_{0}-t_{1})^{2}}+\frac{3r_{0}+r_{1}}{2(t_{0}-t_{1})}\,,$
(35) $\displaystyle A$ $\displaystyle=y_{1}+r_{1}t_{1}\,,$ $\displaystyle C$
$\displaystyle=2\,\frac{a-A}{(t_{0}-t_{1})^{2}}+\frac{3r_{1}+r_{0}}{2(t_{1}-t_{0})}\,.$

As is clear, the shape of the scale $\mu_{\rm{run}}(\tau)$ in (33), and
therefore of all $\tau$-dependent scales in (32), is controlled by the
transition points $t_{i}$, which we choose to be

$\displaystyle t_{0}$ $\displaystyle=\frac{n_{0}}{Q},\,\,\,$ $\displaystyle
t_{2}=0.295\,n_{2}\,,$ $\displaystyle t_{1}$
$\displaystyle=\frac{n_{1}}{Q},\,\,\,$ $\displaystyle t_{3}=0.5\,n_{3}\,.$
(36)

The particular forms in (36) were designed somewhat empirically with
particular purposes: $t_{0,1}$ demarcate the boundary between the fully non-
perturbative and resummation regions; they roughly track the peak of the
differential distribution, and thus scale as the leading power correction
$1/Q$. Meanwhile $t_{2}$ approximates the crossover where singular and non-
singular contributions to the cross section are of equal magnitude, which
indicates that resummation effects become less important and should hence be
turned off. Finally, we set $t_{3}$ to be just below the maximal thrust value
of the spherically symmetric configuration, $\tau^{\text{sph}}=1/2$, so that
our predictions reduce to their fixed-order values slightly below this
kinematic endpoint.

### 2010 Profiles

In their earlier study Abbate  _et al._ (2011), the authors designed the
following set of profile functions:

$\displaystyle\mu_{H}$ $\displaystyle=e_{H}Q\,,$ $\displaystyle\mu_{S}(\tau)$
$\displaystyle=\begin{cases}\mu_{0}+\frac{b}{2t_{1}}\tau^{2},&\tau\leq
t_{1},\\\ b\tau+d,&t_{1}\leq\tau\leq t_{2},\\\
\mu_{H}-\frac{b}{1-2t_{2}}\left(\frac{1}{2}-\tau\right)^{2},&\tau\geq
t_{2}\end{cases},$ $\displaystyle\mu_{J}(\tau)$
$\displaystyle=\left[1+e_{J}\,\left(\frac{1}{2}-\tau\right)^{2}\right]\,\sqrt{\mu_{H}\,\mu_{S}(\tau)},$
$\displaystyle\mu_{ns}(\tau)$ $\displaystyle=\begin{cases}\mu_{H}&n_{s}=1\\\
\mu_{J}(\tau)&n_{s}=0\\\
\frac{1}{2}\big{(}\mu_{J}(\tau)+\mu_{S}(\tau)\big{)}&n_{s}=-1\\\
\end{cases}\,,$ $\displaystyle\mu_{R}(\tau)=R(\tau)$
$\displaystyle\equiv\begin{cases}R_{0}+\mu_{1}\tau+\mu_{2}\tau^{2},&\tau\leq
t_{1}\\\ \mu_{S}(\tau),&\tau\geq t_{1}\end{cases},$ (37)

with parameters $t_{i}$ that are akin to those introduced in (36), whereas the
parameters $b,d$ are fixed by the continuity requirement of the soft profile
and its first derivatives at $\tau=t_{i}$—see Abbate  _et al._ (2011) for
complete details. We also note that these profiles have been somewhat
superseded by the 2015 thrust analysis of Hoang  _et al._ (2015a), in which
the change of profile functions was found to have only a minor impact on the
final $\\{\alpha_{s},\Omega_{1}\\}$ extractions. We nevertheless decided to
compare against the 2010 profiles here, since we will see in the following
that this impact can be greater in different renormalon cancellation schemes.

Besides their overall functional form, the profiles in (32) and (37) depend on
multiple parameters that allow us to vary these scales in order to estimate
unknown perturbative corrections in all sectors of our calculation. The
parameters $e_{H,J,S}$ control the overall magnitude and width of the hard,
jet, and soft profile bands, $R_{0}\lesssim\mu_{0}$ ensures that no large
logarithms between the soft and renormalon scales arise while ensuring a non-
zero subtraction in the peak region Abbate  _et al._ (2011), the adjustable
parameters $n_{i}$ allow a variation about our estimates for the transition
points between the different regions, the parameter $r$ represents a variable
slope in the transition regions, and $n_{s}$ picks out different values of the
non-singular scale $\mu_{ns}$ entering the matching contributions to
$\sigma^{PT}_{c}$ in (II.2), whose variation probes the size of missing
higher-order terms in these fixed-order predictions. Varying between the
choices for $\mu_{ns}$ in (32) and (37) that track the scales $\mu_{J}$ and/or
$\mu_{S}$ also probes the missing effect of resumming sub-leading logarithms
present in the non-singular contribution. The choice of whether $\mu_{H}$ or
$\mu_{J}$ is the default $\mu_{ns}$ scale, and thus whether the range of
variations probe higher or lower values of $\mu_{ns}$, represents the primary
difference between the two sets of profile functions. Both are valid and free
choices based on considerations of improving perturbative convergence and
obtaining reasonable theory uncertainty estimates. As we will see below,
however, this choice can have a measurable effect on the (far-)tail region of
the distribution and therefore the $\alpha_{s}$ extractions that include these
bins.

| 2018 Profiles of Bell  _et al._ (2019) | 2010 Profiles of Abbate  _et al._ (2011)
---|---|---
$e_{H}$ | $0.5\leftrightarrow 2$ | $0.5\leftrightarrow 2$
$e_{J}$ | $-0.75\leftrightarrow 0.75$ | $-1\leftrightarrow 1$
$e_{S}$ | 0 | N.A.
$n_{0}$ |  $1\leftrightarrow 2$ GeV | N.A.
$n_{1}$ |  $8.5\leftrightarrow 11.5$ GeV |  $2\leftrightarrow 8$ GeV
$n_{2}$ | $0.9\leftrightarrow 1.1$ | $0.678\leftrightarrow 1.017$
$n_{3}$ | $0.8\leftrightarrow 0.9$ | N.A.
$\mu_{0}$ |  $0.8\leftrightarrow 1.2$ GeV |  $1.5\leftrightarrow 2.5$ GeV
$R_{0}$ |  $\mu_{0}-0.4$ GeV | $0.85\,\mu_{0}$
$r$ | $0.75\leftrightarrow 1.33$ | N.A.
$\delta c_{\tilde{S}}^{3}$ | $-1\leftrightarrow 1$ | $-1\leftrightarrow 1$
$\delta r^{2}$ | $-1\leftrightarrow 1$ | $-1\leftrightarrow 1$
$n_{s}$ | $\\{-1,0,1\\}$ | $\\{-1,0,1\\}$
Table 1:  Parameter ranges for both sets of profile functions we consider in
this work. In both scenarios, all parameters are chosen randomly within the
ranges shown, and the ‘central’ values are the centers of the given ranges,
except for the parameters $e_{H}$ and $r$, whose central values are given by
$e_{H}=r=1$. Note that we have refined this variation for the R0 scheme—see
Section V.5 for details.

Figure 3: The 2018 and 2010 profiles we implement. The top left figure
compares the central values for $\mu_{H,J,S,R}$, while the top right figure
shows the non-singular scale $\mu_{ns}$. Bottom plots show the variations of
the 2018 and 2010 profiles from a random scan of parameters in the ranges
shown in Tab. 1 (as $\mu_{H}$ controls the overall scale of all profiles, its
variation is indicated by the black arrows). In the bottom plots
$\mu_{S}<\mu_{J}$ is realized for any given set of scales, despite the fact
that the overall bands overlap. The R⋆ scales from (29) are shown in all
figures, and all plots correspond to $Q=m_{Z}$.
---
Figure 4: Top Row: Normalized thrust distribution across the $\tau$ domain
relevant for the $\alpha_{s}$ extractions in the R scheme at different
resummed and matched accuracies as indicated by the colors. The plots are
generated with 64 variations of the embedded profile parameters, with 2018
(2010) profiles shown in the left (right) column. Bottom Row: The same for the
R⋆ scheme.

Apart from these two sets of profile functions, we recall that our purpose
consists in comparing schemes that use different prescriptions to subtract the
leading soft renormalon, as described in Section III.2. Specifically, we
choose $\mu_{R}(\tau)=R(\tau)$ as shown in (32) and (37) in the R and R0
schemes, whereas we impose a cutoff on the subtraction scale $\mu_{R}$ in the
R∗ scheme via the modified profile $\mu_{R}(\tau)=R^{*}(\tau)$ given in (30)
(where the function $R(\tau)$ that enters (30) is the same function that is
used in the R and R0 schemes). Where to impose the cutoff is a free choice,
and we have chosen, for both sets of profile functions, to define

$R_{\text{\rm{max}}}=R(t_{1})\,,$ (38)

i.e. we freeze the growth of $\mu_{R}$ to the value it takes at $t_{1}$, which
marks the transition from the non-perturbative (peak) region of the
distribution to the resummed (tail) region. We freely admit this is an
arbitrary choice—it represents an attempt to limit the growth of the effective
shift $\zeta_{\text{eff}}$ as illustrated in Fig. 2 to a mild but measurable
extent and to study its impact on the $\alpha_{s}$ determination. This will be
a primary focus of our studies in the remainder of this paper. There are
certainly other choices of $R_{\text{max}}$ and/or the shape of the $\mu_{R}$
profile one could choose to study.

In order to estimate the total perturbative uncertainty, we adopt the
procedure used in Abbate  _et al._ (2011) and Bell  _et al._ (2019), where
random values of all parameters that enter the profile functions are scanned
over a predetermined range with each instantiation yielding a different
profile (including values of $t_{1}$ and thus $R_{\text{\rm{max}}}$ per (38)).
The overall envelope of these profiles is then taken as the total theory
uncertainty. The ranges we scan over for both sets of profile functions are
shown in Tab. 1.

The central values for the 2018 and 2010 $\mu_{H,J,S,R}$ profiles are shown in
the top left panel of Fig. 3, while the top right panel gives the central
values of the non-singular scale $\mu_{ns}$ for each choice of
$n_{s}\in\\{-1,0,1\\}$ in both scenarios. The results of 64 random scans of
profile parameters are then illustrated in the bottom two panels of Fig. 3 for
both sets of profile functions.444Our $\\{\alpha_{s},\Omega_{1}\\}$ fits in
Sec. V will actually scan over $\mathcal{O}(500-1000)$ random profile
parameters within the same ranges. Note that in all panels we have plotted the
renormalon scale with the modified R⋆ prescription of (29), as otherwise it
would be visibly indistinguishable from $\mu_{S}$ across the bulk of the
domain. The hard scale parameter $e_{H}$ furthermore controls the overall
scale of all profiles, and for illustration purposes the curves have been
normalised appropriately, whereas its variation is indicated by the black
arrows. All plots in Fig. 3 are shown for $Q=m_{Z}$.

From Fig. 3 one notices that the choices embedded in the two sets of profile
functions in (32) and (37), which are both legitimate and robust frameworks,
lead to a qualitatively different behavior of the various scales throughout
the bulk of the $\tau$ domain. This statement is especially true when
considering the non-singular scale $\mu_{ns}$ in the top right panel of Fig.
3. In Sec. V.3 we will study the impact this has on the $\alpha_{s}$
extractions.

###  IV.2 Predictions for Differential Distributions

In addition to the uncertainties from unknown higher-order corrections, which
we estimate via the procedure described above, we must also account for
systematic errors associated to the numerical extraction of some of the theory
parameters from event generators. In our setup there are two such parameters
of concern, the three-loop soft matching coefficient $c_{\tilde{S}}^{3}$,
which we take from Brüser  _et al._ (2018), and the
$\mathcal{O}(\alpha_{s}^{2})$ remainder functions $r^{2}$, which we have
extracted from EVENT2. Following Bell  _et al._ (2019), we assign an error
function to the central values found for these objects
$X\in\\{c_{\tilde{S}}^{3},r^{2}\\}$ as

$\displaystyle X$ $\displaystyle=X^{\rm{central}}+\begin{cases}\delta
X\;\Delta X^{\rm{upper}}&\left(\delta X>0\right)\\\ \delta X\;\Delta
X^{\rm{lower}}&\left(\delta X<0\right)\end{cases}\,,$ (39)

where $\delta X$ is varied between $\pm 1$ as presented in Tab. 1, and $\Delta
X$ represents the associated uncertainty on the extracted parameters, which is
given in (9) for $c_{\tilde{S}}^{3}$, whereas we follow the strategy described
in Bell  _et al._ (2019) to assign an uncertainty to the remainder function
$r^{2}$. This procedure allows us to account for systematic uncertainties
associated to these quantities, despite the fact that these errors are not the
primary focus of our study.

Given the theory inputs described above, we can use (20) to predict cumulative
and differential thrust distributions with well-defined theory
uncertainties.555We stress that all our predictions are based on the
cumulative cross section in (20). In order to determine uncertainty bands for
the differential distributions, we have calculated the derivative of _all_
cumulative distribution curves coming from individual profile variations and
then maximized/minimized these across the $\tau$ domain, as opposed to simply
taking the derivative of the max/min cumulative curves. We do so in Fig. 4 for
the R and R⋆ schemes in the top and bottom row, respectively, whereas we
relegate a discussion of the R0 scheme to Sec V.5 and App. C. Here the left
(right) column of panels corresponds to predictions made with the 2018 (2010)
profile choices, and for convenience we indicated the scale choice in each
renormalon cancellation scheme with a subscript. Each panel shows the
prediction for the differential cross section, normalized by the total cross
section and multiplied by a factor of $\tau$, across the $\tau$ domain that
will be relevant for the $\alpha_{s}$ fits in Sec. V. We use the
$\mathcal{O}(\alpha_{s}^{3})$ fixed-order hadronic expression from Chetyrkin
_et al._ (1996) for the total cross section normalization. The various colors
in Fig. 4 indicate different resummed and matched accuracies, ranging from
NLL${}^{\prime}+\mathcal{O}(\alpha_{s})$ in gray to our best
N3LL′$+\mathcal{O}(\alpha_{s}^{2})$ prediction in purple. All panels,
including Fig. 5, are produced with $Q=m_{Z}$, $\alpha_{s}(m_{Z})=0.11$,
$\Omega_{1}(R_{\Delta},R_{\Delta})=0.4$ GeV, and $R_{\Delta}=1.5$ GeV.

Broadly speaking, we observe excellent convergence between different resummed
and matched accuracies, using both 2018 and 2010 profile scales and
considering different renormalon cancellation schemes. This convergence is
evident even towards large $\tau$ values where, at least for the 2018
profiles, we observe a slight widening of the uncertainty bands around the
$t_{2}$ transition point, which controls the $\tau$-dependent scales in these
schemes via (33), and which we also vary according to Tab. 1. We also observe
from the bottom panels of Fig. 4 that logarithms of $\mu_{S}/\mu_{\delta}$,
which are non-zero above $t_{1}$ in the R⋆ scheme (but absent in the R
scheme), do not qualitatively impact the perturbative convergence of the
differential distributions—these logarithms, which are of $\mathcal{O}(1)$
across the relevant domain, therefore seem safe. Finally, we have checked that
all schemes in Fig. 4 exhibit excellent qualitative agreement when compared
against data in the relevant fitting regime.

Independent of the renormalon cancellation scheme, we also notice throughout
Fig. 4 that the error bands associated to 2010 profile variations are
generally smaller than those associated to 2018 profile variations; the latter
can thus be considered more conservative than the former. Of course, this
behavior could have been anticipated from the profile variations shown in Fig.
3, where the width of the 2018 profiles is already wider than its 2010
counterpart. This feature is largely due to the desire to predict a global set
of observables (angularities) in Bell  _et al._ (2019), where the 2018
profiles were determined, all with well-behaved convergence and uncertainties,
as opposed to thrust alone. Doing so required a somewhat more conservative set
of variations. We have, however, checked that narrowing the 2018 soft, jet,
and renormalon scale variations in Fig. 3 to widths heuristically similar to
the 2010 profiles does not alter our conclusions below in Sec. V. We have
therefore chosen to leave the 2018 variations in the range used in Bell  _et
al._ (2019).

In Fig. 5 we examine the differences between the differential distributions
shown in Fig. 4 more closely. Specifically, we display
$\tau/\sigma_{\rm{tot}}\,d\sigma_{c}/d\tau$ for the four schemes in in Fig. 4,
normalized to the central profiles of the R2010 scheme. All predictions are
made at N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$ accuracy with $Q=m_{Z}$.
As can be seen, scheme variations can lead to multi-percent effects on the
differential distributions, and these effects are especially pronounced
towards the far tails of the distributions. Later, in Sec. V, we will explore
the possibility of performing fits within a restricted $\tau$ domain isolated
more towards the purely dijet region (i.e. $\tau\leq 0.225$), where the
differences between the various schemes is less pronounced, rather than the
default fitting window that is shown in Fig. 5. One should also keep in mind
that the results in Fig. 5 are computed with fixed values of
$\alpha_{s}(m_{Z})$ and $\Omega_{1}(R_{\Delta},R_{\Delta})$ as reported above,
and therefore do not necessarily translate into overlapping best-fit values
for $\\{\alpha_{s},\Omega_{1}\\}$ for a given theory profile—this merely
provides a guide to a potentially more stable fit window that we will analyse
below.

Figure 5: Comparison of the differential distributions shown in Fig. 4 for
four renormalon-cancellation and profile-variation schemes, normalized to the
central profile of the R2010 scheme.

##  V Global $\\{\alpha_{s},\Omega_{1}\\}$ Extractions

We will now fit the differential distributions that we designed in the
previous sections to the global thrust data in order to extract the strong
coupling $\alpha_{s}$ and the dominant NP parameter $\Omega_{1}$. To do so, we
first collect the data and describe the method we use for this extraction,
before we present our main observations regarding the impact of varying
renormalon schemes and perturbative profile scale choices on the fit results.
We also perform a more detailed comparison of our theory framework to the one
that was used in prior SCET-based thrust analyses, and briefly explore other
sources of systematic uncertainties that may affect the extraction. Throughout
this section, we report values for the fit parameters at standard scales, i.e.
$\alpha_{s}(m_{Z})$ and $\Omega_{1}(R_{\Delta},R_{\Delta})$, where
$R_{\Delta}=1.5$ GeV is an arbitrary (perturbative) reference scale, although
we will suppress this scale dependence in the following for brevity.

###  V.1 Experimental Data

As is already evident from (1) and the related discussion in the introduction,
the fit parameters $\alpha_{s}$ and $\Omega_{1}$ will turn out to be fairly
correlated, since the effects of varying one can to some extent be compensated
by the other. In order to lift this degeneracy, it will be important to
include data at varying c.o.m. energies $Q$. We will therefore perform a
global fit to the available thrust data that includes 52 different datasets
with c.o.m. energies spread over $Q\in\\{35,207\\}$ GeV. Specifically, we
include data from ALEPH Heister  _et al._ (2004) at 91.2, 133, 161, 172, 183,
189, 200 and 206 GeV, DELPHI Abreu  _et al._ (1999, 2000); Abdallah  _et al._
(2003) at 45, 66, 76, 91.2, 133, 161, 172, 183, 189, 192, 196, 200, 202, 205
and 207 GeV, JADE Movilla Fernandez  _et al._ (1998) at 35 and 44 GeV, L3
Achard  _et al._ (2004) at 41.4, 55.3, 65.4, 75.7, 82.3, 85.1, 91.2, 130.1,
136.1, 161.3, 172.3, 182.8, 188.6, 194.4, 200 and 206.2 GeV, OPAL Ackerstaff
_et al._ (1997); Abbiendi  _et al._ (2000, 2005) at 91, 133, 161, 172, 177,
183, 189 and 197 GeV, SLD Abe  _et al._ (1995) at 91.2 GeV and TASSO
Braunschweig  _et al._ (1990) at 35 and 44 GeV. We note that this dataset
largely mimics the one that was used in Abbate  _et al._ (2011), in an effort
to minimize systematic differences between the two extractions that are
associated to the experimental input.

In the following our default fits will include 488 bins whose centers fall
within the $6/Q\leq\tau\leq 0.33$ domain. This fit window corresponds to the
default choice that was also used in Abbate  _et al._ (2011). Roughly
speaking, it corresponds to a domain which starts at the onset of the
resummation-sensitive tail region of the distribution, going up to the
kinematic endpoint of a three-parton final-state, beyond which further
corrections to those captured in the dijet factorization theorem (5) become
important.

###  V.2 Extraction Method

Figure 6: Extractions of $\\{\alpha_{s},\Omega_{1}\\}$ in the R2018 scheme
(top panel) and the R${}^{*}_{2018}$ scheme (bottom panel) for different
resummed and matched accuracies as indicated by the colors. Here $Q=m_{Z}$,
and we note that we have restricted
$\Omega_{1}\geq\Delta(R_{\Delta},R_{\Delta})=0.1$ GeV in these scans.
---
Figure 7: Extractions of $\\{\alpha_{s},\Omega_{1}\\}$ using
N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$ theory predictions fitted to
global thrust data in our default fitting window, $6/Q\leq\tau\leq 0.33$.
Upper panels compare the systematic impact of varying 2010 vs. 2018 profile
scales within a given renormalon cancellation scheme, while bottom panels vary
renormalon schemes within a given set of profile functions. Dashed (solid)
ellipses correspond to 68% (95%) C.L.. Note that $\Omega_{1}$ is formally
defined differently in the R and R∗ schemes.

Given the framework outlined above, we are now in the position to compare
various scheme-dependent theory predictions to the available thrust data, with
the aim of extracting values for the strong coupling constant $\alpha_{s}$ and
the leading NP shift parameter $\Omega_{1}$ in each setup. To assess the
quality of these fits quantitatively, we will perform a $\chi^{2}$ analysis at
the binned level:

$\chi^{2}\equiv\sum_{i,j}\,\Delta_{i}\,V_{ij}^{-1}\Delta_{j}\,,$ (40)

where we have defined $\Delta_{i}$ as the difference between the theoretical
prediction and the data in the $i$-th bin,

$\Delta_{i}\equiv\frac{1}{\sigma}\frac{d\sigma}{d\tau}(\tau_{i})\Bigr{\rvert}^{\text{exp}}-\frac{1}{\sigma}\frac{d\sigma}{d\tau}(\tau_{i})\Bigr{\rvert}^{\text{th}}\,,$
(41)

and where the theory predictions for a bin between
$\left[\tau_{1},\tau_{2}\right]$ (with $\tau_{2}>\tau_{1}$) are calculated
from the difference of the cumulative thrust distribution evaluated at the
endpoints of the bin, but with profile scales evaluated at the bin’s center
$\overline{\tau}\equiv(\tau_{1}+\tau_{2})/2$:

$\frac{1}{\sigma}\frac{d\sigma}{d\tau}(\tau_{i})\Big{|}^{\text{th}}_{\rm{MP}}\equiv\frac{1}{\sigma_{tot}}\frac{\sigma_{c}\big{(}\tau_{2},\mu_{a}(\overline{\tau})\big{)}-\sigma_{c}\big{(}\tau_{1},\mu_{a}(\overline{\tau})\big{)}}{\tau_{2}-\tau_{1}}\,,$
(42)

as advocated in Abbate  _et al._ (2011). Here $\mu_{a}$ refers to any of the
$\tau$-dependent scales present in our framework. We deem this the _midpoint_
(MP) binning procedure.

Then, for a given dataset, $V_{ij}$ incorporates the statistical and
systematic errors quoted by the experimental collaborations. The statistical
errors of each bin, $e^{\rm{stat}}_{i}$, are considered to be independent and
contribute to the diagonal entries of $V_{ij}$. However, the correlated
systematic uncertainties between different bins yield non-zero off-diagonal
entries. As no information on these correlations is available in the
literature, we estimate them using the Minimal Overlap Model (MOM), as
advocated by the LEP QCD working group Abbiendi  _et al._ (2005); Heister  _et
al._ (2004), and which was also employed in the analyses of Abbate  _et al._
(2011); Hoang  _et al._ (2015a). According to the MOM prescription, the off-
diagonal elements of $V_{ij}$ can be estimated as the minimum of the two
systematic errors of the individual bins $i$ and $j$, such that

$V_{ij}\big{|}_{\text{MOM}}=(e^{\rm{stat}}_{i})^{2}\delta_{ij}+\min(e_{i}^{\rm{sys}},e_{j}^{\rm{sys}})^{2}\,.$
(43)

For a given value of the c.o.m. energy $Q$ and profile function parameters, we
may then find the values $\\{\alpha_{s},\Omega_{1}\\}$ that minimize
$\chi^{2}$.

To estimate the theoretical uncertainty on the extracted
$\\{\alpha_{s},\Omega_{1}\\}$ values, we repeat the same procedure for
multiple random draws of the profile function parameters within the ranges
given in Tab. 1. In particular, for scans evaluated at
N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$ accuracy we perform
$\mathcal{O}(1000)$ scans, while we perform $\mathcal{O}(500)$ scans at lower
perturbative accuracies. Upon collecting all of the minimized
$\\{\alpha_{s},\Omega_{1}\\}$ pairs, we then define the theoretical error
estimate as the 68%/95% confidence-level (C.L.) ellipse for the two
parameters, i.e. by the ellipse centered at the coordinates
$\\{\mu_{\alpha},\mu_{\Omega}\\}$, with $\mu_{X}$ the mean value from all
theory draws, and with the standard covariance matrix for the dataset
multiplied by the 68th/95th percentile value for $\chi^{2}$ distributions with
two degrees of freedom ($\sim$2.279/5.991), serving as a weight matrix. We
then recognize the overall error ellipse defined by these two inputs (center
and weight matrix) as an object that can be parameterized as

$K_{\rm{theory}}=\begin{pmatrix}\sigma_{\alpha}^{2}&\rho_{\alpha\Omega}\,\sigma_{\alpha}\sigma_{\Omega}\\\
\rho_{\alpha\Omega}\,\sigma_{\alpha}\sigma_{\Omega}&\sigma_{\Omega}^{2}\end{pmatrix}\,,$
(44)

where we interpret $\sigma_{X}$ as weighted variances for the two parameters
$\alpha_{s}$ and $\Omega_{1}$, and with off-diagonal entries involving the
weighted correlation coefficients $\rho_{\alpha\Omega}$ of the two parameters.
We report these parameters for our global scans in the following sections. To
conclude, we interpret (44) as the area within which, upon randomly drawing
new sets of theory parameters and executing the $\chi^{2}$ minimization
procedure mentioned above, roughly 68%/95% of the extracted
$\\{\alpha_{s},\Omega_{1}\\}$ values will fall.

One could also consider incorporating a further experimental uncertainty as
described in Hoang  _et al._ (2015a), by constructing a $\chi^{2}$
distribution as a function of $\alpha_{s}$ and $\Omega_{1}$ around the minimum
$\chi^{2}$ obtained using the central profiles in Tab. 1. However, this
uncertainty is expected Hoang  _et al._ (2015a) to be much smaller in
comparison to the perturbative theory errors parameterized in (44) and, given
that our motivation consists in highlighting systematic theory uncertainties,
we will only quote the errors embedded in (44).

###  V.3 Results and Discussion

Profiles | Parameters |  $R$ (Default) |  $R^{\star}$ (Default) |  $R$ (Dijet) |  $R^{\star}$ (Dijet)
---|---|---|---|---|---
2018 Profiles | $\\{\alpha_{s},\Omega_{1}\\}$ | $\\{0.1124,0.383\\}$ | $\\{0.1134,0.364\\}$ | $\\{0.1126,0.379\\}$ | $\\{0.1135,0.353\\}$
$\\{\sigma_{\alpha},\sigma_{\Omega}\\}$ | $\\{0.0021,0.102\\}$ | $\\{0.0022,0.107\\}$ | $\\{0.0016,0.098\\}$ | $\\{0.0021,0.118\\}$
$\rho_{\alpha\Omega}$ | $-0.847$ | $-0.881$ | $-0.872$ | $-0.908$
2010 Profiles | $\\{\alpha_{s},\Omega_{1}\\}$ | $\\{0.1132,0.378\\}$ | $\\{0.1161,0.320\\}$ | $\\{0.1133,0.373\\}$ | $\\{0.1161,0.308\\}$
$\\{\sigma_{\alpha},\sigma_{\Omega}\\}$ | $\\{0.0012,0.045\\}$ | $\\{0.0024,0.075\\}$ | $\\{0.0010,0.049\\}$ | $\\{0.0023,0.072\\}$
$\rho_{\alpha\Omega}$ | $-0.633$ | $-0.758$ | $-0.703$ | $-0.731$
Table 2: Central values for the 95% C.L. $\\{\alpha_{s},\Omega_{1}\\}$
ellipses, and their associated weighted variances $\sigma_{\alpha,\Omega}$ and
correlation coefficient $\rho_{\alpha\Omega}$, in different renormalon-
cancellation and profile-variation schemes. Values for $\Omega_{1}$ are given
in units of GeV. The columns ‘Default’ and ‘Dijet’ indicate different fitting
windows—see text for details (the numbers in the ‘Default’ columns correspond
to the ellipses shown in Fig. 7).

We first study the convergence of the $\\{\alpha_{s},\Omega_{1}\\}$
extractions when progressively higher perturbative orders are used in the
theoretical predictions. Specifically, we display fit results for
NLL${}^{\prime}+\mathcal{O}(\alpha_{s})$ up to
N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$ resummed and matched accuracies
for the R2018 scheme (top panel) and the R${}^{\star}_{2018}$ scheme (bottom
panel) for $Q=m_{Z}$ in Fig. 6.666Note that we have restricted
$\Omega_{1}\geq\Delta(R_{\Delta},R_{\Delta})$ in our extraction code, hence
some lower-order scan results are artificially aligned at this boundary, but
we have verified that this does not affect any of the
N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$ fits shown in the following.
Also note that the polygonal, lower-order results represent the bounding
region for all fits in this plot, not just those lying within a certain
confidence level. As expected from Fig. 4, the convergence to the highest
perturbative accuracy is excellent in each scheme. In particular, the trend is
towards lower values of $\alpha_{s}$ and $\Omega_{1}$ for these $n=1$ schemes,
when the perturbative accuracy is increased. This observation is consistent
with prior analyses Abbate  _et al._ (2011); Hoang  _et al._ (2015a). We find
a similar behaviour for the R2010 and R${}^{\star}_{2010}$ schemes. However,
the pattern turns out to be rather different in the R0 schemes, as will be
described in App. C. Because of the instabilities found in these schemes, we
do not consider the R0 schemes any further in this section.

In Fig. 7 we show the result of global $\\{\alpha_{s},\Omega_{1}\\}$
extractions at the highest perturbative accuracy in four different renormalon-
cancellation and profile-variation schemes. For illustration purposes, we only
display two of these schemes in each panel. Although our framework does not
exactly correspond to the one used in Abbate  _et al._ (2011); Hoang  _et al._
(2015a) (a detailed comparison will be given in Sec. V.4 below), the results
in red for the R2010 scheme are certainly closest to these numbers, and they
may hence serve as a proxy for these results.

By comparing the uncertainty ellipses in all four schemes, the one in the
R2010 scheme clearly sticks out to be significantly smaller than the others.
From the upper left panel, we furthermore observe that the impact of switching
profiles is only moderate in the R scheme, which is in line with the
observations made in Hoang  _et al._ (2015a). This effect is, however, more
pronounced in the R∗ scheme, as can be seen in the upper right panel. More
importantly, the lower panels compare different renormalon-cancellation
schemes for the same set of profile functions (2010/2018 in the bottom
left/right panel). It is worth emphasizing that switching between the R and R∗
schemes changes the formal definition of the $\Omega_{1}$ parameter, and one
should therefore be careful in interpreting the vertical axis in these plots.
For what concerns the horizontal axis, on the other hand, $\alpha_{s}$ is of
course independent of the chosen renormalon scheme, and is defined by the
usual $\overline{\text{MS}}$ prescription. The definition of $\Omega_{1}$
_does_ , however, affect the $\alpha_{s}$ extraction in an indirect way, since
the fitting routine determines the best fit values of the correlated
$\\{\alpha_{s},\Omega_{1}\\}$ pair.

From the lower panels in Fig. 7 we observe that, relative to the R scheme, the
error ellipses are noticeably shifted to larger $\alpha_{s}$ values in the R∗
scheme, with the effect being more significant for 2010 profiles mainly
because of its lower value of $t_{1}$. On a more quantitative level, the
results from Fig. 7, which were obtained using the default fitting window with
$6/Q\leq\tau\leq 0.33$, are summarized in the columns labelled ‘Default’ in
Tab. 2. Here $\sigma_{\alpha}$ and $\sigma_{\Omega}$ are the weighted
variances of $\alpha_{s}$ and $\Omega_{1}$, respectively, whereas
$\rho_{\alpha\Omega}$ quantifies their correlation as described above.

Figure 8: Result of a global fit to thrust data using
N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$ theory predictions and a fitting
window $6/Q\leq\tau\leq 0.33$. The plots show the $\alpha_{s}-\Omega_{1}$
plane (Left) and $\alpha_{s}-\chi^{2}_{\rm{dof}}$ plane (Right) for four
different renormalon-cancellation and profile-variation schemes. The ellipses
in the left panel correspond to 95% C.L..

In Fig. 8 we combine all of the 95% C.L. ellipses from Fig. 7 into the left
panel, and show the corresponding $\alpha_{s}$-$\chi^{2}_{\rm{dof}}$ plane in
the right panel. We re-emphasize that the vertical axis should not be
misinterpreted in the left plot, since the definition of the $\Omega_{1}$
parameter depends on the chosen renormalon scheme. From the
$\chi^{2}_{\rm{dof}}$ distribution one observes that all schemes provide good
fits to the data, with the R${}^{*}_{2010}$ scheme, which yields the largest
$\alpha_{s}$ values, being slightly less preferred than the others. Moreover,
the left plot shows that, upon considering all four schemes, the extracted 95%
C.L. values of $\alpha_{s}$ and $\Omega_{1}$, while largely consistent with
one another, span a range that extends well beyond that of any one scheme
alone. This represents the _first core observation_ of our work, namely that
changing between different, well-defined renormalon-cancellation schemes
and/or profile scale variations can lead to noticeably different extractions
of $\\{\alpha_{s},\Omega_{1}\\}$. This may be viewed as an indication of
additional systematic theory uncertainties. On the other hand, a theoretically
motivated reason to stick to a single scheme could, of course, remove this
uncertainty. That being said, the schemes we include in Fig. 8 represent only
a subset of possible schemes exhibiting reasonably good perturbative
convergence and quality of fit—two criteria that could be used to prefer a
particular scheme.

The NP corrections embedded in our framework have been derived from a dijet
soft function that enters the factorization theorem (5). It has recently been
argued Luisoni  _et al._ (2021); Caola  _et al._ (2022a, b) that other sources
of NP corrections become relevant in the far-tail region of the distribution
that may not be negligible over the entire domain $6/Q\leq\tau\leq 0.33$ used
in the previous fits. Related to this, we observe in Fig. 5 that the scheme
dependence we consider has its most prominent effects in the far-tail region
as well. We therefore consider an alternative fit window with $6/Q\leq\tau\leq
0.225$ in the remainder of this section that concentrates more on purely dijet
events. The total number of bins in this setup is then reduced from 488 to
371. As an indication, the perturbative variations of the $Q=m_{Z}$ thrust
distributions in the different schemes shown in Fig. 5 can be reduced from
$\lesssim 12\%$ at the upper boundary of our default fit window ($\tau=0.33$)
to $\lesssim 5\%$ at the upper boundary of the reduced fit window
($\tau=0.225$). We examine next if this has a noticeable imprint on the
$\\{\alpha_{s},\Omega_{1}\\}$ extractions.

Figure 9: The same as in Fig 8 for the default fit window with
$6/Q\leq\tau\leq 0.33$ (solid contours) and a reduced dijet fit window with
$6/Q\leq\tau\leq 0.225$ (dashed contours).

In the two panels of Fig. 9 we compare results that were obtained using this
‘dijet fit window’ (dashed contours) to the ones with the default fit window
(solid contours) that were already shown in Fig. 8. In the left plot one sees
that this change has only a mild effect on the $\\{\alpha_{s},\Omega_{1}\\}$
extractions (as previously noted in Abbate  _et al._ (2011)), with the most
prominent effect being a shift of the R2018 ellipse to slightly larger
$\alpha_{s}$ values. As a result, the overall spread of the fit results among
the four considered schemes is just slightly reduced in this setup. These
observations are also given numerically in the last two columns of Tab. 2. In
the right panel of Fig. 9, on the other hand, one observes that the more
prominent effect of narrowing the fit window is a universal trend towards
lower $\chi^{2}_{\rm{dof}}$ values among all considered schemes, despite the
fact that the number of bins used in these analyses has been reduced
significantly. This improvement is especially noticeable for the
R${}^{\star}_{2010}$ scheme (in blue), where the overwhelming number of fits
drops below the $\chi^{2}_{\rm{dof}}=1$ contour, and which yields $\alpha_{s}$
values that are more compatible with the PDG world average than the other
schemes.

To summarize, we find that fits that are based on a more central dijet-type
$\tau$ domain seem to yield higher-quality results for
$\\{\alpha_{s},\Omega_{1}\\}$ extractions than those including data from the
far-tail region, where multi-jet events start to dominate. This represents the
_second core observation_ of our analysis, which suggests that precision fits
on a more limited dijet window may provide an alternative to the strategy
proposed e.g. in Nason and Zanderighi (2023), in the absence of a model-
independent understanding of NP corrections associated with tri- and multi-jet
events. From Fig. 9, on the other hand, we only see slight evidence for any
substantial qualitative improvement in the agreement between different schemes
when considering more central $\tau$ fits. In other words, over the range of
schemes we consider, using a narrower, higher-quality fit window does not by
itself remove the potential systematic uncertainty on
$\\{\alpha_{s},\Omega_{1}\\}$ coming from this scheme dependence.

###  V.4 Comparison to Prior Results

While there exist a number of dedicated thrust-based $\alpha_{s}$ extractions
in the literature (see e.g. Davison and Webber (2009); Abbate  _et al._
(2011); Gehrmann  _et al._ (2013); Hoang  _et al._ (2015a)), our framework is
particularly close to the one used in Abbate  _et al._ (2011); Hoang  _et al._
(2015a), and we will therefore perform a more detailed comparison to these
analyses in this section. First of all, both setups have in common that they
use methods from SCET to resum large logarithmic corrections in the dijet
limit, and the SCET predictions are matched to fixed-order QCD calculations to
account for non-singular corrections in the far-tail region of the
distributions. Even more importantly, both setups use a gapped dijet shape
function to account for NP corrections in a scheme that is free of the leading
soft renormalon. Higher-order corrections are furthermore estimated using a
scan over profile scale parameters, and the resulting theory predictions are
fitted to global thrust data for c.o.m energies $Q\in\\{35,207\\}$ GeV. Apart
from these similarities, there exist, however, a number of differences between
the two frameworks that we will now explain in detail, but which do not make a
significant impact on our conclusions:

  * •

There are two differences concerning the perturbative treatment described in
Sec. II. First, we only match the N3LL′ resummed prediction to
$\mathcal{O}(\alpha_{s}^{2})$ fixed-order calculations, due to our findings in
Sec. II.2 regarding instabilities of EERAD3 in the small $\tau$-region. In
contrast to this, Abbate  _et al._ (2011); Hoang  _et al._ (2015a) do include
the $\mathcal{O}(\alpha_{s}^{3})$ matching using EERAD3, though with lower
statistics and thus larger errors that overshadow these instabilities. In any
case they do not measurably affect the region used for the fits. The second
point concerns the three-loop soft constant $c_{\tilde{S}}^{3}$, for which
Abbate  _et al._ (2011); Hoang  _et al._ (2015a) used the Padé approximant in
(10), whereas we implemented the estimate in (9) that only became available in
2018 Brüser  _et al._ (2018). As we will show in the following section,
switching to the Padé value brings our predictions into even better
quantitative agreement with the ones from Abbate  _et al._ (2011); Hoang  _et
al._ (2015a).

  * •

Our NP treatment in the R schemes, on the other hand, closely resembles the
one used in Abbate  _et al._ (2011); Hoang  _et al._ (2015a) for what concerns
both the renormalon scheme definition and the profile functions used. The
differences between the two frameworks are absolutely minor in this respect,
and they concern e.g. only slightly different input values for the gap
parameter $\Delta(R_{\Delta},R_{\Delta})=0.05$ GeV at a reference scale
$R_{\Delta}=2$ GeV used in Abbate  _et al._ (2011); Hoang  _et al._ (2015a),
compared to our value $\Delta(R_{\Delta},R_{\Delta})=0.1$ GeV at
$R_{\Delta}=1.5$ GeV. Note that the gap parameter anyway gets absorbed into
$\Omega_{1}$, which we extract from data, according to (18).

  * •

The analyses in Abbate  _et al._ (2011); Hoang  _et al._ (2015a) furthermore
include a number of small effects (bottom and hadron masses, QED effects),
which we did not consider, mainly because we were not aiming at a competitive
$\alpha_{s}$ extraction in this work. In Abbate  _et al._ (2011); Hoang  _et
al._ (2015a) the authors showed that these have only a minor impact on the
fits, and given that they are unrelated to the main concern of our
paper—renormalon schemes and profile variations—we do not expect that these
will change any of our findings in a significant way.

  * •

In addition we use a slightly different method for calculating binned
distributions. Whereas Abbate  _et al._ (2011); Hoang  _et al._ (2015a)
integrate directly the resummed differential cross section, we start from the
cumulative distribution (3), and calculate the bins according to the midpoint
prescription in (42). As explained in Abbate  _et al._ (2011), this avoids the
presence of any spurious contributions, and the difference to the method from
Abbate  _et al._ (2011); Hoang  _et al._ (2015a) is expected to be at the sub-
percent level in the considered fit windows.

  * •

Finally, there are some minor differences in the applied fitting procedure.
While the two approaches use highly similar datasets, the authors of Abbate
_et al._ (2011); Hoang  _et al._ (2015a) calculate the uncertainty ellipses by
calculating the best-fit values for $\\{\alpha_{s},\Omega_{1}\\}$ obtained
from 500 profile variations, drawing the minimum ellipse surrounding the
minimum convex polygon encapsulating all of these points, and centered at the
average of max and min $\alpha_{s}/\Omega_{1}$ values obtained in each
direction. Our method for obtaining $K_{\rm{theory}}$ (and the associated
error ellipse) is instead outlined in Sec. V.2 above. Moreover, the authors of
Abbate  _et al._ (2011); Hoang  _et al._ (2015a) quote an additional
experimental uncertainty in their final results, which is however much smaller
than the theory error we also implement.

Despite these differences, it is a non-trivial result that our final numbers
in the R schemes shown in Fig. 8 and summarized in Tab. 2 are very similar to
the ones obtained in Abbate  _et al._ (2011); Hoang  _et al._ (2015a). To make
this statement more quantitative, we quote the final numbers obtained in
Abbate  _et al._ (2011),

$\displaystyle\alpha_{s}(m_{Z})$ $\displaystyle=0.1135\pm 0.0011\,,$
$\displaystyle\Omega_{1}(R_{\Delta},R_{\Delta})$ $\displaystyle=(0.323\pm
0.051)\text{\leavevmode\nobreak\ GeV}\,,$ (45)

which should be compared to our R2010 numbers, whereas the results of the
thrust analysis in Hoang  _et al._ (2015a),

$\displaystyle\alpha_{s}(m_{Z})$ $\displaystyle=0.1128\pm 0.0012\,,$
$\displaystyle\Omega_{1}(R_{\Delta},R_{\Delta})$ $\displaystyle=(0.322\pm
0.068)\text{\leavevmode\nobreak\ GeV}\,,$ (46)

were derived in a setup that is closer to our R2018 scheme. Note that the NP
parameter is evaluated here at a slightly larger reference scale
$R_{\Delta}=2$ GeV than for the numbers quoted above (the evolution to $1.5$
GeV is only a minor effect that reduces the value of $\Omega_{1}$ by $\sim
4\%$). Roughly speaking, these numbers translate into an error ellipse that is
similar in size to the red one in Fig. 8, but slightly shifted downwards. We
will, in fact, identify one effect that drives the ellipse into this direction
in the following section.

In view of this agreement, we may thus state that our analysis, which is based
on a completely independent set of codes, for the first time confirms the
results of Abbate  _et al._ (2011); Hoang  _et al._ (2015a), regarding the
degree to which their extracted value of $\alpha_{s}$ sits lower than the PDG
world average. We consider this important cross-check to be the _third core
conclusion_ of our work.

###  V.5 On the Three-Loop Soft Constant $c_{\tilde{S}}^{3}$

As mentioned in Sec. II, there exist two approximations for the three-loop
soft matching coefficient $c_{\tilde{S}}^{3}$ in the literature, namely the
EERAD3 extraction from Brüser  _et al._ (2018) we have reported in (9) and
which is used in our analysis, and the Padé approximant given in (10) that was
implemented in prior $\alpha_{s}$ fits Abbate  _et al._ (2011); Hoang  _et
al._ (2015a). A priori the impact of varying this constant should be small,
given that it represents a three-loop effect. In this section we point out at
least two ways that this statement should be qualified.

First, we study the role of $c_{\tilde{S}}^{3}$ in
$\\{\alpha_{s},\Omega_{1}\\}$ fits, and provide a more direct comparison
between our results presented above and those of Abbate  _et al._ (2011);
Hoang  _et al._ (2015a). To this end, we compare in Fig. 10 the 95% C.L.
ellipses for $\\{\alpha_{s},\Omega_{1}\\}$ obtained when
N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$ theory predictions are fitted to
$Q=m_{Z}$ datasets, in the R2010 scheme. In particular, the red ellipse
corresponds to our default choice of $c_{\tilde{S}}^{3}$ given in (9),
including its error, whereas the brown ellipse in Fig. 10 uses instead the
Padé approximant in (10), including again its quoted error. The errors quoted
in (9)-(10) are by themselves actually very minor contributions to the total
error ellipses in Fig. 10. The differences in the two central values, however,
is much larger. As a result, one notices a significant downward shift of our
default $95\%$ C.L. ellipse, which in fact brings our numbers into even better
agreement with the results from Abbate  _et al._ (2011), as can be verified by
comparing to the numbers quoted in (V.4). As the two values of
$c_{\tilde{S}}^{3}$ differ by more than $3\sigma$, on the other hand, the two
ellipses barely overlap. We thus consider the variation of $c_{\tilde{S}}^{3}$
as another systematic theory uncertainty that may be larger than previously
expected (especially due to the instablities in EERAD3 at small $\tau$
described in Sec. II.2), although its impact on the determination of
$\alpha_{s}$ itself may only be limited as the plot in Fig. 10 suggests.

Figure 10: Comparison of two $\\{\alpha_{s},\Omega_{1}\\}$ extractions that
either use the EERAD3 value $c_{\tilde{S}}^{3}=-19988\pm 5440$ (red) or the
Padé approximant $c_{\tilde{S}}^{3}=691\pm 1000$ (brown).

Figure 11: The same as in Fig. 4 for the R0 scheme. The upper plots use the
EERAD3 value for $c_{\tilde{S}}^{3}$ from (9), and the lower plots the Padé
approximant given in (10). Left (right) plots refer to 2018 (2010) profile
functions.

In addition there is a second reason why a better determination of the three-
loop soft constant may be warranted, which is related to the R0 renormalon
scheme that we introduced in Sec. III.2, but which we largely disregarded in
this section because of stability issues. To illustrate these, we show in Fig.
11 the analogous plots to the ones in Fig. 4, but for the R0 schemes. In
particular, the upper panels show the result for 2018 profiles (left) and 2010
profiles (right), when the default value of $c_{\tilde{S}}^{3}$ from (9) is
used. Away from the central $\tau$ domain, one clearly observes that the
theoretical predictions are not improved when increasingly higher perturbative
orders are included. This effect is particularly pronounced for the 2018
profile scans, and it is true despite the fact that we have already tuned the
ranges for some of the parameters in Tab. 1, and only allowed for $0.25\leq
e_{H}\leq 1.25$ and $-0.5\leq e_{J}\leq 0.5$ ($-0.75\leq e_{J}\leq 0.75$)
variations in the R${}^{0}_{2018}$ (R${}^{0}_{2010}$) scenarios. We recall
that the R0 scheme is special, since it is sensitive to one higher power of
logarithms in its subtraction terms and, critical to the present discussion,
the three-loop soft constant $c_{\tilde{S}}^{3}$, which is not yet exactly
known—cf. App. B. We therefore expect that a concrete determination of this
constant, and perhaps a more refined set of profile variations, could
eventually stabilize these curves as well. Regardless, to probe the impact of
different $c_{\tilde{S}}^{3}$ values on the distributions we have included in
the lower panels of Fig. 11 the corresponding distributions when the Padé
approximant in (10) is used. Here one notices that the convergence is somewhat
improved in comparison to the corresponding upper plots. Indeed, we have also
checked more generally that the breadth of the purple N3LL′
$+\mathcal{O}(\alpha_{s}^{2})$ bands in the R(⋆) schemes of Fig. 4 is also
artificially enhanced due to the present uncertainty on $c_{\tilde{S}}^{3}$
given in (9). Hence we can readily conclude that $c_{\tilde{S}}^{3}$ has a
visually noticeable impact on the quality of perturbative convergence across
all the considered schemes.

To summarize, the three-loop soft matching coefficient $c_{\tilde{S}}^{3}$ is
the only N3LL′ ingredient that is currently not known precisely. Despite being
a higher-order coefficient, the spread of the two existing estimates for this
coefficient has a noticeable impact on both the $\\{\alpha_{s},\Omega_{1}\\}$
extraction and the perturbative stabilty of the distributions in the R0
scheme, whose renormalon cancellation terms explicitly depend on
$c_{\tilde{S}}^{3}$. This issue also grows in importance in view of our
observations in Sec. II about the accuracy of EERAD3 calculations for small
$\tau$ values, which the estimate in (9) has relied upon. A resolution could
come from an independent numerical extraction, or, ideally, a complete
analytic calculation of the three-loop thrust soft function, of which partial
results have already been published Baranowski  _et al._ (2022); Chen  _et
al._ (2020). These observations about the importance of a more precise
$c_{\tilde{S}}^{3}$ determination represent the _fourth core conclusion_ of
our work.

##  VI Conclusion

We have revisited extractions of the strong coupling $\alpha_{s}$ from a
global fit to thrust data motivated by the current tension between SCET-based
extractions in Abbate  _et al._ (2011); Hoang  _et al._ (2015a) and the PDG
world average. To this end, we have used N3LL′ resummed theory predictions
that are matched to $\mathcal{O}(\alpha_{s}^{2})$ fixed-order calculations,
and we have implemented a gapped, renormalon-free, dijet shape function to
account for non-perturbative (NP) corrections. In the tail region of the
distributions, where the fits are performed, the dominant NP effects then
manifest as a shift of the distributions that is driven by the first moment of
the shape function $\Omega_{1}$. The fits are therefore formulated as two-
parameter extractions in the $\alpha_{s}-\Omega_{1}$ plane.

Whereas recent analyses with similar motivations have focused on NP effects
from three-jet configurations Luisoni  _et al._ (2021); Caola  _et al._
(2022a, b) that are most relevant in the far-tail region of the distributions,
our study addressed systematic effects in the dijet factorization theorem
itself. Specifically, we examined the impact from changing between different
renormalon cancellation schemes and profiled scale variations. Concerning the
former, we observed that the specific implementation of a renormalon-free gap
parameter that models the minimal energy of a hadronic final state leads to an
effective $\tau$-dependent shift as shown in Fig. 2. We then selected and
studied alternative renormalon schemes, following the approach in Bachu  _et
al._ (2021), that tame the growth of this shift for larger $\tau$ values, and
which amounts to a $\mathcal{O}(10\%)$ variation of the NP corrections in the
considered fit windows. As for scale variations, we implemented two different
sets of profile functions that are very close to the ones used previously in
Abbate  _et al._ (2011) and Hoang  _et al._ (2015a), respectively. While the
choice of profile functions did not make a major difference in those
extractions, we found that it can be magnified if varied together with the
renormalon scheme. In total this defines six different combinations of
renormalon-cancellation and profile-variation schemes.

As our predictions in the R0 schemes, which are more sensitive to higher-order
terms in the subtractions, turned out to be less stable, we focused on the
remaining four schemes for the $\\{\alpha_{s},\Omega_{1}\\}$ extractions. In
order to lift the degeneracy between the fit parameters, we performed a global
fit to thrust data with c.o.m. energies $Q\in\\{35,207\\}$ GeV. Our results
are shown in Fig. 8 and Tab. 2. In particular, we found that they are
consistent with previous extractions when the same renormalon-cancellation and
profile-scale choices are implemented, which can be viewed as an independent
confirmation of the analyses in Abbate  _et al._ (2011); Hoang  _et al._
(2015a). The scheme dependence of our fit results in the
$\alpha_{s}-\Omega_{1}$ plane, on the other hand, turned out to be
unexpectedly large, possible hinting at additional systematic theory
uncertainties. Whereas we do not advance any one particular scheme to be
favored over others on theoretical grounds, we remark that some of our
$\alpha_{s}$ extractions are more consistent with the PDG average. This
provides significant motivation for further scrutinizing the observed scheme
dependence of the dijet predictions in the future.

Since the scheme dependence is more pronounced in the far-tail region of our
predictions, we also performed fits that concentrate more on a central dijet-
like $\tau$ domain. While these fits turned out to be of higher quality (see
Fig. 9), the pattern in the $\alpha_{s}-\Omega_{1}$ plane did not change
significantly in this case. Finally, we explored the impact of the three-loop
soft matching coefficient $c_{\tilde{S}}^{3}$ on both the
$\\{\alpha_{s},\Omega_{1}\\}$ extractions and the perturbative stability of
the R0 schemes, and we argued that an improved determination of this quantity
would help to further improve the theoretical predictions.

As our primary goal consisted in bringing the scheme dependence of the SCET-
based $\alpha_{s}$ extractions to attention, we did not include a number of
sub-leading effects that will not change the main conclusions of our analysis,
but which may nevertheless be included in future precision fits to thrust
data. These include, first and foremost, the $\mathcal{O}(\alpha_{s}^{3})$
remainder function $r_{c}^{3}$ and the aforementioned soft three-loop constant
$c_{\tilde{S}}^{3}$, which we could not extract using EERAD3 because of
instabilities in the small $\tau$ region. Furthermore, one may account for QED
effects or non-zero bottom and hadron masses that were included in the
analysis of Abbate  _et al._ (2011); Hoang  _et al._ (2015a). Finally,
depending on the chosen fit window, it could be helpful to incorporate
resummation effects at the Sudakov shoulder $\tau=1/3$ Bhattacharya  _et al._
(2022).

Our findings also motivate renewed attention to sub-leading power
contributions and the associated resummation of logarithmic effects. As
discussed in Sec. IV.1, the main difference between the considered sets of
profile functions is sourced by the non-singular scale $\mu_{ns}$, which was
chosen to track the other dynamical scales in order to partially account for
missing resummation effects that appear in the non-singular contribution.
Hence a dedicated resummation analysis of next-to-leading power corrections
could eliminate the need to correlate $\mu_{ns}$ to the other scales, and in
this way possibly alleviate the systematic uncertainty due to scale-profile
choices that our study indicates. Resummations at sub-leading power to the
thrust distribution have been studied e.g. in Moult  _et al._ (2018, 2020);
Beneke  _et al._ (2022).

As a final point, we mention that, regardless of the rich physics discussed in
this analysis, a definitive extraction of $\alpha_{s}$ from $e^{+}e^{-}$ event
shapes will have to lift its degeneracy with the NP shift parameter
$\Omega_{1}$. While we have used thrust data from different c.o.m energies $Q$
for this purpose in the current analysis, an additional promising method
consists in extending the class of observables to $e^{+}e^{-}$ angularity
distributions $\tau_{a}$ Berger  _et al._ (2003). These observables generalize
the thrust variable, and their dominant NP effect is controlled by the same
parameter $\Omega_{1}$ (see (1)) with a coefficient Berger and Sterman (2003);
Berger and Magnea (2004); Lee and Sterman (2007); Becher and Bell (2014)

$c_{\tau_{a}}=\frac{2}{1-a}\,\,\,\,\,\text{with}\,\,\,\,\,a<1$ (47)

that depends on the angularity $a$ ($a=0$ corresponds to thrust). This allows
for the possibility of performing $\\{\alpha_{s},\Omega_{1}\\}$ extractions
globally in _both_ $Q$ and the new theory parameter $a$, offering a dual probe
in the decorrelation effort. Since angularities have been resummed and matched
to N2LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$ accuracy Bell  _et al._
(2019), a precision fit to available data Achard  _et al._ (2011) is well
overdue.

## Acknowledgements

We thank Andre Hoang, Vicent Mateu and Iain Stewart for helpful comments and
discussions, and in particular for pointing us to the different renormalon
schemes introduced in Bachu  _et al._ (2021). We thank Thomas Gehrmann for
help with EERAD3, and Boram Yoon for advice and assistance with our large
numerical computations. The work of C.L., Y.M. and B.Y. was supported by the
U.S. Department of Energy, Office of Science, through the Office of Nuclear
Physics, an Early Career Research Award (C.L.) and the LDRD Program at Los
Alamos National Laboratory. The research of G.B. was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762
- TRR 257. J.T. gratefully acknowledges funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No. 101022203, and guest support from the T-2 group at
LANL. This research used resources provided by the Los Alamos National
Laboratory Institutional Computing Program. LANL is operated by Triad National
Security, LLC, for the National Nuclear Security Administration of the U.S.
Department of Energy under Contract No. 89233218CNA000001. J.T. also thanks
the HPC and HEP teams at the DAMTP for their assistance with local cluster
computing. J.T. and C.L. are grateful to the organizers and participants of
the INT workshops ‘New physics searches at the precision frontier’ and
‘Probing QCD at High Energy and Density with Jets’ for their hospitality and
discussions during the completion of this work.

##  Appendix A Perturbative Expansions

Here we expand the prediction (5) for the resummed $\tau$ cross section in
fixed order to $\mathcal{O}(\alpha_{s}^{3})$, which we will use to extract the
non-singular remainder functions from numerical computations. It is expedient
to begin by computing the fixed-order expansion of the Laplace-transformed
cross section,

$\widetilde{\sigma}(\nu)=\frac{1}{\sigma_{0}}\int_{0}^{\infty}d\tau
e^{\nu\tau}\frac{d\sigma}{d\tau}\,,$ (48)

which obeys the factorized form,

$\widetilde{\sigma}(\nu)=H(L_{H},\mu)\tilde{J}^{2}(\tilde{L}_{J},\mu)\widetilde{S}(\tilde{L}_{S},\mu)\,,$
(49)

where $L_{H}=\ln(\mu/Q)$, $L_{J}=\ln(\mu^{j}e^{\gamma_{E}}\nu/Q^{j})$ (where
$j=2$), $L_{S}=\ln(\mu e^{\gamma_{E}}\nu/Q)$, and $\mu$ is generic
factorization scale. Each function $H,\widetilde{J},\widetilde{S}$ obeys a
perturbative expansion predicted by their RG evolution,

$F(L_{F},\mu)=F(0,Q_{F})e^{K_{F}(\mu,\mu_{F})}\,,$ (50)

where $Q_{H}=Q,Q_{J}=Qe^{\gamma_{E}}/\nu^{1/j}$, and
$Q_{S}=Qe^{\gamma_{E}}/\nu$, and $K_{F}$ is given by

$K_{F}(\mu,Q_{F})=-j_{F}\kappa_{F}K_{\Gamma}(\mu,Q_{F})+K_{\gamma_{F}}(\mu,Q_{F})\,,$
(51)

with $\kappa_{H}=4,\kappa_{J}=-2$ and $\kappa_{S}=4$. $K_{\Gamma},K_{\gamma}$
have the fixed-order expansions Almeida  _et al._ (2014)

$\displaystyle K_{\Gamma}(\mu,Q_{F})$
$\displaystyle=\frac{\alpha_{s}(\mu)}{4\pi}\frac{\Gamma_{0}}{2}\ln^{2}\frac{\mu}{Q_{F}}$
(52)
$\displaystyle\quad+\Bigl{(}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{)}^{2}\Bigl{(}\frac{1}{3}\Gamma_{0}\beta_{0}\ln^{3}\frac{\mu}{Q_{F}}+\frac{\Gamma_{1}}{2}\ln^{2}\frac{\mu}{Q_{F}}\Bigr{)}$
$\displaystyle\quad+\Bigl{(}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{)}^{3}\biggl{[}\frac{1}{3}\Gamma_{0}\beta_{0}^{2}\ln^{4}\frac{\mu}{Q_{F}}+\frac{1}{3}(\Gamma_{0}\beta_{1}+2\Gamma_{1}\beta_{0})\ln^{3}\frac{\mu}{Q_{F}}+\frac{\Gamma_{2}}{2}\ln^{2}\frac{\mu}{Q_{F}}\biggr{]}\,,$

and

$\displaystyle K_{\gamma_{F}}(\mu,Q_{F})$
$\displaystyle=\frac{\alpha_{s}(\mu)}{4\pi}\gamma_{F}^{0}\ln\frac{\mu}{Q_{F}}$
(53)
$\displaystyle\quad+\Bigl{(}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{)}^{2}\Bigl{(}\gamma_{F}^{0}\beta_{0}\ln^{2}\frac{\mu}{Q_{F}}+\gamma_{F}^{1}\ln\frac{\mu}{Q_{F}}\Bigr{)}$
$\displaystyle\quad+\Bigl{(}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{)}^{3}\biggl{[}\frac{4}{3}\gamma_{F}^{0}\beta_{0}^{2}\ln^{3}\frac{\mu}{Q_{F}}+\gamma_{F}^{0}\beta_{1}+2\gamma_{F}^{1}\beta_{0})\ln^{2}\frac{\mu}{Q_{F}}+\gamma_{F}^{2}\ln\frac{\mu}{Q_{F}}\biggr{]}\,,$

where $\Gamma_{n},\gamma_{F}^{n}$ are the coefficients in the perturbative
expansions of the anomalous dimensions,

$\Gamma_{\text{cusp}}[\alpha_{s}(\mu)]=\sum_{n=0}^{\infty}\Bigl{(}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{)}^{n+1}\Gamma_{n}\,,\quad\gamma_{F}[\alpha_{s}(\mu)]=\sum_{n=0}^{\infty}\Bigl{(}\frac{\alpha_{s}(\mu)}{4\pi}\Bigr{)}^{n+1}\gamma_{F}^{n}\,.$
(54)

In (50), the functions $F$ evaluated at their “natural” scale $Q_{F}$ have no
explicit logs in their expansions:

$\displaystyle
F(0,Q_{F})=1+\sum_{n=1}^{\infty}\Bigl{(}\frac{\alpha_{s}(Q_{F})}{4\pi}\Bigr{)}^{n}c_{F}^{n}\,.$
(55)

Combining all the pieces in (49), we can express the fixed-order Laplace-space
cross section in the form

$\widetilde{\sigma}(\nu)=\widetilde{C}(L_{\nu},Q)e^{K(L_{\nu},Q)}\,,$ (56)

with $L_{\nu}=\ln(e^{\gamma_{E}}\nu)$, and the ‘constant’ in front has the
form

$\widetilde{C}(L,Q)=1+\frac{\alpha_{s}(Q)}{4\pi}C_{10}+\Bigl{(}\frac{\alpha_{s}(Q)}{4\pi}\Bigr{)}^{2}(C_{20}+C_{21}L)+\Bigl{(}\frac{\alpha_{s}(Q)}{4\pi}\Bigr{)}^{3}(C_{30}+C_{31}L+C_{32}L^{2})+\cdots\,,$
(57)

where the extra log terms come from expanding out $\alpha_{s}(Q_{J,S})$
appearing in $F(Q_{J,S})$ in (55) in powers of $\alpha_{s}(Q)$. These
coefficients are given up to $\mathcal{O}(\alpha_{s}^{3})$ by

$\displaystyle C_{10}$
$\displaystyle=c_{H}^{1}+2c_{\tilde{J}}^{1}+c_{\tilde{S}}^{1}\,,$ (58)
$\displaystyle C_{20}$
$\displaystyle=c_{H}^{2}+2c_{\tilde{J}}^{2}+(c_{\tilde{J}}^{1})^{2}+c_{\tilde{S}}^{2}+c_{H}^{1}(2c_{\tilde{J}}^{1}+c_{\tilde{S}}^{1})+2c_{\tilde{J}}^{1}c_{\tilde{S}}^{1}\,,\quad
C_{21}=2\beta_{0}\Bigl{(}\frac{2c_{\tilde{J}}^{1}}{j}+c_{\tilde{S}}^{1}\Bigr{)}\,,$
$\displaystyle C_{30}$
$\displaystyle=c_{H}^{3}+2(c_{\tilde{J}}^{3}+c_{\tilde{J}}^{2}c_{\tilde{J}}^{1})+c_{\tilde{S}}^{3}+c_{H}^{2}(2c_{\tilde{J}}^{1}+c_{\tilde{S}}^{1})+c_{H}^{1}[2c_{\tilde{J}}^{2}+(c_{\tilde{J}}^{1})^{2}+c_{\tilde{S}}^{2}+2c_{\tilde{J}}^{1}c_{\tilde{S}}^{1}]+2c_{\tilde{J}}^{1}c_{\tilde{S}}^{2}+[2c_{\tilde{J}}^{2}+(c_{\tilde{J}}^{1})^{2}]c_{\tilde{S}}^{1}\,,$
$\displaystyle C_{31}$
$\displaystyle=2\beta_{1}\Bigl{(}\frac{2c_{\tilde{J}}^{1}}{j}+c_{\tilde{S}}^{1}\Bigr{)}+4\beta_{0}\Bigl{(}\frac{2c_{\tilde{J}}^{2}+(c_{\tilde{J}}^{1})^{2}}{j}+c_{\tilde{S}}^{2}\Bigr{)}+2\beta_{0}\Bigl{[}\frac{2c_{\tilde{J}}^{1}}{j}(c_{H}^{1}+c_{\tilde{S}}^{1})+(c_{H}^{1}+2c_{\tilde{J}}^{1})c_{\tilde{S}}^{1}\Bigr{]}\,,\quad
C_{32}=4\beta_{0}^{2}\Bigl{(}\frac{2c_{\tilde{J}}^{1}}{j^{2}}+c_{\tilde{S}}^{1}\Bigr{)}\,.$

Meanwhile, the expansion of the exponent $K$ in (56) is given by

$\displaystyle K(L,Q)$
$\displaystyle=\quad\Bigl{(}\frac{\alpha_{s}(Q)}{4\pi}\Bigr{)}\bigl{(}K_{12}L^{2}+K_{11}L\bigr{)}$
(59)
$\displaystyle\quad+\Bigl{(}\frac{\alpha_{s}(Q)}{4\pi}\Bigr{)}^{2}\bigl{(}K_{23}L^{3}+K_{22}L^{2}+K_{21}L)$
$\displaystyle\quad+\Bigl{(}\frac{\alpha_{s}(Q)}{4\pi}\Bigr{)}^{3}\bigl{(}K_{34}L^{4}+K_{33}L^{3}+K_{32}L^{2}+K_{31}L\bigr{)}+\cdots\,,$

with coefficients given up to $\mathcal{O}(\alpha_{s}^{3})$ by

$\displaystyle K_{12}$
$\displaystyle=-\frac{\Gamma_{0}}{2}\Bigl{(}\frac{2\kappa_{J}}{j}+\kappa_{S}\Bigr{)}\,,\qquad
K_{11}=\frac{2\gamma_{J}^{0}}{j}+\gamma_{S}^{0}\,,$ (60) $\displaystyle
K_{23}$
$\displaystyle=-\frac{1}{3}\Gamma_{0}\beta_{0}\Bigl{(}\frac{2\kappa_{J}}{j^{2}}+\kappa_{S}\Bigr{)}\,,\quad
K_{22}=-\frac{\Gamma_{1}}{2}\Bigl{(}\frac{2\kappa_{J}}{j}+\kappa_{S}\Bigr{)}+\Bigl{(}\frac{2\gamma_{J}^{0}}{j^{2}}+\gamma_{S}^{0}\Bigr{)}\beta_{0}\,,\quad
K_{21}=\frac{2\gamma_{J}^{1}}{j}+\gamma_{S}^{1}\,,$ $\displaystyle K_{34}$
$\displaystyle=-\frac{1}{3}\Gamma_{0}\beta_{0}^{2}\Bigl{(}\frac{2\kappa_{J}}{j^{3}}+\kappa_{S}\Bigr{)}\,,\quad
K_{33}=-\frac{1}{3}(\Gamma_{0}\beta_{1}+2\Gamma_{1}\beta_{0})\Bigl{(}\frac{2\kappa_{J}}{j^{2}}+\kappa_{S}\Bigr{)}+\frac{4}{3}\beta_{0}^{2}\Bigl{(}\frac{2\gamma_{J}^{0}}{j^{3}}+\gamma_{S}^{0}\Bigr{)}\,,$
$\displaystyle K_{32}$
$\displaystyle=-\frac{\Gamma_{2}}{2}\Bigl{(}\frac{2\kappa_{J}}{j}+\kappa_{S}\Bigr{)}+\Bigl{(}\frac{2\gamma_{J}^{0}}{j^{2}}+\gamma_{S}^{0}\Bigr{)}\beta_{1}+2\beta_{0}\Bigl{(}\frac{2\gamma_{J}^{1}}{j^{2}}+\gamma_{S}^{1}\Bigr{)}\,,\quad
K_{31}=\frac{2\gamma_{J}^{2}}{j}+\gamma_{S}^{2}\,,$

recalling $\gamma_{S}^{0}=0$ (though we have kept it to make the patterns more
obvious).

Plugging the expansions (58)-(60) into (56), we obtain the fixed-order
expansion of the Laplace-space cross section,

$\widetilde{\sigma}(\nu)=1+\sum_{n=1}^{\infty}\Bigl{(}\frac{\alpha_{s}(Q)}{4\pi}\Bigr{)}^{n}\sum_{k=0}^{2n}\widetilde{\sigma}_{nk}L_{\nu}^{k}\,,$
(61)

where up to $\mathcal{O}(\alpha_{s}^{3})$,

$\displaystyle\widetilde{\sigma}_{12}$
$\displaystyle=K_{12},\qquad\widetilde{\sigma}_{11}=K_{11},\qquad\widetilde{\sigma}_{10}=C_{10},$
(62) $\displaystyle\widetilde{\sigma}_{24}$
$\displaystyle=\frac{K_{12}^{2}}{2},\,\,\,\,\,\widetilde{\sigma}_{23}=K_{11}K_{12}+K_{23},\,\,\,\,\,\widetilde{\sigma}_{22}=\frac{K_{11}^{2}}{2}+C_{10}K_{12}+K_{22},\,\,\,\,\,\widetilde{\sigma}_{21}=C_{21}+C_{10}K_{11}+K_{21},\,\,\,\,\,\widetilde{\sigma}_{20}=C_{20},$
$\displaystyle\widetilde{\sigma}_{36}$
$\displaystyle=\frac{K_{12}^{3}}{6},\qquad\widetilde{\sigma}_{35}=K_{12}\left(K_{23}+\frac{K_{11}K_{12}}{2}\right),\qquad\widetilde{\sigma}_{34}=K_{12}\left(K_{22}+\frac{K_{11}^{2}}{2}+\frac{C_{10}K_{12}}{2}\right)+K_{11}K_{23}+K_{34}\,,$
$\displaystyle\widetilde{\sigma}_{33}$
$\displaystyle=C_{10}(K_{11}K_{12}+K_{23})+C_{21}K_{12}+K_{12}K_{21}+K_{11}K_{22}+K_{33}+\frac{K_{11}^{3}}{6},$
$\displaystyle\widetilde{\sigma}_{32}$
$\displaystyle=C_{32}+C_{21}K_{11}+C_{20}K_{12}+K_{11}K_{21}+C_{10}\left(K_{22}+\frac{K_{11}^{2}}{2}\right)+K_{32}\,,$
$\displaystyle\widetilde{\sigma}_{31}$
$\displaystyle=C_{31}+C_{20}K_{11}+C_{10}K_{21}+K_{31},\qquad\widetilde{\sigma}_{30}=C_{30}\,.$

It is also useful to know how to transform (61) back to momentum space. A
quick way to do this is from (5), which implies that the fixed-order $\tau$
cross section is given simply by

$\frac{1}{\sigma_{0}}{\sigma_{c}(\tau)}=\widetilde{\sigma}(\nu;L_{\nu}\to\partial_{\Omega}+L_{\tau})\frac{e^{\gamma_{E}\Omega}}{\Gamma(1-\Omega)}\Bigr{\rvert}_{\Omega\to
0}\,,$ (63)

where each log $L_{\nu}$ in (61) is replaced by the differential operator
shown, with $L_{\tau}=\ln(1/\tau)$. Then the cross section is expanded in the
form,

$\frac{1}{\sigma_{0}}{\sigma_{c}(\tau)}=1+\sum_{n=1}^{\infty}\Bigl{(}\frac{\alpha_{s}(Q)}{4\pi}\Bigr{)}^{n}\sum_{k=0}^{2n}\sigma_{nk}L_{\tau}^{k}\,,$
(64)

where the coefficients $\sigma_{nk}$ are given in terms of the Laplace-space
coefficients $\widetilde{\sigma}_{nk}$ in (62) by

$\displaystyle\sigma_{12}$
$\displaystyle=\widetilde{\sigma}_{12},\,\,\,\,\,\sigma_{11}=\widetilde{\sigma}_{11},\,\,\,\,\,\sigma_{10}=\widetilde{\sigma}_{10}-\frac{\pi^{2}}{6}\widetilde{\sigma}_{12},$
(65) $\displaystyle\sigma_{24}$
$\displaystyle=\widetilde{\sigma}_{24},\,\,\,\,\,\sigma_{23}=\widetilde{\sigma}_{23},\,\,\,\,\,\,\sigma_{22}=\widetilde{\sigma}_{22}-\pi^{2}\widetilde{\sigma}_{24},\,\,\,\,\,\sigma_{21}=\widetilde{\sigma}_{21}-\frac{\pi^{2}}{2}\widetilde{\sigma}_{23}-8\zeta_{3}\widetilde{\sigma}_{24},\,\,\,\,\,\sigma_{20}=\widetilde{\sigma}_{20}-\frac{\pi^{2}}{6}\widetilde{\sigma}_{22}-2\zeta_{3}\widetilde{\sigma}_{23}+\frac{\pi^{4}}{60}\widetilde{\sigma}_{24}$
$\displaystyle\sigma_{36}$
$\displaystyle=\widetilde{\sigma}_{36},\,\,\,\,\,\sigma_{35}=\widetilde{\sigma}_{35},\,\,\,\,\,\sigma_{34}=\widetilde{\sigma}_{34}-\frac{5\pi^{2}}{2}\widetilde{\sigma}_{36},\,\,\,\,\,\sigma_{33}=\widetilde{\sigma}_{33}-\frac{5\pi^{2}}{3}\widetilde{\sigma}_{35}-40\zeta_{3}\widetilde{\sigma}_{36},\,\,\,$
$\displaystyle\sigma_{32}$
$\displaystyle=\widetilde{\sigma}_{32}-\pi^{2}\widetilde{\sigma}_{34}-20\zeta_{3}\widetilde{\sigma}_{35}+\frac{\pi^{4}}{4}\widetilde{\sigma}_{36},\qquad\sigma_{31}=\widetilde{\sigma}_{31}-\frac{\pi^{2}}{2}\widetilde{\sigma}_{33}-8\zeta_{3}\widetilde{\sigma}_{34}+\frac{\pi^{4}}{12}\widetilde{\sigma}_{35}+\left(20\pi^{2}\zeta_{3}-144\zeta_{5}\right)\widetilde{\sigma}_{36},$
$\displaystyle\sigma_{30}$
$\displaystyle=\widetilde{\sigma}_{30}-\frac{\pi^{2}}{6}\widetilde{\sigma}_{32}-2\zeta_{3}\widetilde{\sigma}_{33}+\frac{\pi^{4}}{60}\widetilde{\sigma}_{34}+\left(\frac{10\pi^{2}\zeta_{3}}{3}-24\zeta_{5}\right)\widetilde{\sigma}_{35}+\left(40\zeta_{3}^{2}-\frac{5\pi^{6}}{168}\right)\widetilde{\sigma}_{36}\,.$

In obtaining these expressions we made use of the following handy identities:

$\displaystyle\mathcal{G}(\Omega)\equiv\frac{e^{\gamma_{E}\Omega}}{\Gamma(1-\Omega)}\,,\qquad\partial_{\Omega}\mathcal{G}\rvert_{\Omega\to
0}$ $\displaystyle=0\,,\qquad\partial_{\Omega}^{2}\mathcal{G}\rvert_{\Omega\to
0}=-\frac{\pi^{2}}{6}\,,\qquad\partial_{\Omega}^{3}\mathcal{G}\rvert_{\Omega\to
0}=-2\zeta_{3}\,,$ (66)
$\displaystyle\partial_{\Omega}^{4}\mathcal{G}\rvert_{\Omega\to 0}$
$\displaystyle=\frac{\pi^{4}}{60}\,,\quad\partial_{\Omega}^{5}\mathcal{G}\rvert_{\Omega\to
0}=\frac{10\pi^{2}}{3}\zeta_{3}-24\zeta_{5}\,,\quad\partial_{\Omega}^{6}\mathcal{G}\rvert_{\Omega\to
0}=40\zeta_{3}^{2}-\frac{5\pi^{6}}{168}\,.$

##  Appendix B Renormalon Cancellation Formulae

In this appendix we list all of the renormalon-cancellation formulae required
for our N3LL${}^{\prime}+\mathcal{O}(\alpha_{s}^{2})$ thrust calculation.
Beginning at the level of the unexpanded, renormalon-corrected cross section
in (22)-(23), and further defining

$f_{\text{mod}}^{(i)}(k-2\Delta)\equiv\left(\frac{\alpha_{s}(\mu_{S})}{4\pi}\right)^{i}\bar{f}_{\text{mod}}^{(i)}(k-2\Delta)\,,$
(67)

the components of (23) can be written explicitly as

$\displaystyle\bar{f}_{\text{mod}}^{(0)}(k-2\Delta)$
$\displaystyle=f_{\text{mod}}(k-2\Delta)\,,$ (68a)
$\displaystyle\bar{f}_{\text{mod}}^{(1)}(k-2\Delta)$
$\displaystyle=-2\,\delta^{1}\left(\frac{\mu_{R}}{\xi}\right)f^{\prime}_{\text{mod}}(k-2\Delta)\,,$
(68b) $\displaystyle\bar{f}_{\text{mod}}^{(2)}(k-2\Delta)$
$\displaystyle=\Bigl{[}-2\,\delta^{2}\left(\frac{\mu_{R}}{\xi}\right)f^{\prime}_{\text{mod}}(k-2\Delta)+2\left(\delta^{1}\,\frac{\mu_{R}}{\xi}\right)^{2}\,f^{\prime\prime}_{\text{mod}}(k-2\Delta)\Bigr{]}\,,$
(68c) $\displaystyle\bar{f}_{\text{mod}}^{(3)}(k-2\Delta)$
$\displaystyle=\left[-2\,\delta^{3}\left(\frac{\mu_{R}}{\xi}\right)f^{\prime}_{\text{mod}}(k-2\Delta)+4\,\delta^{1}\,\delta^{2}\,\left(\frac{\mu_{R}}{\xi}\right)^{2}f^{\prime\prime}_{\text{mod}}(k-2\Delta)-\frac{4}{3}\left(\delta^{1}\,\frac{\mu_{R}}{\xi}\right)^{3}f^{\prime\prime\prime}_{\text{mod}}(k-2\Delta)\right]\,,$
(68d)

where the subtraction terms $\delta^{n}$ and $\Delta$ have an implicit
functional dependence on the soft and renormalon scales,
$\\{\delta^{n},\Delta\\}\equiv\\{\delta^{n}(\mu_{S},\mu_{R}),\Delta(\mu_{S},\mu_{R})\\}$.
Indeed, the form of $\Delta(\mu_{\delta},\mu_{R})$ for arbitrary reference
$\mu_{\delta}$ and subtraction $\mu_{R}$ scales (i.e. accounting for RGE in
both scales, according to (III.2)) is known to the required three-loop order,
and is given by

$\displaystyle\Delta(\mu_{\delta},\mu_{R})$
$\displaystyle=\Delta(R_{\Delta},R_{\Delta})+2\left(\frac{\mu_{R}}{\xi}\right)\eta_{\Gamma}(\mu_{\delta},\mu_{R})$
$\displaystyle+\frac{R_{\Delta}}{2\beta_{0}}e^{-G\left[\alpha_{s}(R_{\Delta})\right]}\left(\frac{2\pi}{\beta_{0}}e^{i\pi}\right)^{\frac{\beta_{1}}{2\beta_{0}^{2}}}\bigg{\\{}\gamma^{0}_{R}\left[\Gamma\left(-\frac{\beta_{1}}{2\beta_{0}^{2}},-\frac{2\pi}{\beta_{0}\alpha_{s}(\mu_{R})}\right)-\Gamma\left(-\frac{\beta_{1}}{2\beta_{0}^{2}},-\frac{2\pi}{\beta_{0}\alpha_{s}(R_{\Delta})}\right)\right]$
$\displaystyle-\frac{1}{2\beta_{0}}\left[\gamma_{R}^{1}-\frac{\gamma_{R}^{0}}{\beta_{0}}\left(\beta_{1}+\frac{B_{2}}{2}\right)\right]\left[\Gamma\left(-\frac{\beta_{1}}{2\beta_{0}^{2}}-1,-\frac{2\pi}{\beta_{0}\alpha_{s}(\mu_{R})}\right)-\Gamma\left(-\frac{\beta_{1}}{2\beta_{0}^{2}}-1,-\frac{2\pi}{\beta_{0}\alpha_{s}(R_{\Delta})}\right)\right]$
$\displaystyle+\frac{1}{4\beta_{0}^{2}}\left[\gamma_{R}^{2}-\frac{\gamma_{R}^{1}}{\beta_{0}}\left(\beta_{1}+\frac{B_{2}}{2}\right)+\gamma_{R}^{0}\left(B_{2}+\frac{B_{2}\beta_{1}}{2\beta_{0}^{2}}-\frac{B_{3}}{4\beta_{0}}+\frac{B_{2}^{2}}{8\beta_{0}^{2}}\right)\right]$
$\displaystyle\times\left[\Gamma\left(-\frac{\beta_{1}}{2\beta_{0}^{2}}-2,-\frac{2\pi}{\beta_{0}\alpha_{s}(\mu_{R})}\right)-\Gamma\left(-\frac{\beta_{1}}{2\beta_{0}^{2}}-2,-\frac{2\pi}{\beta_{0}\alpha_{s}(R_{\Delta})}\right)\right]\bigg{\\}},$
(69)

where we take the input gap parameter to be
$\Delta(R_{\Delta},R_{\Delta})=0.1$ GeV at an arbitrary reference scale
$R_{\Delta}=1.5$ GeV (the exact value for this reference scale is not expected
to be particularly consequential in the tail regions relevant to our fits
Abbate  _et al._ (2011).) and where $G$ is the anti-derivative of
$1/\beta\left[\alpha\right]$,

$G\left[\alpha\right]=\frac{2\pi}{\beta_{0}}\left[\frac{1}{\alpha}+\frac{\beta_{1}}{4\pi\beta_{0}}\ln\alpha-\frac{B_{2}}{(4\pi)^{2}}\alpha+\frac{B_{3}}{(4\pi)^{3}}\frac{\alpha^{2}}{2}\right]\,.$
(70)

Here $B_{2,3}$ are given as

$\displaystyle B_{2}$
$\displaystyle\equiv-\frac{\beta_{2}}{\beta_{0}}+\frac{\beta_{1}^{2}}{\beta_{0}^{2}}\,,\,\,\,\,\,\,\,\,\,\,B_{3}\equiv-\frac{\beta_{3}}{\beta_{0}}+2\frac{\beta_{1}\beta_{2}}{\beta_{0}^{2}}-\frac{\beta_{1}^{3}}{\beta_{0}^{3}}\,,$
(71)

in terms of the QCD $\beta$-function coefficients defined through

$\mu\frac{d\alpha_{s}(\mu)}{d\mu}=-2\beta_{0}\frac{\alpha_{s}(\mu)^{2}}{4\pi}\left[1+\left(\frac{\alpha_{s}(\mu)}{4\pi}\right)\frac{\beta_{1}}{\beta_{0}}+\left(\frac{\alpha_{s}(\mu)}{4\pi}\right)^{2}\frac{\beta_{2}}{\beta_{0}}+\left(\frac{\alpha_{s}(\mu)}{4\pi}\right)^{3}\frac{\beta_{3}}{\beta_{0}}\right]\,.$
(72)

The solution to (72) for the running coupling up to four loops is

$\displaystyle\frac{1}{\alpha_{s}(\mu)}$
$\displaystyle=\frac{1}{\alpha_{s}(m_{Z})}\biggl{\\{}X+\frac{\alpha_{s}(m_{Z})}{4\pi}\frac{\beta_{1}}{\beta_{0}}\ln
X+\biggl{(}\frac{\alpha_{s}(m_{Z})}{4\pi}\biggr{)}^{2}\biggl{[}\frac{\beta_{2}}{\beta_{0}}\Bigl{(}1-\frac{1}{X}\Bigr{)}+\frac{\beta_{1}^{2}}{\beta_{0}^{2}}\Bigl{(}\frac{\ln
X}{X}+\frac{1}{X}-1\Bigr{)}\biggr{]}$ (73)
$\displaystyle\qquad\qquad+\biggl{(}\frac{\alpha_{s}(m_{Z})}{4\pi}\biggr{)}^{3}\frac{1}{X^{2}}\biggl{[}\frac{\beta_{3}}{2\beta_{0}}(X^{2}-1)+\frac{\beta_{1}\beta_{2}}{\beta_{0}^{2}}(X+\ln
X-X^{2})+\frac{\beta_{1}^{3}}{2\beta_{0}^{3}}[(1-X)^{2}-\ln^{2}X]\biggr{]}\biggr{\\}}\,,$

where

$X\equiv 1+\frac{\alpha_{s}(m_{Z})}{2\pi}\beta_{0}\ln\frac{\mu}{m_{Z}}\,.$
(74)

The actual functional dependence of the subtraction terms $\delta^{n}$ in (68)
on the reference and subtraction scales is scheme-dependent, as is the
structure of the $\mu_{R}$-anomalous dimensions appearing in (69). We now
present the ingredients $\delta^{i}$ and $\gamma_{R}^{i}$ in both the $n=0$
and $n=1$ schemes considered in the main text.

###  B.1 $n=0$ Schemes

From (26) and the fixed-order expansion of the soft function that can be
constructed from (50), one reads off that

$\delta(\mu_{\delta},\mu_{R})=\frac{\mu_{R}}{2\xi}\ln\tilde{S}\Bigl{(}v=\frac{\xi}{\mu_{R}},\mu_{\delta}\Bigr{)}\equiv\frac{\mu_{R}}{2\xi}\sum_{n=1}\left(\frac{\alpha_{s}(\mu_{\delta})}{4\pi}\right)^{n}\,\delta^{n}(\mu_{\delta},\mu_{R})$
(75)

for a generic $\xi$ scheme and reference scale $\mu_{\delta}$. Given the known
three-loop expansion of the Laplace-space soft function $\tilde{S}$, one can
then easily obtain the order-by-order expressions for $\delta$ as follows:

$\displaystyle\delta^{1}(\mu_{\delta},\mu_{R})$
$\displaystyle=\Gamma_{s}^{0}\,L^{2}+c_{\tilde{S}}^{1}\,,$ (76a)
$\displaystyle\delta^{2}(\mu_{\delta},\mu_{R})$
$\displaystyle=\frac{2}{3}\Gamma_{s}^{0}\beta_{0}\,L^{3}+\Gamma_{s}^{1}L^{2}+\left(\gamma_{s}^{1}+2c_{\tilde{S}}^{1}\beta_{0}\right)L+c_{\tilde{S}}^{2}-\frac{1}{2}(c_{\tilde{S}}^{1})^{2}\,,$
(76b) $\displaystyle\delta^{3}(\mu_{\delta},\mu_{R})$
$\displaystyle=\frac{2}{3}\Gamma_{s}^{0}\beta_{0}^{2}\,L^{4}+\frac{2}{3}\left(2\Gamma_{s}^{1}\beta_{0}+\Gamma_{s}^{0}\beta_{1}\right)\,L^{3}+\left(\Gamma_{s}^{2}+2\gamma_{s}^{1}\beta_{0}+4c_{\tilde{S}}^{1}\beta_{0}^{2}\right)\,L^{2}$
(76c)
$\displaystyle\quad+\left(\gamma_{s}^{2}+2c_{\tilde{S}}^{1}\beta_{1}+4c_{\tilde{S}}^{2}\beta_{0}-2(c_{\tilde{S}}^{1})^{2}\beta_{0}\right)\,L-c_{\tilde{S}}^{1}c_{\tilde{S}}^{2}+\frac{1}{3}(c_{\tilde{S}}^{1})^{3}+c_{\tilde{S}}^{3}\,,$

with $L=\ln\frac{\mu_{\delta}e^{\gamma_{E}}\xi}{\mu_{R}}$. Then, from (III.2)
and (72), one immediately obtains the relevant $\mu_{R}$-anomalous dimensions
in this scheme:

$\displaystyle\gamma_{R}^{0}$
$\displaystyle=\frac{1}{2\xi}\left(\Gamma_{s}^{0}\,L^{2}+c_{\tilde{S}}^{1}\right)\,,$
(77a) $\displaystyle\gamma_{R}^{1}$
$\displaystyle=\frac{1}{2\xi}\left(\frac{2}{3}\Gamma_{s}^{0}\beta_{0}\,L^{3}+\Gamma_{s}^{1}\,L^{2}+\left(\gamma_{s}^{1}+2c_{\tilde{S}}^{1}\beta_{0}\right)L+c_{\tilde{S}}^{2}-\frac{1}{2}(c_{\tilde{S}}^{1})^{2}-2\beta_{0}\left(\Gamma_{s}^{0}\,L^{2}+c_{\tilde{S}}^{1}\right)\right)\,,$
(77b) $\displaystyle 2\xi\,\gamma_{R}^{2}$
$\displaystyle=\frac{2}{3}\Gamma_{s}^{0}\beta_{0}^{2}\,L^{4}+\frac{2}{3}\left(2\Gamma_{s}^{1}\beta_{0}+\Gamma_{s}^{0}\beta_{1}\right)L^{3}+\left(\Gamma_{s}^{2}+2\gamma_{s}^{1}\beta_{0}+4c_{\tilde{S}}^{1}\beta_{0}^{2}\right)L^{2}+\left(\gamma_{s}^{2}+2c_{\tilde{S}}^{1}\beta_{1}+4c_{\tilde{S}}^{2}\beta_{0}-2(c_{\tilde{S}}^{1})^{2}\beta_{0}\right)L$
$\displaystyle-
c_{\tilde{S}}^{2}c_{\tilde{S}}^{1}+\frac{1}{3}(c_{\tilde{S}}^{1})^{3}+c_{\tilde{S}}^{3}$
$\displaystyle-2\left(\frac{4}{3}\Gamma_{s}^{0}\beta_{0}^{2}\,L^{3}+\left(\Gamma_{s}^{0}\beta_{1}+2\Gamma_{s}^{1}\beta_{0}\right)\,L^{2}+\beta_{0}\left(4c_{\tilde{S}}^{1}\beta_{0}+2\gamma_{s}^{1}\right)\,L+2c_{\tilde{S}}^{2}\beta_{0}+c_{\tilde{S}}^{1}\beta_{1}-(c_{\tilde{S}}^{1})^{2}\beta_{0}\right)\,,$
(77c)

where now $L=\ln\xi e^{\gamma_{E}}$ ($\mu_{\delta}=\mu_{R}$ in this
calculation, see (III.2)). This log of course vanishes in schemes where
$\xi=e^{-\gamma_{E}}$.

In practice, the subtraction series (75) appears in (20) as part of a