aainstitutetext: Institut für Theoretische Teilchenphysik und Kosmologie, RWTH
Aachen University, D-52056 Aachen, Germanybbinstitutetext: Cavendish
Laboratory, University of Cambridge, Cambridge CB3 0HE, United
Kingdom111Preprint: Cavendish-HEP-22/08, P3H-22-100, TTK-22-31

# NNLO B-fragmentation fits and their application to $t\bar{t}$ production and
decay at the LHC

Michał Czakon a Terry Generet b Alexander Mitov b Rene Poncelet
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this work we derive three sets of non-perturbative fragmentation functions,
with uncertainties, for $B$-hadrons, $J/\psi$’s and muons resulting from
semileptonic $B$ decays. All three sets are with next-to-next-to leading order
accuracy and include next-to-next-to leading logarithmic soft gluon
resummation. The novel feature of these new sets is that they are fully
consistent with our formalism for next-to-next-to leading order (NNLO)
calculations for final states with identified $B$, $J/\psi$ or a $\mu$. We
employ the fragmentation functions derived in this work to make state of the
art predictions for such final states in $t\bar{t}$ events at the LHC. A
special emphasis is placed on observables sensitive to the top quark mass. The
present work opens the door for many LHC applications, like, open $B$
production or $B$ production in association with bosons.

##  1 Introduction

Studies of identified hadrons containing a heavy flavor (i.e. a valence $b$ or
$c$ quark) have a long history. Indeed, calculations with next to leading
order (NLO) QCD accuracy, including resummation of collinear and soft
logarithms at next-to-leading logarithmic accuracy, have been compared to data
from lepton and/or hadron colliders for almost three decades Mele:1990cw ;
Binnewies:1997xq ; Cacciari:1998it ; Cacciari:2001td ; Cacciari:2002pa ;
Cacciari:2003uh ; Cacciari:2005rk ; Cacciari:2005ry ; Kniehl:2005de ;
Cacciari:2005uk ; Kniehl:2006mw ; Kneesch:2007ey ; Cacciari:2012ny ;
Cacciari:2015fta ; Kniehl:2020szu . Increasing the accuracy of theoretical
predictions beyond NLO in QCD has been a longstanding goal, see for example
ref. Cacciari:2002pa . A first step in this direction was the calculation of
the massless coefficient functions for ${e^{+}e^{-}}$ colliders Rijken:1996vr
; Rijken:1996npa ; Rijken:1996ns ; Mitov:2006wy . The resummation of soft and
collinear logarithms has been possible for at least a decade with the
derivation of certain process-independent ingredients like the heavy quark
perturbative fragmentation function (PFF) Melnikov:2004bm ; Mitov:2004du , the
time-like splitting functions Mitov:2006ic ; Moch:2007tx ; Almasy:2011eq  and
soft anomalous dimensions Cacciari:2001cw ; Cacciari:2002xb ; Gardi:2005yi .
Public codes implementing these results exist Bertone:2013vaa ;
Bertone:2015cwa . Combined, these results have made it possible to extract the
fragmentation function (FF) for heavy quarks from ${e^{+}e^{-}}$ collisions
with next-to-next-to leading order (NNLO) accuracy and with resummed soft and
collinear logs to next-to-next-to leading logarithmic (NNLL) accuracy
Salajegheh:2019ach ; Fickinger:2016rfd ; Ridolfi:2019bch . Ideas about heavy
flavor fragmentation that do not involve data fits have also been put forward
Braaten:1994bz ; Aglietti:2006yf . Heavy flavor production in a jet has also
received a lot of attention Anderle:2017cgl ; Dai:2018ywt ; Dai:2021mxb ;
Caletti:2022glq ; Caletti:2022hnc .

Progress towards applying these results at other processes (like the LHC) have
been hampered by the lack of process-dependent coefficient functions with NNLO
accuracy. First step in this direction was undertaken only very recently
Czakon:2021ohs  where a numerical method for the calculation of the
coefficient functions at NNLO in any process was developed and applied to
$b$-fragmentation in ${t\bar{t}}$ production and decay at the LHC.

Interesting parallels can be made with the field of light-flavor fragmentation
which historically has been developed largely independently from that of heavy
flavor fragmentation, despite the many similarities, common features and
ingredients between the two. A number of sets of light flavor fragmentation
functions exist, mostly at NLO but also more recently at NNLO QCD
Binnewies:1994ju ; Kniehl:2000cr ; Kniehl:2000fe ; Kniehl:2000hk ;
Albino:2005me ; deFlorian:2007aj ; deFlorian:2014xna ; Anderle:2015lqa ;
Sato:2016wqj ; Bertone:2017tyb ; Moffat:2021dji ; Khalek:2021gxf ;
Abdolmaleki:2021yjf ; Borsa:2021ran ; Borsa:2022vvp ; Khalek:2022vgy . Just
like in the case of heavy flavor fragmentation, progress towards their
application to other processes has been rather restricted due to the lack of
NNLO-accurate process-dependent coefficient functions. The development of ref.
Czakon:2021ohs  should allow for the calculation of coefficient functions with
full flexibility, broadening the past practice of exclusively using analytical
calculations which, despite their many useful features, are restricted to a
small class of very inclusive observables.

The calculation of ref. Czakon:2021ohs  has utilized non-perturbative
fragmentation functions (NPFF) extracted previously. While this made it
possible to present predictions for specific realistic observables at the LHC,
the predictions were not ideal since, as explained there, the NPFFs used in
ref. Czakon:2021ohs  were not fully consistent with the calculation of the
perturbative part. Although such a mismatch is numerically small, a consistent
program for NNLO predictions for processes with identified $B$ hadrons
requires the use of a consistently extracted set of NPFFs.

The main goal of this article is to extract such a set from the available
${e^{+}e^{-}}$ data. As a first application, and building upon the large body
of literature Corcella:2001hz ; Cacciari:2002re ; Corcella:2005dk ;
Biswas:2010sa ; Kniehl:2012mn ; Agashe:2012bn ; MoosaviNejad:2013xqs ;
MoosaviNejad:2016gcd ; Agashe:2016bok ; Brancaccio:2021gcz , in this work we
also provide NNLO QCD accurate predictions for a set of realistic
$b$-fragmentation observables in ${t\bar{t}}$ events at the LHC. A particular
emphasis is placed on observables sensitive to the value of the top quark
mass. Indeed, the current work allows for the first time to perform
extractions of $m_{t}$ from this class of observables with full NNLO QCD
accuracy. A second novel feature of this work is that we consider observables
constructed not only at the level of $b$-flavored hadrons but also for their
decay products: $J/\psi$ and muons from semileptonic $B$ decays.

This work is organized as follows: in sec. 2 we present the main result of
this work, namely, the extraction of NNLO-accurate NPFFs for $B$-hadrons and
$J/\psi$ and soft muons resulting from $B$-decays. Theoretical aspects of this
extraction related to the treatment of coefficient functions, perturbative
fragmentation functions and soft-gluon resummation are given in, respectively,
sections 2.1, 2.2 and 2.3; the global fit of NPFF for $B$-hadrons from
$e^{+}e^{-}$ data is presented in sec. 2.4 while sec. 2.5 presents the
derivation of NNLO-accurate NPFFs for $J/\psi$ and $\mu$. Sec. 3 is devoted to
phenomenological applications of the newly derived NPFFs, in particular, for
state-of-the-art predictions for precision $m_{t}$ determination at the LHC.
Some lengthier expressions related to soft-gluon resummation are collected in
appendix A while appendix B presents extraction of $B$ NPFF from one dataset
at a time. Our conclusions are summarized in sec. 4.

##  2 Extraction of NNLO $B$-fragmentation functions from $e^{+}e^{-}$ data

The main source of data on $b$-fragmentation at present are measurements of
$b$-flavored hadrons in ${e^{+}e^{-}}$ collisions at $E_{\rm c.m.}=m_{Z}$.
What is measured is the energy distribution of $b$-flavored hadrons $B$.
Throughout this work we denote this energy through the dimensionless variable
$z=E/E_{max}$. The relevant cross-section can be written as

$\frac{d\sigma_{B}(m_{b},z)}{dz}=\sum_{i}\Bigg{\\{}\frac{d\sigma_{i}(\mu_{Fr},z)}{dz}\otimes
D_{i\to B}(\mu_{Fr},m_{b},z)\Bigg{\\}}(z)+{\cal O}(m_{b}^{2})\,,$ (1)

where $\otimes$ denotes the usual integral convolution, $\mu_{Fr}$ is the
factorization scale associated with the fragmentation function, $m_{b}$ is the
pole mass of the $b$ quark and the index $i$ runs over all partons that can
contribute. Typically, one includes the heavy flavor ($b$ and $\bar{b}$ in
this case), all lighter flavors and antiflavors as well as the gluon. Heavier
quarks (top in this case) are not included. Further details can be found in
ref. Czakon:2021ohs .

The process-independent set of fragmentation functions $D_{i\to B}$ can be
further factorized into perturbative and non-perturbative components

$D_{i\to B}(\mu_{Fr},m_{b},z)=\Big{\\{}D_{i\to
b}(\mu_{Fr},\mu_{Fr,0},m_{b},z)\otimes{D^{\rm NP}_{b\to
B}}(\mu_{Fr,0},z)\Big{\\}}(z)\,.$ (2)

The scale $\mu_{Fr,0}$ appearing in the above equation sets the scale at which
the non-perturbative fragmentation process $b\to B$ takes place. It is low, of
the order of the mass $m_{b}$ of the heavy quark. The standard choice for this
scale, which we also adopt in this work, is $\mu_{Fr,0}=m_{b}$. On the other
hand, the scale $\mu_{Fr}$ is large, of the order of the hard scale in the
process i.e. one has $\mu_{Fr}\sim m_{Z}\gg\mu_{Fr,0}$. This hierarchy of
scales introduces large collinear logarithms to all orders of the strong
coupling constant. These logs are resummed with the help of the DGLAP
evolution equation which the functions $D_{i\to B}$ satisfy.

One of the main goals of this work is to provide a fit from ${e^{+}e^{-}}$
data of the non-perturbative fragmentation function ${D^{\rm NP}_{b\to B}}$
for the fragmentation of a $b$-quark into a $B$-hadron. The function ${D^{\rm
NP}_{b\to B}}$ is fully defined only once the perturbative cross-section
$d\sigma_{i}/dz$ and the perturbative fragmentation function $D_{i\to b}$ have
been fully specified, to which we turn next.

###  2.1 Treatment of the coefficient functions

In this work we have utilized the analytically known through NNLO massless
partonic cross-sections $d\sigma_{i}/dz$ Rijken:1996vr ; Rijken:1996npa ;
Rijken:1996ns ; Mitov:2006wy . As indicated in eq. (1), the use of massless
coefficient functions leads to missing power corrections ${\cal O}(m_{b}^{2})$
which are negligible for all kinematics except for very small $z$. On the
other hand, the mass-independent term and associated logarithms
$\ln^{n}(m_{b}^{2}/m_{Z}^{2})$ are correctly accounted for by eq. (1). We have
resummed soft-gluon logarithms in the coefficient functions through NNLL
following the approach of ref. Cacciari:2001cw . The resummation follows the
standard Mellin space approach with $N$ being the Mellin variable conjugated
to $z$. In the soft limit $z\to 1$ (i.e. $N\to\infty$) the cross-section
further factorizes as

$\sigma_{b}(N)=H(\alpha_{S})\exp\left[G(N,\alpha_{S})\right]+{\cal O}(1/N)\,.$
(3)

The explicit expressions for the function $G$ (through NNLL) and for the hard
function $H$ (through NNLO) can be found in appendix A.2. The complete
resummed expression requires the matching of the above result with the exact
fixed order result. Usually this is done through the so-called additive
matching procedure where the resummed result takes the form

$\sigma_{b}(N)=H(\alpha_{S})\exp\left[G(N,\alpha_{S})\right]+M(N,\alpha_{S})\,,$
(4)

where the function $M$ contains only ${\cal O}(1/N)$ terms and is defined in
such a way that the expression in eq. (4), when expanded through NNLO,
coincides with the Mellin transform of the coefficient function
$d\sigma_{b}/dz$.

In this work we use a different, multiplicative, version of the above matching
procedure, where the resummed and matched cross-section reads

$\sigma_{b}(N)=\left(H(\alpha_{S})+\tilde{M}(N,\alpha_{S})\right)\exp\left[G(N,\alpha_{S})\right]\,.$
(5)

The function $\tilde{M}$ above is defined from the same requirement as the
function $M$ appearing in the additive approach. The multiplicative and
additive versions differ only beyond NNLO. The main reason we choose to work
with the multiplicatively matched resummed version is that for $z\to 1$ the
resummed cross-section does not become negative, see sec. 2.3 for details.

Finally, any $N$-space result needs to be inverted back to $z$-space. To this
end we utilize the Minimal Prescription of ref. Catani:1996yz .

###  2.2 Treatment of the perturbative fragmentation function

The perturbative fragmentation function $D_{i\to b}$ can be further factorized
as follows

$D_{i\to b}(\mu_{Fr},\mu_{Fr,0},m,z)=\Bigg{\\{}\sum_{j}E_{i\to
j}(\mu_{Fr},\mu_{Fr,0},z)\otimes D^{\rm ini}_{j\to
b}(\mu_{Fr,0},m,z)\Bigg{\\}}(z)\,,$ (6)

where the index $j$ runs over all partons that contribute above the scale
$\mu_{Fr,0}$.

The evolution factor $E_{i\to j}$ is standard; it satisfies the DGLAP
evolution equation with a trivial boundary condition. In particular, it does
not require soft gluon resummation. In our implementation we take this factor
from the public code APFEL Bertone:2013vaa  where timelike evolution is
implemented through NNLO. We do not restrict our implementation to the flavor
non-singlet combination $(b-\bar{b})$ but include all partonic contributions.

The non-trivial information about the PFF at the scale $\mu_{Fr,0}\sim m_{b}$
is contained within the initial condition $D^{\rm ini}_{j\to b}$ Mele:1990cw .
In the case $j=b$ (i.e. the transition $b\to b$), the function $D^{\rm
ini}_{b\to b}$ contains terms that are enhanced in the soft-limit and thus
require soft-gluon resummation. This has been first discussed in ref.
Cacciari:2001cw  through NLL and extended to NNLL in ref. Gardi:2005yi .
Similarly to the case of the coefficient function discussed above, we
implement the soft-gluon resummation of the initial condition through a
multiplicative matching

$D^{\rm ini}_{b\to b}(m_{b},N)=\left(H^{\rm
ini}(m_{b},\alpha_{S})+\tilde{M}^{\rm
ini}(m_{b},N,\alpha_{S})\right)\exp\left[G^{\rm
ini}(m_{b},N,\alpha_{S})\right]\,.$ (7)

For brevity, we have not shown the dependence on $\mu_{Fr,0}$ in the above
equation. Further discussion of the soft gluon resummation is given in sec.
2.3. The explicit expressions for the functions $G^{\rm ini}$ and $H^{\rm
ini}$, through NNLL and NNLO respectively, can be found in appendix A.1.

###  2.3 Implementation of the soft gluon resummation

There are subtleties in the implementation of the soft gluon resummation of
$D^{\rm ini}$ related to the decoupling of the heavy flavor $b$, which require
a more detailed discussion. Let us start by recalling the “master” resummation
formula for $D^{\rm ini}$ Cacciari:2001cw

$\displaystyle D^{\rm ini}_{b\to b}(\mu_{Fr,0},m_{b},N)$ $\displaystyle\sim$
$\displaystyle\exp\int_{0}^{1}dz\frac{z^{N-1}-1}{1-z}\left(\int_{(1-z)^{2}m_{b}^{2}}^{\mu_{Fr,0}^{2}}\frac{dk^{2}}{k^{2}}A\left(\alpha_{S}(k^{2})\right)+D\left(\alpha_{S}((1-z)^{2}m_{b}^{2})\right)\right)$
(8) $\displaystyle+{\cal O}(1/N)\,.$

The function $D$ is often denoted as $H$ in the literature; we refrain from
using this notation here since it can lead to confusion with the generic
notation $H$ for the hard matching function.

Recalling that $\mu_{Fr,0}\geq m_{b}$, the scale $k^{2}$ of the running
coupling in the function $A$ above takes values both below and above the mass
$m_{b}$ of the heavy quark $b$. The scale $(1-z)^{2}m_{b}^{2}$ that enters the
wide-angle emission function $D$ never exceeds $m_{b}^{2}$. With this in mind,
and following eq. (7), we write the above master expression in the following
way

$\displaystyle D^{\rm ini}_{b\to b}(\mu_{Fr,0},m_{b},N)$ $\displaystyle=$
$\displaystyle H^{\rm ini}\exp\left[G^{\rm ini}\right]C_{\rm dec}(\mu_{\rm
th}=m_{b},m_{b})+{\cal O}(1/N)\,.$ (9)

The LHS of the above equation is evaluated at a scale $\mu_{Fr,0}$ which
implies that it is considered in a theory with $n_{f}=n_{l}+1$ active flavors
where the heavy flavor is considered active. In contrast, some of the terms in
the RHS of eq. (9) are defined in a theory with $n_{f}=n_{l}$ active flavors.
Therefore, care must be taken that the final expression for the RHS of eq. (9)
as a whole is converted to a theory with $n_{f}=n_{l}+1$ active flavors. In
the following we explain how this is done as well as the meaning of each term.

