# Transverse Momentum-Dependent Heavy-Quark Fragmentation at Next-to-Leading
Order

Rebecca von Kuk  Johannes K. L. Michel  and Zhiquan Sun

(April 25, 2024)

###### Abstract

The transverse momentum-dependent fragmentation functions (TMD FFs) of heavy
(bottom and charm) quarks, which we recently introduced, are universal
building blocks that enter predictions for a large number of observables
involving final-state heavy quarks or hadrons. They enable the extension of
fixed-order subtraction schemes to quasi-collinear limits, and are of
particular interest in their own right as probes of the nonperturbative
dynamics of hadronization. In this paper we calculate all TMD FFs involving
heavy quarks and the associated TMD matrix element in heavy-quark effective
theory (HQET) to next-to-leading order in the strong interaction. Our results
confirm the renormalization properties, large-mass, and small-mass consistency
relations predicted in our earlier work. We also derive and confirm a
prediction for the large-$z$ behavior of the heavy-quark TMD FF by extending,
for the first time, the formalism of joint resummation to capture quark mass
effects in heavy-quark fragmentation. Our final results in position space
agree with those of a recent calculation by another group that used a highly
orthogonal organization of singularities in the intermediate momentum-space
steps, providing a strong independent cross check. As an immediate
application, we present the complete quark mass dependence of the energy-
energy correlator (EEC) in the back-to-back limit at
$\mathcal{O}(\alpha_{s})$.

##  1 Introduction

The fragmentation of heavy (bottom or charm) quarks into the experimentally
observed heavy meson and baryon states is of particular interest because the
mass of the quark imprints as a perturbative scale on the otherwise
nonperturbative process of hadronization. In ref. [1], we initiated the study
of transverse momentum-dependent (TMD) fragmentation functions (FFs) for the
formation of heavy hadrons from a parent heavy quark. Our work generalized the
well-studied case where only the collinear momentum fraction of the hadron is
resolved [2, 3, 4, 5, 6] (see refs. [7] for a recent precision study on data),
and is part of a larger, ongoing research effort to understand differential
jet and fragmentation functions involving heavy quarks, heavy hadrons, and
heavy-quark bound states [8, 9, 1, 10, 11, 12, 13, 14]. In this context, TMD
FFs are unique in the wealth of information they can provide on the
hadronization process, essentially offering a full three-dimensional view of
the fragmentation cascade. The main goal of the present paper is to evaluate
all unpolarized heavy-quark TMD FFs at perturbative transverse momenta to
complete next-to-leading order (NLO) in the strong coupling, with the aim of
improving the baseline perturbative precision for studying their rich
nonperturbative structure.

Importantly, the factorization paradigm ensures that the heavy-quark TMD FFs
appear as universal building blocks across predictions for a large number of
processes involving final-state heavy quarks. In ref. [1], we explicitly
considered their phenomenology in heavy-hadron pair production in the back-to-
back limit in $e^{+}e^{-}$ collisions and in (polarized) semi-inclusive DIS at
the future Electron-Ion Collider [15], cf. ref. [16] for recent dedicated
projections for polarized _collinear_ parton distribution functions (PDFs) and
FFs involving heavy quarks at the EIC. However, the applicability of heavy-
quark TMD FFs is much wider: Applications of immediate phenomenological
interest range from $Z$$+$hadron [17] and dihadron [18] azimuthal
decorrelations in $pp$ and $pA$ collisions to hadron-in-jet transverse-
momentum distributions [19, 20], transverse momentum-like event shapes
(extending the calculation of mass effects in thrust-like event shapes [21,
22, 23, 24]), and energy-energy correlators in the back-to-back limit in
$e^{+}e^{-}$ [25], $ep$ and $eA$ [26, 27], and $pp$ collisions [28, 29]. In
all of these cases, factorization formulas involving TMD FFs have been
derived, and the NLO results we present in this paper enable one to fully
account for the effect of quark masses and accurately capture their highly
nontrivial interplay with transverse momenta in all of these processes.

Figure -179:  The probability density
$\mathrm{d}\mathcal{P}_{Q\bar{Q}\to(Qg)\bar{Q}}$ for a gluon to split off a
heavy quark in the quasi-collinear limit, which coincides with the heavy-quark
TMD FF $d_{1\,Q/Q}$ up to a factor of $z^{2}$ for finite collinear gluon
momentum fraction $\bar{z}>0$ and transverse momentum $k_{T}>0$. On the rear
walls of the plot we show the integrals of $d_{1\,Q/Q}$ over one of $k_{T}$ or
$z$ as predicted by the NLO calculation we perform in this paper, showcasing
that the TMD regularization and renormalization procedure renders them finite
despite the intricate singular behavior of the integrand. The same is true for
the total integral over both $z$ and $k_{T}<m$, which is not shown in the
plot, but at the values of the scales chosen here is equal to
$\frac{\alpha_{s}C_{F}}{4\pi}\bigl{(}5-\frac{\pi}{2}-\frac{2\pi^{2}}{3}\bigr{)}\approx-0.067$.
For numerical illustration, we take $\alpha_{s}(\mu=m)=0.2$.

The heavy-quark TMD FF is intimately connected with the universal gluon
splitting probability off a heavy-quark in the quasi-collinear limit [30, 31],
e.g. after heavy-quark pair production has taken place,

$\displaystyle\frac{\mathrm{d}\mathcal{P}_{Q\bar{Q}\to(Qg)\bar{Q}}}{\mathrm{d}z\,\mathrm{d}(k_{T}^{2})}=\frac{\alpha_{s}C_{F}}{2\pi}\,\frac{k_{T}^{2}z^{2}(1+z^{2})+m^{2}(1-z)^{4}}{(1-z)\bigl{[}k_{T}^{2}z^{2}+m^{2}(1-z)^{2}\bigr{]}^{2}}\,,$
(1.1)

where $m$ is the heavy-quark mass, $k_{T}$ is the transverse momentum of the
gluon with respect to the heavy quark, and $z$ is the collinear momentum
fraction retained by the heavy quark. The quasi-collinear splitting amplitude
is shown in a three-dimensional plot as a function of $k_{T}$ and the momentum
fraction $\bar{z}\equiv 1-z$ carried by the gluon in figure -179. Famously,
and in contrast to the behavior of light partons, the splitting amplitude is
suppressed for energetic gluons with large momentum fractions $\bar{z}\gtrsim
0.2$, an effect known as the dead cone [32]. The splitting probability is
strongly peaked at small $\bar{z}$ and small $k_{T}\lesssim m$. Unlike for
light partons, it eventually turns over for any finite value of $\bar{z}>0$
and vanishes (near the front right edge of the plot) as $k_{T}\to 0$, yet the
simultaneous limit $\bar{z},k_{T}\to 0$ diverges. A key challenge we address
in this paper is how to consistently extend this complicated singular behavior
of the quasi-collinear splitting probability into the limits $k_{T}\to
0,\bar{z}\to 0$ where the gluon is unresolved. The right framework for doing
so is – indeed – TMD factorization [33], which provides precise field-
theoretic prescriptions to isolate, regulate, and renormalize the associated
divergences, and it is precisely in this sense that the heavy-quark TMD FF is
the formal generalization of eq. (1.1) to $\bar{z}\to 0$ and/or $k_{T}\to 0$
in a way that is manifestly finite in four spacetime dimensions. As a
demonstration, we include on the rear walls of the three-dimensional plot in
figure -179 the integral of the heavy-quark TMD FF over one of $z$ or
$k_{T}\leq m$, respectively, as predicted by our calculations in this paper,
showcasing that they are indeed individually finite. Similarly, the total
integral over both variables $0\leq z\leq 1$ and $0\leq k_{T}\leq m$ (and thus
any integral over a subset of the plane) is also finite in the usual
distributional sense, and for the indicated values of the scales simply
evaluates to
$\frac{\alpha_{s}C_{F}}{4\pi}\bigl{(}5-\frac{\pi}{2}-\frac{2\pi^{2}}{3}\bigr{)}$.
We stress that isolating this singular behavior at fixed order in terms of a
universal and finite heavy-quark TMD FF is not an academic exercise, but
instead is of immediate interest to extend fixed-order subtraction schemes
like that of $q_{T}$ subtractions [34] and its multi-jet generalizations [35,
36, 37, 38] to quasi-collinear limits near final-state heavy quarks. Such
extensions, which can potentially yield large gains in numerical efficiency,
are likewise an active topic of current research, see e.g. [39, 40, 41, 42].

The remainder of this paper is structured as follows: Section 2 contains a
review of the theoretical framework introduced in ref. [1] to describe heavy-
quark fragmentation in the transverse plane. In section 3, we calculate all
heavy-quark TMD FFs to next-to-leading order in the strong coupling, and
present the results in both momentum and position space. In section 4, we
perform an analogous calculation for the effective TMD matrix elements that
arise in boosted heavy-quark effective theory (bHQET) in the limit where the
transverse momentum is small compared to the mass. We employ these results and
known one-loop results in the literature to verify the expected small and
large-mass consistency relations in section 5. Section 6 is dedicated to an
extension of the joint resummation formalism to massive quark fragmentation,
which explains the surprising simplicity of our results at large $z$. As an
immediate, first application of our results, we present the complete mass
dependence of the energy-energy correlator in the back-to-back limit in
section 7. In section 8, we compare our main results to those of ref. [13],
where an independent calculation of the one-loop heavy-quark TMD FF was
recently performed. We find full agreement at the level of the final position-
space results, but stress the differences in how the intermediate momentum-
space calculations are organized (amounting to a strong independent cross
check) and – on the phenomenological side – the differences in how the $z$
dependence is treated in the bHQET limit. We conclude in section 9.

##  2 Theoretical framework

###  2.1 Heavy-quark TMD FFs

We study the fragmentation of a heavy quark $Q$ with pole mass
$m\gg\Lambda_{\mathrm{QCD}}$ into an unpolarized heavy hadron $H$ and
additional fragmentation products $X$. We work in QCD with $n_{f}=n_{\ell}+1$
flavors, where $n_{\ell}$ is the number of massless quark flavors. We work in
the “hadron frame” for fragmentation [33] where
$P_{H}^{\mu}=(P_{H}^{-},M_{H}^{2}/P_{H}^{-},0)$, with $M_{H}$ the hadron mass,
and $P_{H,\perp}=0$ by definition. (We refer to Appendix A.1 for our lightcone
conventions.)

In ref. [1] we showed that up to a factor of the total probability for $Q$ to
fragment into $H$, the distribution differential in $k_{T}=P_{X,T}\sim M_{H}$
and the lightcone momentum fraction $z_{H}$ of the original quark retained by
$H$ is governed by a new perturbative matching coefficient, the _partonic
heavy-quark TMD FF_ $d_{1\,Q/Q}(z,k_{T},\mu,\zeta)$,

$\displaystyle
D_{1\,H/Q}(z_{H},b_{T},\mu,\zeta)=d_{1\,Q/Q}(z_{H},b_{T},\mu,\zeta)\,\chi_{H}+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{m}\Bigr{)}+\mathcal{O}(\Lambda_{\mathrm{QCD}}b_{T})\,.$
(2.1)

Making all regulators explicit, the formal definition of $d_{1\,Q/Q}$ at the
bare level reads

$\displaystyle
d_{1\,Q/Q}^{\mathrm{bare}}(z,b_{T},\epsilon,\eta,\zeta/\nu^{2})$
$\displaystyle=\frac{1}{2z^{1-2\epsilon}N_{c}}\int\\!\frac{\mathrm{d}b^{+}}{4\pi}\,e^{\mathrm{i}b^{+}(p^{-}/z)/2}$
(2.2)
$\displaystyle\quad\times\operatorname{Tr}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{X}\operatorname{tr}\Bigl{[}\frac{\not{{\bar{n}}}}{2}\,\big{\langle}0\big{\lvert}W_{\eta}^{\dagger}(b)\,\psi_{Q}(b)\big{\rvert}QX\big{\rangle}\big{\langle}QX\big{\lvert}\bar{\psi}_{Q}(0)\,W_{\eta}(0)\big{\rvert}0\big{\rangle}\Bigr{]}\,,$
(2.11)

where we work in $d=4-2\epsilon$ spacetime dimensions, $N_{c}$ is the number
of colors, $\operatorname{Tr}$ ($\operatorname{tr}$) indicates a trace over
color (spin), and
$\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{X}$
indicates a sum over all possible partonic final states combined with an
integral over their phase space. The fields in the first matrix element are
evaluated at a spacetime position $b\equiv(0,b^{+},b_{\perp})$ with
$b_{\perp}$ Fourier conjugate to $k_{\perp}$, and $b_{T}^{2}\equiv-
b_{\perp}^{2}$ and $k_{T}^{2}=-k_{\perp}^{2}$, respectively. The heavy quark
in the external state carries momentum $p^{\mu}=(p^{-},m^{2}/p^{-},0)$, i.e.,
the above definition is equal to the hadron-level definition of the heavy-
quark TMD FF $D_{1\,H/Q}$, but with $H$ replaced by $Q$ itself in the external
state and restricting to partonic final states $X$ [1]. The definition in eq.
(2.2) is given in terms of the heavy-quark field renormalized on shell,

$\displaystyle\psi^{\mathrm{bare}}_{Q}(x)=Z_{\psi,\mathrm{OS}}^{1/2}(m,\mu,\epsilon)\,\psi_{Q}(x)\,,$
(2.12)

and in terms of Wilson lines $W_{\eta}(x)$ defined as anti-path ordered
exponentials of gauge fields extending to positive infinity along the
lightcone direction $\bar{n}^{\mu}$,

$\displaystyle
W_{\eta}(x)=\bar{P}\Bigl{[}\exp\Bigl{(}-\mathrm{i}g\int_{0}^{\infty}\\!\mathrm{d}s\,{\bar{n}}\cdot
A(x+{\bar{n}}s)\Bigr{)}\Bigr{]}_{\eta}\,.$ (2.13)

The subscript $\eta$ indicates the presence of an additional rapidity
regulator. In explicit calculations, we will use the so-called $\eta$
regulator of refs. [43, 44] whose action on a single gluon with momentum
$\ell$ attached to the Wilson line amounts to inserting a factor of
$\lvert\ell^{-}/\nu\rvert^{\eta}$ under phase-space or loop integrals over
$\ell$, with $\nu$ a dimensionful rapidity scale. The final result for the
bare matrix element depends on $\zeta/\nu^{2}$, with
$\zeta\equiv(p^{-}/z)^{2}$ the Collins-Soper scale.

In terms of the bare collinear matrix element in eq. (2.2), the renormalized
partonic heavy-quark TMD FF is given by

$\displaystyle d_{1\,Q/Q}(z,b_{T},\mu,\zeta)$ $\displaystyle=\lim_{\epsilon\to
0}Z_{\mathrm{UV}}(\mu,\zeta,\epsilon)\lim_{\eta\to
0}\Bigl{[}d_{1\,Q/Q}^{\mathrm{bare}}(z,b_{T},\epsilon,\eta,\zeta/\nu^{2})\,\sqrt{S}(b_{T},m,\epsilon,\eta,\nu)\Bigr{]}\,,$
(2.14)

where $S$ is the universal bare TMD soft function for the $\eta$ regulator,
which cancels all poles of $\eta$ and the associated $\nu$ dependence, and is
independent of the heavy-quark mass up to secondary quark mass effects
starting at two loops. By contrast, the $\overline{\text{MS}}$ renormalization
factor $Z_{\mathrm{UV}}$ for quark TMDs is independent of the mass to all
orders by RG consistency with the hard matching coefficient at the hard
scattering energy $Q\sim p^{-}/z\gg m\sim k_{T}\sim 1/b_{T}$. While soft
subtractions (known as zero bins in the SCET literature) generally need to be
accounted for in both virtual and real collinear diagrams when computing the
bare collinear matrix element itself, we have explicitly verified that they
lead to scaleless integrals for our choice of regulator, also in the presence
of the quark mass.

###  2.2 Nonvalence contributions

There are two further ways in which heavy quarks can participate in the TMD
fragmentation process for $m\sim k_{T}\sim 1/b_{T}\gg\Lambda_{\mathrm{QCD}}$;
these were sketched in ref. [1], but we spell them out here explicitly. In one
case, the resolved final state is a heavy hadron $H$, while the parent parton
$i=g,q,\bar{q}$ is light,

$\displaystyle
D_{1\,H/i}(z_{H},b_{T},\mu,\zeta)=d_{1\,Q/i}(z_{H},b_{T},\mu,\zeta)\,\chi_{H}+\mathcal{O}(\Lambda_{\mathrm{QCD}})\,.$
(2.15)

Here $\chi_{H}$ is again the universal inclusive fragmentation probability for
$Q\to H$, while the bare matching coefficient $d_{1\,Q/i}$ is obtained from
eqs. (2.2) and (2.14) after replacing the heavy-quark fields in the bare
correlator by suitable unpolarized combinations of gluon or light-quark
fields.

In the other case, the resolved final-state hadron $h$ is light. In this case
the TMD FF $D_{1\,h/i}$ for $i=g,q,\bar{Q},Q,\bar{Q}$ has to be matched onto
collinear fragmentation functions $D_{h/j}$ at the scale
$\Lambda_{\mathrm{QCD}}$ where all degrees of freedom $j=g,q,\bar{q}$ are
light,

$\displaystyle
D_{1\,h/i}(z_{H},b_{T},\mu,\zeta)=\frac{1}{z_{H}^{2}}\sum_{j}\int\\!\frac{\mathrm{d}z}{z}\,\mathcal{J}_{j/i}(z,b_{T},m,\mu,\zeta)\,D_{h/j}\Bigl{(}\frac{z_{H}}{z},\mu\Bigr{)}+\mathcal{O}(\Lambda_{\mathrm{QCD}})\,.$
(2.16)

This has the standard form of matching TMD FFs onto collinear FFs, but note
that we include the mass as an additional third argument of the matching
coefficient to make explicit that even when $i$ is light, the heavy quark may
in general contribute to the matching at $\mu\sim m\sim k_{T}$ through closed
loops.

###  2.3 bHQET fragmentation factors and the large-mass limit

In ref. [1], we showed that for transverse momenta $k_{T}\ll m$, the TMD
dynamics are governed by new (and in general nonperturbative) matrix elements
defined in boosted HQET [45, 46] that we dubbed TMD fragmentation factors.
Specifically, for $\Lambda_{\mathrm{QCD}}\lesssim k_{T}\ll m$ and counting

$\displaystyle 1-z_{H}\sim 1\,,$ (2.17)

which amounts to integrating over wide bins in $z_{H}$ or taking low $z_{H}$
moments, the TMD FF for producing a heavy hadron $H$ off a heavy quark $Q$
factorizes as [1]

$\displaystyle
D_{1\,H/Q}(z_{H},b_{T},\mu,\zeta)=\delta(1-z_{H})\,C_{m}(m,\mu,{\zeta})\,\chi_{1,H}\Bigl{(}b_{T},\mu,{\frac{\sqrt{\zeta}}{m}}\Bigr{)}+\mathcal{O}\Bigl{(}\frac{1}{m}\Bigr{)}\,,$
(2.18)

where $C_{m}(m,\mu,\zeta)$ is generated by separately matching collinear and
soft modes at the scale $\mu\sim m$ onto HQET and QCD with $n_{\ell}$ flavors,
respectively [47]; explicit expressions for $C_{m}$ and other perturbative
ingredients to the order required for our perturbative checks are given in
appendix B.

In eq. (2.18), the TMD dynamics are encoded in the unpolarized TMD
fragmentation factor $\chi_{1,H}$. For $\Lambda_{\mathrm{QCD}}\ll k_{T}\sim
1/b_{T}$, we showed in ref. [1] that it is given by a product of the total
fragmentation probability $\chi_{H}$ and a perturbative matching coefficient
$C_{1}$,

$\displaystyle\chi_{1,H}(b_{T},\mu,{\rho})=C_{1}(b_{T},\mu,{\rho})\,\chi_{H}+\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})\,,$
(2.19)

where $C_{1}$ will be the target of our calculation in section 4.3. Here we
have also introduced the shorthand $\rho\equiv v^{-}=\sqrt{\zeta}/m$ for the
boost of the hadron; in ref. [1] we argued from consistency with TMD
factorization that $\chi_{1,H}$ develops an anomalous dependence on this
variable governed by the Collins-Soper kernel, and we will make this
expectation more explicit in section 4.1.