The hard matching function contains terms that behave as constant in the limit
$N\to\infty$. It has the following expansion

$\displaystyle H^{\rm ini}$ $\displaystyle=$
$\displaystyle\sum_{i=0}\left(\frac{\alpha_{S}(\mu_{R,0})}{2\pi}\right)^{i}H_{i}^{\rm
ini}\,.$ (10)

The coefficients $H_{i}^{\rm ini},i=0,1,2$ can be found in appendix A.1.

The Sudakov exponent $G^{\rm ini}$ can be written as

$\displaystyle G^{\rm ini}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}dz\frac{z^{N-1}-1}{1-z}\int_{m_{b}^{2}}^{\mu_{Fr,0}^{2}}\frac{dk^{2}}{k^{2}}A^{\prime}\left(\alpha_{S}(k^{2})\right)$
$\displaystyle+\int_{0}^{1}dz\frac{z^{N-1}-1}{1-z}\left(\int_{(1-z)^{2}m_{b}^{2}}^{m_{b}^{2}}\frac{dk^{2}}{k^{2}}A\left(\alpha_{S}^{(n_{l})}(k^{2})\right)+D\left(\alpha_{S}^{(n_{l})}((1-z)^{2}m_{b}^{2})\right)\right)\,,$

with anomalous dimensions $A$ and $D$ defined through the following expansions

$\displaystyle A^{(^{\prime})}(a)$ $\displaystyle=$
$\displaystyle\sum_{i=1}\left(\frac{a}{\pi}\right)^{i}A^{(^{\prime})}_{i}\,,$
(12) $\displaystyle D(a)$ $\displaystyle=$
$\displaystyle\sum_{i=1}\left(\frac{a}{\pi}\right)^{i}D_{i}\,.$ (13)

Throughout this work, all quantities denoted with a prime are defined in a
theory with $n_{f}=n_{l}+1$ active flavors, while quantities without primes
are in a theory with $n_{f}=n_{l}$ active flavors. Similarly,
$\alpha_{S}^{(n_{l})}$ is the coupling in a theory with $n_{f}=n_{l}$ active
flavors while $\alpha_{S}$ is the one in a theory with $n_{f}=n_{l}+1$ active
flavors. The explicit expressions can be found in appendix A.1.

The term in the first line of eq. (2.3) is defined above the flavor threshold
$\mu_{\rm th}=m_{b}$ i.e. in a theory with $n_{f}=n_{l}+1$ active flavors. In
contrast, the term in the second line of that equation contains scales below
$\mu_{\rm th}=m_{b}$ and is therefore defined in a theory with $n_{f}=n_{l}$
active flavors. Upon its evaluation (discussed below) it needs to be converted
to a theory with $n_{l}+1$ active flavors. This is achieved with the help of
the decoupling relations for the strong coupling constant Larin:1994va ;
Chetyrkin:1997un ; Schroder:2005hy ; Chetyrkin:2005ia  and the decoupling
constant $C_{\rm dec}(\mu_{\rm th},m_{b})$ relating the fragmentation function
$D^{\rm ini}_{b\to b}$ between theories with $n_{l}$ and $n_{l}+1$ active
flavors. Schematically, the latter relation takes the following form

$D^{{\rm ini}\,(n_{f}=n_{l}+1)}_{b\to b}(\mu_{Fr,0}=\mu_{\rm th},m_{b})=C_{\rm
dec}(\mu_{\rm th},m_{b})\,D^{{\rm ini}\,(n_{f}=n_{l})}_{b\to
b}(\mu_{Fr,0}=\mu_{\rm th},m_{b})\,.$ (14)

The decoupling coefficient $C_{\rm dec}$ depends on the scales $\mu_{\rm th}$
and $m_{b}$. Since, as indicated in eq. (9), we have chosen the threshold to
be at $\mu_{\rm th}=m_{b}$ in this work we only need $C_{\rm dec}$ for this
value of $\mu_{\rm th}$. The matching coefficient is known through NNLO
Neubert:2007je . In this work we are only interested in its $z\to 1$ limit

$\displaystyle C_{\rm dec}(\mu_{\rm th},m_{b})$ $\displaystyle=$
$\displaystyle
1+\bigg{(}\frac{\alpha_{S}(\mu_{R,0})}{2\pi}\bigg{)}^{2}C_{F}T_{F}$ (15)
$\displaystyle\times\ln(N)\bigg{[}-\frac{2}{3}\ln^{2}\bigg{(}\frac{m_{b}^{2}}{\mu_{\text{th}}^{2}}\bigg{)}-\frac{20}{9}\ln\bigg(\frac{m_{b}^{2}}{\mu_{\text{th}}^{2}}\bigg{missing})-\frac{56}{27}\bigg{]}+{\cal
O}(N^{0})\,.$

We will not be concerned with the explicit expression for the $N$-independent
constant term above since we have absorbed it into the hard matching
coefficient $H_{2}^{\rm ini}$, see eq. (9).

The agreement between the NNLO expansion of the resummation formula eq. (9)
and the direct fixed order calculation of $D^{\rm ini}_{b\to b}$ represents a
non-trivial consistency check between the soft-gluon resummation formalism,
the anomalous dimension $D_{2}$ (which is specific to $D^{\rm ini}$) computed
in ref. Gardi:2005yi  with the help of two independent approaches and the
decoupling coefficient $C_{\rm dec}$ derived in ref. Neubert:2007je . On the
other hand, given $C_{\rm dec}$ has been extracted from $D^{\rm ini}_{b\to b}$
itself, one may wonder how strong this consistency check is. To that end, in
the following we invoke another independent argument confirming the result
(15) for $C_{\rm dec}$.

Since we are only interested in the strict soft limit $z\to 1$, one may
exploit the relationship between parton distribution functions (PDF) and FF in
this limit Gardi:2005yi  to confirm the decoupling relation (15) from the
known decoupling relations for PDF Buza:1996wv . One difference between the
case of heavy flavor FF we study in this work and PDF for heavy flavor
production is that the literature on PDF almost exclusively deals with the
case when “intrinsic” heavy flavor is not present (“intrinsic” heavy flavor
would be the direct analogue of heavy flavor FF). The results of ref.
Buza:1996wv  are derived in this context.

The results of ref. Buza:1996wv  are sufficient to address the question of
interest here. Recall that we are interested in the double limit $z\to 1$
(which corresponds to the limit $x\to 1$ for PDF) and $m_{b}\to 0$ where all
power corrections ${\cal O}(m_{b}^{n})$ are neglected. In this double limit
one expects that the “valence” heavy and light quarks behave the same way,
i.e. any (cut) loop corrections involving heavy quarks affect them equally.
From this one may conclude that the decoupling constant relevant for the heavy
quark is the same as the one for the light quarks. From the results of ref.
Buza:1996wv  one can see that it is the operator matrix element $A^{{\rm
NS},(2)}_{qq,H}$ which describes the decoupling of both the non-singlet and
the singlet quark densities in the soft limit. One can easily see than in the
limit $x\to 1$ the terms $\sim 1/(1-x)_{+}$ in $A^{{\rm NS},(2)}_{qq,H}$ are
identical to $C_{\rm dec}$ in (15). In fact, at NNLO, the argument that the
“valence” heavy and light quarks behave the same way can also be extended to a
“valence” gluon upon the replacement of overall color factors. Indeed, the
terms $\sim 1/(1-x)_{+}$ in the diagonal matrix element $A^{{\rm
S},(2)}_{gg,H}$ coincide with the ones in $A^{{\rm NS},(2)}_{qq,H}$ after the
color factor replacement $C_{A}\to C_{F}$.

Regarding our numerical implementation, we do not resum the $\ln(N)$ term in
$C_{\rm dec}$ since, as it has been pointed out in ref. Neubert:2007je , it is
of UV origin. Similar approach to this term has been taken in ref.
Fickinger:2016rfd . We have checked that the numerical contribution from this
term is very small, less than 1% of the hard matching function $H^{\rm ini}$
for $N<450$.

The evaluation of the integrals in eq. (2.3) is standard: one first performs
the integration over $k^{2}$ via a change of variables from $k^{2}$ to
$\alpha_{S}(k^{2})$ (first line) or $\alpha_{S}^{(n_{l})}(k^{2})$ (second
line), perturbatively expanding the integrand in $\alpha_{S}$. This transforms
the integrands into sums of powers of $\alpha_{S}$, which can be integrated
trivially. The result of the integral over $k^{2}$ on the first line of eq.
(2.3) is independent of $z$, so the integral over $z$ on the same line can
also be performed trivially. The result of the integral over $k^{2}$ on the
second line, however, is expressed in terms of
$\alpha_{S}^{(n_{l})}((1-z)^{2}m_{b}^{2})$. The integral over $z$ is performed
by another change of variables from $z$ to
$\alpha_{S}^{(n_{l})}((1-z)^{2}m_{b}^{2})$, perturbatively expanding the
integrand in $\alpha_{S}^{(n_{l})}((1-z)^{2}m_{b}^{2})$. The handling of the
term $z^{N-1}-1$ through NNLL accuracy can be found in appendix A of ref.
Catani:2003zt . This procedure yields a result for eq. (2.3) which is
expressed in terms of $\alpha_{S}(m_{b}^{2})$, $\alpha_{S}(\mu_{Fr,0}^{2})$,
$\alpha_{S}^{(n_{l})}(m_{b}^{2}/N^{2})$ and $\alpha_{S}^{(n_{l})}(m_{b}^{2})$,
all of which are re-expressed perturbatively in terms of
$\alpha_{S}(\mu_{R,0}^{2})$, using the decoupling relations for the strong
coupling constant Larin:1994va ; Chetyrkin:1997un ; Schroder:2005hy ;
Chetyrkin:2005ia  to replace $\alpha_{S}^{(n_{l})}(m_{b}^{2})$ with
$\alpha_{S}(m_{b}^{2})$. As a result, one arrives at the following explicit
expression for the function $G^{\rm ini}$ defined in eq. (2.3)

$G^{\rm
ini}(\mu_{Fr,0},\mu_{R,0},m_{b},N,\alpha_{S}(\mu_{R,0}))=\sum_{i=1}^{\infty}\alpha_{S}^{i-2}(\mu_{R,0})g_{i}(\lambda)\,,$
(16)

with $\lambda=\alpha_{S}(\mu_{R,0})b_{0}\ln(N)$. The QCD $\beta$-function
coefficients relevant for this work can be found in appendix A. Only the terms
with $i=1,2,3$ in eq. (16) contribute through NNLL. The explicit expressions
for the functions $g_{1},g_{2}$ and $g_{3}$ can be found in appendix A.1.

It is a well-known fact that the functions $g_{i}$ contain terms of the form
$\ln^{j}(1-2\lambda)/(1-2\lambda)^{k}$ which become singular for
$\lambda=1/2$. In terms of the Mellin variable $N$, this corresponds to a
singularity at $N=N_{L}\approx 34$. While such singularities can be
consistently treated within the Minimal Prescription approach Catani:1996yz
for inverting resummed expressions from $N$ to $z$ space, it is known that
such resummed expressions tend to become negative for large values of $z$.
This, in turn, may reduce the usefulness of the soft gluon resummation. Ref.
Cacciari:2005uk  has introduced a way around this issue, based on remapping
the Mellin variable $N$ in the soft-gluon resummed exponent

$N\to N\frac{1+f/N_{L}}{1+fN/N_{L}}\,.$ (17)

The parameter $f$ is chosen to have a value $f=1.25$.

Given the importance of the large $z$ region and the NNLO+NNLL precision of
our work, we have found it important to reanalyse the above issues. In fact,
we propose another prescription, where the function $G^{\rm ini}$ is truncated
in an expansion around $\lambda=0$

$\displaystyle G^{\rm
ini}(N,\alpha_{S})=\phantom{g_{1,2}\alpha_{S}\ln^{2}(N)+g_{1,3}\alpha_{S}^{2}\ln^{3}(N)+g_{1,4}\alpha_{S}^{3}\ln^{4}(N)+g_{1,5}\alpha_{S}^{4}\ln^{5}(N)+g_{1,6}}$
$\displaystyle
g_{1,2}\alpha_{S}\ln^{2}(N)+g_{1,3}\alpha_{S}^{2}\ln^{3}(N)+g_{1,4}\alpha_{S}^{3}\ln^{4}(N)+g_{1,5}\alpha_{S}^{4}\ln^{5}(N)+g_{1,6}\alpha_{S}^{5}\ln^{6}(N)$
$\displaystyle+g_{2,1}\alpha_{S}\ln(N)+g_{2,2}\alpha_{S}^{2}\ln^{2}(N)+g_{2,3}\alpha_{S}^{3}\ln^{3}(N)+g_{2,4}\alpha_{S}^{4}\ln^{4}(N)$
$\displaystyle+g_{3,1}\alpha_{S}^{2}\ln(N)+g_{3,2}\alpha_{S}^{3}\ln^{2}(N)$
$\displaystyle\,,$ (18)

where $g_{i,j}$ is the coefficient of $\ln^{j}(N)$ in the expansion of the
function $g_{i}$ around $\lambda=0$:

$g_{i}(\lambda)=\sum_{j=1}^{\infty}g_{i,j}\alpha_{S}^{j}\ln^{j}(N)\;,\qquad
g_{1,1}=0\;.$ (19)

The depth of the expansion of the functions $g_{1},g_{2}$ and $g_{3}$ in eq.
(18) is sufficient for NNLL accuracy as we now explain.

Define the NnLL accurate truncation

$G^{\rm
ini}(N,\alpha_{S})\Big{|}_{\text{N${}^{n}$LL}}=\sum_{i=1}^{n+1}\alpha_{S}^{i-2}g_{i}(\lambda)\;,$
(20)

with

$\exp[G^{\rm
ini}(N,\alpha_{S})\Big{|}_{\text{N${}^{n}$LL}}]=1+\sum_{l=1}^{\infty}\sum_{k=1}^{2l}c_{l,k}\alpha_{S}^{l}\ln^{k}(N)\;,$
(21)

where $c_{l,k}$ are polynomial in $g_{i,j}$, $i\leq n+1$.

At Nn+1LL accuracy, the coefficients $c_{l,k}$ with $k<2(l-n)$ depend on
$g_{n+2,j}$. Indeed, write

$\exp[G^{\rm ini}(N,\alpha_{S})\Big{|}_{\text{N${}^{n+1}$LL}}]=\exp[G^{\rm
ini}(N,\alpha_{S})\Big{|}_{\text{N${}^{n}$LL}}]\exp[\alpha_{S}^{n}\,g_{n+2}]\;.$
(22)

The LL term $\big{(}g_{1,2}\alpha_{S}\ln^{2}(N)\big{)}^{l-n-1}/(l-n-1)!$ of
the expansion of $\exp[G^{\rm ini}(N,\alpha_{S})\Big{|}_{\text{N${}^{n}$LL}}]$
enhances the Nn+1LL term $\alpha_{S}^{n}\,g_{n+2,1}\alpha_{S}\ln(N)$ of the
expansion of $\exp[\alpha_{S}^{n}\,g_{n+2}]$ to contribute

$g_{1,2}^{l-n-1}g_{n+2,1}\alpha_{S}^{l}\ln^{2(l-n)-1}(N)/(l-n-1)!\;,$ (23)

to $\exp[G^{\rm ini}(N,\alpha_{S})\Big{|}_{\text{N${}^{n+1}$LL}}]$. Hence,
$c_{l,2(l-n)-1}$ contains the term $g_{1,2}^{l-n-1}g_{n+2,1}/(l-n-1)!$.
Similar arguments apply for $c_{l,k}$ with $k<2(l-n)-1$.

Since $c_{l,k}$ with $k<2(l-n)$ requires Nn+1LL resummation, the NnLL resummed
expression is correct at most for coefficients $c_{l,k}$ with $k\geq 2(l-n)$.
Furthermore, two expressions that only differ in the coefficients $c_{l,k}$
with $k<2(l-n)$ have the same NnLL accuracy. Hence, any coefficient $g_{i,j}$
that does not contribute to $c_{l,k}$ with $k\geq 2(l-n)$ but contributes to
$k<2(l-n)$ is irrelevant for NnLL accuracy and can be modified at will. The
maximal $k$ of $c_{l,k}$ to which $g_{i,j}$ contributes, is determined
similarly as above eq. (23) by considering LL enhancement

$\big{(}g_{1,2}\alpha_{S}\ln^{2}(N)\big{)}^{l-i-j+2}g_{i,j}\alpha_{S}^{i-2+j}\ln^{j}(N)=g_{1,2}^{l-i-j+2}g_{i,j}\,\alpha_{S}^{l}\ln^{2(l-i+2)-j}(N)\;.$
(24)

Thus, if

$2(l-i+2)-j<2(l-n)\qquad\Leftrightarrow\qquad j>2(n-i+2)\;,$ (25)

then $g_{i,j}$ is irrelevant for NnLL resummation. In our prescription, we
simply set such $g_{i,j}$ to zero, which results in eq. (18) for $n=2$. Notice
finally that, owing to the negative sign of all coefficients $g_{1,j}$, the
numeric inversion of the resummed $D^{\rm ini}$ within the expanded
prescription is unproblematic.