To evaluate $C_{1}$, we will use that – for fully general
$\Lambda_{\mathrm{QCD}}\lesssim 1/b_{T}$ – the unpolarized TMD fragmentation
factor $\chi_{1,H}$ reduces to a vacuum matrix element of Wilson lines when
summed over all hadrons $H$ containing $Q$ [1]. Since $\sum_{H}\chi_{H}=1$ and
$C_{1}$ in eq. (2.19) is independent of the hadronic state $H$, it is then
straightforward to obtain $C_{1}$ by evaluating the Wilson line correlator at
perturbative $b_{T}$.

At the partonic level, eqs. (2.18) and (2.19) together imply the following
consistency condition for the partonic heavy-quark TMD FF in the limit $m\gg
k_{T}\sim 1/b_{T}$ [1],

$\displaystyle
d_{1\,Q/Q}(z,b_{T},\mu,\zeta)=\delta(1-z)\,C_{m}(m,\mu,{\zeta})\,C_{1}\Bigl{(}b_{T},\mu,{\frac{\sqrt{\zeta}}{m}}\Bigr{)}+\mathcal{O}\Bigl{(}\frac{1}{b_{T}m}\Bigr{)}\,,$
(2.20)

and in section 5.1 we will use this relation as a check of our one-loop
results. By contrast, nonvalence (or disfavored) partonic heavy-quark TMD FFs
become power-suppressed in the heavy-quark limit,

$\displaystyle i\neq Q\,:\quad
d_{1\,Q/i}(z,b_{T},\mu,\zeta)=\mathcal{O}\Bigl{(}\frac{1}{b_{T}m}\Bigr{)}\,.$
(2.21)

The indirect effect ($D_{1\,h/i}$, $i\neq Q,\bar{Q}$) of the heavy quark on
light-hadron production for $m\gg k_{T}$ is leading in $1/m$. It is governed
by virtual contributions from collinear and soft so-called mass modes and thus
follows exactly the TMD PDF case [48]. Finally, the direct contribution
($D_{1\,h/Q}$) to light-hadron production becomes strongly peaked at $z_{h}\to
0$ in the limit $m\gg k_{T}\gtrsim\Lambda_{\mathrm{QCD}}$, but for any finite
value of $z_{h}$ (or argument $z$ of $\mathcal{J}_{j/Q}$ at the partonic
level) is kinematically suppressed as $m\gg k_{T}$; we will verify this latter
behavior in section 5.1, but leave a dedicated analysis of its $z_{h}\to 0$
behavior – which is known to be subtle in fragmentation [49] – to future work.

###  2.4 Consistency conditions in the light-quark limit

In the case when the heavy-quark is light compared to $k_{T}$, the matching at
the scale $\mu\sim k_{T}$ takes exactly the standard form for matching TMD FFs
onto twist-2 collinear FFs [1],

$\displaystyle
D_{1\,\mathfrak{h}/Q}(z_{\mathfrak{h}},b_{T},\mu,\zeta)=\frac{1}{z_{\mathfrak{h}}^{2}}\sum_{i}\int\\!\frac{\mathrm{d}z}{z}\,\mathcal{J}_{i/q}(z,b_{T},\mu,\zeta)\,D_{\mathfrak{h}/i}^{(n_{\ell}+1)}\Bigl{(}\frac{z_{\mathfrak{h}}}{z},\mu\Bigr{)}+\mathcal{O}(m^{2}b_{T}^{2})\,,$
(2.22)

where $\mathfrak{h}=h,H$ may be both a heavy or a light hadron and
$D_{\mathfrak{h}/j}^{(n_{\ell}+1)}$ are collinear FFs in a theory with
$n_{\ell}$ light and one massive flavor, i.e., the mass is the highest IR
scale in the twist-2 matching. Importantly, the TMD FF matching coefficients
$\mathcal{J}_{i/q}$ are independent of these details in the IR and are thus
given by the universal TMD FF matching coefficients in a theory with
$n_{f}=n_{\ell}+1$ _light_ degrees of freedom, which are known to N3LO [50,
51]. The problem of the mass dependence is thus reduced to the well-understood
behavior of the collinear FFs, which differs depending on whether
$\mathfrak{h}$ is heavy or light [52, 53],

$\displaystyle D_{H/i}^{(n_{\ell}+1)}(z_{H},\mu)$
$\displaystyle=d_{Q/i}(z_{H},\mu)\,\chi_{H}+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{m}\Bigr{)}\,,$
$\displaystyle D_{h/i}^{(n_{\ell}+1)}(z_{h},\mu)$
$\displaystyle=\sum_{j}\int\\!\frac{\mathrm{d}z}{z}\,\mathcal{M}_{j/i,T}(z,m,\mu)\,D_{h/j}^{(n)}\Bigl{(}\frac{z_{h}}{z},\mu\Bigr{)}\,.$
(2.23)

Here $d_{Q/i}(z_{H},\mu)$ is the partonic _collinear_ heavy-quark FF and
$\mathcal{M}_{j/i,T}$ is the timelike matching function governing the flavor
decoupling in light-hadron collinear FFs.

Comparing eqs. (2.1), (2.15), and (2.16) to eqs. (2.22) and (2.4), we can read
off the following leading-power behavior of the relevant perturbative heavy-
quark TMD matrix elements in the limit $m\ll k_{T}\sim 1/b_{T}$ [1],

$\displaystyle d_{1\,Q/i}(z,b_{T},\mu,\zeta)$
$\displaystyle=\frac{1}{z^{2}}\sum_{j}\int\\!\frac{\mathrm{d}z^{\prime}}{z^{\prime}}\,\mathcal{J}_{j/i}(z^{\prime},b_{T},\mu,\zeta)\,d_{Q/j}\Bigl{(}\frac{z}{z^{\prime}},\mu\Bigr{)}\,,$
$\displaystyle\mathcal{J}_{k/i}(z,b_{T},m,\mu,\zeta)$
$\displaystyle=\frac{1}{z^{2}}\sum_{j}\int\\!\frac{\mathrm{d}z^{\prime}}{z^{\prime}}\,\mathcal{J}_{j/i}(z^{\prime},b_{T},\mu,\zeta)\,\mathcal{M}_{k/j,T}\Bigl{(}\frac{z}{z^{\prime}},m,\mu\Bigr{)}\,.$
(2.24)

In section 5.2, we will use these expressions, together with the various known
ingredients on the right-hand side, to perform cross checks of our one-loop
results in section 3 in all channels.

##  3 Partonic heavy-quark TMD fragmentation at NLO

For ease of calculation and to make contact with the splitting probability in
eq. (1.1), we will perform our calculation of the bare collinear matrix
element as well as its renormalization in momentum space, i.e., as a function
of $z$ and $k_{T}$, and later compute various integral transforms of the
renormalized object directly. Passing to momentum space, the one-loop
correction to the renormalized heavy-quark TMD FF in eq. (2.14) reads

$\displaystyle d_{1\,Q/Q}^{(1)}(z,k_{T},\mu,\zeta)$
$\displaystyle=\lim_{\epsilon\to
0}\biggl{\\{}\delta(1-z)\,\frac{1}{\pi}\delta(k_{T}^{2})\,Z_{\mathrm{UV}}^{(1)}(\mu,\zeta,\epsilon)$
(3.1) $\displaystyle\quad+\lim_{\eta\to
0}\Bigl{[}d_{1\,Q/Q}^{\mathrm{bare}\,(1)}(z,k_{T},\epsilon,\eta,\zeta/\nu^{2})+\delta(1-z)\,\tfrac{1}{2}S^{(1)}(k_{T},\epsilon,\eta,\nu)\Bigr{]}\biggr{\\}}\,,$

where we have used that
$d_{1\,Q/Q}^{(0)}=\delta(1-z)\,\frac{1}{\pi}\delta(k_{T}^{2})$ at tree level.
(Here and in the following we use the same symbol for functions of $b_{T}$ and
their Fourier transforms, as the meaning is always clear from the argument.)
Note that in our implementation of dimensional regularization in the
transverse plane, we make use of the azimuthal symmetry of the unpolarized TMD
FF and take the curly braces in eq. (3.1) to be a density in $\pi k_{T}^{2}$
for convenience, i.e., they have an integer mass dimension of $-2$, where the
factor of $\pi$ ensures that final results in $d=4$ are properly normalized
azimuthally symmetric densities in vectorial $\vec{k}_{T}$ [54, 55]. With
these preliminaries, the bare collinear matrix element in momentum space is
given by

$\displaystyle
d_{1\,Q/Q}^{\mathrm{bare}}(z,k_{T},\epsilon,\eta,\zeta/\nu^{2})$
$\displaystyle=\frac{1}{2z^{1-2\epsilon}N_{c}}\int\\!\frac{\mathrm{d}\Omega_{2-2\epsilon}}{2\pi}\,k_{T}^{-2\epsilon}\operatorname{Tr}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{X}\operatorname{tr}\Bigl{[}\frac{\not{{\bar{n}}}}{2}\,\big{\langle}0\big{\lvert}W_{\eta}^{\dagger}\,\psi_{Q}\big{\rvert}QX\big{\rangle}$
(3.10)
$\displaystyle\quad\times\big{\langle}QX\big{\lvert}\bigl{[}\delta\bigl{(}p^{-}/z+\mathrm{i}\partial^{-}\bigr{)}\,\delta^{(2-2\epsilon)}\bigl{(}k_{\perp}^{\mu}+\mathrm{i}\partial_{\perp}^{\mu}\bigr{)}\,\bar{\psi}_{Q}\,W_{\eta}\bigr{]}\big{\rvert}0\big{\rangle}\Bigr{]}\,.$

Here $\mathrm{d}\Omega_{2-2\epsilon}$ is the solid angle element for
$k_{\perp}^{\mu}$ in $2-2\epsilon$ dimensions and the $\delta$ functions
involving derivative operators on the second line act on the fields to their
right, i.e., they fix the total minus and perpendicular momentum injected into
the compound operator.

Compared to a direct calculation in position space, which involves expanding
Bessel-like hypergeometric functions in the dimensional regulator, the main
challenge in the momentum-space calculation is the careful distributional
treatment of the singularity structure prior to expanding in $\epsilon$, which
we address in section 3.1. Renormalized results in momentum space are
presented in section 3.2, which we independently verify and further simplify
by moving to cumulant space in section 3.3. Results in $b_{T}$ space are then
obtained by a straightforward integral transform mapping finite functions onto
functions, and are given in section 3.4.

###  3.1 Bare real-emission diagrams and distributional expansion

(a)  (b)  (c)

Figure -178:  Diagrams contributing to the partonic heavy-quark TMD FF in
Feynman gauge. Heavy quark propagators and field insertions are indicated in
red. In diagrams (a) and (c) the gluon attaches to the Wilson lines that are
part of the operators denoted by ${\color[rgb]{1,0,0}\otimes}$. The dashed
line indicates on-shell cuts. Mirror diagrams for (a) and (c) are understood
and are included in expressions given in the text.

The two contributing real-emission diagrams in Feynman gauge are given in
figure -178 (a) and (b), where $X=g$ is a single gluon. Including the mirror
diagram for (a), their contributions to $d_{1\,Q/Q}^{\mathrm{bare}\,(1)}$
evaluate to

$\displaystyle d_{1\,Q/Q}^{\mathrm{bare}\,(a)}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{e^{\epsilon\gamma_{E}}}{\Gamma(1-\epsilon)}\Bigl{(}\frac{\mu}{k_{T}}\Bigr{)}^{2\epsilon}\Bigl{(}\frac{\sqrt{\zeta}}{\nu}\Bigr{)}^{-\eta}\frac{1}{\pi
z^{2-2\epsilon}}\frac{z^{\eta}}{\bar{z}^{1+\eta}}\frac{4z^{3}}{k_{T}^{2}z^{2}+m^{2}\bar{z}^{2}}\,,$
(3.11) $\displaystyle d_{1\,Q/Q}^{\mathrm{bare}\,(b)}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{e^{\epsilon\gamma_{E}}}{\Gamma(1-\epsilon)}\,\Bigl{(}\frac{\mu}{k_{T}}\Bigr{)}^{2\epsilon}\,\frac{1}{\pi
z^{2-2\epsilon}}\,2z^{2}\bar{z}$
$\displaystyle\quad\times\Bigl{[}\frac{k_{T}^{2}z^{2}+m^{2}(1-4z+z^{2})-\epsilon(k_{T}^{2}z^{2}+m^{2}\bar{z}^{2})}{(k_{T}^{2}z^{2}+m^{2}\bar{z}^{2})^{2}}\Bigr{]}\,,$
(3.12)

where we defined the shorthand $\bar{z}\equiv 1-z$ and
$\alpha_{s}=\alpha_{s}(\mu)$ is the renormalized $\overline{\text{MS}}$
coupling. Note that to our working order, all results in this paper are
independent of whether they are expressed in terms of $\alpha_{s}^{(n_{f})}$
or $\alpha_{s}^{(n_{\ell})}$. We suppress overall factors of $\Theta(z)$,
$\Theta(\bar{z})$, or $\Theta(k_{T}^{2})$ from the final-state on-shell
conditions in the following, but stress that they are formally present in all
the distributional identities we use or derive, and also multiply all finite
remainder terms.

The two seemingly simple expressions above feature an intricate interplay of
singularities as $k_{T}\to 0$ and/or $z\to 1$ ($\bar{z}\to 0$), which in
particular arise as the quark propagator goes on shell for
$k_{T}^{2}z^{2}+m^{2}\bar{z}^{2}\to 0$. At this step, the singularities are
regulated by powers of $\bar{z}^{-\eta}$, dimensional regularization, or both.
However, if we restrict to finite $z<1$ and $k_{T}>0$, the limits $\eta\to 0$
and $\epsilon\to 0$ can be taken right away and all other contributions in the
renormalized one-loop formula including the virtual diagrams drop out in eq.
(3.1), leaving behind a finite result from the sum of (a) and (b),

$\displaystyle z<1~{}\text{and}~{}k_{T}>0\,:\quad
d_{1\,Q/Q}^{(1)}(z,k_{T},\mu,\zeta)=\frac{\alpha_{s}C_{F}}{4\pi}\frac{2}{\pi}\frac{k_{T}^{2}z^{2}(1+z^{2})+m^{2}\bar{z}^{4}}{\bar{z}(k_{T}^{2}z^{2}+m^{2}\bar{z}^{2})^{2}}\,.$
(3.13)

This is, of course, exactly the differential splitting probability for
$Q\bar{Q}\to(Qg)\bar{Q}$ in the quasi-collinear limit that we considered in
eq. (1.1), but now the field-theoretic definition of $d_{1\,Q/Q}$ provides us
with explicit regulators controlling the singular limits. It is in this sense
that the heavy-quark TMD FF, after including the virtual diagram and
performing the renormalization, will provide a fully differential extension of
the splitting probability into the unresolved limit(s) that is finite in
$d=4$, universal, and embedded in factorization theorems.

To proceed, we first expand the Wilson line diagram (a) in the rapidity
regulator $\eta$, using the standard identity

$\displaystyle\frac{z^{\eta}}{\bar{z}^{1+\eta}}=-\frac{\delta(\bar{z})}{\eta}+\mathcal{L}_{0}(\bar{z})+\mathcal{O}(\eta)\,.$
(3.14)

This expansion is valid here because the remaining terms in eq. (3.11) are
still dimensionally regulated at $\bar{z}=0$ as $k_{T}\to 0$ i.e., the
remainder of the expression is a good test function for general $\epsilon$.
Since the mass enters diagram (a) only through the combination $m\bar{z}$, the
action of the $\delta(\bar{z})/\eta$ term on the rest of the diagram is to
reduce it to the massless limit. This in particular means that the coefficient
of the $1/\eta$ pole (including its exact $\epsilon$ dependence) is equal to
the one for the massless TMD FF, as it must for it to cancel the pole in the
$\eta$-regulated soft function, which is independent of the mass at this
order. For our implementation of dimensional regularization, the soft function
is given by

$\displaystyle\frac{1}{2}S^{(1)}(k_{T},\epsilon,\eta,\nu)=\frac{\alpha_{s}C_{F}}{4\pi}\frac{e^{\epsilon\gamma_{E}}}{\Gamma(1-\epsilon)}\Bigl{(}\frac{\mu}{k_{T}}\Bigr{)}^{2\epsilon}\frac{1}{\pi
k_{T}^{2}}\Bigl{[}+\frac{4}{\eta}+2\ln\frac{\nu^{2}}{k_{T}^{2}}+\mathcal{O}(\eta)\Bigr{]}\,.$
(3.15)

Combining this with diagram (a), the $\eta\to 0$ limit can be taken, and we
find

$\displaystyle d_{1\,Q/Q}^{\mathrm{bare}\,(a+s)}$
$\displaystyle\equiv\lim_{\eta\to
0}\Bigl{[}d_{1\,Q/Q}^{\mathrm{bare}\,(a)}(z,k_{T},\epsilon,\eta,\zeta/\nu^{2})+\delta(1-z)\,\tfrac{1}{2}S^{(1)}(k_{T},\epsilon,\eta,\nu)\Bigr{]}$
(3.16)
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{e^{\epsilon\gamma_{E}}}{\Gamma(1-\epsilon)}\frac{1}{\pi
z^{2-2\epsilon}}\biggl{\\{}\frac{\mu^{2\epsilon}}{k_{T}^{2+2\epsilon}}\Bigl{[}2\ln\frac{\zeta}{k_{T}^{2}}\delta(\bar{z})+4\mathcal{L}_{0}(\bar{z})\Bigr{]}+\frac{\mu^{2\epsilon}}{m^{2+2\epsilon}}\,f^{(a)}(k_{T}^{2}/m^{2},z,\epsilon)\biggr{\\}}\,.$

Here the first term in curly braces, which is left behind from the
cancellation against the soft function and the same as for the massless case,
is straightforward to expand in $\epsilon$ using identities like eqs. (A.4)
and (A.5). This produces, among others, a term
$\tfrac{1}{\epsilon^{2}}\delta(\bar{z})\delta(k_{T}^{2})$ as well as a term
$\tfrac{1}{\epsilon}\mathcal{L}_{0}(\bar{z})\delta(k_{T}^{2})$. The nontrivial
part is to expand the dimensionless function

$\displaystyle f^{(a)}(x,z,\epsilon)\equiv
x^{-1-\epsilon}\frac{-4+4z-4xz^{2}}{\bar{z}^{2}+xz^{2}}$ (3.17)

in a Laurent series in $\epsilon$ involving plus distributions of $x$ and $z$.
To do so, we recall that plus distributions in one variable $x$ formally arise
in perturbative calculations as

$\displaystyle
f(x,\epsilon)=\bigl{[}f(x,\epsilon)\bigr{]}_{+}+\delta(x)\,F(\epsilon)\,,\qquad
F(\epsilon)=\int_{0}^{1}\\!\mathrm{d}x^{\prime}\,f(x^{\prime},\epsilon)\,,$
(3.18)

where $f(x,\epsilon)\sim x^{-1-\epsilon}+\mathcal{O}(x^{0})$ has a
dimensionally regulated singular limit $x\to 0$ on the first left-hand side.
By contrast, on the right-hand side $f(x,\epsilon)$ may be expanded in
$\epsilon$ since each term in the expansion is now dimensionally regulated,
whereas the explicit poles are typically isolated in the Laurent series of the
total integral $F(\epsilon)$. Of course, this is precisely how relations like
eq. (3.14) are derived.

The identity in eq. (3.18) generalizes to our case of two variables $x$ and
$z$, cf. ref. [56],

$\displaystyle
f(x,z,\epsilon)=\bigl{[}f(x,z,\epsilon)\bigr{]}_{+,+}+\delta(x)\bigl{[}F_{x}(z,\epsilon)\bigr{]}_{+}+\delta(\bar{z})\bigl{[}F_{z}(x,\epsilon)\bigr{]}_{+}+\delta(x)\,\delta(\bar{z})\,F_{xz}(\epsilon)\,,$
(3.19)

where

$\displaystyle F_{x}(z,\epsilon)$
$\displaystyle\equiv\int_{0}^{1}\\!\mathrm{d}z^{\prime}\,f(x,z^{\prime},\epsilon)\,,\qquad
F_{z}(x,\epsilon)\equiv\int_{0}^{1}\\!\mathrm{d}x^{\prime}\,f(x^{\prime},z,\epsilon)\,,$
$\displaystyle\qquad
F_{xz}(\epsilon)\equiv\int_{0}^{1}\\!\mathrm{d}x\,\int_{0}^{1}\\!\mathrm{d}z^{\prime}\,f(x,z^{\prime},\epsilon)$
(3.20)

and the two-dimensional plus bracket $[\dots]_{+,+}$ is defined through its
action on test functions $g(x,z)$ as

$\displaystyle\int_{0}^{1}\\!\mathrm{d}x\int_{0}^{1}\\!\mathrm{d}z\,\bigl{[}f(x,z,\epsilon)\bigr{]}_{+,+}\,g(x,z)$
$\displaystyle\equiv\int_{0}^{1}\\!\mathrm{d}x\int_{0}^{1}\\!\mathrm{d}z\,f(x,z,\epsilon)\,\bigl{[}g(x,z)-g(0,z)-g(x,0)+g(0,0)\bigr{]}\,.$
(3.21)

Note that in our case with $x=k_{T}^{2}/m^{2}$, the integration domain can
also extend to $x_{*}>x>1$, in which case the $[\dots]_{+,+}$ bracket reduces
to a one-dimensional one,

$\displaystyle\int_{1}^{x_{*}}\\!\mathrm{d}x\int_{0}^{1}\\!\mathrm{d}z\,\bigl{[}f(x,z,\epsilon)\bigr{]}_{+,+}\,g(x,z)\equiv\int_{1}^{x_{*}}\\!\mathrm{d}x\int_{0}^{1}\\!\mathrm{d}z\,f(x,z,\epsilon)\,\bigl{[}g(x,z)-g(x,0)\bigr{]}\,,$
(3.22)

while integration domains with other shapes can be implemented through
suitable test functions. So far, eq. (3.19) together with eq. (3.1) is a
tautology; however, as in the one-dimensional case in eq. (3.18), the crucial
step is that the expansion in $\epsilon$ may now be performed under each of
the various plus brackets.

For the case of $f=f^{(a)}$ as defined in eq. (3.17), the various integrals in
eq. (3.1) are all straightforward to evaluate for general $\epsilon$ in terms
of incomplete Beta and digamma functions. If we restrict to their expansions
in $\epsilon$ for brevity, eq. (3.19) evaluates to

$\displaystyle f^{(a)}(x,z,\epsilon)$
$\displaystyle=\bigl{[}f^{(a)}(x,z,0)\bigr{]}_{+,+}+\delta(x)\biggl{\\{}\frac{4}{\epsilon}\mathcal{L}_{0}(\bar{z})+4\biggl{[}\frac{4z\ln\bigl{(}1+z^{2}/\bar{z}^{2}\bigr{)}}{\bar{z}}\biggr{]}_{+}\biggr{\\}}$
$\displaystyle\quad+\delta(\bar{z})\biggl{[}\frac{2(1-x)\ln
x-4(\pi\sqrt{x}+x+x^{2})}{x(1+x)^{2}}\biggr{]}_{+}$
$\displaystyle\quad+\delta(\bar{z})\delta(x)\biggl{\\{}-\frac{2}{\epsilon^{2}}-2\pi-\frac{5\pi^{2}}{6}\biggr{\\}}+\mathcal{O}(\epsilon)\,.$
(3.23)

In this form, the double pole $1/\epsilon^{2}$ is isolated at $\bar{z}=x=0$.
The term proportional to $\tfrac{1}{\epsilon}\mathcal{L}_{0}(\bar{z})$ along
the $x=0$ boundary ends up canceling a corresponding term from the $\epsilon$
expansion of the first term in curly braces on the second line of eq. (3.16),
such that $d_{1\,Q/Q}^{\mathrm{bare}\,(a+s)}$ is finite as $\epsilon\to 0$ for
$z<1$ for any $x\geq 0$. The double pole
$\tfrac{1}{\epsilon}\delta(\bar{z})\delta(k_{T}^{2})$ similarly cancels.

The second diagram in figure -178 does not require rapidity regularization. In
order to expand it in the dimensional regulator, we repeat the above
procedure, where we define

$\displaystyle d_{1\,Q/Q}^{\mathrm{bare}\,(b)}$
$\displaystyle\equiv\frac{\alpha_{s}C_{F}}{4\pi}\frac{e^{\epsilon\gamma_{E}}}{\Gamma(1-\epsilon)}\frac{1}{\pi
z^{2-2\epsilon}}\frac{\mu^{2\epsilon}}{m^{2+2\epsilon}}\,f^{(b)}(k_{T}^{2}/m^{2},z,\epsilon)\,,$
$\displaystyle f^{(b)}(x,z,\epsilon)$
$\displaystyle=2z^{2}\bar{z}\,x^{-2\epsilon}\frac{xz^{2}+(1-4z+z^{2})-\epsilon(xz^{2}+\bar{z}^{2})}{(xz^{2}+\bar{z}^{2})^{2}}\,.$
(3.24)

The integrals in eq. (3.1) for $f=f^{(b)}$ can again be performed including
their exact dependence on $\epsilon$, leading to digamma and ${}_{2}F_{1}$
hypergeometric functions. Expanding the latter in $\epsilon$ using the HypExp
2.0 package [57], we find the following distributional expansion of the
dimensionless function $f^{(b)}$,

$\displaystyle f^{(b)}(x,z,\epsilon)$
$\displaystyle=\bigl{[}f^{(b)}(x,z,0)\bigr{]}_{+,+}+\delta(x)\biggl{[}\frac{4z^{3}}{\bar{z}(2\bar{z}z-1)}+2\bar{z}\ln\Bigl{(}1+\frac{z^{2}}{\bar{z}^{2}}\Bigr{)}\biggr{]}_{+}$
(3.25) $\displaystyle\quad+\delta(\bar{z})\biggl{[}\frac{x(5+x)\ln
x+(-1+\pi\sqrt{x}-x)(2-3x-3x^{2})}{x(1+x)^{3}}\biggr{]}_{+}$ (3.26)
$\displaystyle\quad+\delta(\bar{z})\delta(x)\biggl{\\{}\frac{2}{\epsilon}+1+\frac{3\pi}{2}\biggr{\\}}+\mathcal{O}(\epsilon)\,.$
(3.27)

Combining eqs. (3.1) and (3.25) produces an explicit result for the sum of
real-emission diagrams expanded in $\epsilon$. The finite terms involving the
one and two-dimensional plus distributions are lengthy and reappear in our
renormalized result below, so we only quote the poles for reference:

$\displaystyle
d_{1\,Q/Q}^{\mathrm{bare}\,(a+s)}+d_{1\,Q/Q}^{\mathrm{bare}\,(b)}=\frac{\alpha_{s}C_{F}}{4\pi}\frac{1}{\pi
z^{2}}\,\biggl{\\{}\delta(\bar{z})\delta(k_{T}^{2})\Bigl{(}-\frac{2}{\epsilon}\ln\frac{\zeta}{m^{2}}+\frac{2}{\epsilon}\Bigr{)}+\mathcal{O}(\epsilon^{0})\biggr{\\}}\,.$
(3.28)

###  3.2 Virtual contributions and UV renormalization

The virtual diagram (c) from figure -178 is common to any collinear matrix
element involving massive fermions and can be evaluated in a straightforward
way using integration by contours for the $\ell^{+}$ component of the loop
momentum $\ell$, Feynman parametrization for the $\ell_{\perp}$ integral, and
an analytic $\ell^{-}$ integral at the end. Including the mirror diagram, the
result reads

$\displaystyle
d_{1\,Q/Q}^{\mathrm{bare}\,(c)}=\frac{\alpha_{s}C_{F}}{4\pi}\,\delta(1-z)\frac{1}{\pi}\delta(k_{T}^{2})\,\frac{2e^{\epsilon\gamma_{E}}\Gamma(\epsilon)}{\epsilon(1-2\epsilon)}\frac{\mu^{2\epsilon}}{m^{2\epsilon}}\,,$
(3.29)

and contributes a double $1/\epsilon^{2}$ pole. So far our calculation was
performed in terms of bare quark fields, and we still have to account for eq.
(2.12). Using the standard on-shell renormalization factor for a quark field
with pole mass $m$ in dimensional regularization,

$\displaystyle Z_{\psi,\mathrm{OS}}\equiv
1+Z_{\psi,\mathrm{OS}}^{(1)}+\mathcal{O}(\alpha_{s}^{2})\,,\qquad
Z_{\psi,\mathrm{OS}}^{(1)}=-\,\frac{\alpha_{s}C_{F}}{4\pi}\Big{(}\frac{3}{\epsilon}+3\ln\frac{\mu^{2}}{m^{2}}+4\Big{)}$
(3.30)

the rapidity-renormalized collinear correlator on the second line of eq. (3.1)
is given by

$\displaystyle\lim_{\eta\to
0}\Bigl{[}d_{1\,Q/Q}^{\mathrm{bare}\,(1)}(z,k_{T},\epsilon,\eta,\zeta/\nu^{2})+\delta(1-z)\,\tfrac{1}{2}S^{(1)}(k_{T},\epsilon,\eta,\nu)\Bigr{]}$
$\displaystyle=d_{1\,Q/Q}^{\mathrm{bare}\,(a+s)}+d_{1\,Q/Q}^{\mathrm{bare}\,(b)}+d_{1\,Q/Q}^{\mathrm{bare}\,(c)}+\delta(1-z)\frac{1}{\pi}\delta(k_{T}^{2})Z_{\psi}^{(1)}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{1}{\pi
z^{2}}\,\biggl{\\{}\delta(\bar{z})\delta(k_{T}^{2})\Bigl{(}\frac{2}{\epsilon^{2}}+\frac{3}{\epsilon}-\frac{2}{\epsilon}\ln\frac{\zeta}{\mu^{2}}\Bigr{)}+\mathcal{O}(\epsilon^{0})\biggr{\\}}\,,$
(3.31)

and we continue to only quote the poles for brevity. The
$\overline{\text{MS}}$ UV renormalization factor for quark TMD PDFs and FFs is
likewise well known and independent of the mass [58],

$\displaystyle
Z_{\mathrm{UV}}^{(1)}(\mu,\zeta,\epsilon)=-\frac{\alpha_{s}C_{F}}{4\pi}\biggl{\\{}\frac{2}{\epsilon^{2}}+\frac{3}{\epsilon}-\frac{2}{\epsilon}\ln\frac{\zeta}{\mu^{2}}\biggr{\\}}\,.$
(3.32)

It precisely cancels the poles from eq. (3.2) when inserting both into eq.
(3.1), as expected. We thus confirm by an explicit one-loop calculation that
the TMD FF for heavy quarks obeys standard TMD evolution equations, with
secondary mass effects in the Collins-Soper kernel at higher loop orders
understood [48].

Taking the $\epsilon\to 0$ limit, we obtain our final result for the
$\mathcal{O}(\alpha_{s})$ correction to the renormalized partonic heavy-quark
TMD FF:

$\displaystyle d_{1\,Q/Q}^{(1)}(z,k_{T},\mu,\zeta)$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{1}{\pi
z^{2}}\,\biggl{\\{}\delta(\bar{z})\Bigl{[}2\ln\frac{\zeta}{\mu^{2}}\mathcal{L}_{0}(k_{T}^{2},\mu^{2})+2\ln\frac{\mu^{2}}{m^{2}}\mathcal{L}_{0}(k_{T}^{2},m^{2})-\ln^{2}\frac{\mu^{2}}{m^{2}}\delta(k_{T}^{2})+3\ln\frac{\mu^{2}}{m^{2}}\delta(k_{T}^{2})\Bigr{]}$
$\displaystyle\qquad+\frac{1}{m^{2}}\biggl{[}\frac{2xz^{4}(1+z^{2})+2z^{2}\bar{z}^{4}}{\bar{z}(xz^{2}+\bar{z}^{2})^{2}}\biggr{]}_{+,+}+\delta(k_{T}^{2})\biggl{[}\frac{2(1+z^{2})}{\bar{z}}\ln\Bigl{(}1+\frac{z^{2}}{\bar{z}^{2}}\Bigr{)}-\frac{4z^{3}}{\bar{z}(1-2\bar{z}z)}\biggr{]}_{+}$
$\displaystyle\qquad+\delta(\bar{z})\,\frac{1}{m^{2}}\biggl{[}-\frac{2+3x+4x^{2}+3x^{3}+\pi\sqrt{x}(2+7x+x^{2})+x(1+7x+2x^{2})\ln
x}{x(1+x)^{3}}\biggr{]}_{+}$
$\displaystyle\qquad+\delta(\bar{z})\,\delta(k_{T}^{2})\Bigl{(}5-\frac{\pi}{2}-\frac{2\pi^{2}}{3}\Bigr{)}\biggr{\\}}\,.$
(3.33)

Here we still use the dimensionless variable $x\equiv k_{T}^{2}/m^{2}$. To
simplify the result we have combined terms under $+,+$ and $+$ brackets as
much as possible, exploiting their linearity, reexpressed plus distributions
$\mathcal{L}_{n}(k_{T}^{2},\mu^{2})$ in terms of
$\mathcal{L}_{n}(k_{T}^{2},m^{2})$ and logarithms of $\mu^{2}/m^{2}$, and used
that products of one-dimensional plus distributions can be written as
(trivial) two-dimensional ones, i.e.,

$\displaystyle\mathcal{L}_{0}(\bar{z})\,\mathcal{L}_{0}(k_{T}^{2},m^{2})=\mathcal{L}_{0}(\bar{z})\,\frac{1}{m^{2}}\mathcal{L}_{0}(k_{T}^{2}/m^{2})=\frac{1}{m^{2}}\Bigl{[}\frac{1}{x}\frac{1}{\bar{z}}\Bigr{]}_{+,+}\,.$
(3.34)

Eq. (3.2) is the main result of this section. Up to the terms on the first
line of the right-hand side, which are predicted by the RGE and vanish for
$\zeta=\mu^{2}=m^{2}$, we have cast the heavy-quark TMD FF precisely in the
form of eq. (3.19), i.e., as a two-dimensional plus distribution with
nontrivial functions of $z$ and $k_{T}$ living on the respective boundaries
and an overall boundary contribution at $z=1$ and $k_{T}=0$. This structure
was already illustrated in figure -179, where the partial integrals over one
of the variables can immediately be read off from the respective boundary
terms. Of course, the total contribution in the bulk at $z<1$ and $k_{T}>0$,
which can simply be read off from the content of the $[\dots]_{+,+}$ brackets
accounting for the prefactor of $1/z^{2}$, is still equal to the simple result
we found in eq. (3.13) and thus to the quasi-collinear splitting probability
in eq. (1.1).

###  3.3 Results in cumulant space

Fixed-order subtraction methods, as well as certain $q_{T}$ resummation
formalisms like the RadISH approach [59], require the cumulative distribution
in $k_{T}$ integrated over transverse momenta $\vec{k}_{T}$ with
$\lvert\vec{k}_{T}\rvert\leq{k_{T}^{\mathrm{cut}}}$ as an input,

$\displaystyle
d_{1\,Q/Q}(z,{k_{T}^{\mathrm{cut}}},\mu,\zeta)\equiv\pi\\!\int^{({k_{T}^{\mathrm{cut}}})^{2}}\\!\mathrm{d}(k_{T}^{2})\,d_{1\,Q/Q}(z,k_{T},\mu,\zeta)$
(3.35)

often simply referred to as the cumulant. As for the Fourier transform we will
indicate the cumulative distribution simply through its argument
${k_{T}^{\mathrm{cut}}}$.

The organization of our final result in eq. (3.2) in terms of plus
distributions defined on rectangular domains in $(z,k_{T}^{2})$ makes it very
straightforward to compute the cumulative distribution because the cumulant
constitutes a special case of eq. (3.1) for a unit test function.
Specifically, for plus-bracketed functions $f$ of either a single variable
$x=k_{T}^{2}/m^{2}$ or of both $x$ and $z$, we have the following integrals up
to $x_{\mathrm{cut}}\equiv({k_{T}^{\mathrm{cut}}}/m)^{2}$:

$\displaystyle\int^{x_{\mathrm{cut}}}\\!\mathrm{d}x\,\bigl{[}f(x,z)\bigr{]}_{+,+}$
$\displaystyle=\biggl{[}\int_{1}^{x_{\mathrm{cut}}}\\!f(x,z)\biggr{]}_{+}\,,$
$\displaystyle\int^{x_{\mathrm{cut}}}\\!\mathrm{d}x\,\bigl{[}f(x)\bigr{]}_{+}$
$\displaystyle=\int_{1}^{x_{\mathrm{cut}}}\\!f(x)\,,$ (3.36)

where the remaining one-dimensional plus bracket on the right-hand side of the
first line refers to the $z$ dependence of the integral inside it, while the
latter’s $x_{\mathrm{cut}}$ dependence is function valued.

The integrals in eq. (3.3) are straightforward to evaluate for the respective
terms in eq. (3.2). For purposes of our presentation, we find it convenient to
further decompose the $z$ plus distribution arising from the cumulant of the
$[\dots]_{+,+}$ term as

$\displaystyle\bigl{[}s(x_{\mathrm{cut}},z)+r(x_{\mathrm{cut}},z)\bigr{]}_{+}=\bigl{[}s(x_{\mathrm{cut}},z)\bigr{]}+r(x_{\mathrm{cut}},z)-\delta(\bar{z})\int_{0}^{1}\\!\mathrm{d}z^{\prime}\,r(x_{\mathrm{cut}},z^{\prime})\,,$
(3.37)

where $s(x_{\mathrm{cut}},z)$ is defined as the leading term under the
brackets for $\bar{z}\ll x_{\mathrm{cut}}\sim 1$, with homogeneous $\sim
1/\bar{z}$ scaling, while $r(x_{\mathrm{cut}},z)$ is defined as the remainder
that is at most logarithmically divergent as $z\to 1$. Notably, we find that
in this operation, the $x_{\mathrm{cut}}$ dependence of $s(z)$ is purely
logarithmic, while the power-like $x_{\mathrm{cut}}$ dependence of the
integrated $r(x_{\mathrm{cut}},z^{\prime})$ multiplying $\delta(\bar{z})$
precisely compensates that of the complicated cumulant integral of the one-
dimensional plus distribution in $x$ in eq. (3.2).

Collecting terms, our final result for the NLO correction to the renormalized
heavy-quark TMD FF in cumulant space reads

$\displaystyle d_{1\,Q/Q}^{\,(1)}(z,{k_{T}^{\mathrm{cut}}},\mu,\zeta)$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{1}{z^{2}}\bigg{\\{}\delta(1-z)\bigg{[}-4\ln\frac{\mu}{{k_{T}^{\mathrm{cut}}}}\ln\frac{\zeta}{m^{2}}+4\ln^{2}\frac{\mu}{m}+6\ln\frac{\mu}{m}-4\ln^{2}\frac{{k_{T}^{\mathrm{cut}}}}{m}+4-\frac{\pi^{2}}{6}\bigg{]}$
$\displaystyle\qquad\qquad\qquad-4\Bigl{(}1-2\ln\frac{{k_{T}^{\mathrm{cut}}}}{m}\Bigr{)}\,\mathcal{L}_{0}(1-z)-8\mathcal{L}_{1}(1-z)+\mathcal{R}\Bigl{[}z,\Bigl{(}\frac{{k_{T}^{\mathrm{cut}}}}{m}\Bigr{)}^{\\!2\,}\Bigr{]}\bigg{\\}}\,,$
(3.38)

where we defined

$\displaystyle\mathcal{R}(z,x_{\mathrm{cut}})\equiv$
$\displaystyle\,\,\frac{4-4z(1-x_{\mathrm{cut}}z)}{\bar{z}^{2}+x_{\mathrm{cut}}z^{2}}+4\ln\frac{\bar{z}^{2}}{x_{\mathrm{cut}}}+\frac{1}{\bar{z}}\Bigl{[}2(1+z^{2})\ln\Bigl{(}1+\frac{x_{\mathrm{cut}}z^{2}}{\bar{z}^{2}}\Bigr{)}-4z\ln\frac{x_{\mathrm{cut}}}{\bar{z}^{2}}\Bigr{]}\,.$
(3.39)