The Landau pole discussed above is not specific to $D^{\rm ini}$. From the
results in appendix A.2 one can see that it is also present in the
${e^{+}e^{-}}$ coefficient function, albeit at a different value for $N$, due
to terms $\ln^{j}(1-\lambda)/(1-\lambda)^{k}$. In order to deal with this pole
one can choose to either remap the Mellin variable $N$ or, as we have done in
this work, truncate the corresponding function $G$ for the ${e^{+}e^{-}}$
coefficient function in analogy with eq. (18). Note that unlike the case of
$D^{\rm ini}$, the leading coefficients $g_{i,j}$ of the ${e^{+}e^{-}}$
coefficient functions are positive which precludes the inversion to $z$ space
of the coefficient function. However, the $N$-space product of the resummed
coefficient function and $D^{\rm ini}$ is invertible as long as
$2\alpha_{S}(\mu_{R,0})>\alpha_{S}(\mu_{R})$. This requirement is always
satisfied since $\mu_{R,0}\leq\mu_{R}$. For this reason in our implementation
we first multiply in $N$-space the resummed expressions for $D^{\rm ini}$ and
the ${e^{+}e^{-}}$ coefficient functions and only invert to $z$-space their
product. This combination of $D^{\text{ini}}$ and the coefficient function is
matched to its perturbative expansion in $\alpha_{s}(\mu_{R,0})$ and
$\alpha_{s}(\mu_{R})$ through
$\mathcal{O}(\alpha_{s}(\mu_{R,0})^{2},\alpha_{s}(\mu_{R,0})\alpha_{s}(\mu_{R}),\alpha_{s}(\mu_{R})^{2})$.

The three possible ways of dealing with the Landau pole in the soft gluon
resummed result, namely, no prescription, rescaling of $N$ (17) and truncation
of $G$ (18), are compared in fig. 1. Also shown in that figure is the
numerical impact of the additive (4) versus multiplicative (5) matching. The
$N$ variable is only remapped in the resummed part of the PFF, i.e. in the
term $H^{\rm ini}\exp(G^{\rm ini})$ but not in the matching term $M^{\rm
ini}$. We show separately the PFFs and their convolution with a fixed NPFF,
specifically, the NPFF fitted in ref. Cacciari:2005uk

$D_{\text{NP}}(x)=\frac{\Gamma(a+b+2)}{\Gamma(a+1)\Gamma(b+1)}x^{a}(1-x)^{b},\;\;\;\;a=24\pm
2,\;b=1.5\pm 0.2\;.$ (26)

All PFFs are with NNLO+NNLL accuracy and are evaluated for
$\mu_{Fr,0}=\mu_{R,0}=m_{b}$. The green curve shows the PFF with a truncated
Sudakov exponent $G^{\rm ini}$ that is matched multiplicatively. It
corresponds to the PFF used in this work. The unregulated (in black) and
remapped (in red) PFFs are matched additively.

Figure 1: Different ways for dealing with the Landau pole in the soft gluon
resummed PFF. All curves are at NNLO+NNLL for $\mu_{Fr,0}=\mu_{R,0}=m_{b}$.
Shown are both the PFFs (left) and their convolution with a fixed NPFF eq.
(26) (right). The various curves are: the PFF without any regularization of
the Landau pole (black), with $N$ remapped (red) and the PFF with a truncated
exponent matched to the fixed-order result additively (blue) and
multiplicatively (green). The coloured bands in the plot on the right
correspond to the uncertainty on the NPFF as shown in eq. (26).

From fig. 1 we conclude that the impact of the way the Landau pole in the
resummed result for $D^{\rm ini}$ is handled can be significant. In
particular, the PFF without any regularization becomes negative and a
significant portion of its area is in this region. This negative contribution
is compensated by a much larger value in the peak region. On the other hand,
the two regularized expressions, with remapped $N$ and with truncated $G^{\rm
ini}$, are much closer to each other than the un-regularized result. Still,
their behaviors in the large $z$ region show notable differences. The remapped
PFF does not vanish at large $z$ while the truncated one rapidly approaches
zero as $z\to 1$. Given that both curves are normalized, this difference
translates into differences for both their peak values and the positions of
their maxima. The phenomenological significance of this difference is harder
to assess since although the two curves are similar, they are not consistent
with each other within the NPFF uncertainty. The un-regularized one (in black)
is significantly different from the other two curves. We also note that the
difference between the additive (in blue) and multiplicative (in green)
matching is negligible everywhere except the extremely large $z$ region where
the curve with additive matching becomes slightly negative above $z=0.98$ and
diverges in the limit $z\to 1$, while the multiplicative one stays non-
negative for all values of $z$ and vanishes in the limit $z\to 1$.

Before closing this section, we would like to make one more remark. When
estimating the uncertainty due to scale variation we do not vary the scale
$\mu_{Fr,0}$ but keep it fixed at $m_{b}$. The reason for this is a subtlety
in the soft-gluon resummation of $D^{\rm ini}$. As is well known and as can be
seen from eq. (8), due to the running of the coupling to scales as low as the
Landau pole $\Lambda$, the resummed result breaks down for $z\geq z_{max}$,
where $z_{max}$ is defined by the relation $\mu(1-z_{max})\sim\Lambda$. The
scale $\mu$ is a fixed hard scale specific to the resummed observable. For the
case of $D^{\rm ini}$, one has $\mu=\mu_{R,0}$ which is quite low. Indeed, as
can be seen in fig. 1, this breakdown of the resummed result is consistent
with the behavior of the curve with no regularization which turns negative at
large $z$. Thus, variation of the scale $\mu_{R,0}$ by a factor of two can
noticeably change the value of $z_{max}$ and, therefore, the large-$z$
behavior of the resummed PFF. We also note that varying $\mu_{R,0}$ below
$m_{b}$ would be problematic in the context of the decoupling of the heavy
flavor. For these reasons in this work we keep the scale $\mu_{R,0}$ fixed at
$m_{b}$. We have verified that the difference between the $b$-quark PFF
evaluated at $\mu_{R,0}=m_{b}$ and at $\mu_{R,0}=2m_{b}$ is tiny everywhere
except for large $z$, as expected, due to the shift in the breakdown point
$z_{max}$. We also note that in this context the scales $\mu_{Fr,0}$ and
$\mu_{R,0}$ are typically set equal, meaning that the above comments are also
applicable for the setting of the scale $\mu_{Fr,0}$. We have checked that the
dependence of the PFF on the scale $\mu_{Fr,0}$ is much milder than the one on
$\mu_{R,0}$.

###  2.4 Fits to data

Having specified our implementation of the ${e^{+}e^{-}}$ coefficient
functions and perturbative fragmentation functions, we are ready to proceed
with the fit of the NPFF from ${e^{+}e^{-}}$ data. We perform a combined fit
to datasets from ALEPH ALEPH:2001pfo , DELPHI DELPHI:2011aa , OPAL
OPAL:2002plk  and SLD SLD:2002poq . We note that these datasets contain
weighted averages of various $B$ mesons, and some of the sets include baryons.
This implies that predictions based on NPFF extracted from these data sets
will also be for the same combination of $B$ hadrons.

####  2.4.1 Treatment of the datasets

The data – each dataset, in fact – correspond to the normalized differential
distribution

${\rm Data}_{B}(z)=\frac{1}{\sigma_{B}}\frac{d\sigma_{B}(m_{b},z)}{dz}\,.$
(27)

In practice, the distribution ${\rm Data}_{B}(z)$ is not measured directly.
The experiments’ published data represents the binned version of this
differential distribution, supplemented by a covariance matrix representing
the correlated experimental uncertainty of this binned measurement. To aid our
fits, we would like to represent the binned data, including its uncertainty,
in a differential form like in eq. (27). To that end we proceed as follows. We
interpret the binned data and its uncertainty as a multidimensional normal
random distribution. Since the data is normalized, the dimension of that
distribution is given by the number of bins minus one. We then generate 100
binned “datasets”; one of them corresponds to the measured binned data (i.e.
the “central value”) while the remaining 99 sets are generated randomly, by
sampling this multidimensional random distribution. We then treat each one of
these 100 datasets as if it is the true data and produce a data fit of the
form ${\rm Data}_{B}(z)$ for each one of those binned sets. As a whole, the
100 data fits represent the true differential data distribution and its
uncertainty.

As a last step, the 100 fitted data distributions are used to fit a set of 100
NPFF’s, one being the central member, accompanied by 99 error sets (the
detailed description of this procedure is given in sec. 2.4.2 below). When
predictions based on this set of NPFFs are made, 100 predictions need to be
computed. One represents the central prediction for the given observable while
the other 99 predictions are used to derive the NPFF uncertainty for this
variable. This uncertainty is taken to be the root mean square (r.m.s.)
$1\sigma$ deviation with respect to the mean of the full FF set of 100 curves.
Throughout this work we present this uncertainty in such a way that it is
always centered around the mean and not around the central member.

For our data fits we assume two functional forms. One has 5 parameters (as in
ref. DELPHI:2011aa ), while the other has 8

$\displaystyle{\rm Data}_{B}(z)$ $\displaystyle=$ $\displaystyle
n\frac{\Gamma(a+b+2)}{\Gamma(a+1)\Gamma(b+1)}z^{a}(1-z)^{b}+(1-n)\frac{\Gamma(c+d+2)}{\Gamma(c+1)\Gamma(d+1)}z^{c}(1-z)^{d}\,,$
$\displaystyle{\rm Data}_{B}(z)$ $\displaystyle=$ $\displaystyle
m\bigg{[}n\frac{\Gamma(a+b+2)}{\Gamma(a+1)\Gamma(b+1)}z^{a}(1-z)^{b}+(1-n)\frac{\Gamma(c+d+2)}{\Gamma(c+1)\Gamma(d+1)}z^{c}(1-z)^{d}\bigg{]}$
(28)
$\displaystyle+(1-m)\frac{\Gamma(e+f+2)}{\Gamma(e+1)\Gamma(f+1)}z^{e}(1-z)^{f}\,.$

We perform fits for each individual measurement as well as to their
combination. The following procedure was followed for both models. First, the
statistical and systematic covariance matrices of the data were added. The
OPAL data has asymmetrical systematic uncertainties, meaning four systematic
covariance matrices are needed, but only two are made available. Because of
this, the two provided systematic covariance matrices were averaged before
being added to the statistical one. SLD only provides the diagonal elements of
the covariance matrices, so the off-diagonal elements were set to zero. The
SLD data also has two systematic uncertainties, which were added in
quadrature. All covariance matrices were then diagonalized by determining
their eigenvalues and eigenvectors. This procedure is necessary for both the
ALEPH and OPAL datasets, since the provided covariance matrices are not
positive-definite, likely due to rounding errors. By diagonalizing the
covariance matrices it is possible to reject all negative eigenvalues and only
use the positive ones in the fit.

Dataset | ALEPH | DELPHI | OPAL | SLD
---|---|---|---|---
Number of bins | $19$ | $9$ | $20$ | $22$
Number of positive eigenvalues | $13$ | $9$ | $13$ | $22$
Number of used eigenvalues ($5$ parameters) | $7$ | $8$ | $8$ | $22$
$\chi^{2}/n_{\text{dof}}$ ($5$ parameters) | $3.41/2$ | $1.76/3$ | $4.97/3$ | $5.31/17$
Number of used eigenvalues ($8$ parameters) | $13$ | $9$ | $12$ | $22$
$\chi^{2}/n_{\text{dof}}$ ($8$ parameters) | $5.39/5$ | $0.574/1$ | $4.09/4$ | $1.06/14$
Table 1: The number of bins, number of positive eigenvalues of the covariance
matrix, number of used eigenvalues for both central fits and the
$\chi^{2}/n_{\text{dof}}$ for both central fits for all four datasets.

It is still possible, however, that some of the positive eigenvalues are
artificially small due to the rounding errors. Because of this, only the $n$
largest eigenvalues were included in the fit, where $n$ is chosen such that
$n+1$ is the smallest number of used eigenvalues that leads to an unacceptably
large $\chi^{2}/n_{\text{dof}}$ if the data model is fitted to the respective
dataset only. Since, typically, the difference in $\chi^{2}$ between $n$ and
$n+1$ is at least a factor of ten, this choice is quite unambiguous in
practice, despite the seemingly subjective criterion. The results of the fits
to each individual dataset are given in table 1.

Figure 2: Comparison of the 5-parameter (blue) and 8-parameter (orange) fits
for the distribution ${\rm Data}_{B}(z)$. Shown are the absolute distributions
(top left) and their ratio (bottom left) as well as their absolute (top right)
and relative (bottom right) uncertainties. The uncertainty bands correspond to
the root mean square (r.m.s.) of the set of 100 data fits.

Finally, a fit is performed on the combined data. For the $5$-parameter model,
this yields $\chi^{2}/n_{\text{dof}}=85.4/40$, whereas the $8$-parameter fit
has $\chi^{2}/n_{\text{dof}}=115/48$. Since both the 5- and 8-parameter models
can describe each individual dataset very well by construction, these
unacceptably large values of $\chi^{2}/n_{\text{dof}}$ must be due to tensions
between the datasets. We resolve this by rescaling all covariance matrices by
$2.14$ for the $5$-parameter case and by $2.40$ for the $8$-parameter case,
thus enforcing $\chi^{2}/n_{\text{dof}}=1$ for the combined fit. This
procedure is analogous to the one used in ref. DELPHI:2011aa . The additional
99 variations of the datasets are then generated using the rescaled covariance
matrices.

Table 1 shows that more data points are used for the $8$-parameter fit than
for the $5$-parameter one. Normally, this should make the $8$-parameter fit
more precise. However, as discussed above, the tensions between the datasets
are larger for the $8$-parameter case, so the errors are rescaled by a larger
factor. As a result the final uncertainties of the $8$\- and $5$-parameter
fits are very similar. This can be seen in fig. 2 which shows a comparison of
both fits, including their uncertainties. The uncertainty bands correspond to
the r.m.s. of the complete set of data fits.

Due to its more restrictive functional form, the $5$-parameter fits exhibit a
somewhat sharp turn near $x\sim 0.6$. The fits based on the $8$-parameter
model do not have this behavior. Apart from this feature, the two fit models
are essentially equivalent.

It is interesting to compare the results from our combined 5-parameter fit to
the combined fit performed in ref. DELPHI:2011aa , since both fits are based
on the same datasets and utilize the same functional forms. This comparison
can be found in fig. 3. The two uncertainty bands are modeled differently: our
uncertainty band (in blue) corresponds to the r.m.s. of the 100 data fits
while the DELPHI error band (in brown) corresponds to the propagated errors on
the fitted parameters taken from table 12 of ref. DELPHI:2011aa . The error
bands of our new fits are significantly smaller. This is partly due to the
fact our new fits retain more eigenvalues of the covariance matrices of the
data sets, in particular of the one from OPAL. We have checked that even when
using the same eigenvalues as in ref. DELPHI:2011aa , the error bands obtained
using our method are still significantly smaller. We suspect the reason for
this difference is related to an inconsistency between the published results
for the data fit and the first 40 moments of the same fit also given in that
reference. Assuming for a moment that the published moments are correct, the
size of the uncertainties of the fit parameters must be smaller by a factor of
about 5 which would make the uncertainties roughly equal. As an argument that
the uncertainty estimate of our data fit is of the correct magnitude, we
compare in the following section the size of the uncertainty of the derived
NPFF and find it to be of similar size to other extractions in the literature.
We have also verified that when using the same number of eigenvalues as ref.
DELPHI:2011aa  our fit reproduces the first 10 moments given explicitly in
table 14 of that reference within 1 permile, while the moments of our fit
differ from the moments of the fit published in ref. DELPHI:2011aa  by about 4
percent.

Figure 3: Comparison of the $5$-parameter fit (blue) with the fit from ref.
DELPHI:2011aa  (brown). Absolute distributions (left) and their ratio (right).
The blue error band corresponds to the r.m.s. of the 100 data fits while the
brown error band corresponds to the propagated errors on the fitted parameters
taken from table 12 of ref. DELPHI:2011aa .

####  2.4.2 NPFF fitting procedure

The NPFF is extracted by comparing theory prediction eqns. (1) and (2) with
the combined data as in eq. (27). Since the data fits themselves already
encode the experimental uncertainties, the NPFF should not be fitted to the
data fits, but rather obtained by exactly deconvolving the data fits using the
perturbative prediction. However, since this cannot be achieved analytically,
the NPFF is actually fitted to the data fits, but with the ultimate goal to
achieve $\chi^{2}=0$, unlike a proper fit in the usual sense which should have
$\chi^{2}/n_{\text{dof}}=1$. Because of this, there is no reason to limit the
number of degrees of freedom. In fact, it is desirable to use as many degrees
of freedom as computationally acceptable, since this will enable the most
precise approximation of the exact deconvolution. We have found it best to use
a $z$-space parametrization consisting of a high-degree polynomial
supplemented with a $\delta$-function

${D^{\rm NP}_{b\to B}}(z)=a_{0}\delta(1-z)+\sum_{i=3}^{31}a_{i}z^{i}(1-z)\,.$
(29)

The coefficients of the polynomial in eq. (29) are fixed by minimizing the
squared difference between the data fit eq. (27) and the theory prediction,
integrated over the range $0.1\leq z\leq 0.95$. This range is chosen in order
to avoid the low $z$ region which is sensitive to missing mass power
corrections and the extremely large $z$ region where soft-gluon resummation
starts to be important. We repeat this procedure for each one of the 100 data
fits, thus producing 100 NPFFs. In this we utilize the 8 parameter data fit
discussed in sec. 2.4.