We note that $\mathcal{R}(z,x_{\mathrm{cut}})=\mathcal{O}(\bar{z}^{0})$ for
$z\to 1$ and also $\mathcal{R}(z,x_{\mathrm{cut}})=\mathcal{O}(z^{0})$ for
$z\to 0$, so it is integrable in all limits. For later reference we note its
total integral,

$\displaystyle\int_{0}^{1}\\!\mathrm{d}z\,\mathcal{R}(z,x_{\mathrm{cut}})$
$\displaystyle=-4\mathrm{Li}_{2}\Bigl{(}1+\frac{1}{\mathrm{i}\sqrt{x_{\mathrm{cut}}}}\Bigr{)}-4\mathrm{Li}_{2}\Bigl{(}1-\frac{1}{\mathrm{i}\sqrt{x_{\mathrm{cut}}}}\Bigr{)}$
$\displaystyle\quad-\frac{3x_{\mathrm{cut}}^{2}+3x_{\mathrm{cut}}+2}{(1+x_{\mathrm{cut}})^{2}}\ln
x_{\mathrm{cut}}+\frac{2x_{\mathrm{cut}}(x_{\mathrm{cut}}-\pi\sqrt{x_{\mathrm{cut}}}+1)}{(1+x_{\mathrm{cut}})^{2}}\,.$
(3.40)

To check eq. (3.3), we have performed an independent calculation of the bare
collinear correlator in dimensional regularization and renormalized it
directly in cumulant space, This approach involves more complicated
hypergeometric functions of combinations of $x_{\mathrm{cut}}$ and $z$ in
intermediate steps; in particular, the challenge in this case is that their
expansion in $\epsilon$ produces terms that are singular as $z\to 1$, even
after subtracting terms that scale as $1/\bar{z}^{1+\epsilon}$ for general
$\epsilon$. This means that the full hypergeometric result in dimensional
regularization needs to be treated as in eq. (3.18) before expanding in
$\epsilon$. The final result is in full agreement with eq. (3.3).

Eq. (3.3) has the – at first sight – remarkable property that the dependence
of the coefficient of $\delta(1-z)$ on ${k_{T}^{\mathrm{cut}}}/m$ is purely
logarithmic, while a power-like dependence on ${k_{T}^{\mathrm{cut}}}/m$ would
in general be allowed by our power counting assumption
${k_{T}^{\mathrm{cut}}}\sim m$ for the heavy-quark TMD FF. Instead, after the
cancellations discussed below eq. (3.37) have taken place, the nontrivial
power-like dependence on ${k_{T}^{\mathrm{cut}}}/m$ arises solely in the
$\mathcal{R}(z,x_{\mathrm{cut}})$ term, which is regular as $z\to 1$. In
section 6 we will explain this fact from the perspective of joint threshold
and $k_{T}$ factorization, and give all-order expressions for the $z\to 1$
limit of $d_{1\,Q/Q}$ at fixed $k_{T}/m\sim b_{T}m\sim 1$.

###  3.4 Results in position space

We now move on to position-space results in terms of the $\vec{b}_{T}$
variable Fourier-conjugate to $\vec{k}_{T}$, which are the key input for
producing resummed predictions from solving multiplicative $b_{T}$-space RGEs
and Fourier-transforming the final result. As is standard, the two-dimensional
Fourier transform of the azimuthally symmetric renormalized TMD FF at hand is
simply given by a Fourier-Bessel (or Hankel) transform of zeroth order,

$\displaystyle
d_{1\,Q/Q}(z,b_{T},\mu,\zeta)=2\pi\int_{0}^{\infty}\\!\mathrm{d}k_{T}\,k_{T}\,J_{0}(b_{T}k_{T})\,d_{1\,Q/Q}(z,k_{T},\mu,\zeta)\,,$
(3.41)

where $J_{n}$ is the $n$th-order Bessel function of the first kind. It is
convenient to integrate the above relation by parts,

$\displaystyle
d_{1\,Q/Q}(z,b_{T},\mu,\zeta)=b_{T}\int_{0}^{\infty}\\!\mathrm{d}{k_{T}^{\mathrm{cut}}}\,J_{1}(b_{T}{k_{T}^{\mathrm{cut}}})\,d_{1\,Q/Q}(z,{k_{T}^{\mathrm{cut}}},\mu,\zeta)\,.$
(3.42)

This lets us compute the Fourier-space result directly from the cumulant in
eq. (3.3) obtained in the previous section, which already had function-valued
(and simpler) ${k_{T}^{\mathrm{cut}}}$ dependence, cf. the more complicated
original momentum-space result in eq. (3.2).

Evaluating the relevant Hankel transforms, our result for the NLO correction
to the renormalized heavy-quark TMD FF in position space reads

$\displaystyle d_{1\,Q/Q}^{\,(1)}(z,b_{T},\mu,\zeta)$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{1}{z^{2}}\bigg{\\{}\delta(1-z)\bigg{[}-2L_{b}\ln\frac{\zeta}{m^{2}}+4\ln^{2}\frac{\mu}{m}+6\ln\frac{\mu}{m}-L_{y}^{2}+4-\frac{\pi^{2}}{6}\bigg{]}$
$\displaystyle\qquad\qquad-4(1+L_{y})\,\mathcal{L}_{0}(1-z)-8\mathcal{L}_{1}(1-z)+\tilde{\mathcal{R}}(z,b_{T}m)\bigg{\\}}\,,$
(3.43)

where we have used the shorthands

$\displaystyle L_{b}$ $\displaystyle\equiv
2\ln\frac{b_{T}\mu}{2e^{-\gamma_{E}}}\,,\qquad y\equiv b_{T}m\,,\quad
L_{y}\equiv 2\ln\frac{y}{2e^{-\gamma_{E}}}\,,$ (3.44)

and the regular term $\tilde{\mathcal{R}}(y)$ is again integrable for $z\to
0,1$ and given by

$\displaystyle\tilde{\mathcal{R}}(z,y)$ $\displaystyle\equiv
y\int_{0}^{\infty}\\!\mathrm{d}t\,J_{1}(ty)\,\mathcal{R}(z,t^{2})$
$\displaystyle=\frac{4}{\bar{z}}\Bigl{[}1+L_{y}+(1+z^{2})K_{0}\Bigl{(}\frac{y\bar{z}}{z}\Bigr{)}-y\bar{z}K_{1}\Bigl{(}\frac{y\bar{z}}{z}\Bigr{)}+2\ln\bar{z}\Bigr{]}\,.$
(3.45)

While we were not able to find a closed-form expression for the total $z$
integral of $\tilde{\mathcal{R}}$, it can readily be obtained numerically or
analytically in asymptotic limits through the Hankel transform of eq. (3.3).

###  3.5 Nonvalence channels at $\mathcal{O}(\alpha_{s})$

(a)  (b)  (c)

Figure -177:  Diagrams contributing to the gluon-initiated partonic heavy-
quark TMD FF and the mass-dependent TMD FF matching coefficients at
$\mathcal{O}(\alpha_{s})$. Heavy quark propagators and field insertions are
indicated in red. The mirror diagrams for (a) and (b) are understood and are
included in expressions given in the text.

For completeness we also compute the partonic TMD FFs (or TMD FF matching
coefficients onto collinear distributions) for all remaining partonic channels
involving heavy quarks. The relevant diagrams at $\mathcal{O}(\alpha_{s})$ are
given in figure -177, where we labeled the identified final-state “hadron” by
the momentum $p$. Note that for diagrams (b) and (c) compared to the otherwise
identical diagrams in figure -178, this also changes the definition of
$k_{\perp}$, as the coordinate system needs to be adjusted such that the gluon
now has vanishing transverse momentum. Results for antiquarks are identical at
this order. For the $g\to Q$ contribution, diagram (a) evaluates to

$\displaystyle
d_{1\,Q/g}^{(1)}(z,k_{T})=\frac{\alpha_{s}}{4\pi}T_{F}\frac{2}{\pi
z^{2}}\frac{k_{T}^{2}z^{4}\bigl{(}2z^{2}-2z+1\bigr{)}+m^{2}z^{2}}{\bigl{(}k_{T}^{2}z^{2}+m^{2}\bigr{)}^{2}}\,.$
(3.46)

For gluon production off the heavy quark, we find after summing over diagrams
(b) and (c) and including the mirror diagram for (b),

$\displaystyle\mathcal{J}_{g/Q}^{(1)}(z,k_{T},m)=\frac{\alpha_{s}}{4\pi}C_{F}\frac{2}{\pi}\frac{k_{T}^{2}\bigl{(}z^{2}-2z+2\bigr{)}+m^{2}z^{2}}{z\bigl{(}k_{T}^{2}+m^{2}\bigr{)}^{2}}\,.$
(3.47)

It is straightforward to verify from eq. (3.13) that at $0<z<1$ and $k_{T}>0$
the above expression satisfies
$1/z^{2}\mathcal{J}_{g/Q}^{(1)}(z,k_{T},m)=d_{1Q/Q}^{(1)}(1-z,k_{T}\frac{z}{1-z})$,
as expected from $z\leftrightarrow 1-z$ and the simultaneous change of
coordinate system; similarly,
$d_{1\,Q/g}^{(1)}(z,k_{T})=d_{1Q/g}^{(1)}(1-z,k_{T}\frac{z}{1-z})$. Note that
the bare results in eqs. (3.46) and (3.47) are finite at $z>0$ in $d=4$ and
without a rapidity regulator on their own, as expected because there is no
cross term with the soft function or the UV renormalization at this order.
This is in particular the case as $k_{T}\to 0$, where all singularities in
$k_{T}$ are now cut off by the mass for any value of $z>0$.

The effect of the heavy quark on the gluon-onto-gluon matching coefficient at
one loop can be written as

$\displaystyle\mathcal{J}_{g/g}^{(1)}(z,k_{T},m,\mu,\zeta)\equiv\mathcal{J}_{g/g}^{(1,n_{\ell})}(z,k_{T},\mu,\zeta)+\Delta\mathcal{J}_{g/g}^{(1,h)}(z,k_{T},m)\,,$
(3.48)

and arises purely from its contribution to the gluon wave function
renormalization; here $\mathcal{J}_{g/g}^{(1,n_{\ell})}(z,k_{T})$ is the well-
known NLO matching coefficient for the gluon TMD onto the collinear gluon FF
in the case where both are defined in a purely light theory with $n_{\ell}$
flavors, see eq. (B). After $\overline{\text{MS}}$-renormalizing the matching
coefficient, the net result is

$\displaystyle\Delta\mathcal{J}_{g/g}^{(1)}(z,k_{T},m,\mu)=\frac{\alpha_{s}T_{F}}{4\pi}\,\delta(1-z)\frac{1}{\pi}\delta(k_{T}^{2})\,\frac{4}{3}\ln\frac{m^{2}}{\mu^{2}}\,.$
(3.49)

This agrees with ref. [60], where the same contribution was recently revisited
as part of computing all secondary heavy-quark effects in the two-loop gluon
TMD PDF; we refer to that reference for a detailed discussion of how it
ensures the proper TMD renormalization properties when connecting the theories
with $n_{\ell}$ and $n_{\ell}+1$ quark degrees of freedom.

Taking cumulant integrals, we find

$\displaystyle d_{1\,Q/g}^{(1)}(z,{k_{T}^{\mathrm{cut}}})$
$\displaystyle=\frac{\alpha_{s}}{4\pi}T_{F}\frac{2}{z^{2}}\biggl{\\{}(2z^{2}-2z+1)\ln(x_{\mathrm{cut}}z^{2}+1)+\frac{2x_{\mathrm{cut}}(1-z)z^{3}}{{x_{\mathrm{cut}}}z^{2}+1}\biggr{\\}}\,,$
$\displaystyle\mathcal{J}_{g/Q}^{(1)}(z,{k_{T}^{\mathrm{cut}}},m)$
$\displaystyle=\frac{\alpha_{s}}{4\pi}C_{F}\frac{2}{z}\biggl{\\{}(z^{2}-2z+2)\ln(x_{\mathrm{cut}}+1)+\frac{2x_{\mathrm{cut}}(z-1)}{1+x_{\mathrm{cut}}}\biggr{\\}}\,,$
(3.50)

where $x_{\mathrm{cut}}\equiv({k_{T}^{\mathrm{cut}}}/m)^{2}$. Finally, taking
a $J_{1}$ Hankel transform of the cumulant results as in eq. (3.42), we arrive
at the following $b_{T}$-space results:

$\displaystyle d_{1\,Q/g}^{(1)}(z,b_{T})=$
$\displaystyle\frac{\alpha_{s}}{4\pi}T_{F}\frac{4}{z^{2}}\biggl{\\{}(2z^{2}-2z+1)K_{0}\Bigl{(}\frac{b_{T}m}{z}\Bigr{)}+b_{T}m(1-z)K_{1}\Bigl{(}\frac{b_{T}m}{z}\Bigr{)}\biggr{\\}}\,,$
$\displaystyle\mathcal{J}_{g/Q}^{(1)}(z,b_{T},m)$
$\displaystyle=\frac{\alpha_{s}}{4\pi}C_{F}\frac{4}{z}\bigl{\\{}(z^{2}-2z+2)K_{0}(b_{T}m)-b_{T}m(1-z)K_{1}(b_{T}m)\bigr{\\}}\,.$
(3.51)

##  4 bHQET TMD fragmentation factors at NLO

###  4.1 All-order renormalization properties

In ref. [1] we showed that up to a normalization factor $N_{H}$ involving
Clebsch-Gordan coefficients, the unpolarized bHQET TMD fragmentation factor
$\chi_{1,H}$ in eq. (2.18) for a heavy hadron state $H$ is given by the trace
of the spin density matrix $\rho_{\ell}$ encoding the dynamics of the light
constituents $\ell$ of the hadron,

$\displaystyle\chi_{1,H}(b_{T},\mu,\rho)=N_{H}\sum_{h_{\ell}}\rho_{\ell,h_{\ell}h_{\ell}}(b_{\perp},\mu,\rho)\,.$
(4.1)

Here we use the shorthand $\rho\equiv\bar{n}\cdot v$ as introduced below eq.
(2.19), which can be interpreted as the total boost of the fragmenting heavy-
quark system with velocity $v$ in the frame where the Wilson line direction
$\bar{n}$ has unit spatial component. A priori, the matrix element is
invariant under boosts – or reparametrization invariance type-III
transformations (RPI-III) in SCET [61, 62] – that would take
$\bar{n}\to\alpha\,\bar{n}$ for an arbitrary factor $\alpha$, but similarly to
the case of relativistic TMD FFs, the presence of rapidity divergences breaks
this invariance such that $\chi_{1,H}$ and $\rho_{\ell}$ develop an anomalous
dependence on $\rho$. In this subsection, we make the discussion of this point
in ref. [1] fully explicit, allowing us to cross check our predictions for the
renormalization properties at one loop in section 4.3.

Spelling out all regulators and its open light helicity indices
$h_{\ell},h_{\ell}^{\prime}$, the renormalized light-spin density matrix is
defined as

$\displaystyle\rho_{\ell,h_{\ell}h_{\ell}^{\prime}}(b_{\perp},\mu,\rho)=\lim_{\epsilon\to
0}Z_{\chi_{1}}^{-1}(\mu,\rho,\epsilon)\lim_{\eta\to
0}\Bigl{[}\rho_{\ell,h_{\ell}h_{\ell}^{\prime}}^{\mathrm{bare}}\Bigl{(}b_{\perp},\epsilon,\eta,\frac{\rho}{\nu}\Bigr{)}\,\sqrt{S^{(n_{\ell})}}(b_{T},\epsilon,\eta,\nu)\Bigr{]}\,,$
(4.2)

where $S^{(n_{\ell})}$ is the bare TMD soft function in the theory where only
the $n_{\ell}$ light quark degrees of freedom are active. This is the same
soft function as appears in light-quark TMD PDFs or FFs at transverse momenta
below the heavy-quark threshold: The fact that it also arises here and cancels
the rapidity divergences as $\eta\to 0$ in the bare light-spin density matrix
follows from consistency of the TMD factorization for flavor-changing hard
processes, e.g. charged-current DIS $e^{-}d\to\nu_{e}\,c$.111 We assume as in
section 2.1 that a rapidity regulator $\eta$ is used for which soft (zero-bin)
subtractions in the bare collinear matrix element are scaleless. If this is
not the case, it is crucial that zero-bin subtractions are defined by modes
that are soft in a frame where the heavy quark is highly boosted, i.e., their
interactions with the heavy-quark sector are encoded in lightlike Wilson lines
and the multiplicative zero-bin subtraction factor is obtained by taking
$\bar{n}\cdot v\to\infty$.  The UV $\overline{\text{MS}}$ renormalization
factor $Z_{\chi_{1}}$ likewise follows from consistency and is given by

$\displaystyle Z_{\chi_{1}}^{-1}\bigl{(}\mu,\sqrt{\zeta}/m,\epsilon\bigr{)}$
$\displaystyle=Z_{\mathrm{UV}}(\mu,\zeta,\epsilon)\,Z_{C_{m}}(m,\mu,\zeta,\epsilon)$
$\displaystyle=1-\frac{\alpha_{s}C_{F}}{4\pi}\frac{1}{\epsilon}\Bigl{(}-4\ln\frac{\sqrt{\zeta}}{m}+2\Bigr{)}+\mathcal{O}(\alpha_{s}^{2})\,,$
(4.3)

where $Z_{\mathrm{UV}}$ is the standard UV TMD renormalization factor in eqs.
(2.14) and (3.32), and $Z_{C_{m}}$ is the $\overline{\text{MS}}$
renormalization factor of the Wilson coefficient in eq. (2.18) generated by
integrating out the heavy quark [47],

$\displaystyle C_{m}^{\mathrm{bare}}(m,\zeta,\epsilon)$
$\displaystyle=Z_{C_{m}}(m,\mu,\zeta,\epsilon)\,C_{m}(m,\mu,\zeta)\,.$ (4.4)

Both the soft function and the UV renormalization are independent of the
hadronic states, and thus universal for all possible traces of $\rho_{\ell}$,
i.e., they are the same for all possible TMD fragmentation factors
$\chi_{1,H}$ (including polarized ones). In particular, they are the same also
for the total TMD fragmentation factor $\chi_{1}$ that we consider in section
4.2. Note that while $Z_{\mathrm{UV}}$ and $Z_{C_{m}}$ individually depend on
the dimensionful Collins-Soper scale $\zeta=(\bar{n}\cdot P_{H}/z)^{2}$,
consistency requires that only the dimensionless combination
$\rho=\sqrt{\zeta}/m$ appears in $Z_{\chi_{1}}$ on the left-hand side of eq.
(4.1). The fact that $\rho$ is dimensionless also implies that
$Z_{\chi_{1}}(\mu,\rho)$ can only subtract single poles in $\epsilon$, but
this is indeed precisely what one expects for the renormalization of the
product of Wilson lines with a timelike-lightlike cusp [63].

Finally, the bare light-spin density matrix itself is defined as

$\displaystyle\rho_{\ell,h_{\ell}h_{\ell}^{\prime}}^{\mathrm{bare}}\Bigl{(}b_{\perp},\epsilon,\eta,\frac{\bar{n}\cdot
v}{\nu}\Bigr{)}$
$\displaystyle=\frac{1}{N_{c}}\operatorname{Tr}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{X}\big{\langle}0\big{\lvert}W_{\eta}^{\dagger}(b_{\perp})\,Y_{v}(b_{\perp})\big{\lvert}s_{\ell},h_{\ell},f_{\ell};X\big{\rangle}$
(4.13)
$\displaystyle\qquad\qquad\quad\times\big{\langle}s_{\ell},h_{\ell}^{\prime},f_{\ell};X\big{\rvert}Y^{\dagger}_{v}(0)\,W_{\eta}(0)\big{\rvert}0\big{\rangle}\,,$
(4.14)