The coefficient of the $\delta$-function is fixed by normalization. Its value
for the mean of the NNLO NPFF fits reads

$a_{0}=0.159\pm 0.010\,.$ (30)

The reason we have introduced a $\delta$-function in our fit is to increase
the hardness of the NPFF since we have not been able to achieve this
satisfactorily by means of a pure polynomial fit. Indeed, as can be concluded
from $a_{0}$’s fit value in eq. (30), the coefficient of the $\delta$-function
is incompatible with zero. Note that a $\delta$-function contribution in the
NPFF has already been used in the literature in the past Cacciari:2005uk .

The new NPFF set derived in this work, including its uncertainty, is shown in
fig. 4. Also shown there is a comparison to existing NPFF sets. One should
keep in mind that due to the $\delta$-function term – which is not plotted –
the normalization of our new set appears different. As can be seen from eq.
(30) the coefficient of this term is about 0.16. Normalization aside, we note
that the shape of the new NPFF agrees with the one from ref. Cacciari:2005uk
but not with the one from ref. Fickinger:2016rfd . That later set was used
previously in ref. Czakon:2021ohs . The difference between the new set derived
in this work and the one from ref. Fickinger:2016rfd  illustrates the need for
a NPFF which is consistently derived within our formalism.

It is interesting to clarify the reasons behind the differences in shapes and
uncertainties between the three sets shown in fig. 4. First there are
differences in the perturbative parts: the set of ref. Cacciari:2005uk  is at
NLO while the one of ref. Fickinger:2016rfd  is with NNLO accuracy but is
obtained in a slightly different formalism (a detailed discussion can be found
in ref. Czakon:2021ohs ). There are differences in the datasets that have been
used in the fits: ref. Cacciari:2005uk  only fits the ALEPH and SLD datasets.
Since the ALEPH data is not very constraining, the error bands of ref.
Cacciari:2005uk  originate, essentially, from the SLD data. Furthermore, ref.
Cacciari:2005uk  uses a two-parameter model for the NPFF and does not provide
a correlation matrix for the uncertainties on the fit parameters. For this
reason the error bands for that set that can be seen in fig. 4 are likely an
overestimate. Both of those factors combined could explain why our new NPFF
appears to have much smaller uncertainty. The NPFF of ref. Fickinger:2016rfd
does not include in their fit the OPAL data, while using less eigenvalues from
the other experiments than our new analysis. This explains the reduced
uncertainty of our new set relative to the one of ref. Fickinger:2016rfd .

Figure 4: Left: a comparison of the new NPFF fit (blue) with the one from ref.
Fickinger:2016rfd  (orange) and ref. Cacciari:2005uk  (green). The error bands
correspond to one standard deviation. Center: the ratio plot of the full new
NPFF set to its central member, and (right) the relative r.m.s. of the new
NPFF set.

One may wonder how the individual datasets that enter our combination affect
the determination of the NPFF. To clarify this point in appendix B we discuss
the extraction of NPFF from each individual dataset and compare the result
with the combined extraction.

Before closing this section, we would like to briefly discuss alternative
fitting models. We have also tried to fit the NPFF directly in Mellin space.
To that end we have utilized the first 30 or so moments of the ${e^{+}e^{-}}$
coefficient functions and PFF. We have assumed that in $z$-space the NPFF is
built out of an exponential (like the one used in ref. Fickinger:2016rfd  )
times a polynomial of a high degree and have inverted this expression
analytically to $N$-space. We have been able to fit all of the first 30
moments with very high precision, typically better than few permil. Despite
that, the predicted $z$-space distribution at the level of $B$ hadrons does
not fit well the data it was derived from in the very large $z$ region and
close to the peak of the distribution. This seem to suggest that this approach
is limited by biases in NPFF’s assumed functional form.

The results for all fits derived in this work are attached as convenient to
use LHAPDF Buckley:2014ana  interpolation grids for the set of fragmentation
functions $D_{i\to B}$, evolved to a scale $\mu_{Fr}$. Note that, as in ref.
Czakon:2021ohs , these grids have been symmetrized with respect to charge.
This means that any results generated with these grids cannot distinguish
between hadrons containing a $b$-quark and hadrons containing a
$\bar{b}$-quark; this restriction also applies to the calculations presented
in sec. 3 below. The $b\leftrightarrow\bar{b}$ symmetrization is not an
intrinsic limitation of the calculational approach, however. If desirable,
future calculations can instead employ unsymmetrized FF sets, allowing a
distinction of B-hadrons based on $b$-quark charge to be made.

###  2.5 NPFF for $B$ decays to a $J/\psi$ or to a soft muon

Not all measurements of $B$-flavored hadrons measure the $B$ directly. What is
often measured are $B$ decay products like $J/\psi$ or the $\mu$ originating
in semileptonic $B$ decays. Due to their very clean signatures, these final
states have been advocated in the context of precise top quark mass
determination, see refs. Kharchilava:1999yj ; CDF:2009mbf ; CMS:2016ixg ;
ATLAS:2022jbw . In the following we explain how these decays can be
incorporated into the $B$-fragmentation formalism discussed in this work.

Consider the decay $B\to d+X$ of a spin-0 $B$-meson at rest, where $d$ is a
descendent particle with mass $m_{d}$. Keeping the decay rate fully
differential in the momentum of $d$, the decay rate can be written as

$d\Gamma(B\to
d+X)=\frac{1}{4\pi}f(E_{d}^{\text{rest}})dE_{d}^{\text{rest}}d\cos(\theta)d\phi\,,$
(31)

where $E_{d}^{\text{rest}}$ is the energy of the descendant in the $B$-hadron
rest frame and the normalization is such that

$\frac{d\Gamma(B\to d+X)}{dE_{d}^{\text{rest}}}=f(E_{d}^{\text{rest}})\,.$
(32)

Introducing the dimensionless variable $y=E_{d}^{\text{rest}}/m_{B}$, with
$0\leq y\leq 1$, eq. (31) becomes

$d\Gamma(B\to d+X)=\frac{m_{B}}{4\pi}f(y\>m_{B})dyd\cos(\theta)d\phi\,.$ (33)

This phase space can next be boosted from the $B$-hadron rest frame to a frame
with $E_{B}\gg m_{B}$. Without loss of generality the boost can be taken to be
along the $z$-axis. The goal is to express the energy of $d$ in this boosted
frame as a fraction $z$ of $E_{B}$, i.e. $E_{d}=zE_{B}$. Using

$\displaystyle\delta\bigg{(}z-\frac{E_{d}}{E_{B}}\bigg{)}$
$\displaystyle=\delta\bigg{(}z-\gamma_{B}\frac{E_{d}^{\text{rest}}+\beta_{B}\cos(\theta)\sqrt{(E_{d}^{\text{rest}})^{2}-m_{d}^{2}}}{E_{B}}\bigg{)}$
$\displaystyle\approx\delta\bigg{(}z-\frac{E_{d}^{\text{rest}}+\cos(\theta)\sqrt{(E_{d}^{\text{rest}})^{2}-m_{d}^{2}}}{m_{B}}\bigg{)}$
$\displaystyle=\delta\bigg{(}z-y-\cos(\theta)\sqrt{y^{2}-\frac{m_{d}^{2}}{m_{B}^{2}}}\bigg{)}\;,$
(34)

where on the second line $\gamma_{B}=E_{B}/m_{B}$ and $\beta_{B}\approx 1$
were used, this yields

$\displaystyle\frac{d\Gamma(B\to d+X)}{dz}$
$\displaystyle=\int_{0}^{2\pi}\int_{-1}^{1}\int_{0}^{1}\frac{m_{B}}{4\pi}f(y\>m_{B})\delta\bigg{(}z-y-\cos(\theta)\sqrt{y^{2}-\frac{m_{d}^{2}}{m_{B}^{2}}}\bigg{)}dyd\cos(\theta)d\phi$
$\displaystyle=\int_{0}^{1}\frac{m_{B}}{2\sqrt{y^{2}-\frac{m_{d}^{2}}{m_{B}^{2}}}}f(y\>m_{B})\theta\bigg{(}1-\frac{(z-y)^{2}}{y^{2}-\frac{m_{d}^{2}}{m_{B}^{2}}})\bigg{)}dy$
$\displaystyle=\int_{\frac{z}{2}+\frac{m_{d}^{2}}{2z\>m_{B}^{2}}}^{1}\frac{m_{B}}{2\sqrt{y^{2}-\frac{m_{d}^{2}}{m_{B}^{2}}}}f(y\>m_{B})dy\equiv\Gamma_{B}\>D_{B\to
d}(z)\;.$ (35)

The function $D_{B\to d}(z)$ can be interpreted as a fragmentation function
for the transition $B\to d$. Since the consecutive fragmentations ${D^{\rm
NP}_{b\to B}}$ and $D_{B\to d}$ can be combined via convolution, the decay
$B\to d+X$ can be included in the calculation by simply convolving ${D^{\rm
NP}_{b\to B}}$ and $D_{B\to d}$ and then running the calculation as one
normally would for the production of $B$-hadrons. The only ingredient
necessary in order to incorporate $B$-hadron decays in a fragmentation
calculation is $f(E_{d}^{\text{rest}})$, i.e. $d\Gamma(B\to
d+X)/dE_{d}^{\text{rest}}$. Such energy spectra are highly constrained by very
precise measurements at $B$-factories, which means that the total theory error
is essentially unaffected by the further inclusion of the $B$ decay.

The decay energy spectra $B\to J/\psi+X$ and $B\to\mu+X$ needed in this work
are obtained from ver. 2.1.1 of the program EvtGen Lange:2001uf . We study
$B$-hadron decays with $J/\psi$ or muons for several different types of
$B$-hadrons: $B^{+}$, $B^{0}$, $B_{s}^{0}$ and $\overline{\Lambda}_{b}^{0}$.
The mixture of $B$-hadrons produced in high-energy collisions corresponds to
$40.8\%$ $B^{\pm}$, $40.8\%$
$\stackon[.1pt]{$B$}{\raisebox{-0.6pt}{\scalebox{0.2}{(}}\raisebox{-1.9pt}{\scalebox{0.7}{$-$}}\raisebox{-0.6pt}{\scalebox{0.2}{)}}}^{0}$,
$10.0\%$
$\stackon[.1pt]{$B$}{\raisebox{-0.6pt}{\scalebox{0.2}{(}}\raisebox{-1.9pt}{\scalebox{0.7}{$-$}}\raisebox{-0.6pt}{\scalebox{0.2}{)}}}_{s}^{0}$
and $8.4\%$ $B$-baryons ParticleDataGroup:2020ssz . In the following, it is
assumed that the vast majority of produced B-baryons are
$\stackon[.1pt]{$\Lambda$}{\raisebox{-0.6pt}{\scalebox{0.2}{(}}\raisebox{-1.9pt}{\scalebox{0.7}{$-$}}\raisebox{-0.6pt}{\scalebox{0.2}{)}}}_{b}^{0}$’s,
such that the contributions from other $B$-baryons can be neglected.

####  2.5.1 Muon spectra

In $B\to\mu+X$ decays one needs to distinguish between same-sign decays
($\overline{b}\to\mu^{+}$) and opposite-sign decays
($\overline{b}\to\mu^{-}$), as well as between prompt muons
($\overline{b}\to\mu^{+}$) and non-prompt muons (mostly
$\overline{b}\to\overline{c}\to\mu^{-}$). There are also two types of spectra,
corresponding to different ways of treating final states with multiple muons:
multiplicities (one counts all muons) and hardest-only (one counts only the
hardest muon). Additionally, we could study the spectra in the rest frame or
in the boosted one. This results in a total of $144$ different spectra (a few
of which are identically equal to zero), all of which have been generated
using EvtGen. We note that eq. (35) assumes isotropic decays, which is a valid
assumption for spin-0 $B$-mesons, but not for $\overline{\Lambda}_{b}^{0}$.
This does not appear to be an issue in practice since, we have checked, the
analytically boosted spectra and the MC histograms revealed no systematic
difference.

We fit the total (the sum of same-sign, opposite-sign, prompt and non-prompt)
hardest-only spectrum of each individual $B$-hadron species. This is the type
of spectrum needed to describe the measurement presented in ref. ATLAS:2022jbw
. The choice to include both same-sign and opposite-sign decays is based on
the fact that, as explained in sec. 2.4.2, the $B$-hadron FF sets used in this
work have been charge-symmetrized. All spectra are normalized to unity for the
fit and then rescaled to their physical normalization according to EvtGen. The
fit model is a cubic spline, continuous up to the second derivative, which
splits the range $[y_{\text{min}},y_{\text{max}}]$ into 15 regions, where
$y_{\text{min}}=m_{\mu}/m_{B}$ and
$y_{\text{max}}=1/2(1+m_{\mu}^{2}/m_{B}^{2})$. The exact positions of the
splits is chosen in such a way that all spectra can be fitted simultaneously.
The fit is performed by binning the decay events in $32$ bins of equal width,
covering the interval $[y_{\text{min}},y_{\text{max}}]$, and then minimizing
$\chi^{2}$ with uncertainties from the Monte Carlo integration only.

The final spectrum is the linear combination of the spectra of the four
$B$-hadron species according to the hadronization fractions given above. This
spectrum is shown in fig. 5. The error bands are estimated in the following
way, in analogy with ref. ATLAS:2022jbw . The total error is given by the
combination of the errors on the $B$-hadron hadronization fractions and those
on the branching ratios. First, the errors on the $B$-hadron hadronization
fractions are propagated. Then, these errors are rescaled by the ratio of the
total systematic error coming from $B$-hadron hadronization fractions and
branching ratios to the systematic error coming from $B$-hadron hadronization
fractions only, as can be extracted from fig. 9 in ref. ATLAS:2022jbw . This
method gives a reasonable estimate of the size of the uncertainties, keeping
in mind the obtained errors are so small compared to the errors on the NPFF,
that they can be safely neglected even if they are underestimated by a factor
of 3.

Figure 5: Left: the muon decay fragmentation function for the high-energy
B-hadron mixture. Right: the corresponding uncertainty plot. The error bands
are obtained as described in sec. 2.5.1.

####  2.5.2 $J/\psi$ spectra

We follow a strategy very similar to the one for decays to a muon. There are,
however, a few differences. First, since it is kinematically impossible for
any of the $B$-hadrons considered here to decay to more than one $J/\psi$ and
because the $J/\psi$ is charge-neutral, there are only $24$ different spectra
of interest. Only the combination of prompt and non-prompt $J/\psi$’s will be
studied here. Note that in EvtGen 2.1.1 the decay $B^{0}\to J/\psi\;K^{0}$ is
included twice: once as just $B^{0}\to J/\psi\;K^{0}$ and once distinguishing
between $B^{0}\to J/\psi\;K_{S}^{0}$ and $B^{0}\to J/\psi\;K_{L}^{0}$. We have
made sure to avoid double counting, which otherwise could lead to a visible
difference in the spectrum.

Second, each one of the $B$-hadrons considered here has at least one two-body
decay to a $J/\psi$, leading to sharp peaks in the rest-frame energy spectrum.
Some of these peaks appear delta-function-like (due to decays involving
particles with very small decay widths) while others, despite being
resolvable, are hard to fit with a pure spline. Since the spectra of interest
are the boosted ones, we have fitted such decays directly in the boosted frame
using the fact that boosting turns delta-distributions into step functions and
sharp Breit-Wigner peaks into smooth steps. The latter can be easily fitted
with a cubic spline, whereas the former have been included analytically with
the same masses and branching ratios as in EvtGen.

The fit model is again a cubic spline, continuous up to the second derivative.
However, the fitting of boosted spectra requires extra effort, since they
satisfy the following symmetry relation

$D(x)=D\bigg{(}\frac{y_{\text{min}}^{2}}{x}\bigg{)}\;,$ (36)

which is taken into account by covering the region $[y_{\text{min}},1]$ with
the spline and mapping it to $[y_{\text{min}}^{2},y_{\text{min}}]$ using eq.
(36). Additionally, a much larger number of bins – 256 – was used in order to
resolve the smooth steps. The spline splits the range $[y_{\text{min}},1]$
into $25$ regions, the exact positions of the splits having been manually
chosen to ensure a good fit for all spectra simultaneously. The exception is
the $\overline{\Lambda}_{b}^{0}$ spectrum, which differs considerably from the
others in its shape. It was fitted using a spline with $15$ regions.

The spectra of all four $B$-hadron species and the $B$-hadron mixture are
shown in fig. 6. There is a sizable ($\sim 8\%$) uncertainty on the
normalization but the uncertainty on the shape is again small enough and can
be neglected.

Figure 6: The $J/\psi$ decay fragmentation functions for various $B$ decays
(left) and corresponding uncertainty plot (right).

##  3 Application: NNLO-accurate $B,J/\psi$ and soft-$\mu$ observables in
$t\bar{t}$ events

Top quark decays at the LHC are a major source of $b$ quarks. For this reason
$t\bar{t}$ events at the LHC allow one to investigate $b$-fragmentation in a
high-purity, high-statistics environment. Additionally, a number of
$b$-fragmentation observables are sensitive to the value of the top quark mass
and, since long ago, have been advocated in the context of precision
determination of $m_{t}$. With such applications in mind, we apply the NPFFs
extracted in the present work to $t\bar{t}$ production and decay at the LHC
and make NNLO accurate predictions for a number of observables involving an
identified $B$-hadron, $J/\psi$ or a soft muon.