where $\lvert s_{\ell},h_{\ell},f_{\ell}\rangle$ indicates a light constituent
state $\ell$ of total spin $s_{\ell}$, spin component $h_{\ell}$ along the
direction $z^{\mu}\equiv v^{\mu}-\frac{{\bar{n}}^{\mu}}{{\bar{n}}\cdot v}$,
and valence flavor content $f_{\ell}$. The Wilson lines $Y_{v}$ along the
timelike direction encode the static interactions with the heavy quark,

$\displaystyle Y_{v}(x)\equiv
P\Bigl{[}\exp\Bigl{(}\mathrm{i}g\int_{0}^{\infty}\\!\mathrm{d}s\,v\cdot
A(x+vs)\Bigr{)}\Bigr{]}\,.$ (4.15)

The rapidity regulator only acts on the Wilson lines $W_{\eta}$ along the
lightlike direction $\bar{n}$, which are still defined as in eq. (2.13). Here
the rapidity regulator is the only source of boost invariance (RPI-III)
breaking, which implies that the bare spin density matrix can only depend on
the combination $\rho/\nu=\bar{n}\cdot v/\nu$. A priori, it is not obvious
that $\rho=\sqrt{\zeta}/m$, as it appears in the renormalization factor in eq.
(4.1), should coincide with the $\rho=\bar{n}\cdot v$ that appears in the bare
matrix element, but this is in fact guaranteed by the crucial relation

$\displaystyle\sqrt{\zeta}=\bar{n}\cdot P_{H}/z=m\,\bar{n}\cdot v$ (4.16)

that arises from consistently expanding to leading power in $1/m$ [1].

We conclude this section by summarizing the renormalization group equations
for the renormalized $\rho_{\ell}$ in $d=4$ dimensions that follow from the
above renormalization procedure:

$\displaystyle\mu\frac{\mathrm{d}}{\mathrm{d}\mu}\rho_{\ell,h_{\ell}h_{\ell}^{\prime}}(b_{\perp},\mu,\rho)$
$\displaystyle=\gamma_{\chi_{1}}\bigl{[}\alpha_{s}(\mu),\rho\bigr{]}\,\rho_{\ell,h_{\ell}h_{\ell}^{\prime}}(b_{\perp},\mu,\rho)\,,$
$\displaystyle\rho\frac{\mathrm{d}}{\mathrm{d}\rho}\rho_{\ell,h_{\ell}h_{\ell}^{\prime}}(b_{\perp},\mu,\rho)$
$\displaystyle=\gamma_{\zeta}^{(n_{\ell})}(b_{T},\mu)\,\rho_{\ell,h_{\ell}h_{\ell}^{\prime}}(b_{\perp},\mu,\rho)\,,$
(4.17)

where the $\mu$ anomalous dimension is given by the difference of the standard
TMD UV anomalous dimension $\gamma_{\mu}^{q}$ and the anomalous dimension of
$C_{m}$, evaluating both at $\zeta=\rho^{2}m^{2}$,

$\displaystyle\gamma_{\chi_{1}}\bigl{[}\alpha_{s}(\mu),\rho\bigr{]}=\gamma_{\mu}^{q}(\mu,\rho^{2}m^{2})-\gamma_{C_{m}}(\mu,m^{2},\rho^{2}m^{2})=\frac{\alpha_{s}C_{F}}{4\pi}\Bigl{(}-8\ln\rho+4\Bigr{)}+\mathcal{O}(\alpha_{s}^{2})\,.$
(4.18)

We stress that when extracting the anomalous dimension at higher orders, care
must be taken to perturbatively reexpress $\alpha_{s}^{(n_{\ell}+1)}$ in terms
of $\alpha_{s}(\mu)\equiv\alpha_{s}^{(n_{\ell})}(\mu)$.

From the second line of eq. (4.1), we see that $\rho_{\ell}$ develops an
anomalous dependence on the boost parameter $\rho$ governed by the Collins-
Soper CS kernel $\gamma_{\zeta}^{(n_{\ell})}$ of the theory where only
$n_{\ell}$ light quarks are active. For perturbative $b_{T}$, it is given in
our normalization by

$\displaystyle\gamma_{\zeta}^{(n_{\ell})}(b_{T},\mu)=\frac{\alpha_{s}C_{F}}{4\pi}\,\biggl{(}-8\ln\frac{b_{T}\mu}{2e^{-\gamma_{E}}}\biggr{)}+\mathcal{O}(\alpha_{s}^{2})+\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})\,.$
(4.19)

A natural boundary condition for $\rho_{\ell}$, where it is free of large
rapidity logarithms, is given by $\rho\sim 1$, i.e., the rest frame of the
heavy quark. By contrast, in physical cross sections, the boost is fixed to
$\rho=Q/m\gg 1$, i.e., the lab frame in the $e^{+}e^{-}$ case, or the Breit
frame for SIDIS. The $\rho$-evolution between those scales resums the
associated large rapidity logarithms. As in the relativistic case, the two-
dimensional evolution equations satisfy an integrability condition,

$\displaystyle\rho\frac{\mathrm{d}}{\rho}\gamma_{\chi_{1}}\bigl{[}\alpha_{s}(\mu),\rho\bigr{]}=\mu\frac{\mathrm{d}}{\mu}\gamma_{\zeta}^{(n_{\ell})}(\mu,b_{T})=-2\Gamma_{\mathrm{cusp}}\bigl{[}\alpha_{s}(\mu)\bigr{]}$
(4.20)

which involves the lightlike-lightlike cusp anomalous dimension
$\Gamma_{\mathrm{cusp}}(\alpha_{s})=\tfrac{\alpha_{s}}{4\pi}\,4C_{F}+\mathcal{O}(\alpha_{s}^{2})$
and can readily be verified at one loop.

###  4.2 Sum over hadronic states

An important property of the unpolarized bHQET fragmentation factors is that
their description greatly simplifies when summing over all hadrons $H$
containing the heavy quark,

$\displaystyle\chi_{1}(b_{T},\mu,\rho)=\sum_{H}\chi_{1,H}(b_{T},\mu,\rho)=\sum_{\ell}\sum_{h_{\ell}}\rho_{\ell}(b_{\perp},\mu,\rho)\,.$
(4.21)

As indicated, this sum in turn is equal to the sum of traces of the light-spin
density matrices for all the possible light states $\ell$. Unlike the case of
inclusive fragmentation, $\chi_{1}$ does not simply evaluate to unity by
probability conservation. However, _as_ in the inclusive case, the complete
sum over states can be performed also here, leaving behind a simple vacuum
Wilson line correlator [1],

$\displaystyle\chi_{1}^{\mathrm{bare}}(b_{T},\epsilon,\eta,\bar{n}\cdot
v)\equiv\frac{1}{N_{c}}\operatorname{Tr}\big{\langle}0\big{\lvert}W_{\eta}^{\dagger}(b_{\perp})\,Y_{v}(b_{\perp})\,Y^{\dagger}_{v}(0)\,W_{\eta}(0)\big{\rvert}0\big{\rangle}\,.$
(4.22)

Its renormalization proceeds exactly as for the whole $\rho_{\ell}$ in eq.
(4.2),

$\displaystyle\chi_{1}(b_{T},\mu,\rho)=\lim_{\epsilon\to
0}Z_{\chi_{1}}^{-1}(\mu,\rho,\epsilon)\lim_{\eta\to
0}\Bigl{[}\chi_{1}^{\mathrm{bare}}\Bigl{(}b_{T},\epsilon,\eta,\frac{\rho}{\nu}\Bigr{)}\,\sqrt{S^{(n_{\ell})}}(b_{T},\epsilon,\eta,\nu)\Bigr{]}\,.$
(4.23)

While eqs. (4.22) and (4.23) are valid nonperturbatively and likely exhibit
interesting properties at long distances $b_{T}\sim 1/\Lambda_{\mathrm{QCD}}$,
an immediate application of their simple operator structure is to evaluate the
correlator perturbatively at $b_{T}\ll 1/\Lambda_{\mathrm{QCD}}$, which we do
in the next section. As a byproduct of this, we may then read off the
universal matching coefficient $C_{1}$ in eq. (2.19) by summing eq. (2.19)
over all hadrons $H$ such that $\sum\chi_{H}=1$ drops out,

$\displaystyle\chi_{1}(b_{T},\mu,\rho)=C_{1}(b_{T},\mu,\rho)+\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})\,.$
(4.24)

Of course, one can also determine $C_{1}$ by a standard matching calculation,
e.g. using non-decoupled heavy-quark fields $h_{v}$ in the operators and
heavy-quark external states in both the theory at the scale $\mu\sim 1/b_{T}$
and at the low scale $\mu\sim\Lambda_{\mathrm{QCD}}$ to read off the matching
coefficient, and the procedure (in pure dimensional regularization) involves
equivalent diagrams and leads to the same result.

###  4.3 One-loop calculation

(a)  (b)  Figure -176:  Left: Diagrams contributing to the total bHQET TMD
fragmentation factor in the perturbative domain. Wilson lines along the
lightlike ${\bar{n}}$ direction are indicated by the slanted double lines. The
horizontal double lines indicate the Wilson lines along the timelike $v$
direction. The dashed line indicates on-shell cuts after introducing a
complete set of states. The mirror diagram for (a) is understood and included
in expressions given in the text. Right: Phase-space integrands $\mathcal{A}$
for the two real-emissions diagrams as a function of the gluon emission
rapidity $y$.

To evaluate $\chi_{1}$ perturbatively, it is convenient to reinstate a
complete sum over states and transition to momentum space to evaluate the bare
matrix element,

$\displaystyle\chi_{1}^{\mathrm{bare}}\Bigl{(}k_{T},\epsilon,\eta,\frac{\bar{n}\cdot
v}{\nu}\Bigr{)}=\frac{1}{N_{c}}\operatorname{Tr}\int\\!\frac{\mathrm{d}\Omega_{2-2\epsilon}}{2\pi}\,k_{T}^{-2\epsilon}\big{\langle}0\big{\lvert}W_{\eta}^{\dagger}\,Y_{v}\big{\rvert}X\big{\rangle}\big{\langle}X\big{\lvert}\bigl{[}\delta^{(2-2\epsilon)}\bigl{(}k_{\perp}^{\mu}+\mathrm{i}\partial_{\perp}^{\mu}\bigr{)}\,Y^{\dagger}_{v}\,W_{\eta}\bigr{]}\big{\rvert}0\big{\rangle}\,,$
(4.25)

in analogy to eq. (3.10). Virtual diagrams with $X$ the vacuum state are
scaleless. The two contributing real emission diagrams at
$\mathcal{O}(\alpha_{s})$ in Feynman gauge are shown in figure -176b. The one-
gluon Feynman rule for the $Y_{v}$ Wilson line can e.g. be read off from ref.
[46], which uses a covariant derivative convention consistent with the one-
gluon Feynman rule we use for $W_{\eta}$. Including the mirror diagram for
diagram (a), the two diagrams are both given by simple integrals over the
rapidity $y$ of the emitted gluon in the heavy-quark rest frame,

$\displaystyle\chi_{1}^{(a)}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{e^{\epsilon\gamma_{E}}}{\Gamma(1-\epsilon)}\frac{\mu^{2\epsilon}}{\pi
k_{T}^{2+2\epsilon+\eta}}\Bigl{(}\frac{\bar{n}\cdot
v}{\nu}\Bigr{)}^{-\eta}\int\\!\mathrm{d}y\,\mathcal{A}^{(a)}(y,\eta)\,,$
$\displaystyle\chi_{1}^{(b)}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{e^{\epsilon\gamma_{E}}}{\Gamma(1-\epsilon)}\frac{\mu^{2\epsilon}}{\pi
k_{T}^{2+2\epsilon}}\int\\!\mathrm{d}y\,\mathcal{A}^{(b)}(y)\,,$ (4.26)

where the relevant integrands are given by

$\displaystyle\mathcal{A}^{(a)}(y,\eta)=\frac{2e^{-\eta y}}{e^{y}\,\cosh
y}\,,\qquad\mathcal{A}^{(b)}(y)=-\frac{1}{\raisebox{-1.00006pt}{$\cosh^{2}y$}}\,,$
(4.27)

and their integrals evaluate to

$\displaystyle\int\\!\mathrm{d}y\,\mathcal{A}^{(a)}(y,\eta)=-\frac{2\pi}{\sin(\eta\pi/2)}=-\frac{4}{\eta}+\mathcal{O}(\eta)\,,\qquad\int\\!\mathrm{d}y\,\mathcal{A}^{(b)}(y)=-2\,.$
(4.28)

The origin of the rapidity divergence, like the structure of the integrands as
a whole, has a rather intuitive interpretation, which we illustrate in the
right panel of figure -176b: Emissions at large positive gluon rapidity $y$ in
the rest frame of the heavy quark are strongly suppressed for either
contributing diagram, which is the well-known dead-cone effect [32]. For
diagram (b), where the static heavy quark $Y_{v}$ both emits and absorbs the
gluon, the contribution to the emission probability (dot-dashed blue) peaks at
$y=0$ and is suppressed also at large negative rapidities because the
contribution must be even in the rest frame. In contrast to this, the
contribution from diagram (a), where the emission couples to the $W_{\eta}$
that arises from integrating out the interactions with the anticollinear and
soft sectors and is required for gauge invariance, tends to a large constant
as $y\to-\infty$ (dashed red), leading to a rapidity divergence. Of course,
this precisely corresponds to the central soft momentum region, and thus after
introducing the regulator $e^{-\eta y}$ with $\eta<0$ to dampen the behavior
as $y\to\infty$ (solid red), the associated pole can cancel against the soft
function as usual.

Adding the one-loop soft function in eq. (3.15), which at this order is
independent of the matter content of the theory, we find that the $\eta$ poles
cancel for arbitrary $\epsilon$ as they must. Taking $\eta\to 0$, we find,
using the tree-level result
$\chi_{1}^{(0)}(k_{T})=\tfrac{1}{\pi}\delta(k_{T}^{2})$,

$\displaystyle\lim_{\eta\to
0}\Bigl{[}\chi_{1}^{\mathrm{bare}\,(1)}(k_{T},\epsilon,\eta,\frac{\bar{n}\cdot
v}{\nu}\Bigr{)}+\tfrac{1}{2}S^{(n_{\ell})\,(1)}(k_{T},\epsilon,\eta,\nu)\Bigr{]}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{e^{\epsilon\gamma_{E}}}{\Gamma(1-\epsilon)}\frac{\mu^{2\epsilon}}{\pi
k_{T}^{2+2\epsilon}}\Bigl{(}4\ln\bar{n}\cdot v-2\Bigr{)}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\,\frac{1}{\pi}\Bigl{[}-\frac{\delta(k_{T}^{2})}{\epsilon}+\mathcal{L}_{0}(k_{T}^{2},\mu)+\mathcal{O}(\epsilon)\Bigr{]}\Bigl{(}4\ln\bar{n}\cdot
v-2\Bigr{)}\,.$ (4.29)

The remaining pole in $\epsilon$ is precisely canceled by the predicted
counter term in eq. (4.1), confirming the renormalization structure we
discussed in section 4.1 at one loop. Taking a Fourier transform of the finite
remainder, our final result for the renormalized total TMD fragmentation
factor in the perturbative domain, and thus for the matching coefficient
$C_{1}$, reads

$\displaystyle\chi_{1}(b_{T},\mu,\rho)$
$\displaystyle=C_{1}(b_{T},\mu,\rho)+\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})$
$\displaystyle=1+\frac{\alpha_{s}C_{F}}{4\pi}(-L_{b})\bigl{(}4\ln\rho-2\bigr{)}+\mathcal{O}(\alpha_{s}^{2})+\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})\,,$
(4.30)

where $L_{b}$ was defined in eq. (3.44).

##  5 Consistency checks in the large and small-mass limits

###  5.1 Large-mass limit

To check the behavior of our final position-space results for $d_{1\,Q/Q}$ in
eq. (3.4), we first note that the result is exponentially suppressed at $z<1$
for $y\equiv mb_{T}\gg 1$ as a consequence of eq. (2.17) and the large-
argument behavior of the Bessel $K$ functions, i.e., it is given by
$\delta(1-z)$ up to a proportionality factor. To verify the prediction from
the heavy-quark limit in eq. (2.20), it is then sufficient to compare the
total $z$ integral (or any other $z$ moment) of both sides of eq. (2.20). The
total $z$ integral yields

$\displaystyle\int\\!\mathrm{d}z\,d_{1\,Q/Q}^{(1)}(z,b_{T},\mu,\zeta)$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\biggl{\\{}-2L_{b}\ln\frac{\zeta}{m^{2}}+4\ln^{2}\frac{\mu}{m}+6\ln\frac{\mu}{m}+4-L_{y}^{2}-\frac{\pi^{2}}{6}$
$\displaystyle\qquad\qquad-\frac{4\pi}{y}+4\mathcal{L}_{0}(y^{2},b_{0}^{2})+12\delta(y^{2})+\mathcal{O}\Bigl{(}\frac{1}{y^{3}}\Bigr{)}\biggr{\\}}\,.$
(5.1)

Here we used eq. (3.3) for the total $z$ integral of the cumulant-space
remainder term $\mathcal{R}(z,{k_{T}^{\mathrm{cut}}}/m)$, which we then power
expanded for small ${k_{T}^{\mathrm{cut}}}/m$, letting us evaluate the total
$z$ integral of eq. (3.4) power by power in $1/m$. We find that the leading
$\mathcal{O}(1/m^{0})$ term on the first line is in full agreement with the
bHQET prediction from eq. (2.20), using our one-loop result for $C_{1}$ in eq.
(4.3) and the known result for $C_{m}$ in eq. (B.2). We point out the
interesting observation that the first subleading $\mathcal{O}(1/m)$
correction in eq. (5.1) is free of large logarithms $L_{y}\equiv
2\ln\tfrac{y}{b_{0}}$, $b_{0}\equiv 2e^{-\gamma_{E}}$, suggesting a simple
structure if factorized in terms of subleading $\mathcal{O}(1/m)$ TMD matrix
elements in bHQET. On the other hand, the terms at the next order require plus
regularization in $b_{T}$ space due to their $1/y^{2}$ scaling and do feature
a single-logarithmic term $\mathcal{L}_{0}(y^{2},b_{0}^{2})$.

For completeness, we may also expand the nonvalence results in eq. (3.5) for
$y\gg 1$ at finite values of $0<z<1$ in a similar way, which simply yields

$\displaystyle d_{1\,Q/g}^{(1)}(z,b_{T})$
$\displaystyle=\frac{\alpha_{s}T_{F}}{4\pi}\frac{1}{z^{2}}\biggl{\\{}8z^{2}\,\delta(y^{2})+\mathcal{O}\Bigl{(}\frac{1}{y^{3}}\Bigr{)}\biggr{\\}}\,,$
$\displaystyle\mathcal{J}_{g/Q}^{(1)}(z,b_{T},m)$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\biggl{\\{}8z\,\delta(y^{2})+\mathcal{O}\Bigl{(}\frac{1}{y^{3}}\Bigr{)}\biggr{\\}}\,,$
(5.2)

i.e., both channels are in fact suppressed by $1/m^{2}$ at finite $z$.

###  5.2 Small-mass limit

Expanding eqs. (3.4) and (3.5) for small values of the mass, $y\equiv
mb_{T}\ll 1$, we obtain:

$\displaystyle d^{(1)}_{1\,Q/Q}(z_{H},b_{T},\mu,\zeta)$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\frac{1}{z^{2}}\bigg{\\{}\delta(1-z)\bigg{[}-2L_{b}\ln\frac{\zeta}{m^{2}}+4\ln^{2}\frac{\mu}{m}+6\ln\frac{\mu}{m}-L_{y}^{2}+4-\frac{\pi^{2}}{6}\bigg{]}$
$\displaystyle\qquad\qquad\qquad-4(1+L_{y})\,\mathcal{L}_{0}(1-z)-8\mathcal{L}_{1}(1-z)$
$\displaystyle\qquad\qquad\qquad+4+2L_{y}(1+z)+\frac{8\ln\bar{z}-4(1+z^{2})\ln\tfrac{\bar{z}}{z}}{\bar{z}}$
$\displaystyle\qquad\qquad\qquad+y^{2}\biggl{(}\frac{1}{z^{2}}-z-\frac{L_{y}(1+z)^{2}\bar{z}}{2z^{2}}-\frac{(1+z)^{2}\bar{z}\ln\tfrac{\bar{z}}{z}}{z^{2}}\biggr{)}+\mathcal{O}(y^{4})\bigg{\\}}\,,$
$\displaystyle d^{(1)}_{1\,Q/g}(z_{H},b_{T},\mu,\zeta)$
$\displaystyle=\frac{\alpha_{s}T_{F}}{4\pi}\biggl{\\{}-2L_{y}[2(z-1)z+1]-4(z-1)z+4[2(z-1)z+1]\ln
z$
$\displaystyle\qquad\qquad\qquad+y^{2}\biggl{(}-\frac{L_{y}(1-2z)^{2}}{2z^{2}}+\frac{3(z-1)z+1}{z^{2}}+\frac{(1-2z)^{2}\ln
z}{z^{2}}\biggr{)}+\mathcal{O}(y^{4})\biggr{\\}}\,,$
$\displaystyle\mathcal{J}^{(1)}_{g/Q}(z_{H},b_{T},\mu,\zeta)$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\biggl{\\{}-\frac{2L_{y}[(z-2)z+2]}{z}+\frac{4(z-1)}{z}$
$\displaystyle\qquad\qquad\qquad+y^{2}\left(-\frac{L_{y}(z-2)^{2}}{2z}+z+\frac{3}{z}-3\right)+\mathcal{O}(y^{4})\biggr{\\}}\,.$
(5.3)

It is straightforward to verify that the leading $\mathcal{O}(y^{0})$ terms
are in full agreement with the prediction in eq. (2.4) using the one-loop
ingredients collected in our notation in Appendix B. At leading power, the
logarithms $L_{y}\equiv 2\ln\tfrac{y}{2e^{-\gamma_{E}}}$ in all channels are
predicted by the timelike DGLAP evolution between the TMD matching
coefficients at $\mu\sim 1/b_{T}$ and the matching kernels and partonic
collinear FFs at $\mu\sim m$. The $\mathcal{O}(y^{2})$ terms, which we here
give explicitly, likewise feature single logarithms $L_{y}$ in all partonic
channels, and can be used to investigate the $\mathcal{O}(b_{T}^{2})$ twist-4
matching of TMD FFs in a physics setup where the low scale matrix elements at
$\mu\sim m$ are perturbatively calculable.

##  6 Joint threshold and $q_{T}$ resummation for massive quark fragmentation

(a)  (b)  (c)

Figure -175:  Relevant modes contributing to the heavy-quark TMD FF as the
collinear momentum fraction $z_{H}\to 1$ while counting $m\sim k_{T}$. See the
text for a detailed discussion.

The remarkable simplicity of our final $b_{T}$-space result in eq. (3.4) as
$z\to 1$ warrants further investigation. Specifically, the leading terms at
$\mathcal{O}[(1-z)^{-1}]$ feature only a logarithmic dependence on the mass
and $b_{T}$, while in general, a power-like dependence on $mb_{T}\sim 1$ would
be allowed, suggesting that the heavy-quark TMD FF can be factorized further
in this regime such that the logarithms become predicted by RGEs. To
investigate this question, we modify our previous assumption on the power
counting of $\bar{z}_{H}\equiv 1-z_{H}$, eq. (2.17), and instead count

$\displaystyle\text{\hyperref@@ii[sec:joint_res]{Section~{}\ref*{sec:joint_res}}:}\quad\frac{k_{T}}{Q}\sim\frac{m}{Q}\ll\bar{z}_{H}\ll
1\,.$ (6.1)

We stress that we _continue_ to count $m\sim k_{T}\sim 1/b_{T}$. In figure
-175 we illustrate the contributing modes in the EFT for this scenario. Figure
-175 (a) shows the original, unmodified mode setup that led to eq. (2.1). Here
the partonic heavy-quark TMD FF arises from matching collinear modes (green),
accounting for the usual contribution of the central soft modes (orange), onto
the nonperturbative modes (gray) describing the formation of the brown muck,
which is isotropic in the rest frame of the heavy quark indicated by the
dashed black line. The collinear modes are kinematically unconstrained in the
longitudinal direction because $\bar{z}_{H}\sim 1$.

Once the modified assumption eq. (6.1) is made, energetic collinear modes at
the scale $\mu\sim k_{T}\sim m$ can no longer be radiated into the final
state, see figure -175 (b) and (c). Their effect becomes purely virtual and is
encoded in a matching coefficient $H_{m,n}^{Q\to h_{v}}$ (green cross) from
matching onto bHQET modes at lower virtuality. Real emissions into the
$n$-collinear direction at the scale $\mu\sim k_{T}\sim m$ can instead only
have collinear-soft (pink) scaling

$\displaystyle
p_{\mathrm{cs}}\sim\Bigl{(}\bar{z}_{H}Q,\tfrac{k_{T}^{2}}{\bar{z}_{H}Q},k_{T}\Bigr{)}\,,$
(6.2)

which is the characteristic mode of double-differential $z\to 1$ threshold and
$q_{T}$ resummation, often simply dubbed joint resummation [64, 65, 66, 67,
68, 69, 70, 71, 72].222 We note that while ref. [68] also dealt with final-
state heavy quarks, that work considered the joint resummation for the total
transverse momentum of a heavy quark _pair_ with respect to the beam
direction, which is a distinct problem from the one we consider here.

The physical picture at scales $\mu\ll k_{T}\sim m$ in this scenario, which is
governed by one or several bHQET modes at their characteristic lab-frame
rapidity $p^{-}/p^{+}\sim Q^{2}/m^{2}$, further depends on the relative size
of $\Lambda_{\mathrm{QCD}}/m$ and $\bar{z}_{H}$, specifically whether
$\Lambda_{\mathrm{QCD}}/m\sim\bar{z}_{H}$ or
$\Lambda_{\mathrm{QCD}}/m\ll\bar{z}_{H}$. (Since $\Lambda_{\mathrm{QCD}}$
introduces a fourth scale into the problem, i.e., a third independent
dimensionless power counting parameter, either case may hold in general and
independent of $k_{T}/Q$ and $m/Q$.) The relevant modes for the case
$\Lambda_{\mathrm{QCD}}/m\sim\bar{z}_{H}$ are illustrated in figure -175 (c),
and are in fact analogous to the well-understood case of single-inclusive
heavy-hadron production near the endpoint in $e^{+}e^{-}$ collisions [5, 6].
Specifically, the nonperturbative fragmentation process at the scale
$\mu\sim\Lambda_{\mathrm{QCD}}\sim\bar{z}_{H}m$ provides a leading
contribution to the measured longitudinal hadron momentum in this case (but
not to its transverse momentum, where the contribution is power suppressed!).
This effect is encoded in the well-known collinear bHQET fragmentation shape
function $F_{H}(k^{-}/v^{-},\mu)$ [5, 6].333We use the symbol $F_{H}$ here to
distinguish the shape function from the TMD soft function; in ref. [6], the
symbol $S_{H/Q}$ is used instead, and our partonic $\hat{F}$ corresponds to
their $S_{Q/Q}$. In total, the factorization of the hadron-level heavy-quark
TMD FF in this regime reads, with $\zeta=\omega^{2}$ the Collins-Soper scale,

$\displaystyle D_{1\,H/Q}(z_{H},b_{T},\mu,\omega^{2})$
$\displaystyle=H_{m,n}^{Q\to
h_{v}}\Bigl{(}m,\mu,\frac{\nu}{\omega}\Bigr{)}\,\int\\!\mathrm{d}\ell^{-}\,\mathcal{S}(\ell^{-},b_{T},m,\mu,\nu)\,\sqrt{S}(b_{T},m,\mu,\nu)$
$\displaystyle\quad\times\int\\!\frac{\mathrm{d}k^{-}}{v^{-}}\,F_{H}\Bigl{(}\frac{k^{-}}{v^{-}},\mu\Bigr{)}\,\delta\Bigl{(}\bar{z}_{H}-\frac{\ell^{-}+k^{-}}{\omega}\Bigr{)}$
$\displaystyle\quad\times\biggl{[}1+\mathcal{O}(\bar{z}_{H})+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}}\Bigr{)}+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{m}\Bigr{)}\biggr{]}\,,$
(6.3)

where $\mathcal{S}$ and $\sqrt{S}$ are the separately rapidity-renormalized
collinear-soft function from joint resummation [72] and the square root of the
standard TMD soft function, respectively. Similarly, $H_{m,n}^{Q\to h_{v}}$
(defined in ref. [47] without the superscript, which we here include for
clarity) is the rapidity-renormalized matching coefficient encoding the
virtual contributions of massive collinear quark modes $Q$ with label momentum
$\omega$ matched onto bHQET fields $h_{v}$. As in ref. [47], the effect of
these modes is purely virtual due to the kinematic constraint $z_{H}\to 1$. We
note that $H_{m,n}^{Q\to h_{v}}$ also includes the effect of secondary heavy
quark loops not directly connected to the operator starting at
$\mathcal{O}(\alpha_{s}^{2})$, which are responsible for its rapidity
renormalization. By evolving all terms in eq. (6) from their canonical scales
(e.g. in Mellin space $N\leftrightarrow z_{H}$, where
$\nu_{\mathcal{S}}\sim\mu_{F_{H}}v^{-}\sim\omega/N$ and the evolution is
multiplicative) to the common scale $\mu\sim\nu\sim Q$ with the hard function,
all large threshold and transverse momentum logarithms present in the heavy-
hadron TMD cross section in this regime can be resummed simultaneously.

It is instructive to compare the above to the refactorization of the light-
parton TMD FF in the joint resummation regime $k_{T}\ll\bar{z}_{h}Q\ll Q$,
which follows from crossing the light-parton TMD PDF refactorization relations
[72], accounting for the additional effect of secondary heavy quarks starting
at two loops,

$\displaystyle D_{1\,h/q}(z_{h},b_{T},\mu,\omega^{2})$
$\displaystyle=H_{m,n}^{q\to
q}\Bigl{(}m,\mu,\frac{\nu}{\omega}\Bigr{)}\int\\!\mathrm{d}\ell^{-}\,\mathcal{S}(\ell^{-},b_{T},m,\mu,\nu)\,\sqrt{S}(b_{T},m,\mu,\nu)$
$\displaystyle\quad\times\int\\!\mathrm{d}k^{-}\,D^{\mathrm{thr}}_{h/q}\Bigl{(}\frac{k^{-}}{P_{h}^{-}},\mu\Bigr{)}\,\delta\Bigl{(}\bar{z}_{h}-\frac{\ell^{-}+k^{-}}{\omega}\Bigr{)}$
$\displaystyle\quad\times\biggl{[}1+\mathcal{O}(\bar{z}_{H})+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{k_{T}}\Bigr{)}+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}}{m}\Bigr{)}\biggr{]}\,,$
(6.4)

where $D^{\mathrm{thr}}_{h/q}$ is the collinear (threshold) fragmentation
function, whose renormalization is given by the timelike DGLAP kernels
expanded to leading power in $1-z$. The matching coefficient $H_{m,n}^{q\to
q}$ (defined as $H_{c}$ in ref. [48]) encodes the virtual contributions of
secondary heavy-quark loops to the matching of light-quark collinear modes $q$
at the scale $\mu\sim m$ onto their respective counterparts in the theory at
lower scales. By comparing eqs. (6) and (6) and making use of known results
for the heavy-quark threshold FF [6],444 Specifically, ref. [6] features a
combined mass-mode matching coefficient
$C_{m}^{\text{ref.\,\cite[cite]{[\@@bibref{Number}{Fickinger:2016rfd}{}{}]}}}(m,\mu)$
appearing as
$D_{H/Q}(z_{H}=1-k^{-}/P_{H}^{-},\mu)=C_{m}^{\text{ref.\,\cite[cite]{[\@@bibref{Number}{Fickinger:2016rfd}{}{}]}}}(m,\mu)\,F_{H}(k^{-}/v^{-},\mu)\bigl{[}1+\mathcal{O}(\bar{z}_{H})\bigr{]}$,
not to be confused with our $C_{m}(m,\mu,\zeta)$ in eq. (2.18). This relation,
which defines
$C_{m}^{\text{ref.\,\cite[cite]{[\@@bibref{Number}{Fickinger:2016rfd}{}{}]}}}$
as the finite remainder of the collinear FF after dividing out the shape
function, was then used to extract it to two loops from the partonic heavy-
quark FF [73]. We have verified to two loops that $X\equiv
C_{m}^{\text{ref.\,\cite[cite]{[\@@bibref{Number}{Fickinger:2016rfd}{}{}]}}}(m,\mu)/[H_{m,n}^{Q\to
h_{v}}(m,\mu,\omega/\nu)/H_{m,n}^{q\to q}(m,\mu,\omega/\nu)]$ is RG invariant,
which together with eq. (6) guarantees RG invariance of eq. (6). In
particular, the $\nu$ dependence cancels between $H_{m,n}^{Q\to h_{v}}$ and
$H_{m,n}^{q\to q}$, as expected because it arises from the common, secondary
mass effects. We stress, however, that the defining relation for
$C_{m}^{\text{ref.\,\cite[cite]{[\@@bibref{Number}{Fickinger:2016rfd}{}{}]}}}$
cannot simply be inserted back into eq. (6) because it only holds when all
other scales in the problem are large compared to the mass. Accordingly, $X$
does not need to be equal to one, and we indeed find a nontrivial constant
term at $\mathcal{O}(\alpha_{s}^{2})$.  we can conclude that eq. (6) as the
main new result of this section is indeed RGE consistent. We note that both
eqs. (6) and (6) also generically feature a dependence on the heavy quark mass
in the soft and collinear-soft function since $mb_{T}\sim 1$, to which we
return towards the end of this section.

To make contact with the _partonic_ heavy-quark TMD FF, we next take
$\Lambda_{\mathrm{QCD}}\to 0$, which in particular means that now
$\Lambda_{\mathrm{QCD}}/m\ll\bar{z}_{H}$. The mode setup for this regime is
illustrated in figure -175 (b). The only change compared to figure -175 (c) is
that the fragmentation shape function now factorizes into the perturbative
bHQET shape function $\hat{F}$ and the total fragmentation probability
$\chi_{H}$ [6],

$\displaystyle
F_{H}\Bigl{(}\frac{k^{-}}{v^{-}},\mu\Bigr{)}=\hat{F}\Bigl{(}\frac{k^{-}}{v^{-}},\mu\Bigr{)}\,\chi_{H}+\mathcal{O}(\Lambda_{\mathrm{QCD}})\,.$
(6.5)

This implies that the joint resummation of the partonic heavy-quark TMD FF is
governed by the following refactorization relation,

$\displaystyle d_{1\,Q/Q}(z,b_{T},\mu,\omega^{2})$
$\displaystyle=H_{m,n}^{Q\to
h_{v}}\Bigl{(}m,\mu,\frac{\nu}{\omega}\Bigr{)}\,\int\\!\mathrm{d}\ell^{-}\,\mathcal{S}(\ell^{-},b_{T},m,\mu,\nu)\,\sqrt{S}(b_{T},m,\mu,\nu)$
$\displaystyle\quad\times\int\\!\frac{\mathrm{d}k^{-}}{v^{-}}\,\hat{F}\Bigl{(}\frac{k^{-}}{v^{-}},\mu\Bigr{)}\,\delta\Bigl{(}\bar{z}-\frac{\ell^{-}+k^{-}}{\omega}\Bigr{)}+\mathcal{O}\bigl{[}(1-z)^{0}\bigr{]}\,.$
(6.6)

This is precisely the relation governing the leading $1/\bar{z}$ terms in the
partonic heavy-quark TMD FF that we were seeking out, and it is
straightforward to verify eq. (6) at one-loop order against eq. (3.4) using
the fact that $\tilde{\mathcal{R}}(z,y)=\mathcal{O}[(1-z)^{0}]$, the fact that
$v^{-}=\omega/m$, and the ingredients given in appendix B. In particular,
since the soft and soft-collinear function are defined as Wilson line
correlators, they can only acquire mass dependence through secondary mass
effects starting at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$. In turn, this
implies that the mass dependence at $\mathcal{O}(\alpha_{s})$ arises only from
that of $C_{m}$, which is logarithmic to all orders, explaining why eq. (3.4)
has to have the simple logarithmic mass dependence that we found for
$\bar{z}\to 1$. We stress that starting at $\mathcal{O}(\alpha_{s}^{2})$, the
soft and soft-collinear functions _do_ develop a highly nontrivial dependence
on $mb_{T}$, which is thus inherited by the coefficients of
$\mathcal{L}_{0}(\bar{z})$ and $\delta(\bar{z})$ in $d_{1\,Q/Q}^{(2)}$ at two
loops. The required ingredients are readily available in the literature [74,
48],555The secondary mass effects in $\sqrt{S}$ have been calculated in ref.
[48]. For the collinear-soft function one may use that for the exponential
regulator of ref. [71], the collinear-soft function is simply given by an
exponential of the Collins-Soper kernel to all orders [75]. This remains true
in the presence of mass effects in the CS kernel, which have likewise been
calculated in ref. [48]. One potential subtlety is that the results of ref.
[48] for $\sqrt{S}$ were obtained using the $\eta$ regulator of refs. [43, 44]
instead. However, it is easy to see that the one-loop soft function with a
massive gluon (which underlies the dispersive calculation of the secondary
quark mass effects in ref. [48]) leads to identical results for both
regulators, meaning the final renormalized result for the mass-dependent soft
function in ref. [48] is also valid for either regulator.  which can provide a
strong cross check on a future two-loop calculation of the partonic heavy-
quark TMD FF, or simplify it by allowing to restrict to finite $z<1$ from the
start.

##  7 Application: EEC with heavy quarks in the back-to-back limit

###  7.1 Mass-dependent EEC jet functions

In this section we use our results to produce explicit results for the effect
of heavy-quarks on the energy-energy correlator (EEC) observable [76] as
measured in $e^{+}e^{-}$ collisions,

$\displaystyle\mathrm{EEC}(\chi)\equiv\frac{1}{\sigma_{\mathrm{total}}}\,\frac{\mathrm{d}\sigma}{d\cos\chi}=\frac{1}{\sigma_{\mathrm{total}}}\sum_{a,b\,\in
X}\int\\!\mathrm{d}\sigma_{ee\to
X}\frac{E_{a}E_{b}}{Q^{2}}\delta\bigl{[}\cos(\theta_{ab})-\cos(\chi)\bigr{]}\,,$
(7.1)

where the sum runs over all pairs of hadrons $a$ and $b$ in the event with
energies $E_{a},E_{b}$ and $Q$ is the center-of-mass energy. We are
specifically interested in the so-called back-to-back limit $\chi\to\pi$,
where a parton-level factorization theorem has been proven [25] and can be
directly related to TMD factorization. This precisely complements the study of
ref. [77], where the NLO mass effects in the opposite (collinear) limit
$\chi\to 0$ of the EEC were calculated, and can serve as data for a future
calculation of the heavy-quark EEC for general angles.

To do so, it is convenient to rewrite the EEC as a differential distribution
in $z_{\chi}=\bigl{[}1-\cos(\theta_{ab})\bigr{]}/2\to 1$, i.e., the average
pairwise angle of energy-weighted hadron pairs,
$\mathrm{EEC}(\chi)=(1/\sigma_{\mathrm{total}})\,\mathrm{d}\sigma/\mathrm{d}z_{\chi}$.
We are interested in the modifications of the fixed-order EEC from the
presence of heavy quarks in this limit, counting $(1-z_{\chi})Q^{2}\sim
1/b_{T}^{2}\sim m^{2}$, where $b_{T}$ – precisely as in TMD factorization – is
the typical transverse spacetime separation of fields in the relevant
correlators. The $z_{\chi}$-differential distribution is given by
$z_{\mathfrak{h}}$-moments of the TMD FFs for identified hadrons [25],