Our setup is as follows. We consider the production and decay of a top-quark
pair at the LHC at 13 TeV in the Narrow Width Approximation

$pp\to t\bar{t}+X\to\ell\bar{\ell}\nu\bar{\nu}b\bar{b}+X\;.$ (37)

We consider inclusive final states with a single identified particle $F$, with
$F$ being a $B$, a $J/\psi$ or a $\mu$. Unless stated otherwise, the top-quark
mass is set to $m_{t}=172.5$ GeV. We utilize the NNPDF3.1 NNLO pdf set
Ball:2017nwa  and the following two combinations for the renormalization
($\mu_{R}$), factorization ($\mu_{F}$) and fragmentation ($\mu_{Fr}$) scales

$\displaystyle\mu_{R}=\mu_{F}=\mu_{Fr}=m_{t}/2\;,$ (38)
$\displaystyle\mu_{R}=\mu_{F}=H_{T}/4~{},~{}\mu_{Fr}=m_{t}/2\;.$ (39)

The baseline scale variation for the fixed scale choice eq. (38) is the
15-point scale variation where all three scales are independently varied by a
factor of 2 around their central values, subject to the constraints

$\displaystyle 1/2\leq\mu_{R}/\mu_{F}\leq 2\,,$ $\displaystyle
1/2\leq\mu_{R}/\mu_{Fr}\leq 2\,,$ (40) $\displaystyle
1/2\leq\mu_{F}/\mu_{Fr}\leq 2\,.$

For the dynamic scale choice eq. (39) we use the standard 7-point variation of
$\mu_{R}$ and $\mu_{F}$ with fixed $\mu_{Fr}=m_{t}/2$.

As in ref. Czakon:2021ohs , the electroweak parameters are defined within the
$G_{F}$ scheme

$\displaystyle m_{W}$ $\displaystyle=$ $\displaystyle 80.385\;\text{GeV}\,,$
$\displaystyle\Gamma_{W}$ $\displaystyle=$ $\displaystyle
2.0928\;\text{GeV}\,,$ $\displaystyle m_{Z}$ $\displaystyle=$ $\displaystyle
91.1876\;\text{GeV}\,,$ $\displaystyle\Gamma_{Z}$ $\displaystyle=$
$\displaystyle 2.4952\;\text{GeV}\,,$ $\displaystyle G_{F}$ $\displaystyle=$
$\displaystyle 1.166379\cdot 10^{-5}\;\text{GeV}^{-2}\,,$
$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}G_{F}}{\pi}m_{W}^{2}\left(1-(m_{W}/m_{Z})^{2}\right)\,.$
(41)

The top quark width at leading order reads

$\Gamma_{t}^{(0)}=G_{F}\frac{m_{t}^{3}}{8\pi\sqrt{2}}(1-\xi)^{2}(1+2\xi)=1.48063\;\text{GeV}~{}({\rm
for}~{}m_{t}=172.5\;\text{GeV})\,,$ (42)

where $\xi=(m_{W}/m_{t})^{2}$.

###  3.1 Differential distributions of a fragmenting $B$, $J/\psi$ or $\mu$

We consider the following event selection requirements

  * •

$p_{T}(\ell)>25$ GeV, $|\eta(\ell)|<2.5$ ,

  * •

at least 2 anti-kT jets ($R=0.4$) with $p_{T}(j)>25$ GeV and $|\eta(j)|<2.5$ ,

  * •

$\Delta R(\ell,j)>0.4$ ,

  * •

$p_{T}(F)>8$ GeV and $|\eta(F)|<2.5$, $F$ must be part of one jet .

The following differential observables have been considered:
$p_{T}(F),|\eta(F)|,m(F\ell),E(F)$ and $p_{T}(F)/p_{T}(J_{F})$ where $J_{F}$
denotes the jet containing $F$. These observables are computed for all values
of $F=\\{B,J/\psi,\mu\\}$.

Figure 7: The absolute $E(F)$ (left column) and pseudorapidity $|\eta(F)|$
(right column) distributions for $F=B$ (top), $F=J/\psi$ (middle) and $F=\mu$
(bottom). The prediction uses the fixed scale (38). Shown are the 15 point
scale variation bands for LO, NLO and NNLO as well as the NPFF r.m.s.
uncertainty band (in yellow, shown with respect to the NLO).

Figure 8: As in fig. 7 but for $m(F\ell^{-})$ (left column) and $m(F\ell)_{\rm
min}$ (right column).

First we consider the $E(F)$ distribution which can be found in fig. 7. This
distribution is of particular interest in the context of $m_{t}$ determination
and we discuss it extensively in sec. 3.2. By comparing the $E(B)$
distribution presented in this section with the one in sec. 3.2 one can
immediately note the potentially large impact of selection cuts. The
calculation of sec. 3.2 is very inclusive, with minimal selection cuts, unlike
the observable computed in this section which is strongly impacted by cuts. As
can be seen in fig. 7 the region of low $E(F)$ is not reliably described and
higher order corrections are large. However, once $E(F)$ increases
sufficiently (above about 50 GeV for $E(B)$, 30 GeV for $E(J/\psi)$ and 10 GeV
for $E(\mu)$) the behavior of higher order corrections changes dramatically.
The perturbative convergence becomes evident and the improvement due to the
inclusion of NNLO corrections leads to a decrease of the theoretical
uncertainty from scale variation by a factor of about 4 relative to NLO.
Furthermore, in this region, the size of NNLO corrections is only few percent
above NLO further illustrating the reliability of the theory prediction once
NNLO QCD corrections have been included. We also note that the FF uncertainty
is tiny relative to the scale uncertainty even at NNLO.

The pseudorapidity distributions $|\eta(F)|$ for the three final states are
shown in fig. 7. As it may be expected, all three are very well behaved in the
full kinematic range and show no noticeable shape difference between NLO and
NNLO. The NNLO uncertainty bands are much smaller than the NLO ones and the
NNLO/NLO $K$-factor is very small. We conclude that this observable is very
precisely and reliably predicted at NNLO in QCD.

We next turn our attention to $m(F\ell)$, shown in fig. 8, which is another
kinematic variable studied in the context of $m_{t}$ in sec. 3.2. Specifically
we show $m(F\ell^{-})$ and $m(F\ell)_{\rm min}$, the latter being defined by
the requirement that $\ell$ is chosen in such a way that its invariant mass
with $F$ is minimized. We observe a pattern of higher order corrections
roughly in line with $E(F)$ discussed above. The most prominent feature is the
non-overlap of NLO and NNLO uncertainty bands for $m(B\ell)$ below, roughly,
50 GeV. This feature is likely driven by the selection cuts since, as can be
seen in fig. 11, once more inclusive selection cuts are applied the normalized
distribution does not exhibit significant shape change (in this region)
between NLO and NNLO. Outside of this range, the NNLO prediction for
$m(F\ell)$ is consistent with the NLO one. The $m(B\ell)$ and $m(J/\psi\ell)$
distributions at NNLO are affected by a shape difference, while the
$m(\mu\ell)$ one is not, within the NNLO scale uncertainty. In general, the
NNLO uncertainty band is much smaller than the NLO one, especially for
$m(\mu\ell)$, which demonstrates the important role of NNLO QCD corrections in
the precision description of this observable.

Figure 9: As in fig. 7 but for $p_{T}(F)$ (left column) and
$p_{T}(F)/p_{T}(J_{F})$ (right column).

In fig. 9 we show the predictions for $p_{T}(F)$. We observe that for all
final states the NNLO uncertainty band is much smaller than the NLO one. The
pattern of higher order corrections is, however, different for the different
final states. Excluding the region of small $p_{T}(B)$ which is dominated by
missing $m_{b}$ corrections and selection cuts, the $p_{T}(B)$ spectrum shows
remarkable stability and lack of shape correction up to very high $p_{T}$. The
lighter final states, $J/\psi$ and especially $\mu$, on the other hand show
consistent change in shape between NLO and NNLO within the same $p_{T}$
interval. This difference between $B$, $J/\psi$ and $\mu$ is consistent with
the change in the mean $z$ of the FFs for the three final states. The NLO and
NNLO bands overlap in the full kinematic range considered here.

Finally, in fig. 9 we show the predictions for $p_{T}(F)/p_{T}(J_{F})$ for all
three final states. The importance of this observable is related to the fact
it is sensitive to the FF, and that it can be used to extract NPFF from LHC
data. It has already been extensively discussed in ref. Czakon:2021ohs  and we
refer the reader to that reference for further details.

Before concluding this section we would like to make a comment about the size
of the MC uncertainty of the NNLO predictions. Overall it is small;
appreciable fluctuations can only be observed in the $E(F)$ distribution and
in the tails of several $\mu$ distributions. As it turns out, in order to
achieve high quality predictions for all three final states
$F=\\{B,J/\psi,\mu\\}$ considered here, one needs to use dedicated
calculations with dedicated optimization for each final state. The reason for
this is that the fragmentation functions for each final state is different and
as a result different kinematical ranges are sampled.

###  3.2 Top quark mass extraction

In the following we consider the observables $E(F)$ and $m(B\ell)$ due to
their sensitivity to $m_{t}$. We also compute the first moment of $m(B\ell)$
as a function of $m_{t}$. This is the first time this important observable has
been derived with NNLO precision. All predictions shown in this section are
subject to the following event selection requirements

  * •

$p_{T}(\ell)>25$ GeV, $|\eta(\ell)|<2.5$ ,

  * •

$p_{T}(F)>8$ GeV and $|\eta(F)|<2.5$ .

Unlike the selection requirements in sec. 3.1, here we do not impose any jet
selection requirements.

We first consider the $E(B)$ distribution. It is shown in fig. 10 where we
compare predictions with LO, NLO and NNLO accuracy. For reasons to be
explained shortly, besides the predictions with full top quark decay (fig. 10,
left), we also present predictions with top decay kept at LO (fig. 10, right).

Figure 10: The normalized $E(B)$ distribution for fixed scale (38). Shown are
the 15 point scale variation bands for LO, NLO and NNLO as well as the NPFF
r.m.s. uncertainty band. The plot to the right is like the one to the left but
with LO top decay.

The $E(B)$ distribution has been proposed Agashe:2012bn ; Agashe:2016bok  as a
means of measuring the top quark mass relatively free of contamination due to
possible physics beyond the Standard Model (BSM). The main idea is that the
position of the maximum of the $E(B)$ distribution is independent of the
$t\bar{t}$ dynamics. In particular, it is not affected by potential BSM
contributions in $t\bar{t}$ production. Therefore, if one calculates the
dependence of the position of the maximum as a function of $m_{t}$, one can
infer the value of $m_{t}$ from a direct measurement of $E_{\rm max}(B)$. In
practice, for values of $m_{t}$ close to the world average, it is sufficient
to consider a linear fit

$E_{\rm max}(m_{t})=am_{t}+b\,.$ (43)

Our predictions for the coefficients $(a,b)$ at LO, NLO and NNLO can be found
in table 2 (with full top quark decay) and in table 3 (for LO top quark
decay). The uncertainties on the fit parameters are explained in the
following.

An important feature of the method of refs. Agashe:2012bn ; Agashe:2016bok  is
that the top quark decay is a two-body one, i.e. LO-like. Once one includes
higher order QCD corrections the top decay is modified by an additional
radiation which alters the top decay kinematics. This radiation affects the
position of the peak of the $E(B)$ distribution and needs to be accounted for
in any realistic application of this method. The NLO QCD corrections were
first included in ref. Agashe:2016bok  where a significant shift in the peak
position was observed. In this work we also extend, for the first time, this
observable to NNLO in QCD and observe that the NNLO QCD corrections lead to a
further shift in the position of $E_{\rm max}(B)$. As can be seen in table 2,
the inclusion of NLO corrections in top decay leads to a roughly 4 GeV shift
in $E_{\rm max}(B)$ while the NNLO ones lead to a further 2 GeV shift.
Clearly, the inclusion of the NNLO QCD corrections is mandatory for a reliable
description of this observable.

As a calibration of the method, in table 3 we have shown the same predictions
but where the top quark decay has consistently been kept at LO. As can be seen
there, within the uncertainty of the calculation, the position of the maximum
remains stable through NNLO, as expected.

Finally, we would like to mention that while we only show predictions for the
final state with a $B$-hadron, similar predictions can also be made for final
states with a $J/\psi$ or a muon.

$m_{t}$ | LO | NLO | NNLO
---|---|---|---
$171.5$ GeV | $37.553~{}(\pm 0.106)~{}(^{+0.050}_{-0.061})$ | $40.994~{}(\pm 0.147)~{}(^{+1.178}_{-0.710})$ | $42.957~{}(\pm 0.329)~{}(^{+1.087}_{-0.818})$
$172.5$ GeV | $37.816~{}(\pm 0.109)~{}(^{+0.051}_{-0.062})$ | $41.277~{}(\pm 0.158)~{}(^{+1.196}_{-0.717})$ | $43.263~{}(\pm 0.332)~{}(^{+1.073}_{-0.825})$
$173.5$ GeV | $38.093~{}(\pm 0.113)~{}(^{+0.051}_{-0.061})$ | $41.657~{}(\pm 0.168)~{}(^{+1.250}_{-0.745})$ | $43.528~{}(\pm 0.222)~{}(^{+1.010}_{-0.778})$
Lin. fit | LO | NLO | NNLO
$a=$ | $0.270~{}(\pm 0.004)$ | $0.329~{}(\pm 0.028)$ | $0.284~{}(\pm 0.011)$
$b=$ |  $-8.755~{}(\pm 0.708)$ GeV |  $-15.429~{}(\pm 4.820)$ GeV |  $-5.666~{}(\pm 1.816)$ GeV
Table 2: Values of $E_{\rm max}(B)$ for the absolute differential cross section with full top quark decay at LO, NLO and NNLO and for three different values of $m_{t}$. Positions are fit using eq. (45) and 5 GeV bins. Also given are the parameters of the linear fit eq. (43) at LO, NLO and NNLO. $m_{t}$ | LO | NLO | NNLO
---|---|---|---
$171.5$ GeV | $37.553~{}(\pm 0.106)~{}(^{+0.050}_{-0.061})$ | $36.744~{}(\pm 0.169)~{}(^{+0.213}_{-0.313})$ | $36.737~{}(\pm 0.311)~{}(^{+0.081}_{-0.021})$
$172.5$ GeV | $37.816~{}(\pm 0.109)~{}(^{+0.051}_{-0.062})$ | $36.981~{}(\pm 0.182)~{}(^{+0.223}_{-0.330})$ | $37.010~{}(\pm 0.227)~{}(^{+0.109}_{-0.019})$
$173.5$ GeV | $38.093~{}(\pm 0.113)~{}(^{+0.051}_{-0.061})$ | $37.319~{}(\pm 0.193)~{}(^{+0.206}_{-0.296})$ | $37.292~{}(\pm 0.255)~{}(^{+0.113}_{-0.056})$
Lin. fit | LO | NLO | NNLO
$a=$ | $0.270~{}(\pm 0.004)$ | $0.286~{}(\pm 0.029)$ | $0.278~{}(\pm 0.003)$
$b=$ |  $-8.755~{}(\pm 0.708)$ GeV |  $-12.237~{}(\pm 4.962)$ GeV |  $-10.913~{}(\pm 0.556)$ GeV
Table 3: As in table 2 but for LO top quark decay.

Before we conclude our discussion of the $E(B)$ distribution we would like to
explain how the position of the maximum $E_{\rm max}(B)$ has been determined.
This is a non-trivial task since our calculation produces not the differential
spectrum but its binned version, i.e. the position of the maximum is obscured
by finite bin size effects. To this end we extract the position of the maximum
with the help of a fit

$f(E)=n\exp\left(-w\left(\frac{\hat{E}}{E}+\frac{E}{\hat{E}}\right)\right)\,,$
(44)

performed at each perturbative order and for each scale combination. The fit
parameters are $n,w$ and $\hat{E}$, the latter giving the curve’s maximum
position. The functional form in eq. (44) is the same one used in ref.
Agashe:2016bok . We fit the range $E(B)\in[25,100]$ GeV using 5 GeV wide bins.
The $\chi^{2}$ for the fit (44) is computed as

$\chi^{2}=\sum_{i\in\text{Bins}}\left(f(E(B)_{i})-\frac{\differential\sigma}{\differential
E}\Big{|}_{E=E(B)_{i}}\right)^{2}\left(\frac{\differential\sigma^{\rm
MC\,err}}{\differential E}\Big{|}_{E=E(B)_{i}}\right)^{-2}\,,$ (45)

where $E(B)_{i}$ is the middle of each $E(B)$ bin and $\sigma^{\rm MC\,err}$
indicates that only the MC uncertainty of the calculation in each bin is used.
As usual, the uncertainty on the fit parameters is given by the interval
defined by $\Delta\chi^{2}=1$. The uncertainties in tables 2 and 3 are derived
as follows. The first one refers to the uncertainty due to the parameter
$\hat{E}$, while the second bracket gives an estimate of the scale uncertainty
which is derived by repeating the fit for each scale choice and taking the
envelope of all 15 values. The uncertainties of the $a$ and $b$ parameters are
the ones from the linear fit to the central scale values. We have checked that
fits to absolute or normalized distributions produce the same peak position,
as expected.

We next consider the $m(B\ell^{-})$ distribution. We have checked that the
$m(B\ell^{+})$ distribution is equivalent to it within the MC error of the
calculation. The prediction for the normalized $m(B\ell^{-})$ distribution can
be found in fig. 11. We only show the predictions for final states with a
$B$-hadron although predictions for a $J/\psi$ or a muon can be provided. We
observe the much improved perturbative description at higher orders. We also
note that 1 GeV shift in $m_{t}$ has almost identical effect on the normalized
distribution as the inclusion of the NNLO correction relative to the NLO one.
This demonstrates the tight interplay between the inclusion of higher order
corrections to the so-called indirect top quark mass observables and the
ultimate precision of the extracted $m_{t}$.

Figure 11: (Left): the normalized $m(B\ell^{-})$ distribution for fixed scale
(38). Shown are the 15 point scale variation bands for LO, NLO and NNLO as
well as the NPFF r.m.s. uncertainty band. (Right): the NNLO prediction for 3
different values of $m_{t}$.

A feature of the observable $m(B\ell^{-})$ is that it combines $(B,\ell)$
pairs that may originate from either the same top quark or from two different
top quarks. For the purpose of $m_{t}$ determination, one would like to
predominantly have pairs that originate from the same top quark. One way to
achieve this is to consider the modified observable $m(B\ell)_{\rm min}$ where
the $B$ and $\ell$ are paired not based on lepton’s charge but on the
requirement that their invariant mass is minimized. In fig. 12 (left) we show
the predicted $m_{t}$ dependence for the first moment of the normalized
$m(B\ell)_{\rm min}$ distribution. Note that this moment is sensitive to the
selection requirements listed in the beginning of this section. We have
verified that for different selection requirements, for example the ones
listed in sec. 3.1, the behavior of this observable can change significantly.

We find the NNLO correction shifts the dependence relative to NLO, but much
less than the NLO one does relative to LO. The perturbative uncertainty also
decreases at higher orders which means that the inclusion of the NNLO
correction will have an important impact on this observable. We also note that
at NNLO the slope of the $m_{t}$ dependence very slightly decreases which
means a slight decrease with higher orders to the $m_{t}$ sensitivity of this
observable. This effect is nonetheless small relative to the overall theory
uncertainty. While, as expected, the behavior of our prediction is broadly in
line with the one in ref. Biswas:2010sa , a direct comparison at LO and NLO
with that reference is not possible due to the differences between the NPFFs
used and the event selections.

Figure 12: The first moment of the normalized $m(F\ell)_{\rm min}$
distribution for $F=B$ (upper left), $F=J/\psi$ (upper right) and $F=\mu$
(lower left). Shown are the LO, NLO and NNLO predictions as a function of
$m_{t}$ (only the LO and NLO ones for $F=J/\psi,\mu$). Bands at each order
denote 15-point scale variation for the fixed scale (38). The three NLO
curves, divided by their respective central scale values at $m_{t}=172.5$ GeV,
are shown together in the lower-right panel.

There is a longstanding interest Kharchilava:1999yj ; CDF:2009mbf ;
CMS:2016ixg ; ATLAS:2022jbw  in the $m(F\ell)$ distribution with final states
other than $B$, especially $J/\psi$ or a muon. Utilizing the fragmentation
functions to these final states derived in the present work, in the upper-
right and lower-left panels of fig. 12 we show the $m_{t}$ dependence of
$\langle m(J/\psi\ell)_{\rm min}\rangle$ and $\langle m(\mu\ell)_{\rm
min}\rangle$, respectively. Due to the high computational cost we have not
computed the NNLO corrections to these observables. The LO and NLO corrections
to these final states, especially when comparing with the pattern of higher
order corrections for $\langle m(B\ell)_{\rm min}\rangle$, allow one to draw
conclusions about the behavior of these observables. In particular, we see
that through NLO $\langle m(J/\psi\ell)_{\rm min}\rangle$ has a pattern of
higher order corrections very similar to $\langle m(B\ell)_{\rm min}\rangle$.
We expect this trend to continue at NNLO. The $K$-factor for $\langle
m(\mu\ell)_{\rm min}\rangle$ is slightly larger relative to its reduced scale
uncertainties, since LO and NLO uncertainty bands barely overlap. This is due
to the different shape of the muon FF which probes different kinematics
relative to $B$ or $J/\psi$.

Another feature of the three functions $\langle m(F\ell)_{\rm
min}\rangle(m_{t})$ is the change in slope around $m_{t}=172.5$ GeV as the
final state changes. The slope is largest for $F=B$ and decreases towards
$F=J/\psi$ and especially for $F=\mu$. However, the relative change around
$m_{t}=172.5$ GeV given by $\langle m(F\ell)_{\rm min}\rangle(m_{t})/\langle
m(F\ell)_{\rm min}\rangle(172.5\text{ GeV})$, shown at NLO in the lower-right
panel of fig. 12, is the same for $F=B$ and $F=J/\psi$ and only 20% smaller
for $F=\mu$. In fig. 12, the denominator is always evaluated at the central
scale choice in order to preserve the information on the scale uncertainties.
The theoretical uncertainty on a hypothetically extracted value of $m_{t}$ is
determined by the ratio of the width of the scale uncertainty band and the
slope. From fig. 12 one can conclude that the theory uncertainty is slightly
smaller for $F=J/\psi$ than for $F=B$ and almost a factor 2 smaller for
$F=\mu$ than for $F=B$. Additionally, a corresponding measurement for $F=\mu$
would most likely also be the most precise, followed by $F=J/\psi$ and then
$F=B$.

We have also investigated the sensitivity of this important observable to the
choice of scale. To that end in fig. 13 we have compared predictions with the
fixed scale (38) (left) and with the dynamic one (39) (right). Since the
calculation with dynamic scale eq. (39) is performed with a 7-point scale
variation (for fixed scale $\mu_{Fr}$), for a proper comparison of the two
scales in fig. 13 we show the prediction with the dynamic scale (right) and
the one with the fixed scale (38) but with 7-point scale variation (left).
Comparing fig. 13 (left) with fig. 13 (right) we conclude that the two scales
lead to very similar predictions that are fully consistent within
uncertainties, especially at NNLO.

Figure 13: The first moment of the normalized $m(B\ell)_{\rm min}$
distribution computed with two different scales: fixed scale (38) (left) and
the dynamic scale (39) (right). In both cases bands correspond to 7-point
scale variation with fixed $\mu_{Fr}$.

##  4 Conclusions

Recently, ref. Czakon:2021ohs  presented a novel technical approach for the
numerical calculation of NNLO QCD corrections to coefficient functions for
identified hadron production at hadron and/or lepton colliders. The approach
is applicable to both heavy and light hadron production. It utilizes numerical
Monte Carlo integration techniques within the sector-improved residue
subtraction scheme Czakon:2010td ; Czakon:2014oma  which implies that generic
NNLO observables with complicated phase space cuts can be implemented. This
sets the technique of ref. Czakon:2021ohs  apart from the so far pursued,
decades-old approach based on analytical coefficient functions applicable only
to a small number of inclusive observables.

As a first application of this computational technique, ref. Czakon:2021ohs
presented the NNLO QCD corrections to $B$-hadron production in $t\bar{t}$
events at the LHC. This state of the art calculation was supplemented by an
existing non-perturbative fragmentation function for $B$-hadrons. The
significant improvement in the precision of theoretical predictions due to the
inclusion of NNLO QCD corrections, however, is calling for a similar scrutiny
at the precision and quality of the non-perturbative fragmentation
contribution for heavy flavor production. We observe that, at this level of
precision, the NPFF affects the ultimate precision that can be achieved in
theoretical predictions of single heavy identified hadrons.

The main goal of the present paper is to address this gap. It provides a NNLO
NPFF for $B$ hadrons which is extracted from $e^{+}e^{-}$ data in a way which
is consistent with the formalism used by us for state-of-the-art NNLO
predictions for LHC observables.

As a first application of the newly extracted $B$ fragmentation function, in
this work we present NNLO accurate predictions for $B$-production in
$t\bar{t}$ events at the LHC. Such application is of strong immediate interest
due to the large LHC statistics on $t\bar{t}$ events. We observe a significant
improvement from the inclusion of NNLO QCD corrections for a range of generic
observables (by a factor of about four relative to NLO QCD). We have placed a
particular emphasis on predictions for observables that allow precision
determination of the top quark mass. Future applications of the fragmentation
function extracted in this work, and of the computational formalism developed
in ref. Czakon:2021ohs , include open $B$ production at the LHC as well as $B$
production in association with bosons. Our fragmentation function fitting
procedure can also be adapted to charm production which is another application
of strong interest.

The extraction of the $B$ NPFF is based on a global fit of four $e^{+}e^{-}$
measurements, although we have studied its sensitivity to any single dataset.
We have produced a set of 100 NPFFs which allows a detailed FF uncertainty
analysis in any observable of interest. Our results are publicly available as
LHAPDF grids.

The extraction of the $B$ NPFF is based on NNLO perturbative inputs with soft-
gluon resummation through NNLL. The effect of the resummation in the
perturbative initial condition for the $B$ fragmentation function is
important. In analyzing the soft-gluon resummation of the FF initial
condition, we have found complete consistency between its NNLO fixed order
calculation, the anomalous dimensions governing that resummation and the
decoupling through NNLO of the $b$-quark fragmentation function across the
heavy flavor threshold.

Another novelty of the present work is that it also derives sets of FF for
single inclusive production of $J/\psi$ or soft muons originating in $B$
hadron decays. To that end we have convoluted the $B$-hadron NPFF extracted by
us with precisely measured $B$ meson decays spectra. These two sets of NPFFs
are, too, made public with this publication. They allow NNLO accurate
predictions for observables based not only on $B$-hadrons but on their decay
products. Such observables are measured very accurately at the LHC and are
needed, among others, in measurements of $m_{t}$ with high precision or for
prediction of $J/\psi$ and muons in semileptonic $B$ decays in the context of
open $B$ production or in association with bosons.

The NPFFs for all final states considered in the present work are available as
LHAPDF grids through this paper’s electronic submission on arXiv.

###### Acknowledgements.

We would like to thank Gennaro Corcella for a collaboration at an early stage
of this work and for numerous useful comments and suggestions. We would also
like to thank Valerio Bertone for a helpful discussion. The work of M.C. was
supported by the Deutsche Forschungsgemeinschaft under grant 396021762 – TRR
257. The work of T.G. was supported by the DFG under grant 400140256 - GRK
2497: The physics of the heaviest particles at the Large Hadron Collider. The
research of A.M. and R.P. has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 Research and Innovation
Programme (grant agreement no. 683211). A.M. was also supported by the UK STFC
grants ST/L002760/1 and ST/K004883/1. R.P. acknowledges support from the
Leverhulme Trust and the Isaac Newton Trust.

##  Appendix A Results for soft-gluon resummation

In this appendix we collect various results used in the soft-gluon resummation
of the ${e^{+}e^{-}}$ coefficient functions and the initial condition of the
perturbative fragmentation function. Please refer to secs. 2.1 and 2.2 for
definitions.

The relevant coefficients of the QCD $\beta$-function in a theory with $n_{f}$
active flavors read

$\displaystyle b_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\bigg{[}\frac{11}{3}C_{A}-\frac{4}{3}T_{F}n_{f}\bigg{]}\,,$
$\displaystyle b_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{(4\pi)^{2}}\bigg{[}\frac{34}{3}C_{A}^{2}-\frac{20}{3}C_{A}T_{F}n_{f}-4C_{F}T_{F}n_{f}\bigg{]}\,,$
$\displaystyle b_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{(4\pi)^{3}}\bigg{[}-\frac{1415}{27}C_{A}^{2}T_{F}n_{f}+\frac{2857}{54}C_{A}^{3}-\frac{205}{9}C_{A}C_{F}T_{F}n_{f}+2C_{F}^{2}T_{F}n_{f}$
(46)
$\displaystyle+\frac{158}{27}C_{A}T_{F}^{2}n_{f}^{2}+\frac{44}{9}C_{F}T_{F}^{2}n_{f}^{2}\bigg{]}\,,$

where $C_{F}=4/3,C_{A}=3$ and $T_{F}=1/2$ are the usual QCD color factors.

###  A.1 $D^{\rm ini}_{b\to b}$

The cusp anomalous dimension (12) in a theory with $n_{f}$ active flavors
reads though NNLL

$\displaystyle A_{1}$ $\displaystyle=$ $\displaystyle C_{F}\,,$ $\displaystyle
A_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}C_{F}\bigg{[}\bigg{(}\frac{67}{18}-\frac{\pi^{2}}{6}\bigg{)}C_{A}-\frac{10}{9}T_{F}n_{f}\bigg{]}\,,$
$\displaystyle A_{3}$ $\displaystyle=$ $\displaystyle
C_{F}\bigg{[}C_{A}^{2}\bigg{(}\frac{11\zeta_{3}}{24}+\frac{245}{96}-\frac{67\pi^{2}}{216}+\frac{11\pi^{4}}{720}\bigg{)}-C_{A}T_{F}n_{f}\bigg{(}\frac{7\zeta_{3}}{6}+\frac{209}{216}-\frac{5\pi^{2}}{54}\bigg{)}$
(47) $\displaystyle-
C_{F}T_{F}n_{f}\bigg{(}\frac{55}{48}-\zeta_{3}\bigg{)}-\frac{1}{27}T_{F}^{2}n_{f}^{2}\bigg{]}\,.$

The NNLL expansion of the wide-angle anomalous dimension $D$ in eq. (13) reads

$\displaystyle D_{1}$ $\displaystyle=$ $\displaystyle-C_{F}\,,$ $\displaystyle
D_{2}$ $\displaystyle=$ $\displaystyle
C_{F}\bigg{[}C_{A}\bigg{(}-\frac{9\zeta_{3}}{4}+\frac{55}{108}+\frac{\pi^{2}}{12}\bigg{)}+\frac{1}{27}T_{F}n_{l}\bigg{]}\;.$
(48)

The hard matching coefficient (10), expanded through NNLO, reads

$\displaystyle H_{0}^{\rm ini}$ $\displaystyle=$ $\displaystyle 1\,,$
$\displaystyle H_{1}^{\rm ini}$ $\displaystyle=$ $\displaystyle
C_{F}\bigg{[}\bigg{(}\frac{3}{2}-2\gamma\bigg{)}\ln\bigg(\frac{\mu_{Fr,0}^{2}}{m_{b}^{2}}\bigg{missing})-\frac{\pi^{2}}{3}-2\gamma^{2}+2\gamma+2\bigg{]}\,,$
$\displaystyle H_{2}^{\rm ini}$ $\displaystyle=$ $\displaystyle
C_{F}\bigg{\\{}C_{F}\bigg{[}\frac{1}{8}(3-4\gamma)^{2}\ln^{2}\bigg{(}\frac{\mu_{Fr,0}^{2}}{m_{b}^{2}}\bigg{)}+\ln(\frac{\mu_{Fr,0}^{2}}{m_{b}^{2}})\bigg{(}6\zeta_{3}+\frac{27}{8}-7\gamma^{2}+4\gamma^{3}-\pi^{2}$
(49)
$\displaystyle+\gamma\bigg{(}\frac{2\pi^{2}}{3}-1\bigg{)}\bigg{)}-\frac{3\zeta_{3}}{2}+\frac{2}{3}\gamma^{2}(\pi^{2}-3)-\frac{2}{3}\gamma(\pi^{2}-6)-\frac{11\pi^{4}}{180}+\frac{\pi^{2}}{4}+2\gamma^{4}-4\gamma^{3}$
$\displaystyle+\frac{241}{32}-\pi^{2}\ln(4)\bigg{]}+C_{A}\bigg{[}-\frac{11}{24}(4\gamma-3)\ln^{2}\bigg{(}\frac{\mu_{Fr,0}^{2}}{m_{b}^{2}}\bigg{)}+\ln(\frac{\mu_{Fr,0}^{2}}{m_{b}^{2}})\bigg{(}\frac{11}{12}(4\gamma$
$\displaystyle-3)\ln(\frac{\mu_{Fr,0}^{2}}{\mu_{R,0}^{2}})+\frac{1}{72}(-216\zeta_{3}+315-272\gamma-264\gamma^{2}+24\gamma\pi^{2})\bigg{)}+\frac{11}{18}(-6-6\gamma$
$\displaystyle+6\gamma^{2}+\pi^{2})\ln(\frac{\mu_{Fr,0}^{2}}{\mu_{R,0}^{2}})-\frac{97\zeta_{3}}{18}+\gamma\bigg{(}9\zeta_{3}-\frac{55}{27}-\frac{14\pi^{2}}{9}\bigg{)}+\frac{1141}{288}-\frac{22\gamma^{3}}{9}+\frac{\pi^{4}}{12}$
$\displaystyle+\frac{1}{9}\gamma^{2}(3\pi^{2}-34)+\pi^{2}\bigg{(}\frac{7}{54}+\ln(2)\bigg{)}\bigg{]}+T_{F}n_{l}\bigg{[}\frac{1}{6}(4\gamma-3)\ln^{2}\bigg{(}\frac{\mu_{Fr,0}^{2}}{m_{b}^{2}}\bigg{)}$
$\displaystyle+\bigg{(}\bigg{(}1-\frac{4\gamma}{3}\bigg{)}\ln(\frac{\mu_{Fr,0}^{2}}{\mu_{R,0}^{2}})+\frac{4\gamma^{2}}{3}+\frac{8\gamma}{9}-\frac{3}{2}\bigg{)}\ln(\frac{\mu_{Fr,0}^{2}}{m_{b}^{2}})-\frac{2}{9}(-6-6\gamma+6\gamma^{2}$
$\displaystyle+\pi^{2})\ln(\frac{\mu_{Fr,0}^{2}}{\mu_{R,0}^{2}})-\frac{2\zeta_{3}}{9}-\frac{173}{72}-\frac{4\gamma}{27}+\frac{8\gamma^{2}}{9}+\frac{8\gamma^{3}}{9}-\frac{4\pi^{2}}{27}+\frac{4\gamma\pi^{2}}{9}\bigg{]}$
$\displaystyle+T_{F}\bigg{[}\frac{1}{6}(4\gamma-3)\ln^{2}\bigg{(}\frac{\mu_{Fr,0}^{2}}{m_{b}^{2}}\bigg{)}+\ln(\frac{\mu_{Fr,0}^{2}}{m_{b}^{2}})\bigg{(}\bigg{(}1-\frac{4\gamma}{3}\bigg{)}\ln(\frac{\mu_{Fr,0}^{2}}{\mu_{R,0}^{2}})+\frac{4\gamma^{2}}{3}+\frac{8\gamma}{9}$
$\displaystyle-\frac{3}{2}\bigg{)}-\frac{2}{9}(-6-6\gamma+6\gamma^{2}+\pi^{2})\ln(\frac{\mu_{Fr,0}^{2}}{\mu_{R,0}^{2}})+\frac{2\zeta_{3}}{3}-\frac{\pi^{2}}{3}-\frac{56\gamma}{27}+\frac{3139}{648}\bigg{]}\bigg{\\}}\,.$