$\displaystyle\frac{1}{\sigma_{0}}\frac{\mathrm{d}\sigma}{\mathrm{d}z_{\chi}}$
$\displaystyle=\sum_{i,j}\mathcal{H}_{ee\to
ij}(Q^{2},\mu)~{}Q^{2}\\!\\!\int_{0}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}\,J_{0}\bigl{(}\sqrt{1-z_{\chi}}\,b_{T}Q\bigr{)}$
(7.2)
$\displaystyle\quad\times\sum_{\mathfrak{h},\mathfrak{h}^{\prime}}\int\\!\mathrm{d}z_{\mathfrak{h}}\,z_{\mathfrak{h}}^{3}D_{1\,\mathfrak{h}/i}(z_{\mathfrak{h}},b_{T},\mu,Q^{2})\int\\!\mathrm{d}z_{\mathfrak{h}^{\prime}}\,z_{\mathfrak{h}^{\prime}}^{3}D_{1\,\mathfrak{h}^{\prime}/j}(z_{\mathfrak{h}^{\prime}},b_{T},\mu,Q^{2})$
$\displaystyle\equiv\sum_{i,j}\mathcal{H}_{ee\to
ij}(Q^{2},\mu)~{}Q^{2}\\!\\!\int_{0}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}\,J_{0}\bigl{(}\sqrt{1-z_{\chi}}\,b_{T}Q\bigr{)}$
$\displaystyle\quad\times
J_{i}(b_{T},m,\mu,Q/\nu)\,J_{j}(b_{T},m,\mu,Q/\nu)\,S(b_{T},m,\mu,\nu)\,,$

where the sum runs over all possible pairs of partons $i,j$, with nonzero
contributions for $i,j=q,\bar{q}$ a light quark pair or $i,j=Q,\bar{Q}$ a
heavy quark pair, and it is understood that $\mathcal{H}_{ee\to gg}=0$. Here
our normalization is such that
$\mathrm{d}\sigma/\mathrm{d}z_{\chi}=\sigma_{0}\,\delta(1-z_{\chi})+\mathcal{O}(\alpha_{s})$
in the limit where the process is mediated only by virtual photons, and the
hard function $\mathcal{H}_{ee\to ij}$ is given in appendix B. On the second
equality we defined the EEC back-to-back jet function $J_{i}(b_{T},m,\mu,\nu)$
for a – possibly massive – parent parton $i$, and we followed the usual
convention in the EEC literature to factor off a rapidity and UV-renormalized
TMD soft function $S(b_{T},m,\mu,\nu)$. Both $J_{i}$ and $S$ acquire mass
dependence for $1/b_{T}^{2}\sim m^{2}$ from primary and secondary heavy-quark
effects, which we now classify and evaluate further. Note that while the
energy fractions entering eq. (7.1) in general differ from the lightcone
momentum fractions $z_{\mathfrak{h}}$ that the TMD FFs depend on, in
particular for massive particles, the difference is power-suppressed in the
back-to-back limit where $k_{T},m,m_{\mathfrak{h}}\ll z_{\mathfrak{h}}Q$.

To evaluate the sums over hadrons, we break each of them down into two terms,
a sum over light hadrons $\mathfrak{h}=h$ and one over heavy hadrons
$\mathfrak{h}=H$. For light-hadron final states, the standard manipulations of
the collinear FF using its momentum sum rule apply [25], such that also in the
presence of intermediate heavy quarks and for all of $i=q,Q,g$,

$\displaystyle\sum_{h}\int\\!\mathrm{d}z_{h}\,z_{h}^{3}\,D_{1\,h/i}(z_{h},b_{T},\mu,Q^{2})$
$\displaystyle=\sum_{k=q,\bar{q},g}\int\\!\mathrm{d}z\,z\,\mathcal{J}_{k/i}(z,b_{T},m,\mu,Q^{2})+\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})\,,$
(7.3)

where $\mathcal{J}_{1\,k/i}$ is the mass-dependent matching coefficient
introduced in eq. (2.16). Importantly, at the scale
$\mu\sim\Lambda_{\mathrm{QCD}}$ light partons can only fragment to light
hadrons in the theory with $n_{\ell}$ light quark flavors, which means that
the collinear FF momentum sum rule holds for the set of light hadrons
$\\{h\\}$ on their own within that theory.

On the other hand, for heavy-hadron final states $H$ containing either $Q$ or
$\bar{Q}$, we have

$\displaystyle\sum_{H}\int\\!\mathrm{d}z_{H}\,z_{H}^{3}\,D_{1\,H/i}(z_{H},b_{T},\mu,Q^{2})$
$\displaystyle=\sum_{K=Q,\bar{Q}}\int\\!\mathrm{d}z\,z^{3}\,d_{1\,K/i}(z,b_{T},\mu,Q^{2})+\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})\,,$
(7.4)

which again holds for all of $i=q,Q,g$. While the end result looks similar to
eq. (7.1), the moment of the perturbative coefficient arises in a different
(and much simpler) way, namely because the fragmentation probabilities
$\chi_{H}$ from eq. (2.1) sum to unity, while $d_{1\,Q/i}$ and
$d_{1\,\bar{Q}/i}$ are independent of the state $H$. We note in passing that
bound states containing multiple heavy quarks also in general contribute to
the sum over hadrons in eq. (7.2) and e.g. for a $Q\bar{Q}$ bound state can be
evaluated in terms of the TMD FF matching coefficients onto NRQCD long-
distance matrix elements calculated in refs. [9, 10, 11], but this
contribution to the jet function $J_{i}$ only starts at
$\mathcal{O}(\alpha_{s}^{2})$.

The expressions above are particular ($N=1$) Mellin moments of the partonic
heavy-quark TMD FF and the TMD FF matching coefficients, which we in general
define as

$\displaystyle d_{1\,Q/i}^{[N]}(b_{T},\mu,\zeta)$
$\displaystyle\equiv\int\\!\mathrm{d}z\,z^{N+2}d_{1\,Q/i}(z,b_{T},\mu,\zeta)\,,$
$\displaystyle\mathcal{J}_{k/i}^{[N]}(b_{T},m,\mu,\zeta)$
$\displaystyle\equiv\int\\!\mathrm{d}z\,z^{N}\mathcal{J}_{k/i}(z,b_{T},m,\mu,\zeta)\,,$
(7.5)

where for $d_{1\,Q/i}$ we explicitly compensate for the usual measure factor
of $z^{2}$ when taking moments at fixed $b_{T}$ (or $k_{T}$). In summary, we
have

$\displaystyle
J_{i}(b_{T},m,\mu,Q/\nu)\bigl{[}S(b_{T},m,\mu,\nu)\bigr{]}^{1/2}$
$\displaystyle=\sum_{K=Q,\bar{Q}}d_{1\,K/i}^{[1]}(b_{T},\mu,Q^{2})+\sum_{k=q,\bar{q},g}\mathcal{J}_{k/i}^{[1]}(b_{T},m,\mu,Q^{2})\,,$
(7.6)

for the quark, massive quark, or gluon EEC back-to-back jet functions, which
generalizes the defining relation for the massless jet function in ref. [25].
We collect the Mellin moments of our new results that enter the right-hand
side in appendix C. We point out that the organization of our final one-loop
result for the heavy-quark TMD FF in momentum space in terms of plus
distributions on rectangular domains in $(z,k_{T})$, see eq. (3.2), makes it
particularly straightforward to take moments with respect to $z$.

###  7.2 Jet function consistency in the small and large-mass limits

The all-order consistency relation for the mass-dependent EEC jet functions in
the small-mass limit follows from eq. (2.4) and reads

$\displaystyle
J_{i}(b_{T},m,\mu,Q/\nu)=J_{i}(b_{T},\mu,Q/\nu)+\mathcal{O}(mb_{T})\,,$ (7.7)

i.e., they simply reduce to their massless counterparts because the partonic
heavy-quark collinear FFs and the FF flavor matching kernels preserve the
momentum sum rule. (In other words, the measurement becomes fully inclusive
over the details of the IR also in this case, which now includes the mass.)

On the other hand, in the large-mass limit we can infer from the discussion in
section 2.3 that the heavy-quark jet function to all orders satisfies

$\displaystyle
J_{Q}(b_{T},m,\mu,Q/\nu)\bigl{[}S(b_{T},m,\mu,\nu)\bigr{]}^{1/2}=C_{m}(m,\mu,Q^{2})\,C_{1}\Bigl{(}b_{T},\mu,\frac{Q}{m}\Bigr{)}\,.$
(7.8)

Note that it is not power suppressed in this regime because even at scales
below the quark mass, the heavy quark continues to radiate light modes and
pick up recoil, leading to a nontrivial angular dependence, and thus a nonzero
EEC. As for the TMD FFs, the mass effects in the light-parton jet functions
$J_{q}$ and $J_{g}$ for $(1-z_{\chi})Q^{2}\ll m^{2}$ are of purely virtual
origin and encoded in collinear and soft mass matching coefficients onto their
counterparts in the $n_{\ell}$ (purely light) theory, see refs. [48, 60].

###  7.3 Complete mass dependence of the EEC at $\mathcal{O}(\alpha_{s})$

The mass effects in the light-quark mediated cross section $i=q,\bar{q}$ only
enter through closed heavy-quark loops or pair production starting at
$\mathcal{O}(\alpha_{s}^{2})$, such that at $\mathcal{O}(\alpha_{s})$ the
total mass effect simply reads

$\displaystyle\frac{1}{\sigma_{0}}\frac{\mathrm{d}\sigma^{(1)}}{\mathrm{d}z_{\chi}}$
$\displaystyle=\frac{1}{\sigma_{0}}\frac{\mathrm{d}\sigma^{(1)}_{ee\to
q\bar{q}}}{\mathrm{d}z_{\chi}}+2\mathcal{H}_{ee\to
Q\bar{Q}}^{(1)}(Q^{2},\mu)\,\delta(1-z_{\chi})$ (7.9)
$\displaystyle\quad+4\mathcal{H}_{ee\to Q\bar{Q}}^{(0)}\,\pi
Q^{2}\Bigl{\\{}d_{1\,Q/Q}^{[1]\,(1)}\bigl{(}\sqrt{1-z_{\chi}}\,Q,\mu,Q^{2}\bigr{)}+\mathcal{J}_{g/Q}^{[1]\,(1)}\bigl{(}\sqrt{1-z_{\chi}}\,Q,\mu,Q^{2}\bigr{)}\Bigr{\\}}\,,$

where $\mathrm{d}\sigma_{ee\to q\bar{q}}$ refers to the light quark-mediated
cross section, which is known to $\mathcal{O}(\alpha_{s}^{3})$ [78], and the
expression for the hard function $\mathcal{H}_{ee\to Q\bar{Q}}$ is given in
appendix B. We also exploited that the convolution (or the inverse Fourier
transform) is trivial at fixed $\mathcal{O}(\alpha_{s})$, and that the
contributing partonic TMD FFs are all equal. Using the results of appendix C
for $N=1$ and adding the contribution from the one-loop hard-function, we find

$\displaystyle\frac{1}{\sigma_{0}}\frac{\mathrm{d}\sigma^{(1)}}{\mathrm{d}z_{\chi}}$
$\displaystyle=\frac{1}{\sigma_{0}}\frac{\mathrm{d}\sigma^{(1)}_{ee\to
q\bar{q}}}{\mathrm{d}z_{\chi}}+2\mathcal{H}_{ee\to
Q\bar{Q}}^{(0)}\,\frac{\alpha_{s}C_{F}}{4\pi}\biggl{\\{}8(2\ln\rho-1)\,\mathcal{L}_{0}(\bar{z}_{\chi})+8\Bigl{(}2\ln^{2}\rho+\ln\rho+\frac{2\pi^{2}}{3}-2\Bigr{)}\,\delta(\bar{z}_{\chi})$
$\displaystyle\qquad+4\rho^{2}\frac{\pi(-3-18x+x^{2})+\sqrt{x}(9+x-9x^{2}-x^{3})-\sqrt{x}(1+24x+9x^{2}+2x^{3})\ln
x}{\sqrt{x}(1+x)^{4}}\biggr{\\}}$
$\displaystyle\quad+\mathcal{O}(\alpha_{s}^{2})+\mathcal{O}\bigl{(}\bar{z}_{\chi}^{0}\bigr{)}+\mathcal{O}\Bigl{(}\frac{1}{\rho}\Bigr{)}\,,$
(7.10)

where we used the shorthands $\rho\equiv Q/m$ for the boost of the heavy
quarks, $\bar{z}_{\chi}\equiv 1-z_{\chi}$, and
$x\equiv\rho^{2}\bar{z}_{\chi}=k_{T}^{2}/m^{2}$. The result is valid up to
power corrections in $1/Q$, i.e., $\mathcal{O}(\bar{z}_{\chi}^{0})$ or
$\mathcal{O}(1/\rho)$, but captures the exact dependence on
$\rho^{2}\bar{z}_{\chi}\sim 1$.

###  7.4 Extension to nonperturbative
$1-z_{\chi,H}\sim\Lambda_{\mathrm{QCD}}^{2}/Q^{2}$

To close this section, we note that an intriguing modification to the
observable can be made where the sum in the _observable definition_ in eq.
(7.1) is restricted to heavy hadrons $a,b\in\\{H\\}$ containing valence heavy
quarks. For this modified observable, which we may call $\mathrm{EEC}_{H}$
with an associated differential variable $z_{\chi,H}$, it is in fact
straightforward to extend the factorization all the way down into the
nonperturbative region $1-z_{\chi,H}\sim\Lambda_{\mathrm{QCD}}^{2}/Q^{2}$:

$\displaystyle\frac{1}{\sigma_{0}}\frac{\mathrm{d}\sigma}{\mathrm{d}z_{\chi,H}}$
$\displaystyle=\mathcal{H}_{ee\to
Q\bar{Q}}(Q^{2},\mu)\,C_{m}^{2}(m,\mu,Q^{2})$ $\displaystyle\quad\times
Q^{2}\\!\\!\int_{0}^{\infty}\\!\mathrm{d}b_{T}\,b_{T}\,J_{0}\bigl{(}\sqrt{1-z_{\chi,H}}\,b_{T}Q\bigr{)}\Bigl{[}\chi_{1}\Bigl{(}b_{T},\mu,\frac{Q}{m}\Bigr{)}\Bigr{]}^{2}\,,$
(7.11)

where $\chi_{1}$ is the renormalized Wilson-line correlator defined in eq.
(4.23). The simplicity of eq. (7.4), which is valid for
$\Lambda_{\mathrm{QCD}}\lesssim\sqrt{1-z_{\chi,H}}\,Q\ll m\ll Q$ and coincides
with eq. (7.8) for $\Lambda_{\mathrm{QCD}}\ll\sqrt{1-z_{\chi,H}}\,Q$, comes
about exactly as for the individual $\chi_{1}$ in section 4: Summing over
heavy hadrons $H$ allows one to complete the sum over states in the EFT and
reduces the final result to (a square of) vacuum correlators of Wilson lines.
We note that in ref. [1], we already mentioned the factorization in eq. (7.4)
as a factorization theorem for the total heavy-hadron TMD cross section in
$e^{+}e^{-}$ collisions, but that – in hindsight – of course coincides with
the EEC in the back-to-back limit.

##  8 Comparison to 2310.19207 by Dai, Kim, and Leibovich

Recently, an independent NLO calculation of the heavy-quark TMD FFs
$d_{1\,Q/Q}$ and $d_{1\,Q/g}$ was performed in ref. [13]. Before commenting on
the differences between the two calculational approaches, we first want to
stress that after accounting for a measure factor of $1/z^{2}$ and the
renormalized TMD soft function in eq. (B.4), we find complete agreement of our
eqs. (3.4) and (3.5) with the final $b_{T}$-space results given in their eqs.
(2.36) and (4.26). In order to gauge the strength of this cross check, it is
relevant to note that while the calculation of the real-emission diagrams in
ref. [13], like ours in section 3.1, was initially performed in momentum space
and only Fourier-transformed after renormalization, ref. [13] used a
completely orthogonal method to isolate the singularities in
$d^{(1)}_{1\,Q/Q}$, i.e., to address the key source of subtleties.
Specifically, in ref. [13] the characteristic denominator $1/X\equiv
1/(k_{T}^{2}+m^{2}\bar{z}^{2}/z^{2})$ appearing in the real-emission diagrams
was directly expanded into one-dimensional plus distributions of $X$ and in
addition was multiplied by one-dimensional distributions of $1-z$, which
effectively amounts to performing a change of variables to $(z,X)$. This
procedure is completely orthogonal to our approach, namely, we expand in terms
of one-dimensional plus distributions along the boundaries $\bar{z}=0$ and
$k_{T}=0$ and a two-dimensional plus distribution whose defining domain
_aligns_ with the main variables $(z,k_{T})$. The two approaches in particular
differ (and have to differ) for the boundary contribution proportional to
$\delta(1-z)$. Indeed, the process in ref. [13] generates highly nontrivial
twofold integrals $F(m^{2}/\Lambda^{2})$ and $G(m^{2}/\Lambda^{2})$
multiplying $\delta(1-z)$, where $\Lambda$ is the boundary condition of the
$X$ plus distributions, with no closed form for $F$ and $G$ away from
asymptotic limits in $m$.

In view of these substantial differences, it is highly nontrivial that the
final results in $b_{T}$ space are indeed in full agreement, in particular
regarding the latter’s $\delta(1-z)$ boundary condition. In practical terms,
we believe that our final momentum-space result in eq. (3.2) compared to the
one in eq. (2.30) of ref. [13] nevertheless has significant advantages. A
first, practical advantage is that eq. (3.2) is expressed in closed form in
terms of elementary functions. A perhaps more important second advantage is
that aligning the domains of the plus distributions with the main variables of
interest makes it very straightforward to take integral transforms with
respect to either $k_{T}$ or $z$, see the discussion above eq. (3.3) and eq.
(C.2).

Ref. [13] also considered how the heavy-quark TMD FF factorizes in the bHQET
limit $k_{T}\ll m$. At variance with our assumption of $1-z_{H}\sim 1$ in ref.
[1] and in eq. (2.17), which amounts to considering wide bins or moments of
$z_{H}$, ref. [13] counts $1-z_{H}\sim k_{T}/m\ll 1$ in the bHQET regime,
i.e., staying differential in $z_{H}$ but restricting to the peak region,
which leads to a TMD bHQET shape function that depends on $z_{H}$, cf. the
comment below eq. (2.6) in our ref. [1]. Again, before commenting on the
differences between our approaches further, we want to stress that for the
power counting of $1-z_{H}\sim k_{T}/m$ assumed by the authors, their
factorization of the hadron-level TMD FF in terms of $C_{m}$ and their TMD
shape function (which may further be matched onto a collinear shape function,
depending on the hierarchy between $\Lambda_{\mathrm{QCD}}$ and $k_{T}$) is
entirely correct.

One difference between the two power-counting assumptions on $1-z_{H}$ regards
the kind of cross check the respective bHQET calculation provides on the
large-mass behavior of the TMD FF. The approach of ref. [13] checks the
differential $z$ dependence of $d_{1\,Q/Q}^{(1)}$, but in intermediate steps
an identical plus distribution expansion with respect to a very similar,
complicated denominator $1/X^{\prime}=1/(k_{T}^{2}+\bar{z}^{2}m^{2})$ is used,
which is subject to the same subtleties as the full calculation. In this
sense, our calculation provides a strong complementary check of the large-mass
behavior since it directly probes the total integral, i.e., a defining
property of the distributional result.666Ref. [13] also performed a small-mass
consistency check, using the one-loop TMD PDF matching coefficients and a
reciprocity relation. This provides a sufficient check at one loop, but
reciprocity relations between TMD PDF and FF matching coefficients are known
to be subtle at higher orders [51, 50]. Here we stress that the correct all-
order relation for the small-mass limit is the one we originally gave in ref.
[1], i.e., it has to involve the genuine TMD _FF_ matching coefficients onto
twist-2 collinear FFs because the UV matching coefficients are independent of
the details of the IR, that is, whether one matches onto a theory at
$\mu\sim\Lambda_{\mathrm{QCD}}$ or $\mu\sim m$.