The coefficients of the non-truncated function $G^{\rm ini}(m,N,\alpha_{S})$
read through NNLL

$\displaystyle g_{1}(\lambda)$ $\displaystyle=$ $\displaystyle-\frac{1}{2\pi
b_{0}^{2}}\bigg{\\{}2A_{1}\lambda+A_{1}(1-2\lambda)\ln(1-2\lambda)\bigg{\\}}\,,$
$\displaystyle g_{2}(\lambda)$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi^{2}b_{0}^{3}}\bigg{\\{}4\lambda\bigg{[}\pi
b_{0}\bigg{(}(A_{1}b^{\prime}_{0}-A^{\prime}_{1}b_{0})\ln(\frac{\mu_{R,0}^{2}}{m_{b}^{2}})-A^{\prime}_{1}b_{0}\ln(\frac{\mu_{Fr,0}^{2}}{\mu_{R,0}^{2}})\bigg{)}+A_{2}b_{0}-\pi
A_{1}b_{1}\bigg{]}$ $\displaystyle+2\ln(1-2\lambda)\bigg{[}\pi
A_{1}b_{0}b^{\prime}_{0}\ln(\frac{\mu_{R,0}^{2}}{m_{b}^{2}})+A_{2}b_{0}+\pi(2\gamma
A_{1}b_{0}^{2}-A_{1}b_{1}+b_{0}^{2}B_{1})\bigg{]}$ $\displaystyle-\pi
A_{1}b_{1}\ln^{2}(1-2\lambda)\bigg{\\}}\,,$ $\displaystyle g_{3}(\lambda)$
$\displaystyle=$
$\displaystyle\frac{1}{144\pi^{3}b_{0}^{4}(1-2\lambda)}\bigg{\\{}72\pi
b_{0}\lambda\bigg{[}-\pi
b_{0}b^{\prime}_{0}\ln^{2}\bigg{(}\frac{\mu_{R,0}^{2}}{m_{b}^{2}}\bigg{)}\big{(}2A_{1}b^{\prime}_{0}\lambda+A^{\prime}_{1}b_{0}(1-2\lambda)\big{)}$
(50)
$\displaystyle-2\ln(\frac{\mu_{R,0}^{2}}{m_{b}^{2}})\big{(}b_{0}(1-2\lambda)(A^{\prime}_{2}b_{0}-\pi
A_{1}b^{\prime}_{1})+b^{\prime}_{0}(2\lambda(A_{2}b_{0}-\pi A_{1}b_{1})$
$\displaystyle+\pi b_{0}^{2}(2\gamma
A_{1}+B_{1}))\big{)}+b_{0}^{2}(1-2\lambda)\ln(\frac{\mu_{Fr,0}^{2}}{\mu_{R,0}^{2}})\bigg{(}\pi
A^{\prime}_{1}b^{\prime}_{0}\ln(\frac{\mu_{Fr,0}^{2}}{\mu_{R,0}^{2}})-2A^{\prime}_{2}\bigg{)}\bigg{]}$
$\displaystyle-6\lambda\bigg{[}A_{1}\big{(}24\pi^{2}b_{2}b_{0}(1-\lambda)+b_{0}^{2}(-48\gamma\pi^{2}b_{1}-14\lambda+7)+24\pi^{2}b_{1}^{2}\lambda$
$\displaystyle+8(6\gamma^{2}\pi^{2}+\pi^{4})b_{0}^{4}\big{)}+24b_{0}\big{(}\lambda(A_{3}b_{0}-\pi
A_{2}b_{1})+\pi(2\gamma b_{0}^{2}-b_{1})(A_{2}+\pi b_{0}B_{1})$
$\displaystyle+\pi
b_{0}^{2}B_{2}\big{)}\bigg{]}+3\ln(1-2\lambda)\bigg{[}24\pi^{2}A_{1}b_{0}\ln(\frac{\mu_{R,0}^{2}}{m_{b}^{2}})(-2b_{0}b^{\prime}_{1}\lambda+2b_{1}b^{\prime}_{0}\lambda+b_{0}b^{\prime}_{1})$
$\displaystyle+A_{1}\big{(}b_{0}^{2}(48\gamma\pi^{2}b_{1}+14\lambda-7)-24\pi^{2}b_{2}b_{0}(1-2\lambda)-48\pi^{2}b_{1}^{2}\lambda\big{)}$
$\displaystyle+24\pi b_{0}b_{1}(A_{2}+\pi
b_{0}B_{1})\bigg{]}-36\pi^{2}A_{1}b_{1}^{2}\ln^{2}(1-2\lambda)\bigg{\\}}\,.$

Please recall that, as specified in the text following eq. (13), a prime in
the above equations denotes a quantity computed with $n_{f}=n_{l}+1$ while the
ones without a prime are for $n_{f}=n_{l}$.

The functions $g_{1}$ and $g_{2}$ were previously calculated in ref.
Cacciari:2001cw  for $n_{f}=n_{l}$. Substituting primed coefficients with
unprimed ones in our results above should reproduce the results of ref.
Cacciari:2001cw . This is indeed the case. The difference between the results
above and the results of ref. Cacciari:2001cw  is of
$\mathcal{O}(\alpha_{S}^{2})$, i.e. the NLO matching implemented in ref.
Cacciari:2001cw  is insensitive to it.

The function $g_{3}$ was recently derived in ref. Maltoni:2022bpy  for
$n_{f}=n_{l}$. Therefore, in order to compare our result for $g_{3}$ with that
of ref. Maltoni:2022bpy , we convert our expression to a coupling with $n_{l}$
active flavors. In practice, the primes in our result for $g_{3}$ need to be
dropped. Once this replacement is made, the two results still disagree by the
following term

$-\frac{7A_{1}}{48\pi^{3}b_{0}^{2}}(2\lambda+\ln(1-2\lambda))\,,$ (51)

which originates in the matching of $\alpha_{S}$. After subtracting this term
from our result for $g_{3}$, it then exactly matches the $g_{3}$ of ref.
Maltoni:2022bpy .

The result for $H^{\rm ini}$ is in full agreement with those of refs.
Cacciari:2001cw ; Maltoni:2022bpy  after correcting a trivial typo for the
incorrect overall sign of $H_{1}^{\rm ini}$ in ref. Maltoni:2022bpy .

Finally, the term $\sim
C_{F}T_{F}\bigg{(}-\frac{2}{3}\ln^{2}\bigg{(}\frac{m_{b}^{2}}{\mu_{\text{th}}^{2}}\bigg{)}-\frac{20}{9}\ln\bigg(\frac{m_{b}^{2}}{\mu_{\text{th}}^{2}}\bigg{missing})-\frac{56}{27}\bigg{)}$
which appears in the decoupling constant $C_{\rm dec}$ (15) for the
fragmentation function has already been discussed in ref. Maltoni:2022bpy  but
has been given a different interpretation there.

###  A.2 ${e^{+}e^{-}}$ coefficient functions

The results for the cusp anomalous dimension $A_{i}$ can be found in appendix
A.1. The anomalous dimension $B$ reads

$B(a)=\sum_{i=1}\left(\frac{a}{\pi}\right)^{i}B_{i}\,,$ (52)

with coefficients

$\displaystyle B_{1}$ $\displaystyle=$ $\displaystyle-\frac{3}{2}C_{F}\,,$
(53) $\displaystyle B_{2}$ $\displaystyle=$ $\displaystyle
C_{F}\bigg{(}C_{A}\bigg{(}5\zeta_{3}-\frac{3155}{432}+\frac{11\pi^{2}}{36}\bigg{)}-C_{F}\bigg{(}3\zeta_{3}+\frac{3}{16}-\frac{\pi^{2}}{4}\bigg{)}+\bigg{(}\frac{247}{108}-\frac{\pi^{2}}{9}\bigg{)}T_{F}n_{f}\bigg{)}\,.$

The resummation function $G$ defined in eq. (5), with expansion as in eq.
(16), reads

$\displaystyle g_{1}(\lambda)$ $\displaystyle=$ $\displaystyle\frac{1}{\pi
b_{0}^{2}}\bigg{\\{}A_{1}(\lambda+(1-\lambda)\ln(1-\lambda))\bigg{\\}}\,,$
$\displaystyle g_{2}(\lambda)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{2}b_{0}^{3}}\bigg{\\{}2\lambda\bigg{[}\pi
A_{1}b_{0}^{2}\ln(\frac{\mu_{F}^{2}}{\mu_{R}^{2}})-A_{2}b_{0}+\pi
A_{1}b_{1}\bigg{]}+\ln(1-\lambda)\bigg{[}2\pi
A_{1}b_{0}^{2}\ln(\frac{Q^{2}}{\mu_{R}^{2}})$ $\displaystyle+\pi(-2\gamma
A_{1}b_{0}^{2}+2A_{1}b_{1}+b_{0}^{2}B_{1})-2A_{2}b_{0}\bigg{]}+\pi
A_{1}b_{1}\ln^{2}(1-\lambda)\bigg{\\}}\,,$ $\displaystyle g_{3}(\lambda)$
$\displaystyle=$
$\displaystyle\frac{1}{12\pi^{3}b_{0}^{4}(1-\lambda)}\bigg{\\{}\lambda\bigg{[}6\pi
b_{0}^{2}\bigg{(}\big{(}\pi(-2\gamma
A_{1}b_{0}^{2}+2A_{1}b_{1}+b_{0}^{2}B_{1})-2A_{2}b_{0}\big{)}\ln(\frac{Q^{2}}{\mu_{R}^{2}})$
(54)
$\displaystyle+b_{0}(1-\lambda)\ln(\frac{\mu_{F}^{2}}{\mu_{R}^{2}})\bigg{(}2A_{2}-\pi
A_{1}b_{0}\ln(\frac{\mu_{F}^{2}}{\mu_{R}^{2}})\bigg{)}+\pi
A_{1}b_{0}^{2}\ln^{2}\bigg{(}\frac{Q^{2}}{\mu_{R}^{2}}\bigg{)}\bigg{)}$
$\displaystyle+\pi b_{0}\big{(}-6\gamma b_{0}(-2A_{2}b_{0}+2\pi A_{1}b_{1}+\pi
b_{0}^{2}B_{1})+6\pi(2A_{1}b_{2}+b_{0}b_{1}B_{1})-6(2A_{2}b_{1}$
$\displaystyle+b_{0}^{2}B_{2})+\pi^{3}A_{1}b_{0}^{3}+6\gamma^{2}\pi
A_{1}b_{0}^{3}\big{)}+6\lambda\big{(}A_{3}b_{0}^{2}-\pi(A_{2}b_{0}b_{1}+\pi
A_{1}(b_{0}b_{2}-b_{1}^{2}))\big{)}\bigg{]}$
$\displaystyle+6\pi\ln(1-\lambda)\bigg{[}2\pi
A_{1}(b_{2}b_{0}(1-\lambda)+b_{1}^{2}\lambda-\gamma b_{1}b_{0}^{2})+2\pi
A_{1}b_{1}b_{0}^{2}\ln(\frac{Q^{2}}{\mu_{R}^{2}})$
$\displaystyle-2A_{2}b_{1}b_{0}+\pi
b_{1}b_{0}^{2}B_{1}\bigg{]}+6\pi^{2}A_{1}b_{1}^{2}\ln^{2}(1-\lambda)\bigg{\\}}\,.$

The coefficients of the hard matching function in eq. (5), in an expansion in
$\alpha_{S}/(2\pi)$, read

$\displaystyle H_{0}$ $\displaystyle=$ $\displaystyle 1\,,$ $\displaystyle
H_{1}$ $\displaystyle=$ $\displaystyle
C_{F}\bigg{[}\bigg{(}\frac{3}{2}-2\gamma\bigg{)}\ln(\frac{Q^{2}}{\mu_{F}^{2}})+\frac{5\pi^{2}}{6}+\gamma^{2}+\frac{3\gamma}{2}-\frac{9}{2}\bigg{]}\,,$
$\displaystyle H_{2}$ $\displaystyle=$ $\displaystyle
C_{F}\bigg{\\{}C_{F}\bigg{[}\bigg{(}\frac{9}{8}-3\gamma+2\gamma^{2}\bigg{)}\ln^{2}\bigg{(}\frac{Q^{2}}{\mu_{F}^{2}}\bigg{)}+\bigg{(}6\zeta_{3}-\frac{51}{8}+\frac{45\gamma}{4}-\frac{3\gamma^{2}}{2}-2\gamma^{3}+\frac{3\pi^{2}}{4}$
(55)
$\displaystyle-\frac{5\gamma\pi^{2}}{3}\bigg{)}\ln(\frac{Q^{2}}{\mu_{F}^{2}})+6\gamma\zeta_{3}-\frac{33\zeta_{3}}{2}+\frac{61\pi^{4}}{180}+\frac{5\gamma^{2}\pi^{2}}{6}+\frac{3\gamma\pi^{2}}{4}-\frac{35\pi^{2}}{16}+\frac{\gamma^{4}}{2}+\frac{3\gamma^{3}}{2}$
$\displaystyle-\frac{27\gamma^{2}}{8}-\frac{51\gamma}{8}+\frac{331}{32}\bigg{]}+C_{A}\bigg{[}\bigg{(}\frac{11\gamma}{6}-\frac{11}{8}\bigg{)}\ln^{2}\bigg{(}\frac{Q^{2}}{\mu_{F}^{2}}\bigg{)}+\bigg{(}\bigg{(}\frac{11}{4}-\frac{11\gamma}{3}\bigg{)}\ln(\frac{Q^{2}}{\mu_{F}^{2}})$
$\displaystyle+\frac{55\pi^{2}}{36}+\frac{11\gamma^{2}}{6}+\frac{11\gamma}{4}-\frac{33}{4}\bigg{)}\ln(\frac{\mu_{R}^{2}}{\mu_{F}^{2}})+\bigg{(}-3\zeta_{3}+\frac{215}{24}-\frac{367\gamma}{36}-\frac{11\gamma^{2}}{6}-\frac{11\pi^{2}}{12}$
$\displaystyle+\frac{\gamma\pi^{2}}{3}\bigg{)}\ln(\frac{Q^{2}}{\mu_{F}^{2}})-10\gamma\zeta_{3}+\frac{116\zeta_{3}}{9}-\frac{23\pi^{4}}{240}-\frac{\gamma^{2}\pi^{2}}{6}-\frac{11\gamma\pi^{2}}{36}+\frac{1657\pi^{2}}{432}+\frac{11\gamma^{3}}{18}$
$\displaystyle+\frac{367\gamma^{2}}{72}+\frac{3155\gamma}{216}-\frac{5465}{288}\bigg{]}+T_{F}n_{f}\bigg{[}\bigg{(}\frac{1}{2}-\frac{2\gamma}{3}\bigg{)}\ln^{2}\bigg{(}\frac{Q^{2}}{\mu_{F}^{2}}\bigg{)}+\bigg{(}-\frac{19}{6}+\frac{29\gamma}{9}$
$\displaystyle+\frac{2\gamma^{2}}{3}+\frac{\pi^{2}}{3}\bigg{)}\ln(\frac{Q^{2}}{\mu_{F}^{2}})+\bigg{(}\bigg{(}\frac{4\gamma}{3}-1\bigg{)}\ln(\frac{Q^{2}}{\mu_{F}^{2}})-\frac{5\pi^{2}}{9}-\frac{2\gamma^{2}}{3}-\gamma+3\bigg{)}\ln(\frac{\mu_{R}^{2}}{\mu_{F}^{2}})$
$\displaystyle+\frac{2\zeta_{3}}{9}+\frac{\gamma\pi^{2}}{9}-\frac{143\pi^{2}}{108}-\frac{2\gamma^{3}}{9}-\frac{29\gamma^{2}}{18}-\frac{247\gamma}{54}+\frac{457}{72}\bigg{]}\bigg{\\}}\,.$

The above results for $g_{1}$, $g_{2}$ and $g_{3}$ match those of ref.
Vogt:2000ci .

##  Appendix B Fitting NPFF to one dataset at a time

In the following we present the results for fitting NPFF to one dataset at a
time. We use the four datasets ALEPH ALEPH:2001pfo , DELPHI DELPHI:2011aa ,
OPAL OPAL:2002plk  and SLD SLD:2002poq  that enter our combined fit. The
extracted NPFF sets, compared to the combined extraction, can be seen in fig.
14.

Of particular interest is the fit to ALEPH data since it contains only
$B$-mesons, rather than the meson-baryon mixture measured by the other
collaborations. In a ALEPH-only fit one would like to check if the
contribution from baryons has a visible effect on the fits. As can be seen in
fig. 14, when only ALEPH data is considered the uncertainty of the fit is
considerably larger than the one of the combined fit. The fit to all data lies
entirely within the uncertainty band of the ALEPH-only fit which suggests that
it may not be possible to convincingly disentangle the effect of $B$-baryons
using this measurement.

In general, as can be seen in fig. 14, there is agreement between the four
individual fits. The fit showing the largest deviation with respect to the
combined fit is obtained from the DELPHI dataset. Nevertheless, it is not
inconsistent with the combined fit. The OPAL data appears to be the most
precise. For most values of $z$ the combined fit has slightly smaller
uncertainties than the OPAL-only fit. Since there are certain tensions between
the four datasets one may wonder if an improvement from combining them can be
expected. Upon closer inspection of fig. 14 it turns out that for almost all
values of $z$ where the errors on the combined fit are smaller than all the
individual fits, the four datasets actually do agree, which implies the
consistency of the reduced uncertainties of the combined fit.

Figure 14: NPFF fits from one dataset at a time (blue) versus the combined fit
(orange). The four datasets are (from left to right): ALEPH, DELPHI, OPAL and
SLD.

## References

  * (1) B. Mele and P. Nason, Nucl. Phys. B 361, 626-644 (1991) [erratum: Nucl. Phys. B 921, 841-842 (2017)]
  * (2) J. Binnewies, B. A. Kniehl and G. Kramer, Phys. Rev. D 58, 014014 (1998) [arXiv:hep-ph/9712482 [hep-ph]].
  * (3) M. Cacciari, M. Greco and P. Nason, JHEP 05, 007 (1998) [arXiv:hep-ph/9803400 [hep-ph]].
  * (4) M. Cacciari, S. Frixione and P. Nason, JHEP 03, 006 (2001) [arXiv:hep-ph/0102134 [hep-ph]].
  * (5) M. Cacciari and P. Nason, Phys. Rev. Lett. 89, 122003 (2002) [arXiv:hep-ph/0204025 [hep-ph]].
  * (6) M. Cacciari, S. Frixione, M. L. Mangano, P. Nason and G. Ridolfi, JHEP 07, 033 (2004) [arXiv:hep-ph/0312132 [hep-ph]].
  * (7) M. Cacciari, P. Nason and R. Vogt, Phys. Rev. Lett. 95, 122001 (2005) [arXiv:hep-ph/0502203 [hep-ph]].
  * (8) M. Cacciari, P. Nason and C. Oleari, JHEP 10, 034 (2005) [arXiv:hep-ph/0504192 [hep-ph]].
  * (9) B. A. Kniehl and G. Kramer, Phys. Rev. D 71, 094013 (2005) [arXiv:hep-ph/0504058 [hep-ph]].
  * (10) M. Cacciari, P. Nason and C. Oleari, JHEP 04, 006 (2006) [arXiv:hep-ph/0510032 [hep-ph]].
  * (11) B. A. Kniehl and G. Kramer, Phys. Rev. D 74, 037502 (2006) [arXiv:hep-ph/0607306 [hep-ph]].
  * (12) T. Kneesch, B. A. Kniehl, G. Kramer and I. Schienbein, Nucl. Phys. B 799, 34-59 (2008) [arXiv:0712.0481 [hep-ph]].
  * (13) M. Cacciari, S. Frixione, N. Houdeau, M. L. Mangano, P. Nason and G. Ridolfi, JHEP 10, 137 (2012) [arXiv:1205.6344 [hep-ph]].
  * (14) M. Cacciari, M. L. Mangano and P. Nason, Eur. Phys. J. C 75, no.12, 610 (2015) [arXiv:1507.06197 [hep-ph]].
  * (15) B. A. Kniehl, G. Kramer, I. Schienbein and H. Spiesberger, Phys. Rev. D 101, no.11, 114021 (2020) [arXiv:2004.04213 [hep-ph]].
  * (16) P. J. Rijken and W. L. van Neerven, Phys. Lett. B 386, 422-428 (1996) [arXiv:hep-ph/9604436 [hep-ph]].
  * (17) P. J. Rijken and W. L. van Neerven, Phys. Lett. B 392, 207-215 (1997) [arXiv:hep-ph/9609379 [hep-ph]].
  * (18) P. J. Rijken and W. L. van Neerven, Nucl. Phys. B 487, 233-282 (1997) [arXiv:hep-ph/9609377 [hep-ph]].
  * (19) A. Mitov and S. O. Moch, Nucl. Phys. B 751, 18-52 (2006) [arXiv:hep-ph/0604160 [hep-ph]].
  * (20) K. Melnikov and A. Mitov, Phys. Rev. D 70, 034027 (2004) [arXiv:hep-ph/0404143 [hep-ph]].
  * (21) A. Mitov, Phys. Rev. D 71, 054021 (2005) [arXiv:hep-ph/0410205 [hep-ph]].
  * (22) A. Mitov, S. Moch and A. Vogt, Phys. Lett. B 638, 61-67 (2006) [arXiv:hep-ph/0604053 [hep-ph]].
  * (23) S. Moch and A. Vogt, Phys. Lett. B 659, 290-296 (2008) [arXiv:0709.3899 [hep-ph]].
  * (24) A. A. Almasy, S. Moch and A. Vogt, Nucl. Phys. B 854, 133-152 (2012) [arXiv:1107.2263 [hep-ph]].
  * (25) M. Cacciari and S. Catani, Nucl. Phys. B 617, 253-290 (2001) [arXiv:hep-ph/0107138 [hep-ph]].
  * (26) M. Cacciari and E. Gardi, Nucl. Phys. B 664, 299-340 (2003) [arXiv:hep-ph/0301047 [hep-ph]].
  * (27) E. Gardi, JHEP 02, 053 (2005) [arXiv:hep-ph/0501257 [hep-ph]].
  * (28) V. Bertone, S. Carrazza and J. Rojo, Comput. Phys. Commun. 185, 1647-1668 (2014) [arXiv:1310.1394 [hep-ph]].
  * (29) V. Bertone, S. Carrazza and E. R. Nocera, JHEP 03, 046 (2015) [arXiv:1501.00494 [hep-ph]].
  * (30) M. Salajegheh, S. M. Moosavi Nejad, H. Khanpour, B. A. Kniehl and M. Soleymaninia, Phys. Rev. D 99, no.11, 114001 (2019) [arXiv:1904.08718 [hep-ph]].
  * (31) M. Fickinger, S. Fleming, C. Kim and E. Mereghetti, JHEP 11, 095 (2016) [arXiv:1606.07737 [hep-ph]].
  * (32) G. Ridolfi, M. Ubiali and M. Zaro, JHEP 01, 196 (2020) [arXiv:1911.01975 [hep-ph]].
  * (33) E. Braaten, K. m. Cheung, S. Fleming and T. C. Yuan, Phys. Rev. D 51, 4819-4829 (1995) [arXiv:hep-ph/9409316 [hep-ph]].
  * (34) U. Aglietti, G. Corcella and G. Ferrera, Nucl. Phys. B 775, 162-201 (2007) [arXiv:hep-ph/0610035 [hep-ph]].
  * (35) D. P. Anderle, T. Kaufmann, M. Stratmann, F. Ringer and I. Vitev, Phys. Rev. D 96, no.3, 034028 (2017) [arXiv:1706.09857 [hep-ph]].
  * (36) L. Dai, C. Kim and A. K. Leibovich, JHEP 09, 109 (2018) [arXiv:1805.06014 [hep-ph]].
  * (37) L. Dai, C. Kim and A. K. Leibovich, JHEP 09, 148 (2021) [arXiv:2104.14707 [hep-ph]].
  * (38) S. Caletti, A. J. Larkoski, S. Marzani and D. Reichelt, [arXiv:2205.01117 [hep-ph]].
  * (39) S. Caletti, A. J. Larkoski, S. Marzani and D. Reichelt, [arXiv:2205.01109 [hep-ph]].
  * (40) M. Czakon, T. Generet, A. Mitov and R. Poncelet, JHEP 10, 216 (2021) [arXiv:2102.08267 [hep-ph]].
  * (41) J. Binnewies, B. A. Kniehl and G. Kramer, Z. Phys. C 65, 471-480 (1995) [arXiv:hep-ph/9407347 [hep-ph]].
  * (42) B. A. Kniehl, G. Kramer and B. Potter, Phys. Rev. Lett. 85, 5288-5291 (2000) [arXiv:hep-ph/0003297 [hep-ph]].
  * (43) B. A. Kniehl, G. Kramer and B. Potter, Nucl. Phys. B 582, 514-536 (2000) [arXiv:hep-ph/0010289 [hep-ph]].
  * (44) B. A. Kniehl, G. Kramer and B. Potter, Nucl. Phys. B 597, 337-369 (2001) [arXiv:hep-ph/0011155 [hep-ph]].
  * (45) S. Albino, B. A. Kniehl and G. Kramer, Nucl. Phys. B 725, 181-206 (2005) [arXiv:hep-ph/0502188 [hep-ph]].
  * (46) D. de Florian, R. Sassot and M. Stratmann, Phys. Rev. D 75, 114010 (2007) [arXiv:hep-ph/0703242 [hep-ph]].
  * (47) D. de Florian, R. Sassot, M. Epele, R. J. Hernández-Pinto and M. Stratmann, Phys. Rev. D 91, no.1, 014035 (2015) [arXiv:1410.6027 [hep-ph]].
  * (48) D. P. Anderle, F. Ringer and M. Stratmann, Phys. Rev. D 92, no.11, 114017 (2015) [arXiv:1510.05845 [hep-ph]].
  * (49) N. Sato, J. J. Ethier, W. Melnitchouk, M. Hirai, S. Kumano and A. Accardi, Phys. Rev. D 94, no.11, 114004 (2016) [arXiv:1609.00899 [hep-ph]].
  * (50) V. Bertone et al. [NNPDF], Eur. Phys. J. C 77, no.8, 516 (2017) [arXiv:1706.07049 [hep-ph]].
  * (51) E. Moffat et al. [Jefferson Lab Angular Momentum (JAM)], Phys. Rev. D 104, no.1, 016015 (2021) [arXiv:2101.04664 [hep-ph]].
  * (52) R. A. Khalek, V. Bertone and E. R. Nocera, Phys. Rev. D 104, no.3, 034007 (2021) [arXiv:2105.08725 [hep-ph]].
  * (53) H. Abdolmaleki et al. [xfitter Developers’ Team], Phys. Rev. D 104, no.5, 056019 (2021) [arXiv:2105.11306 [hep-ph]].
  * (54) I. Borsa, D. de Florian, R. Sassot and M. Stratmann, Phys. Rev. D 105, no.3, L031502 (2022) [arXiv:2110.14015 [hep-ph]].
  * (55) I. Borsa, R. Sassot, D. de Florian, M. Stratmann and W. Vogelsang, Phys. Rev. Lett. 129, no.1, 012002 (2022) [arXiv:2202.05060 [hep-ph]].
  * (56) R. A. Khalek, V. Bertone, A. Khoudli and E. R. Nocera, [arXiv:2204.10331 [hep-ph]].
  * (57) G. Corcella and A. D. Mitov, Nucl. Phys. B 623, 247-270 (2002) [arXiv:hep-ph/0110319 [hep-ph]].
  * (58) M. Cacciari, G. Corcella and A. D. Mitov, JHEP 12, 015 (2002) [arXiv:hep-ph/0209204 [hep-ph]].
  * (59) G. Corcella and V. Drollinger, Nucl. Phys. B 730, 82-102 (2005) [arXiv:hep-ph/0508013 [hep-ph]].
  * (60) S. Biswas, K. Melnikov and M. Schulze, JHEP 08, 048 (2010) [arXiv:1006.0910 [hep-ph]].
  * (61) B. A. Kniehl, G. Kramer and S. M. Moosavi Nejad, Nucl. Phys. B 862, 720-736 (2012) [arXiv:1205.2528 [hep-ph]].
  * (62) K. Agashe, R. Franceschini and D. Kim, Phys. Rev. D 88, no.5, 057701 (2013) [arXiv:1209.0772 [hep-ph]].
  * (63) S. M. Moosavi Nejad, Phys. Rev. D 88, no.9, 094011 (2013) [arXiv:1310.5686 [hep-ph]].
  * (64) S. M. Moosavi Nejad and M. Balali, Eur. Phys. J. C 76, no.3, 173 (2016) [arXiv:1602.05322 [hep-ph]].
  * (65) K. Agashe, R. Franceschini, D. Kim and M. Schulze, Eur. Phys. J. C 76, no.11, 636 (2016) [arXiv:1603.03445 [hep-ph]].
  * (66) C. Brancaccio, M. Czakon, T. Generet and M. Krämer, JHEP 08, 145 (2021) [arXiv:2106.06516 [hep-ph]].
  * (67) S. Catani, M. L. Mangano, P. Nason and L. Trentadue, Nucl. Phys. B 478, 273-310 (1996) [arXiv:hep-ph/9604351 [hep-ph]].
  * (68) S. A. Larin, T. van Ritbergen and J. A. M. Vermaseren, Nucl. Phys. B 438, 278-306 (1995) [arXiv:hep-ph/9411260 [hep-ph]].
  * (69) K. G. Chetyrkin, B. A. Kniehl and M. Steinhauser, Nucl. Phys. B 510, 61-87 (1998) [arXiv:hep-ph/9708255 [hep-ph]].
  * (70) Y. Schroder and M. Steinhauser, JHEP 01, 051 (2006) [arXiv:hep-ph/0512058 [hep-ph]].
  * (71) K. G. Chetyrkin, J. H. Kuhn and C. Sturm, Nucl. Phys. B 744, 121-135 (2006) [arXiv:hep-ph/0512060 [hep-ph]].
  * (72) M. Neubert, [arXiv:0706.2136 [hep-ph]].
  * (73) M. Buza, Y. Matiounine, J. Smith and W. L. van Neerven, Eur. Phys. J. C 1, 301-320 (1998) [arXiv:hep-ph/9612398 [hep-ph]].
  * (74) S. Catani, D. de Florian, M. Grazzini and P. Nason, JHEP 07, 028 (2003) [arXiv:hep-ph/0306211 [hep-ph]].
  * (75) A. Heister et al. [ALEPH], Phys. Lett. B 512, 30-48 (2001) [arXiv:hep-ex/0106051 [hep-ex]].
  * (76) J. Abdallah et al. [DELPHI], Eur. Phys. J. C 71, 1557 (2011) [arXiv:1102.4748 [hep-ex]].
  * (77) G. Abbiendi et al. [OPAL], Eur. Phys. J. C 29, 463-478 (2003) [arXiv:hep-ex/0210031 [hep-ex]].
  * (78) K. Abe et al. [SLD], Phys. Rev. D 65, 092006 (2002) [erratum: Phys. Rev. D 66, 079905 (2002)] [arXiv:hep-ex/0202031 [hep-ex]].
  * (79) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr and G. Watt, Eur. Phys. J. C 75, 132 (2015) [arXiv:1412.7420 [hep-ph]].
  * (80) A. Kharchilava, Phys. Lett. B 476, 73-78 (2000) [arXiv:hep-ph/9912320 [hep-ph]].
  * (81) T. Aaltonen et al. [CDF], Phys. Rev. D 80, 051104 (2009) [arXiv:0906.5371 [hep-ex]].
  * (82) V. Khachatryan et al. [CMS], JHEP 12, 123 (2016) [arXiv:1608.03560 [hep-ex]].
  * (83) [ATLAS], [arXiv:2209.00583 [hep-ex]].
  * (84) D. J. Lange, Nucl. Instrum. Meth. A 462, 152-155 (2001)
  * (85) P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020)
  * (86) R. D. Ball et al. [NNPDF], Eur. Phys. J. C 77, no.10, 663 (2017) [arXiv:1706.00428 [hep-ph]].
  * (87) M. Czakon, Phys. Lett. B 693, 259-268 (2010) [arXiv:1005.0274 [hep-ph]].
  * (88) M. Czakon and D. Heymes, Nucl. Phys. B 890, 152-227 (2014) [arXiv:1408.2500 [hep-ph]].
  * (89) F. Maltoni, G. Ridolfi, M. Ubiali and M. Zaro, JHEP 10, 027 (2022) [arXiv:2207.10038 [hep-ph]].
  * (90) A. Vogt, Phys. Lett. B 497, 228-234 (2001) [arXiv:hep-ph/0010146 [hep-ph]].