A second and final difference between the power counting assumptions is
phenomenological and regards the kinds of differential distributions one may
describe in either approach. Our stance here is that while – again – the
assumption of ref. [13] is sufficient and entirely correct to describe
specific regions of the multi-differential $(z_{H},k_{T})$ distribution, it is
not the unique power-counting assumption on $1-z_{H}$, $k_{T}$, and $m$ that
would yield a leading-power factorization theorem involving bHQET modes, in
particular not when different possible scalings of $\Lambda_{\mathrm{QCD}}$
are also considered. As one example, we point to our analysis in section 6,
where we continued to count $m\sim k_{T}$, but still took $1-z_{H}$ small, and
likewise obtained a leading-power factorization theorem, in that case also
involving a purely collinear bHQET shape function, but nontrivial
$k_{T}$-dependent correlators of lightlike Wilson lines at the scale $\mu\sim
k_{T}$. In light of this, we believe that capturing the full $(z_{H},k_{T})$
dependence will require dedicated future studies combining all of these
results, while for the time being and within the scope of exploratory theory
studies (or, necessarily less differential, initial experimental measurements)
the restriction to wide bins or low-moment integrals over $z_{H}$ is
reasonable.

##  9 Conclusions

In this paper, we calculated all heavy-quark transverse-momentum dependent
fragmentation functions (TMD FFs) to next-to-leading order (NLO) in the strong
coupling. We provided explicit results in transverse momentum ($k_{T}$) space
– both for the differential and the cumulative distribution – and in the
conjugate position ($b_{T}$) space, retaining the exact dependence on
$k_{T}/m$ and $b_{T}m$, respectively, with $m$ the mass of the heavy quark.
Our calculations provide the last missing, key ingredient for a complete next-
to-next-to-leading logarithmic (NNLL) description of the transverse momentum
distributions of observed heavy hadrons, capturing all quark mass effects. Our
results also enable the extension of fixed-order subtraction methods to quasi-
collinear limits involving final-state heavy quarks.

We further studied the corresponding TMD fragmentation factors and the TMD
light-spin density matrix that arise in boosted Heavy-Quark Effective Theory
(bHQET) when the transverse momentum is small compared to the mass. We
demonstrated explicitly that they satisfy a two-dimensional set of
renormalization group equations (RGEs) with respect to the standard virtuality
scale $\mu$ and a dimensionless analog $\rho=\sqrt{\zeta}/m$ of the usual
Collins-Soper scale $\zeta$. Here $\rho$ has an immediate physical
interpretation as the boost of the heavy hadron in the lab frame. We then
proceeded to also evaluate the unpolarized TMD fragmentation factor at NLO,
allowing us to cross check our results for the heavy-quark TMD FF in the
large-mass and – using known results from the literature – the small-mass
limit, finding full agreement. We further showed how to extend the theory of
joint threshold and $k_{T}$ resummation to fully account for quark mass
effects of $\mathcal{O}(m^{n}b_{T}^{n})$, and used these novel factorization
theorems for the heavy-quark TMD FF to obtain a third independent cross check
on our NLO result in an asymptotic limit. We also compared our final results
for the exact heavy-quark TMD FF in $b_{T}$ space to a recent independent
calculation by another group. We found full agreement despite the highly
orthogonal ways in which singularities were isolated in the two intermediate
momentum-space calculations, which provides a strong independent cross check
of this key universal ingredient for heavy-quark TMD predictions. Finally, and
as a first application, we leveraged our results to obtain the exact quark-
mass dependence of the energy-energy correlator (EEC) in $e^{+}e^{-}$
collisions in the back-to-back limit, where the associated jet functions are
given by Mellin moments of our results. As a side note, we pointed out that if
the EEC is restricted to heavy hadrons ($\mathrm{EEC}_{H}$), substantial
simplifications occur in the back-to-back limit such that the
$\mathrm{EEC}_{H}$ at values below the quark mass is given by a simple product
of Wilson-line vacuum expectation values, even extending down to
nonperturbative values of the observable.

Our results strengthen the theoretical underpinnings of heavy-quark TMD FFs
and pave the way towards precision studies of their nonperturbative dynamics,
as well as for improving the description of heavy flavor differential
distributions in general-purpose Monte-Carlo programs. We look forward to
future applications of our results in precision predictions for transverse
momentum distributions of heavy quarks and hadrons, and to comparing them to
future experimental measurements to shed light on the three-dimensional
dynamics of heavy-quark hadronization.

## Acknowledgments

We thank Iain Stewart, Kyle Lee, and Maximilian Stahlhofen for helpful
discussions. We thank the Erwin Schrödinger Institute for hospitality while
parts of this work were performed. This work was supported in part by the
Office of Nuclear Physics of the U.S. Department of Energy under Contract No.
DE-SC0011090. Z.S. was also supported by a fellowship from the MIT Department
of Physics. R.v.K. was supported by the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (Grant
agreement No. 101002090 COLORFREE). J.M. was supported by the D-ITP
consortium, a program of NWO that is funded by the Dutch Ministry of
Education, Culture and Science (OCW).

##  Appendix A Notation and conventions

###  A.1 Lightcone coordinates

We decompose four-vectors $p^{\mu}$ in terms of lightlike vectors $n^{\mu}$,
${\bar{n}}^{\mu}$ with $n^{2}=\bar{n}^{2}=0$ and $n\cdot{\bar{n}}=2$,

$\displaystyle p^{\mu}=\bar{n}\cdot p\,\frac{n^{\mu}}{2}+n\cdot
p\,\frac{\bar{n}^{\mu}}{2}+p_{\perp}^{\mu}\equiv
p^{-}\frac{n^{\mu}}{2}+p^{+}\frac{\bar{n}^{\mu}}{2}+p_{\perp}^{\mu}\equiv(p^{-},p^{+},p_{\perp})\,,$
(A.1)

such that e.g. $p^{2}=p^{-}p^{+}+p_{\perp}^{2}$. We take collinear momenta to
have large components $p^{-}\gg p_{\perp}\gg p^{+}$. We always take transverse
vectors with subscript $\perp$ to be Minkowskian, $p_{\perp}^{2}\equiv
p_{\perp}\cdot p_{\perp}\leq 0$, and denote their magnitude by
$p_{T}=\sqrt{-p_{\perp}^{2}}$ such that $p_{T}^{2}\geq 0$ is the square of a
real number.

###  A.2 Plus distributions

One-dimensional plus distributions are defined through their action on test
functions $g(x)$,

$\displaystyle\int\\!\mathrm{d}x\bigl{[}f(x)\bigr{]}_{+}\,g(x)\equiv\int\\!\mathrm{d}x\,f(x)\bigl{[}g(x)-g(0)\,\Theta(1-x)\bigr{]}\,,$
(A.2)

where $f(x)$ diverges at most as $x^{-1-\alpha}$ with $\alpha<1$ for $x\to 1$
and $g(x)$ is assumed to be differentiable as $x\to 0$. Integrals of
$[f(x)]_{+}$ against the constant test function over the domain $[0,1]$ vanish
by construction. The generalization of eq. (A.2) to two dimensions is given in
eq. (3.1) in the main text. For logarithmic plus distributions with
homogeneous power counting $\sim 1/x$ and their analogs for dimensionful
variables, we further define the shorthands

$\displaystyle\mathcal{L}_{n}(x)\equiv\Bigl{[}\frac{\Theta(x)\,\ln^{n}x}{x}\Bigr{]}_{+}\,,\quad\mathcal{L}_{n}(k,\mu)\equiv\frac{1}{\mu}\mathcal{L}_{n}\Bigl{(}\frac{k}{\mu}\Bigr{)}\,,\quad\mathcal{L}_{n}(k_{T}^{2},\mu^{2})\equiv\frac{1}{\mu^{2}}\mathcal{L}_{n}\Bigl{(}\frac{k_{T}^{2}}{\mu^{2}}\Bigr{)}\,.$
(A.3)

The definition in eq. (A.2) enables expansions in a regulator
$\alpha\sim\epsilon,\eta$ as follows,

$\displaystyle\frac{\Theta(x)}{x^{1+\alpha}}$
$\displaystyle=-\frac{\delta(x)}{\alpha}+\Bigl{[}\frac{\Theta(x)}{x^{1+\alpha}}\Bigr{]}_{+}$
$\displaystyle=-\frac{\delta(x)}{\alpha}+\mathcal{L}_{0}(x)-\alpha\mathcal{L}_{1}(x)+\mathcal{O}(\alpha^{2})\,,$
(A.4) $\displaystyle\frac{\Theta(x)\,\ln x}{x^{1+\alpha}}$
$\displaystyle=-\frac{\delta(x)}{\alpha^{2}}+\Bigl{[}\frac{\Theta(x)\,\ln
x}{x^{1+\alpha}}\Bigr{]}_{+}$
$\displaystyle=-\frac{\delta(x)}{\alpha^{2}}+\mathcal{L}_{1}(x)+\mathcal{O}(\alpha)\,.$
(A.5)

To relate plus distributions $\mathcal{L}_{n}(k,\mu)$ with different boundary
conditions $\mu_{1,2}$, it is useful to shift one of the boundary conditions
by making use of identities like

$\displaystyle\mathcal{L}_{0}(k,\mu_{2})$
$\displaystyle=\mathcal{L}_{0}(k,\mu_{1})+\delta(k)\,\ln\frac{\mu_{1}}{\mu_{2}}\,,$
$\displaystyle\mathcal{L}_{1}(k,\mu_{2})$
$\displaystyle=\mathcal{L}_{1}(k,\mu_{1})+\mathcal{L}_{0}(k,\mu_{1})\,\ln\frac{\mu_{1}}{\mu_{2}}+\delta(k)\,\frac{1}{2}\ln^{2}\frac{\mu_{1}}{\mu_{2}}\,.$
(A.6)

##  Appendix B Perturbative ingredients

In this appendix we collect various perturbative ingredients from elsewhere in
the literature in our notation, as used for cross-checks of our results in the
main text. The following expressions are all accurate to one-loop order.

#### Hard function.

The hard function describing the pair production of quarks through a virtual
photon or $Z$ boson in $e^{+}e^{-}$ collisions is given by

$\displaystyle\mathcal{H}_{ee\to ij}(Q^{2},\mu)$
$\displaystyle=\delta_{i\bar{j}}\,\Bigl{\\{}e_{i}^{2}-2v_{e}v_{i}e_{i}\operatorname{Re}\bigl{[}P_{Z}(Q^{2})\bigr{]}+(v_{e}^{2}+a_{e}^{2})(v_{i}^{2}+a_{i}^{2})\bigl{\lvert}P_{Z}(Q^{2})\bigr{\rvert}^{2}\Bigr{\\}}\bigl{\lvert}C_{q}\bigr{\rvert}^{2}\,,$
$\displaystyle C_{q}(Q^{2},\mu)$
$\displaystyle=1+\frac{\alpha_{s}C_{F}}{4\pi}\Bigl{[}-\ln^{2}\frac{-Q^{2}-\mathrm{i}0}{\mu^{2}}+3\ln\frac{-Q^{2}-\mathrm{i}0}{\mu^{2}}-8+\frac{\pi^{2}}{6}\Bigr{]}\,.$
(B.1)

Here we have kept the contribution from $Z$ boson exchange and $Z$-photon
interference, as relevant for measurements on the $Z$ pole, where
$P_{Z}(Q^{2})=Q^{2}/(Q^{2}-m_{Z}^{2}+\mathrm{i}\Gamma_{Z}m_{Z})$ is the
reduced $Z$ propagator and $e_{f}$ ($v_{f}$, $a_{f}$) are the electromagnetic
charge (vector, axial coupling to the $Z$) of a fermion $f$.

#### Large-mass and joint-resummation limits.

The combined collinear and soft mass matching coefficient $C_{m}$ see also
ref. [1] for our notation, is given by [47]

$\displaystyle
C_{m}(m,\mu,\zeta)=1+\frac{\alpha_{s}C_{F}}{4\pi}\biggl{(}4\ln^{2}\frac{\mu}{m}+2\ln\frac{\mu}{m}+4+\frac{\pi^{2}}{6}\biggr{)}+\mathcal{O}\bigl{(}\alpha_{s}^{2}\bigr{)}\,.$
(B.2)

At this order, it is also equal to the matching coefficient $H_{m,n}^{Q\to
h_{v}}(m,\mu,\nu/\sqrt{\zeta})$ (defined in ref. [47], without the
superscript) for matching massive collinear quark modes with label momentum
$\sqrt{\zeta}$ onto bHQET on their own, whereas the additional contribution
$\sqrt{H_{m,s}}(m,\mu,\nu)$ from matching secondary soft heavy quarks onto
light soft modes becomes nontrivial at $\mathcal{O}(\alpha_{s}^{2})$. Two-loop
expressions for both are given in ref. [47], where they were extracted from
the NNLO heavy-quark form factor [79, 80]. At that order, there is also a
nonzero $\zeta$ dependence that arises from the rapidity renormalization of
$H_{m,n}^{Q\to h_{v}}$ and $H_{m,s}$. The two-loop expression for secondary
heavy-quark effects on the light-to-light collinear matching coefficient
$H_{m,n}^{q\to q}$ is given in ref. [48] as $H_{c}$, where it was extracted
from the results of refs. [81, 47]. The perturbative HQET shape function is
given by [5, 6]

$\displaystyle\hat{F}\Bigl{(}\frac{k^{-}}{v^{-}},\mu\Bigr{)}=\delta\Bigl{(}\frac{k^{-}}{v^{-}}\Bigr{)}+\frac{\alpha_{s}C_{F}}{4\pi}\Bigl{[}-8\mathcal{L}_{1}\Bigl{(}\frac{k^{-}}{v^{-}},\mu\Bigr{)}-4\mathcal{L}_{0}\Bigl{(}\frac{k^{-}}{v^{-}},\mu\Bigr{)}-\frac{\pi^{2}}{6}\delta\Bigl{(}\frac{k^{-}}{v^{-}}\Bigr{)}\Bigr{]}$
(B.3)

The renormalized TMD soft function for both the $\eta$ regulator of ref. [43,
44] and the exponential regulator of ref. [71] reads

$\displaystyle\sqrt{S}(b_{T},m,\mu,\nu)=1+\frac{\alpha_{s}C_{F}}{4\pi}\Bigl{[}-L_{b}^{2}+4L_{b}\ln\frac{\mu}{\nu}-\frac{\pi^{2}}{6}\Bigr{]}\,.$
(B.4)

The joint-resummation collinear-soft function is given by [72]

$\displaystyle\mathcal{S}(\ell^{-},b_{T},m,\mu,\nu)=\delta(\ell^{-})+\frac{\alpha_{s}C_{F}}{4\pi}\Bigl{[}-8L_{b}\,\mathcal{L}_{0}(\ell^{-},\nu)\Bigr{]}\,.$
(B.5)

#### Small-mass limit.

The massless TMD FF matching coefficients are given by [82, 33, 51]

$\displaystyle\mathcal{J}_{q/q}(z,b_{T},\mu,\zeta)$
$\displaystyle=\delta(1-z)+\frac{\alpha_{s}C_{F}}{4\pi}\biggl{[}\frac{4(1+z^{2})\ln(z)}{1-z}-{4L_{b}\mathcal{L}_{0}(1-z)}$
$\displaystyle\quad+\delta(1-z)\Bigl{(}-{L_{b}^{2}}-2L_{b}\ln\frac{\zeta}{\mu^{2}}-\frac{\pi^{2}}{6}\Bigr{)}+{2L_{b}(z+1)}+2(1-z)\biggr{]}\,,$
$\displaystyle\mathcal{J}_{g/q}(z,b_{T},\mu,\zeta)$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\biggl{[}\frac{4(z^{2}-2z+2)\ln
z}{z}-\frac{2L_{b}(z^{2}-2z+2)}{z}+{2z}\biggr{]}\,,$
$\displaystyle\mathcal{J}_{g/g}(z,b_{T},\mu,\zeta)$
$\displaystyle=\delta(1-z)+\frac{\alpha_{s}C_{A}}{4\pi}\biggl{[}\delta(1-z)\Bigl{(}-{L_{b}^{2}}-2L_{b}\ln\frac{\zeta}{\mu^{2}}-\frac{\pi^{2}}{6}\Bigr{)}-{4L_{b}\mathcal{L}_{0}(1-z)}$
$\displaystyle\quad+\frac{8(z^{2}-z+1)^{2}\ln
z}{(1-z)z}+\frac{4L_{b}(z^{3}-z^{2}+2z-1)}{z}\biggr{]}$
$\displaystyle\mathcal{J}_{q/g}(z,b_{T},\mu,\zeta)$
$\displaystyle=\frac{\alpha_{s}T_{F}}{4\pi}\bigl{[}2(4z^{2}-4z+2)\ln
z-2L_{b}(2z^{2}-2z+1)+4(1-z)z\bigr{]}\,.$ (B.6)

The collinear heavy-quark FFs $d_{Q/Q}(z,\mu)$ and $d_{Q/g}(z,\mu)$ are given
by [52]

$\displaystyle d_{Q/Q}(z,\mu)$
$\displaystyle=\delta(1-z)+\frac{\alpha_{s}}{4\pi}C_{F}\Bigl{\\{}\ln\frac{\mu^{2}}{m^{2}}\bigl{[}4\mathcal{L}_{0}(1-z)+3\delta(1-z)-2(z+1)\bigr{]}$
$\displaystyle\quad-4\mathcal{L}_{0}(1-z)-8\mathcal{L}_{1}(1-z)+4\delta(1-z)+2(z+1)\bigl{[}2\ln(1-z)+1\bigr{]}\Bigr{\\}}\,,$
$\displaystyle d_{Q/g}(z,\mu)=$
$\displaystyle\frac{\alpha_{s}}{4\pi}T_{F}\bigl{[}2z^{2}+2(1-z)^{2}\bigr{]}\ln\frac{\mu^{2}}{m^{2}}\,,\,.$
(B.7)

The decoupling kernels relating collinear FFs in theories with
$n_{f}=n_{\ell}+1$ and $n_{\ell}$ active flavors are given by [53]

$\displaystyle\mathcal{M}_{g/g}(z,\mu)$
$\displaystyle=\delta(1-z)+\frac{\alpha_{s}T_{F}}{4\pi}\,\delta(1-z)\Bigl{(}-\frac{4}{3}\Bigr{)}\ln\frac{\mu^{2}}{m^{2}}$
$\displaystyle\mathcal{M}_{g/Q}(z,\mu)$
$\displaystyle=\frac{\alpha_{s}}{4\pi}C_{F}\frac{1+(1-z)^{2}}{z}\Bigl{(}2\ln\frac{\mu^{2}}{m^{2}}-2-4\ln
z\Bigr{)}\,.$ (B.8)

In both eqs. (B) and (B) we have rewritten the original results in terms of
the minimal distributional basis of $\delta(1-z)$ and $\mathcal{L}_{n}(1-z)$
using eq. (3.37).

##  Appendix C Results for Mellin moments

We first consider the Mellin moments of the momentum-space TMD FFs and
matching coefficients, which are defined in analogy to eq. (7.1) as

$\displaystyle d_{1\,Q/i}^{[N]}(k_{T},\mu,\zeta)$
$\displaystyle\equiv\int\\!\mathrm{d}z\,z^{N+2}d_{1\,Q/i}(z,k_{T},\mu,\zeta)\,,$
$\displaystyle\mathcal{J}_{k/i}^{[N]}(k_{T},m,\mu,\zeta)$
$\displaystyle\equiv\int\\!\mathrm{d}z\,z^{N}\mathcal{J}_{k/i}(z,k_{T},m,\mu,\zeta)\,.$
(C.1)

For $d_{1\,Q/Q}$, the relevant two-dimensional plus distribution simply turns
into

$\displaystyle\int_{0}^{1}\\!\mathrm{d}z\,z^{N}\,\bigl{[}f(x,z)\bigr{]}_{+,+}=\Bigl{[}\int_{0}^{1}\\!\mathrm{d}z\,\bigl{(}z^{N}-1\bigr{)}f(x,z)\Bigr{]}\,,$
(C.2)

which like the moments of the one-dimensional plus distribution in $z$ can
readily be evaluated for general $N$. The remaining one-dimensional plus
brackets in eq. (C.2) act on the $x=k_{T}^{2}/m^{2}$ dependence and can
further be decomposed in terms of singular $1/x$ terms and regular terms of
$\mathcal{O}(x^{0})$ in analogy to eq. (3.37), which is compensated by an
additional $\delta(x)$ boundary condition. As an explicit result, we find for
$N=1$ at one loop as relevant for section 7:

$\displaystyle d_{1\,Q/Q}^{[1]\,(1)}(k_{T},\mu,\zeta)$