values) that can be fixed by taking the (massless) limit $y\to 0$, where the
integrals are known Becher:2006qw . We give the explicit solutions for
$M_{1}$-$M_{7}$ in appendix D. We checked all of them numerically to the
required order in $\varepsilon$ using the Sector Decomposition codes FIESTA4
Smirnov:2015mct  and pySecDec Borowka:2017esm .

In order to take the imaginary part of the planar contributions we proceed as
follows. Defining $f(s)$ as the contribution of the planar diagrams at the
correlator level we want to determine

$\displaystyle\mathrm{Im}[f(s+i0)]=\frac{1}{2i}\mathrm{Disc}[f(s)]=\frac{1}{2i}\lim_{\beta\to
0}\big{[}f(s+i\beta)-f(s-i\beta)\big{]}\,.$ (76)

The result involves the distributions $\delta(s)$ and $\mathcal{L}_{n}(s)$ and
their derivatives prior to pole mass renormalization. (These derivatives only
arise from diagrams that contain quark self-energy subdiagrams.) In practice
we compute eq. (76) for $s>0$ first. We write the HPLs from the solutions of
the master integrals with the help of the HPL package in a form, where the
branch cuts are made explicit as $\log^{n}(-s/\mu^{2})$ terms. We can now
easily take the imaginary part of these terms according to eq. (76). The
remaining terms are regular in $s$ and thus do not contribute. In the next
step we promote the terms $\log^{n}(s/\mu^{2})/s$ and
$[n\log^{n-1}(s/\mu^{2})-\log^{n}(s/\mu^{2})]/s^{2}$ for $s>0$ to the
corresponding plus distributions $\mathcal{L}_{n}(s)$ and their derivatives,
respectively. We fix the coefficient of the $\delta(s)$ term by evaluating the
imaginary part of the integral of $f(s+i0)$ over a small region around $s=0$:

$\displaystyle\mathrm{Im}\bigg{[}\int_{-\Lambda_{1}}^{\Lambda_{2}}\\!\\!\mathrm{d}s\;f(s+i0)\bigg{]}\quad\text{with}\quad\Lambda_{1},\Lambda_{2}>0\,.$
(77)

The result is independent of $\Lambda_{1}$, as the imaginary part of $f(s+i0)$
does not have support for $s<0$. For $\Lambda_{2}\to 0$ the result involves a
constant and terms $\log^{k}(\Lambda_{2})/\Lambda_{2}^{m}$ which match the
structure of the distributions $\mathcal{L}_{n}(s)$ and their derivatives. The
constant term determines the coefficient of $\delta(s)$ in
$\mathrm{Im}[f(s+i0)]$. The coefficients of first and second derivatives of
$\delta(s)$ are determined in an analogous way considering eq. (77) with
additional weight functions $s$ and $s^{2}$, respectively. Upon imposing pole
mass renormalization all derivatives of $\delta(s)$ and $\mathcal{L}_{n}(s)$
consistently cancel, which provides an important cross check.

####  5.2.2 Nonplanar Diagrams

For the calculation of the nonplanar diagrams (a) and (d) we proceed in the
same way as described for the planar diagrams in the previous subsection.
Regarding the IBP reductions, no rapidity regulator is required for diagram
(a), because it does not involve any Wilson line propagators. In the reduction
of diagram (d) issues related to spurious rapidity divergences may in
principle arise. However, we were able to exclude by hand IBP relations such
as eq. (72), which generate (possibly ill-defined) rapidity divergent
integrals at intermediate steps. Interestingly, diagram (d) then happens to
reduce to already known family A master integrals. We explicitly checked the
IBP reduction numerically using pySecDec Borowka:2017esm  together with our
analytical results for the master integrals.

For diagram (a) family B master integrals cannot be avoided and a suitable
basis for them is depicted in figure 3. It contains so-called sunrise and kite
integrals for which generic results can be found in the literature. All three
master integrals admit a triple massive particle cut, which is reflected by
the $\Theta(p^{2}-(3m)^{2})=\Theta(1-9y)$ terms in their imaginary part, and
evaluate to (integrals of) elliptic functions. The master integrals contribute
to the jet field correlator with prefactors proportional to
$(-s-i0)^{-2\varepsilon}$ and upon $\varepsilon$-expansion in practice we only
need the imaginary part of the combination $(-s-i0)^{-2\varepsilon}M_{i}$ when
$s>0$. To determine the coefficients of the $\delta(s)$ and its derivatives we
additionally need the first few terms in the expansion of the master integrals
in $s=p^{2}\bar{y}=p^{2}(1-y)$. All these results are given in appendix D.
They are straightforwardly determined or directly taken from refs.
Laporta:2004rb ; Bauberger:1994by ; Remiddi:2016gno .

$M_{8}$

$M_{9}$

$M_{10}$

Figure 3: Master integrals of family B contributing to diagram (a). Bold solid
lines represent massive propagators, thin solid lines represent massless
propagators. A dot on a line means that the associated propagator is squared.
The solutions for these offshell sunrise and kite type integrals can be found
in the literature. We give the results for their imaginary parts and small $s$
expansion necessary for our calculation in appendix D.

To compute diagrams (b) and (c) we are forced to take a different approach.
The reason is that diagram (b) and (c) already involve rapidity divergent
family B integrals before any reduction. These rapidity divergences cancel in
the sum for each diagram. An IBP reduction method without additional
regulator, as used for diagram (d) therefore seems impractical. Using a
rapidity regulator, on the other hand, results in new types of master
integrals which involve three massive propagators, several massless
propagators, as well as regularized Wilson line propagators. Such integrals
appear too difficult to be solved with the tools at hand. We therefore decided
to refrain from using IBP reductions and compute them directly using numerical
methods and exploiting consistency conditions on their analytical structure as
explained in the following.

Numerical calculation of the contributions from diagrams (b) and (c)

In a first step we use the sector decomposition code pySecDec Borowka:2017esm
to obtain numerical results for the real and imaginary parts of the
contributions of diagrams (b) and (c) to the jet field correlator expanded in
$\varepsilon$ up to ${\mathcal{O}}(\varepsilon^{0})$. To this end we set
$\mu^{2}=s$, separately write the integrands of the two diagrams in Feynman
parameter representation and solve as many integrals analytically as possible.
For both diagrams four integrations over rational functions of the parameters
are remaining and are carried out with pySecDec. We furthermore multiply the
integrands by $p^{2}(1-y)^{2}=s^{2}/p^{2}$ to render them dimensionless and to
damp the singular behavior in the (bHQET) limit $y\to 1$. We then obtained
numerical results at $100$ points in the physical region $0\leq y<1$. From
these numerical data points we verified that the coefficients of the
$1/\varepsilon$ poles together with the divergences of the already computed
diagrams correctly reproduce within numerical uncertainties (which are
typically at the sub-percent level) the divergence structure related to the
renormalization $Z$-factor in eq. (60). This provides an important cross check
for our approach, in particular confirms the correctness concerning the $\mu$
dependence.

Due to the multiplication of the factor $p^{2}(1-y)^{2}=s^{2}/p^{2}$ we do not
have direct numerical access to the coefficient of the $\delta(s)$ term. We
will come back to this point at the end of this subsection. However, analytic
information on the contribution of diagrams (b) and (c) to the jet-function ME
is known in various kinematic limits. Using IBP reduction121212Note that in
the massless case no rapidity divergences arise. and the master integrals
given in Becher:2006qw , it is straightforward to compute the diagrams in the
massless limit. Upon taking the imaginary part the result, including mirror
diagrams, reads

$\displaystyle\\!\mathscr{J}^{\mathrm{b}+\mathrm{c}}_{\mathrm{fin}}(p^{2},0)={}$
$\displaystyle
a_{s}^{2}\left(C_{F}^{2}-\frac{C_{F}C_{A}}{2}\right)\Bigg{[}\left(-\frac{10\zeta_{3}}{3}-68+2\pi^{2}+\frac{\pi^{4}}{6}\right)\delta(p^{2})$
$\displaystyle+\left(-40\zeta_{3}+24+\frac{4\pi^{2}}{3}\right)\mathcal{L}_{0}(p^{2})+\frac{16}{3}\pi^{2}\mathcal{L}_{1}(p^{2})-8\mathcal{L}_{2}(p^{2})\Bigg{]}\,,$
(78)

where the $1/\varepsilon$ divergences have been subtracted. We can also
extract the corresponding expression in the bHQET limit ($y\to 1$) from eq.
(46). Combining the results for the two-loop bHQET jet function
$J_{B}^{(n_{\ell})}$ Jain:2008gb  and the collinear mass mode matching
coefficient $H_{m,n}$ Hoang:2015vua  with our results for the other diagrams,
we find

$\displaystyle\\!\mathscr{J}^{\mathrm{b}+\mathrm{c}}_{\mathrm{fin}}(p^{2},m^{2})={}$
$\displaystyle
a_{s}^{2}\left(C_{F}^{2}-\frac{C_{F}C_{A}}{2}\right)\Bigg{\\{}\delta(s)\Bigg{[}\left(-72\zeta_{3}+24+18\pi^{2}\right)L_{m}$
$\displaystyle\quad+\left(\frac{8\pi^{2}}{3}-16\right)L_{m}^{2}-\frac{154\zeta_{3}}{3}-\frac{\pi^{4}}{45}+\frac{50\pi^{2}}{3}-68-24\pi^{2}\log
2\Bigg{]}$
$\displaystyle+\left[16L_{m}^{2}+\left(32-\frac{16\pi^{2}}{3}\right)L_{m}+64\zeta_{3}-18\pi^{2}\right]\mathcal{L}_{0}(s)$
$\displaystyle+\left[-64L_{m}+\frac{32\pi^{2}}{3}-32\right]\mathcal{L}_{1}(s)+40\mathcal{L}_{2}(s)+\mathcal{O}\\!\left(m^{-2}\right)\Bigg{\\}}\,.$
(79)

Note that this expression contains all distributional terms in
$\\!\mathscr{J}^{\mathrm{b}+\mathrm{c}}(p^{2},m^{2})$.

Figure 4: Left: Contribution of diagrams (b) and (c) to the jet-function ME
finite part multiplied with the singularity-damping factor $p^{2}(1-y)^{2}$ as
a function of $y\equiv m^{2}/p^{2}$ for $\mu^{2}=s$. The numerical data from
pySecDec is shown as the solid blue line. The associated numerical uncertainty
is too small to be visible in the plot. The corresponding heavy quark (bHQET)
limit, according to eq. (79), is shown as dashed green line. The dotted orange
line represents the value of the jet-function ME times damping factor in the
massless limit ($y=0$), i.e. to the $\mathcal{L}_{0}(p^{2})$ term in eq. (78).
Right: Difference between numerical data points and analytic bHQET limit,
normalized to the massless limit, for $\mu^{2}=s$. The error bars indicate the
uncertainty of the numerical integration using pySecDec.

In the left panel of figure 4 we plot our (interpolated) numerical result for
$p^{2}(1-y)^{2}\\!\mathscr{J}^{\mathrm{b}+\mathrm{c}}_{\mathrm{fin}}$ (solid
blue line) and the corresponding analytical expression for $s\ll m^{2}$
according to eq. (79) (dashed green line) for $\mu^{2}=s$ in the range
$0<y<1$. Note that for $s>0$ and $\mu^{2}=s$ only the $\mathcal{L}_{0}(s)$
term in eq. (79) survives. With the factor $p^{2}(1-y)^{2}$ and for
$\mu^{2}=s\to p^{2}>0$ the massless result in eq. (78) collapses to a single
value, namely the coefficient of the $\mathcal{L}_{0}(p^{2})$ distribution.
This is indicated by the dotted orange line which agrees with the numerical
data point at $y=0$ within numerical uncertainty. We also observe that in the
limit $y\to 1$ the numerical data and the bHQET curve correctly approach each
other nicely. To better illustrate this we show in the right panel of figure 4
the difference between numerical data points and bHQET result normalized to
the $y=0$ value (dotted orange line in the left panel) in the range $0.95<y<1$
including the numerical uncertainties as quoted by pySecDec. We clearly see
that for $y\gtrsim 0.996$ the difference is zero within the numerical
uncertainties. This represents a strong cross check of our combined numerical
and analytical evaluation.

In order to obtain a practical fit function that parametrizes our numerical
results respecting at the same time the analytic constraints from the bHQET
and massless limits we proceed as follows. We start with an ansatz for
$\\!\mathscr{J}^{\mathrm{b}+\mathrm{c}}_{\mathrm{fin}}(p^{2},m^{2},\mu^{2}=s)$
valid for $s>0$ of the form

$\displaystyle
s\,\\!\mathscr{J}^{\mathrm{b}+\mathrm{c}}_{\mathrm{fin,ansatz}}(p^{2},m^{2},\mu^{2}\\!=\\!s)=a_{s}^{2}\left(C_{F}^{2}-\frac{C_{F}C_{A}}{2}\right)\sum_{m=0}^{2}\sum_{n=0}^{3}c_{mn}(1-y)^{m}\log^{n}(1-y)\,.$
(80)

Note that we multiplied
$\\!\mathscr{J}^{\mathrm{b}+\mathrm{c}}_{\mathrm{fin}}(p^{2},m^{2})$ on the
LHS with the factor $s$. We can now fix131313We note that a numerical fit of
the coefficients in eq. (81) from the numerical data points provides results
that agree with the values shown within numerical uncertainties.

$c_{00}=-18\pi^{2}+64\zeta_{3},\qquad c_{01}=\frac{16}{3}(\pi^{2}-6),\qquad
c_{02}=16,\qquad c_{03}=0,$ (81)

by comparing eq. (80) in the (bHQET) limit $y\to 1$ to the coefficient of
$\mathcal{L}_{0}(s)$ in eq. (79). Similarly, we can take the massless limit
($y\to 0$) of eq. (80) and compare it to the coefficient of
$\mathcal{L}_{0}(p^{2})$ in eq. (78), which fixes

$c_{10}=-c_{00}-c_{20}+24+\frac{4}{3}\pi^{2}-40\zeta_{3}\,.$ (82)

It is straightforward to generalize the ansatz function to arbitrary values of
$\mu^{2}\neq s$ using the renormalization $Z$-factor in eq. (60) and our
results from the other diagrams. The full distributional ansatz
$\\!\mathscr{J}^{\mathrm{b}+\mathrm{c}}_{\mathrm{fin,ansatz}}(p^{2},m^{2},\mu^{2})$
for arbitrary $\mu$ and $s\geq 0$ is then obtained by promoting the
$\log^{n}(s/\mu^{2})/s$ terms to $\mathcal{L}_{n}(s)$ distributions and adding
the $\delta(s)$ term of eq. (79). To also include the nontrivial analytic
information from the $\delta(p^{2})$ term in the massless quark result of eq.
(78), we integrate the full distributional ansatz over a finite interval
around $s=0$, take the limit $m\to 0$ and compare to the corresponding
integral of eq. (78). In this way we find the additional constraint

$\displaystyle c_{20}={}$
$\displaystyle-\frac{\pi^{2}}{6}c_{01}+2\zeta_{3}c_{02}-\frac{\pi^{4}}{15}c_{03}-\frac{\pi^{2}}{6}c_{11}+2\zeta_{3}c_{12}-\frac{\pi^{4}}{15}c_{13}+\left(1-\frac{\pi^{2}}{6}\right)c_{21}+2(\zeta_{3}-1)c_{22}$
$\displaystyle+\left(6-\frac{\pi^{4}}{15}\right)c_{23}-56\zeta_{3}-\frac{17\pi^{4}}{90}-8\pi^{2}(3\log
2-2)\,.$ (83)

From the same computation (using arbitrary integration limits) one obtains
analogous constraints for $c_{01}$, $c_{02}$, and $c_{03}$ which fully agree
with eq. (81). This represents a nontrivial check of our values for $c_{01}$,
$c_{02}$, and $c_{03}$ and implies consistency with the $L_{m}$ dependent
terms in the $\delta(s)$ coefficient of eq. (79).

Figure 5:  Fit function $G_{\mathrm{fit}}$ (dashed orange line) as given in
eq. (28) compared to the (interpolated) numerical data from pySecDec (solid
blue line). The numerical uncertainties are too small to be visible in the
plot.

Furthermore, one can argue on physical grounds that

$\displaystyle c_{13}=0\,,$ (84)

because we expect the first subleading bHQET power correction at two loops to
have a similar logarithmic structure as eq. (79), i.e. is should not contain
any cubic logarithmic terms.141414Analogously, one could argue that there is
no $c_{23}$ term. Unlike the $c_{13}$ term it is however regular for $y\to 1$
and we will use it to model other regular terms in our fit. A dedicated
numerical analysis in the bHQET regime confirms eq. (84) and moreover yields

$c_{12}=1.0\,.$ (85)

The values for the remaining coefficients in eq. (80) are obtained by a
weighted fit to the numerical data from pySecDec. We find

$c_{11}={}-41.9008,\qquad c_{21}={}-20.2744,\qquad c_{22}={}-10.5870,\qquad
c_{23}={}-7.1277\,.$ (86)

With these values the ansatz approximates our numerical data points to a
precision better than $0.5\%$, while the uncertainty of the points themselves
is at most $1\%$.

To check consistency of our numerical computation with the coefficient of
$\delta(s)$ term given in eq. (79), we performed the same fit as described
above, but without including the additional constraint of eq. (5.2.2). We find
that the fitted value for $c_{20}$ satisfies eq. (5.2.2) within the numerical
uncertainties. Since the massless limit for diagrams (b) and (c) was
explicitly computed, this consistency check reliably verifies the analytic
input for the $\delta(s)$ term in eq. (79).

Finally, we separate the analytical part in eq. (80) from the purely numerical
part given by the terms associated with the coefficients in eqs. (85) and
(86), which we call $G_{\mathrm{fit}}$, see eq. (28). We plot the fit function
$G_{\mathrm{fit}}(y)$ together with the associated numerical results in figure
5. Since the relative contribution of $G_{\mathrm{fit}}(y)$ to the entire
${\mathcal{O}}(\alpha_{s}^{2})$ corrections of the jet-function ME is smaller
then $0.5\%$ in the whole physical range, we do not bother to quote
uncertainties in eqs. (85) and (86).

###  5.3 Secondary Massive Quark Corrections

In this section we compute the $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$
secondary massive quark corrections to the jet-function ME
$\\!\mathscr{J}_{f}(p^{2},m^{2})$. The relevant Feynman graphs are diagrams
(o), (p), (q) and (r) in figure 1, where the vacuum polarization subdiagram is
a massive quark $Q$ bubble. We adopt the approach of refs.
Gritschacher:2013pha ; Pietrulewicz:2014qza  where the corrections are
calculated starting from the ${\cal O}(\alpha_{s})$ graphs describing
radiation of a gluon with mass $M$. The corresponding one-loop diagrams are
shown in figure 6. In a second step a dispersion relation to account for the
gluon splitting into a pair of a quark and antiquark with mass $m$ is applied.
Using the dispersion relation the gluon propagator with the insertion of the
bare massive quark vacuum polarization function in Feynman gauge and including
the (infinitesimal) gluon mass regulator $\Lambda$ can be written as

$\displaystyle\frac{-ig^{\mu\rho}}{p^{2}-\Lambda^{2}+i0}\;\Pi_{\rho\sigma}(p^{2},m^{2})\,\frac{-ig^{\sigma\nu}}{p^{2}-\Lambda^{2}+i0}={}$
$\displaystyle\frac{1}{\pi}\int\limits_{4m^{2}}^{\infty}\\!\frac{\mathrm{d}{M}^{2}}{{M}^{2}}\,\frac{-i(g^{\mu\nu}-\frac{p^{\mu}p^{\nu}}{p^{2}})}{p^{2}-{M}^{2}+i0}\;\mathrm{Im}\left[\Pi({M}^{2},m^{2})\right]$
$\displaystyle-\frac{-i\left(g^{\mu\nu}-\frac{p^{\mu}p^{\nu}}{p^{2}}\right)}{p^{2}-\Lambda^{2}+i0}\;\Pi(0,m^{2})\,.$
(87)

Here $\Pi(p^{2},m^{2})$ is the gluon vacuum polarization function arising from
a massive quark bubble defined by

$\displaystyle\Pi_{\mu\nu}^{AB}(p^{2},m^{2})=-i(p^{2}g_{\mu\nu}-p_{\mu}p_{\nu})\Pi(p^{2},m^{2})\delta^{AB}\equiv\int\mathrm{d}x^{4}\,\mathrm{e}^{ipx}\langle
0|\mathrm{T}J_{\mu}^{A}(x)J_{\nu}^{B}(0)|0\rangle\,,$ (88)

with the vector current $J^{A}_{\mu}\equiv ig{\bar{Q}}T^{A}\gamma_{\mu}Q$. The
imaginary part of the one-loop coefficient $\Pi^{(1)}$ in $d=4-2\varepsilon$
dimensions evaluates to
($\tilde{\mu}\equiv\mu\,e^{\gamma_{E}/2}(4\pi)^{-1/2}$)

$\displaystyle\mathrm{Im}\big{[}\Pi^{(1)}(p^{2},m^{2})\big{]}=\Theta\big{(}p^{2}\\!-\\!4m^{2}\big{)}T_{F}\left(\frac{p^{2}}{\tilde{\mu}^{2}}\right)^{-\varepsilon}\frac{2^{4\varepsilon}\pi^{\frac{3}{2}+\varepsilon}}{\Gamma(\frac{5}{2}-\varepsilon)}\bigg{(}\frac{2m^{2}}{p^{2}}+1-\varepsilon\bigg{)}\bigg{(}1-\frac{4m^{2}}{p^{2}}\bigg{)}^{\\!\\!\frac{1}{2}-\varepsilon}.$
(89)

Note that the Heaviside $\Theta$-function restricts the gluon mass integration
in the dispersion relation to $M>2m$. We could therefore safely set
$\Lambda=0$ in the first term of eq. (87).

(o′)(p′)(q′)(r′)

Figure 6: One-loop Feynman diagrams with massive gluons required for the
calculation of the secondary mass effects, i.e. the
$\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ terms, in the jet-function ME using
the dispersion relation in eq. (87). They are associated with the diagrams
(o), (p), (q) and (r) in figure 1 with a massive quark $Q$ bubble. Left-right
mirror graphs are understood.

Using eq. (87) for the gluon propagator with a heavy quark loop insertion in
figure 1 (o), (p), (q) and (r) we see that we effectively need to compute the
${\cal O}(\alpha_{s})$ jet-function ME $\\!\mathscr{J}_{f}^{(1)}$ with finite
gluon mass $M$, i.e. the corresponding diagrams in figure 6. This is needed
for the contribution from the first term on the RHS of eq. (87) (before
integration over $M^{2}$). To include the contribution from the second term we
also need the result expanded for $M\to\Lambda\ll M,p^{2}$.

The $p^{\mu}p^{\nu}$ terms on the RHS of eq. (87) can be dropped in our
computation due to gauge invariance of quantum electrodynamics with a massive
photon field Stueckelberg:1957zz ; Feldman:1963zz  and the fact that the
secondary $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ massive quark corrections
are up to color factors identical in an Abelian gauge theory. Dropping the
$p^{\mu}p^{\nu}$ terms in eq. (87) we immediately see that diagrams (q′) and
(r′) are proportional to $\bar{n}\cdot\bar{n}=0$ and thus vanish.

The contribution of diagram (o′) to eq. (4) before taking the imaginary part
reads ($s\equiv p^{2}-m^{2}+i0$)

(o′)
$\displaystyle=-i\frac{4\,\alpha_{s}C_{F}\tilde{\mu}^{2\varepsilon}}{(2\pi)^{d}(\bar{n}\cdot
p)s^{2}}\int\mathrm{d}^{d}k\,\frac{(\bar{n}\cdot k)(d-2)s+(\bar{n}\cdot
p)(4m^{2}-(d-2)(2p\cdot k+s))}{[-(p+k)^{2}+m^{2}][-k^{2}+M^{2}]}$ (90)
$\displaystyle=\alpha_{s}C_{F}\,\frac{\tilde{\mu}^{2\varepsilon}\,\Gamma(\varepsilon)}{2^{1-2\varepsilon}\pi^{2-\varepsilon}}\frac{1}{s^{2}}\,\int_{0}^{1}\mathrm{d}u\,\frac{(\varepsilon-1)s(u-1)+2m^{2}((\varepsilon-1)u-1)}{({M}^{2}(1-u)+u(s(u-1)+m^{2}u))^{\varepsilon}}\,,$

where the integration variable $u$ represents a Feynman parameter. This
expression is well-behaved in the large $M$ limit. It is convenient to
separate the leading term for $M\to\infty$ in $d=4-2\varepsilon$ dimensions,

$(\text{o}^{\prime})_{M\to\infty}=\alpha_{s}C_{F}\Gamma(\varepsilon)2^{-1+2\varepsilon}\pi^{-2+\varepsilon}\,\left(\frac{{M}}{\tilde{\mu}}\right)^{-2\varepsilon}\frac{2(2\varepsilon-3)m^{2}+(-1+\varepsilon)^{2}s}{(\varepsilon-2)(\varepsilon-1)s^{2}}\,.$
(91)

In the remaining contribution, which vanishes for $M^{2}\to\infty$, we can
safely set $d=4$:

$\displaystyle(\text{o}^{\prime})-(\text{o}^{\prime})_{M\to\infty}$
$\displaystyle={}-\frac{C_{F}\alpha_{s}}{2\pi^{2}}\,\frac{1}{s^{2}}\int_{0}^{1}\mathrm{d}u\,\left[s(u-1)+2m^{2}(1+u)\right]$
(92)
$\displaystyle\,\,\,\times\left[\log\Big{(}{M}^{2}(1-u)\Big{)}-\log\Big{(}{M}^{2}(1-u)+u(s(u-1)+m^{2}u)\Big{)}\right]+\mathcal{O}(\varepsilon)\,.$

We now take the imaginary part according to eq. (76). For the non-logarithmic
$1/s$ and $1/s^{2}$ terms we use $\mathrm{Im}[(s+i0)^{-1}]=-\pi\delta(s)$ and
$\mathrm{Im}[(s+i0)^{-2}]=\pi\delta^{\prime}(s)$, respectively. In addition,
we get a nonzero contribution from the second logarithm in eq. (92) for
$p^{2}>({M}+m)^{2}$ which also restricts the relevant integration range of the
Feynman parameter $u$. Adding the ${\mathcal{O}}(\alpha_{s})$ counterterm from
the renormalization of $m$ in the pole mass scheme at tree level we obtain
($\xi\equiv\sqrt{1-4m^{2}/M^{2}}$)

$\displaystyle\mathrm{Im}\big{[}\text{(o${}^{\prime}$)}\big{]}-2\,a_{s}m\,\delta
m_{M}^{(1)}\,\delta^{\prime}(s)={}a_{s}C_{F}$
$\displaystyle\times\Bigg{[}-\Theta\left(p^{2}-(m+M)^{2}\right)\frac{\sqrt{(M^{2}-s)^{2}-4m^{2}M^{2}}\left[4m^{4}+2m^{2}(M^{2}+2s)+s(M^{2}-s)\right]}{s^{2}(m^{2}+s)^{2}}$
$\displaystyle\quad-2\,\delta(s)\,\mathrm{e}^{\gamma_{E}\varepsilon}\left(\frac{M}{\tilde{\mu}}\right)^{-2\varepsilon}\Gamma(\varepsilon)\frac{\varepsilon-1}{\varepsilon-2}+\delta(s)\frac{3\,\xi\left(2m^{4}-M^{4}\right)\log\left(\frac{m^{2}}{M^{2}}\right)}{2m^{4}\xi}$
$\displaystyle\quad+\delta(s)\frac{3\left[\left(-4m^{4}-2m^{2}M^{2}+M^{4}\right)\log\left(\frac{1-\xi}{1+\xi}\right)-m^{2}\xi(3m^{2}+2M^{2})\right]}{2m^{4}\xi}\Bigg{]}+\mathcal{O}(\varepsilon)\,.$
(93)

Note that the $\delta^{\prime}(s)$ terms from diagram (o′) and the
contribution due to the quark pole mass counterterm $\delta m_{M}$ (where the
subscript indicates the nontrivial dependence on the gauge boson mass $M$)
cancel in the sum.

Diagram (p′) is rapidity divergent. With the symmetric $\eta$ regulator
Chiu:2011qc ; Chiu:2012ir  we have

(p′)
$\displaystyle=-i\frac{8\,\alpha_{s}C_{F}\tilde{\mu}^{2\varepsilon}\nu^{\eta}}{(2\pi)^{d}}\frac{1}{s}\int\mathrm{d}^{d}k\,\frac{\bar{n}\cdot(p+k)}{[-(p+k)^{2}+m^{2}][-k^{2}+M^{2}][-\bar{n}\cdot
k]^{1+\eta}}$
$\displaystyle=\alpha_{s}C_{F}\frac{\tilde{\mu}^{2\varepsilon}\Gamma(\varepsilon)}{\pi^{2-\varepsilon}2^{1-2\varepsilon}}\left(\frac{\bar{n}\cdot
p}{\nu}\right)^{-\eta}\frac{1}{s}\int_{0}^{\infty}\mathrm{d}u\,\mathrm{d}v\,\frac{vu^{-\eta-1}(u+v)^{2\varepsilon+\eta-2}\delta(1-u-v)}{\left(m^{2}u^{2}+v\left(M^{2}(u+v)-su\right)\right)^{\varepsilon}}\,.$
(94)

To compute the integral we use the same approach as for diagram (o′). The
$M\to\infty$ piece is given by

$(\text{p}^{\prime})_{M\to\infty}=\alpha_{s}C_{F}2^{-1+2\varepsilon}\pi^{-2+\varepsilon}\,\left(\frac{\bar{n}\cdot
p}{\nu}\right)^{-\eta}\left(\frac{{M}}{\tilde{\mu}}\right)^{-2\varepsilon}\frac{\Gamma(\varepsilon)\Gamma(2-\varepsilon)\Gamma(-\eta)}{\Gamma(2-\varepsilon-\eta)}\frac{1}{s}\,,$
(95)

and the remainder term reads

$\displaystyle(\text{p}^{\prime})-(\text{p}^{\prime})_{M\to\infty}={}$
$\displaystyle\frac{C_{F}\alpha_{s}}{2\pi^{2}}\frac{1}{s}\int_{0}^{1}\mathrm{d}u\frac{1-u}{u}$
(96)
$\displaystyle\times\Big{[}\log\Big{(}M^{2}(1-u)\Big{)}-\log\Big{(}m^{2}u^{2}+(1-u)\left(M^{2}-su\right)\Big{)}\Big{]}+{\mathcal{O}}(\varepsilon,\eta)\,.$

At this point the rapidity regularization has done its job and is in
particular not needed anymore for the dispersion integration. One may
therefore expand in $\eta$ keeping only the divergent and finite terms. To
take the imaginary part we proceed in the same way as for diagram (o′).

Accounting for the equivalent contribution from the mirror diagram of (p′) the
complete result for the bare primary massive quark jet-function ME with gluon
mass $M$ to ${\cal O}(\alpha_{s})$ reads

$\displaystyle\\!\mathscr{J}_{Q}^{\mathrm{bare}}(p^{2},m^{2},M^{2})=\delta(s)$
$\displaystyle+a_{s}C_{F}\Bigg{\\{}\Theta\left(p^{2}-(m+M)^{2}\right)\Bigg{[}\frac{4\log\left(\frac{\sqrt{\left(M^{2}-s\right)^{2}-4m^{2}M^{2}}+M^{2}+s}{-\sqrt{\left(M^{2}-s\right)^{2}-4m^{2}M^{2}}+M^{2}+s}\right)}{s}$
$\displaystyle\hskip
142.26378pt-\frac{\left(2m^{2}+s\right)\sqrt{\left(M^{2}-s\right)^{2}-4m^{2}M^{2}}\left(2m^{2}+M^{2}+3s\right)}{s^{2}\left(m^{2}+s\right)^{2}}\Bigg{]}$
$\displaystyle+\delta(s)\Bigg{[}2\mathrm{e}^{\gamma\varepsilon}\Gamma(\varepsilon)\left(\frac{M^{2}}{\tilde{\mu}^{2}}\right)^{-\varepsilon}\left(\frac{2}{\eta}-2\log\left(\frac{Q}{\nu}\right)+2H_{1-\varepsilon}-\frac{\varepsilon-1}{\varepsilon-2}\right)$
$\displaystyle+\frac{\left(4m^{4}-10m^{2}M^{2}+3M^{4}\right)\log\left(\frac{1-\xi}{1+\xi}\right)}{2m^{4}\xi}$
$\displaystyle-\frac{m^{2}\left(-8m^{2}\log\left(\frac{1-\xi}{2}\right)\log\left(\frac{\xi+1}{2}\right)+m^{2}+6M^{2}\right)+\left(2m^{4}-4m^{2}M^{2}+3M^{4}\right)\log\left(\frac{m^{2}}{M^{2}}\right)}{2m^{4}}\Bigg{]}\Bigg{\\}}$
$\displaystyle+\mathcal{O}(\varepsilon,a_{s}^{2})\,,$ (97)

where $m$ is the quark pole mass and $H_{n}$ denotes the $n$-th harmonic
number. For the expression shown here we have kept the exact dependence on
$d=4-2\varepsilon$ for the terms that need regularization in the
$M^{2}\to\infty$ limit of the dispersion integration. Otherwise the expansions
in $\varepsilon$ and $\eta$ are carried out.

From eq. (97) we can determine the ${\cal O}(\alpha_{s})$ bare primary massive
quark jet-function ME with the infinitesimal gluon mass regulator $\Lambda$,

$\displaystyle\\!\mathscr{J}_{Q}^{\mathrm{bare,(1)}}(p^{2},m^{2},\Lambda^{2})={}C_{F}\Bigg{\\{}\left[\frac{4}{\varepsilon}-4\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)+{\mathcal{O}}(\varepsilon)\right]\delta(s)\frac{1}{\eta}+\left[3-4\log\left(\frac{Q}{\nu}\right)\right]\delta(s)\frac{1}{\varepsilon}$
$\displaystyle\quad+\delta(s)\left[2L_{m}^{2}+L_{m}+2\log^{2}\left(\frac{\Lambda^{2}}{\mu^{2}}\right)+4\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)\log\left(\frac{Q}{\nu}\right)+8\right]$
(98)
$\displaystyle\quad+\mathcal{L}_{0}(s)\left[-4L_{m}-4\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)-4\right]+8\mathcal{L}_{1}(s)+\Theta(s)\left[\frac{s}{\left(m^{2}+s\right)^{2}}-\frac{4\log\left(\frac{s}{m^{2}}+1\right)}{s}\right]\Bigg{\\}}\,.$

We stress once again that a dimensionful regularization parameter such as the
gluon mass is mandatory to fully separate the rapidity and IR singularities
such that all remaining $1/\varepsilon$ terms are UV.

Using the dispersion identity eq. (87) the unrenormalized secondary massive
quark ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ correction to the jet-function ME
$\\!\mathscr{J}_{Q}(p^{2},m^{2})$ with an infinitesimal gluon mass regulator
reads

$\displaystyle\left.\\!\mathscr{J}_{Q}^{\mathrm{bare,(2)}}(p^{2},m^{2},\Lambda^{2})\right|_{C_{F}T_{F}}$
$\displaystyle={}\frac{1}{\pi}\int\limits_{4m^{2}}^{\infty}\frac{\mathrm{d}{M}^{2}}{{M}^{2}}\\!\mathscr{J}_{Q}^{\mathrm{bare},(1)}(p^{2},m^{2},{M}^{2})\,\mathrm{Im}\left[\Pi^{(1)}({M}^{2},m^{2})\right]$
$\displaystyle\qquad-\left(\Pi^{(1)}(m^{2},0)-\frac{4T_{F}}{3\varepsilon}\right)\\!\mathscr{J}_{Q}^{\mathrm{bare,(1)}}(p^{2},m^{2},\Lambda^{2})$
(99)
$\displaystyle={}\left.Z_{\\!\mathscr{J}}^{(2)}(s,m^{2},\Lambda^{2},\mu,Q/\nu)\right|_{C_{F}T_{F}}+\left.\\!\mathscr{J}_{Q}^{(2)}(p^{2},m^{2},\Lambda^{2},\mu,Q/\nu)\right|_{C_{F}T_{F}}\,.$

The $4T_{F}/(3\varepsilon)$ term in the second line of eq. (5.3) is associated
with the counterterm coming from the heavy quark loop and is required to
implement the ${\overline{\mathrm{MS}}}$ strong coupling in the $(n_{\ell}+1)$
flavor scheme. The dispersion integral in the first line of eq. (5.3) can be
solved analytically in terms of elementary functions and standard complete
elliptic integrals using Mathematica with appropriate changes of integration
variables, except for the logarithmic real radiation term shown in the second
line of eq. (97). The latter can be integrated numerically very efficiently.
It contributes to the corrections related to the three-particle
$QQ\bar{Q}$-cut and corresponds to eq. (36). The final result for the finite
(renormalized) contribution
$\left.\\!\mathscr{J}_{Q}^{(2)}\right|_{C_{F}T_{F}}$ is displayed in eq. (64)
and the divergent terms $\left.Z_{\\!\mathscr{J}}^{(2)}\right|_{C_{F}T_{F}}$
are given in eq. (65).

##  6 Numerical Analysis

In this section we provide a brief numerical study concerning the structure
and importance of the ${\mathcal{O}}(\alpha_{s}^{2})$ corrections to the
primary massive quark SCET jet function. To be definite we use the universal
jet function $J_{Q}^{\mathrm{uf}(n_{\ell}+1)}(p^{2},m^{2},\mu)$ of eq. (39)
for the analysis since in practice it is sufficient to formulate factorization
theorems for all possible scale hierarchies, see the discussion in section 4.
We remind the reader that the universal SCET jet function agrees with the
massless quark SCET jet function in the limit $m\to 0$, called “massless
limit” in the following. However, its relation to the bHQET jet function in
the limit $s=p^{2}-m^{2}\to 0$ involves the collinear-soft function $S_{c}$
that is singular by itself in this kinematic regime, see eqs. (11), (41) and
(46). Since the bHQET factorization theorem (43) and the factorization theorem
for $p^{2}-m^{2}\sim p^{2}\lesssim m^{2}$ (4) (which is identical to (4)) can
be smoothly connected, a comparison of the universal function
$J_{Q}^{\mathrm{uf}(n_{\ell}+1)}$ with its leading singular limit for
$s=p^{2}-m^{2}\to 0$, called “heavy quark limit” in the following, correctly
visualizes the size of the mass corrections that are not captured by the
leading power bHQET factorization theorem in the region
$0<s=p^{2}-m^{2}\lesssim m^{2}$.

Figure 7: Left: Size of contributions to the $\mathcal{O}(\alpha_{s}^{2})$
coefficient of the universal SCET massive primary quark jet function from
different color structures as a function of $s/m_{t}^{2}$ for
$\mu=\sqrt{s}+5\,\mathrm{GeV}$, $n_{\ell}=5$ and $m=m_{t}=173\,\mathrm{GeV}$
multiplied by $(s/m_{t}^{2})^{1.3}$. Right: Double logarithmic plot of the
modulus of the full ${\mathcal{O}}(\alpha_{s}^{2})$ corrections and the
corresponding contribution coming from $G_{\mathrm{fit}}$. We also show the
sum of all distributional ${\mathcal{O}}(\alpha_{s}^{2})$ corrections and the
non-distributional ones.

In the left panel of figure 7 the size of the contributions of the two-loop
coefficient
$J_{Q}^{(2)}(p^{2},m^{2},\mu)+J_{Q,\mathrm{sec}}^{(2),\mathrm{uf}}(p^{2},m^{2},\mu)$
coming from the different color structures is displayed as a function of
$s/m^{2}=(p^{2}-m^{2})/m^{2}$ for $\mu=\sqrt{s}+5\,\mathrm{GeV}$. The terms
proportional to $C_{F}^{2}$ (green), $C_{F}C_{A}$ (orange),
$C_{F}T_{F}n_{\ell}$ (blue) and $C_{F}T_{F}$ (red) are displayed individually
where we multiplied the factor $(s/m^{2})^{1.3}$ to reduce the variation of
the values over the whole kinematic regime.151515 This convention is also
adopted for figure 8. For all numerical discussions we adopt the
renormalization scale $\mu=\sqrt{s}+5\,\mathrm{GeV}$ and
$m=m_{t}=173\,\mathrm{GeV}$. The total coefficient is displayed in black. We
see that the $C_{F}^{2}$, $C_{F}C_{A}$ and $C_{F}T_{F}n_{\ell}$ terms, not
unexpectedly, provide the bulk of the ${\mathcal{O}}(\alpha_{s}^{2})$
corrections, and that the secondary heavy quark corrections are typically more
than one order of magnitude smaller. Interestingly, due to the behavior of the
$C_{F}^{2}$ and $C_{F}C_{A}$ terms, the ${\mathcal{O}}(\alpha_{s}^{2})$
correction changes sign at $\sqrt{p^{2}}\approx 1.1m$ which is well within the
bHQET regime.

In the left panel of figure 7 we also show the total contribution of the fit
function $G_{\mathrm{fit}}$ (purple dashed) which is part of the non-
distributional terms $G_{C_{F}}^{(1Q)}$ and $G_{C_{A}}^{(1Q)}$ contained in
the $C_{F}^{2}$ and $C_{F}C_{A}$ corrections, see eqs. (24) and (25)
respectively. Typically, it makes up for less than $1\%$ of the total
${\mathcal{O}}(\alpha_{s}^{2})$ correction, except close to the zero at
$\sqrt{p^{2}}\approx 1.1m$. This is also shown in more detail in the right
(double logarithmic) panel of figure 7 where the modulus of the full
${\mathcal{O}}(\alpha_{s}^{2})$ corrections (black) and the corresponding
contribution coming from the function $G_{\mathrm{fit}}$ (purple) are
displayed. Since the ${\mathcal{O}}(\alpha_{s})$ corrections to the jet
function do not have a zero (see below), the overall contribution of
$G_{\mathrm{fit}}$ to the entire jet function never exceeds the per mille
level. In the right panel of figure 7 we have also shown the sum of all
distributional ${\mathcal{O}}(\alpha_{s}^{2})$ corrections (brown) and the
non-distributional ones (light blue) which are constituted by the sum of all
$G$-functions quoted in section 3. We see that the non-distributional
contributions are strongly suppressed compared to the distributional terms in
the bHQET region and still at least an order of magnitude smaller as long as
$s\lesssim 10m^{2}$. In the high energy limit, however, there are strong
cancellations between them because the non-distributional corrections develop
distributional terms in the limit $m\to 0$, see appendix B. So the
contributions of the non-distributional corrections are particularly important
for $s/m^{2}\gtrsim 10$ where they are not at all negligible.

Figure 8: Left: ${\mathcal{O}}(\alpha_{s})$ coefficient of the SCET massive
primary quark jet function (black solid) in comparison to the leading power
results of the massless (black dotted) and bHQET (black dashed) limit. Also
shown are the curves including the respective next terms in the small mass and
bHQET expansions (gray dotted and dashed, respectively). Right:  Analogous
plot for the $\mathcal{O}(\alpha_{s}^{2})$ coefficient of the universal SCET
massive primary quark jet function.

We also would like to show the numerical impact of the genuine SCET massive
primary quark jet function in comparison to the previously known jet functions
for massless quarks (i.e. its high energy limit $s/m^{2}\to\infty$) and for
the bHQET limit $s/m^{2}\to 0$. In the left panel of figure 8 the full
${\mathcal{O}}(\alpha_{s})$ coefficient $J_{Q}^{(1)}(p^{2},m^{2},\mu)$ is
displayed (black solid) in comparison to the singular results of the massless
(black dotted) and bHQET (black dashed) limit. We see that both expansions
work well in the respective limits, but fail to provide reasonable
approximations in the intermediate regime $0.3\lesssim s/m^{2}\lesssim 10$
which corresponds to $1.1\lesssim\sqrt{p^{2}}/m\lesssim 3$. In the right panel
of figure 8 the analogous results are shown for the
${\mathcal{O}}(\alpha_{s}^{2})$ coefficient
$J_{Q}^{(2)}(p^{2},m^{2},\mu)+J_{Q,\mathrm{sec}}^{(2),\mathrm{uf}}(p^{2},m^{2},\mu)$.
Again, the massless as well as the bHQET expansion work well in the respective
limits. But, interestingly, both also approximate the exact result within
$15-20\%$ in the intermediate regime $0.3\lesssim s/m^{2}\lesssim 10$. That
the ${\mathcal{O}}(\alpha_{s})$ and ${\mathcal{O}}(\alpha_{s}^{2})$ correction
of the universal SCET massive primary quark jet function differ in this
respect should be considered as an accidental feature. This can be seen by
also including the respective next terms in the small mass and bHQET
expansions (which represent terms belonging to subleading power expansions in
massless quark SCET and bHQET, respectively). These are shown as the gray
dotted and dashed lines, respectively, in figure 8. At
${\mathcal{O}}(\alpha_{s})$ as well as ${\mathcal{O}}(\alpha_{s}^{2})$ these
subleading terms improve the approximation of the exact results substantially
into the intermediate regime $0.3\lesssim s/m^{2}\lesssim 10$. But we see
that, while at ${\mathcal{O}}(\alpha_{s})$ they also improve the approximation
to $s/m^{2}\approx 1$, the same is not true at
${\mathcal{O}}(\alpha_{s}^{2})$.

##  7 Conclusions

In this work we have calculated the ${\mathcal{O}}(\alpha_{s}^{2})$
corrections to the SCET jet functions for primary massive quarks which are
relevant for the inclusive description of jets initiated by massive quarks
produced by hard interactions with high energy. The result provides important
information for the intermediate region $p^{2}\sim p^{2}-m^{2}=s\sim m^{2}$,
where the previously known ${\mathcal{O}}(\alpha_{s}^{2})$ jet function
results in the massless limit $m^{2}/p^{2}\to 0$ and in the bHQET limit $s\to
0$ do not provide accurate approximations. Thus they represent the essential
input to achieve a coherent description that smoothly interpolates these two
extreme limits.

The interesting conceptual feature of the primary massive quark SCET jet
function is that it can be defined in two different ways depending on whether
- in addition to the zero-bin subtractions - one imposes soft mass mode bin
subtractions or not for diagrams containing secondary massive quark effects,
i.e. a massive quark-antiquark vacuum polarization subdiagram. In the context
of the view that jet functions require so-called zero-bin subtractions
Manohar:2006nz  to constitute infrared finite factorization functions, the
relevance of the soft mass mode bin subtraction Chiu:2009yx ;
Gritschacher:2013pha  depends on the assumed power counting of the quark mass
w.r.t. the jet invariant mass. An alternative view, that leads to equivalent
results, but does not rely on a zero-bin expansion supplemented by power
counting arguments, is to define the jet function with a subtraction defined
from a collinear-soft matrix element Bauer:2011uc ; Procura:2014cba  where the
effects of secondary massive quarks can be optionally included or not. The
collinear-soft matrix elements provide unambiguous prescriptions of the
subtractions eliminating the need to impose an additional power counting. The
definitions and the results for the collinear-soft matrix elements are
identical for massless and massive primary quarks and even apply for bHQET jet
functions (to deal with corrections arising from secondary massive quarks that
are lighter than the super-heavy primary quark).

If the secondary massive quark is included in the collinear-soft matrix
element the resulting jet function is called universal jet function. The
universal jet function reduces to the massless quark jet function in the limit
of zero quark mass and also has the same virtuality ($\mu$) anomalous
dimension. It can be used for the kinematic region ranging from $m^{2}\sim
p^{2}\sim p^{2}-m^{2}$ down to the massless limit $p^{2}\gg m^{2}$, where
$p^{2}$ is the squared jet function invariant mass and $m$ the mass of the
primary quark. If the secondary massive quarks are not included in the
collinear-soft matrix element the jet function is called mass mode jet
function. It is infrared divergent for $m\to 0$ and has a virtuality and a
rapidity anomalous dimension. It can be used for the kinematic region
$m^{2}\sim p^{2}\sim p^{2}-m^{2}$. The two types of jet functions are related
via the so-called collinear-soft function Pietrulewicz:2017gxc , and we have
provided the ${\mathcal{O}}(\alpha_{s}^{2})$ results for both. We have also
demonstrated how both jet functions are used in the factorization approaches
of refs. Gritschacher:2013pha ; Pietrulewicz:2014qza ; Hoang:2015vua  and ref.
Pietrulewicz:2017gxc  which were developed to deal with the effects of
secondary massive quarks in the context of factorization. While the approach
of ref. Pietrulewicz:2017gxc  provides a more transparent view on the modular
structure of the collinear, soft and collinear-soft mass modes, the approach
of refs. Gritschacher:2013pha ; Pietrulewicz:2014qza ; Hoang:2015vua  provides
a method with a smooth quark mass dependence that allows to treat the cases
$p^{2}-m^{2}\ll m^{2}\sim p^{2}$ and $p^{2}-m^{2}\sim m^{2}\sim p^{2}$ in a
single factorization theorem. Taking a jet mass factorization theorem with the
primary massive quark SCET jet function as an example we have shown how both
approaches are related and in which way they are or can be rendered
equivalent. We emphasize that, considered together, both approaches provide a
thorough view on the field theoretic treatment of secondary massive quark
effects within factorization approaches that separate soft and collinear
quantum effects.

Regarding the calculation, we have treated the two-loop secondary mass
corrections due to heavy quark loops separately from the corrections due to
massless quark and gluon loops. In fact, the correct treatment of the
secondary mass corrections requires to carefully disentangle UV, IR and
rapidity divergences present in this case. This is achieved by employing a
dimensionful IR regulator (in our case a small gluon mass). The loop
integrations are than carried out using the dispersive method of refs.
Gritschacher:2013pha ; Pietrulewicz:2014qza ; Hoang:2015vua .

The remaining ${\mathcal{O}}(\alpha_{s}^{2})$ corrections (related to diagrams
without a closed quark-antiquark loop) have been carried out in pure
dimensional regularization. We distinguished two types of relevant diagrams:
planar and nonplanar. While the contributions from the planar diagrams can be
obtained using standard modern multi-loop techniques, we faced interesting
technical issues in the calculation of the nonplanar diagrams. In particular,
for two nonplanar diagrams we have not been able to employ (standard
automated) IBP reductions without generating spurious rapidity divergences
that require an additional rapidity regulator. The latter made the solution of
the corresponding master integrals unfeasible. We therefore treated these two
diagrams in a semi-numerical approach exploiting analytic information from the
massless and bHQET limits. The corresponding numerical results are very
accurate and only represent a very small contribution to the full
${\mathcal{O}}(\alpha_{s}^{2})$ corrections to the jet functions. We have also
provided precise numerical approximations of our analytic results (involving
elliptic functions) for practical implementations.

## Acknowledgments

This work was supported by FWF Austrian Science Fund under the Project No.
P28535-N27 and in part by a GFK fellowship and by the clusters of excellence
PRISMA (EXC 1098) and PRISMA+ (EXC 2118) at JGU Mainz. We acknowledge partial
support by the FWF Austrian Science Fund under the Doctoral Program “Particles
and Interactions” No. W1252-N27. We also thank the Erwin-Schrödinger
International Institute for Mathematics and Physics for partial support. CL
thanks the JGU Mainz and the University of Salamanca for hospitality during
parts of this work. We thank A. Carmona, D. Lechner, M. Löschner, P.
Pietrulewicz, M. Preisser, and D. Samitz for helpful discussions.

##  Appendix A Other Jet Functions

The mass mode jet function for massless primary quarks at
$\mathcal{O}(\alpha_{s}^{2})$ has the form

$\displaystyle
J_{q}^{\mathrm{mf},(n_{\ell}+1)}\left(p^{2},m^{2},\mu,\frac{Q}{\nu}\right)={}$
$\displaystyle\delta(p^{2})+a_{s}^{(n_{\ell})}J_{q}^{(1)}(p^{2},\mu)$ (100)
$\displaystyle+\left(a_{s}^{(n_{\ell})}\right)^{2}\left[J_{q}^{(2)}(p^{2},\mu)+J_{q,\mathrm{sec}}^{(2),\mathrm{mf}}\left(p^{2},m^{2},\mu,\frac{Q}{\nu}\right)\right]+\mathcal{O}(\alpha_{s}^{3})\,,$

with

$\displaystyle
J_{q}^{(1)}(p^{2},\mu)={}C_{F}\Bigg{\\{}\left(7-\pi^{2}\right)\delta(p^{2})-3\mathcal{L}_{0}(p^{2})+4\mathcal{L}_{1}(p^{2})$
$\displaystyle\quad+\varepsilon\Bigg{[}\left(-\frac{28\zeta_{3}}{3}+14-\frac{3\pi^{2}}{4}\right)\delta(p^{2})+\left(\pi^{2}-7\right)\mathcal{L}_{0}(p^{2})+3\mathcal{L}_{1}(p^{2})-2\mathcal{L}_{2}(p^{2})\Bigg{]}$
$\displaystyle\quad+\varepsilon^{2}\Bigg{[}\left(-7\zeta_{3}+28-\frac{7\pi^{2}}{4}-\frac{\pi^{4}}{24}\right)\delta(p^{2})+\left(\frac{28\zeta_{3}}{3}-14+\frac{3\pi^{2}}{4}\right)\mathcal{L}_{0}(p^{2})$
$\displaystyle\qquad+\left(7-\pi^{2}\right)\mathcal{L}_{1}(p^{2})-\frac{3}{2}\mathcal{L}_{2}(p^{2})+\frac{2}{3}\mathcal{L}_{3}(p^{2})\Bigg{]}\Bigg{\\}}\,,$
(101) $\displaystyle
J_{q}^{(2)}(p^{2},\mu)={}C_{F}^{2}\Bigg{[}\left(-18\zeta_{3}+\frac{205}{8}-\frac{67\pi^{2}}{6}+\frac{14\pi^{4}}{15}\right)\delta(p^{2})+\left(-8\zeta_{3}-\frac{45}{2}+7\pi^{2}\right)\mathcal{L}_{0}(p^{2})$
$\displaystyle\quad+\left(37-\frac{20\pi^{2}}{3}\right)\mathcal{L}_{1}(p^{2})-18\mathcal{L}_{2}(p^{2})+8\mathcal{L}_{3}(p^{2})\Bigg{]}$
$\displaystyle\quad+C_{F}C_{A}\Bigg{[}\left(-\frac{206\zeta_{3}}{9}+\frac{53129}{648}-\frac{208\pi^{2}}{27}-\frac{17\pi^{4}}{180}\right)\delta(p^{2})+\left(40\zeta_{3}-\frac{3155}{54}+\frac{22\pi^{2}}{9}\right)\mathcal{L}_{0}(p^{2})$
$\displaystyle\qquad+\left(\frac{367}{9}-\frac{4\pi^{2}}{3}\right)\mathcal{L}_{1}(p^{2})-\frac{22}{3}\mathcal{L}_{2}(p^{2})\Bigg{]}$
$\displaystyle\quad+C_{F}T_{F}n_{\ell}\Bigg{[}\left(\frac{494}{27}-\frac{8\pi^{2}}{9}\right)\mathcal{L}_{0}(p^{2})+\left(\frac{16\zeta_{3}}{9}-\frac{4057}{162}+\frac{68\pi^{2}}{27}\right)\delta(p^{2})$
$\displaystyle\qquad-\frac{116}{9}\mathcal{L}_{1}(p^{2})+\frac{8}{3}\mathcal{L}_{2}(p^{2})\Bigg{]}\,,$
(102) $\displaystyle
J_{q,\mathrm{sec}}^{(2),\mathrm{mf}}\left(p^{2},m^{2},\mu,\frac{Q}{\nu}\right)={}C_{F}T_{F}\Bigg{\\{}\Bigg{[}\left(-\frac{8}{3}L_{m}^{2}-\frac{80}{9}L_{m}-\frac{224}{27}\right)\log\left(\frac{Q}{\nu}\right)+2L_{m}^{2}$
$\displaystyle\qquad+\frac{2}{9}\left(3+4\pi^{2}\right)L_{m}-\frac{8\zeta_{3}}{3}+\frac{20\pi^{2}}{27}+\frac{73}{18}\Bigg{]}\delta(p^{2})+G_{\mathrm{sec}}\Theta\left(p^{2}-(2m)^{2}\right)\Bigg{\\}},$
(103) $\displaystyle
G_{\mathrm{sec}}={}\frac{1}{p^{2}}\Bigg{[}-\frac{32}{3}\mathrm{Li}_{2}\\!\left(\frac{b-1}{1+b}\right)+\frac{16}{3}\log\left(\frac{1-b^{2}}{4}\right)\log\left(\frac{1-b}{1+b}\right)-\frac{8}{3}\log^{2}\left(\frac{1-b}{1+b}\right)$
$\displaystyle\quad+\left(\frac{1}{2}b^{4}-b^{2}+\frac{241}{18}\right)\log\left(\frac{1-b}{1+b}\right)-\frac{5}{27}b^{3}+\frac{241}{9}b-\frac{8\pi^{2}}{9}\Bigg{]}$
$\displaystyle\hphantom{G_{\mathrm{sec}}}=\frac{(p^{2}-(2m)^{2})^{7/2}}{560m^{9}}+\mathcal{O}\\!\left(\frac{(p^{2}-(2m)^{2})^{9/2}}{m^{11}}\right),$
(104)

where $b\equiv\sqrt{1-4m^{2}/p^{2}}$.

The universal jet function for massless primary quarks at
$\mathcal{O}(\alpha_{s}^{2})$ has the form Pietrulewicz:2014qza

$\displaystyle J_{q}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},m^{2},\mu)={}$
$\displaystyle\delta(p^{2})+a_{s}^{(n_{\ell}+1)}J_{q}^{(1)}(p,\mu)$
$\displaystyle+\left(a_{s}^{(n_{\ell}+1)}\right)^{2}\left[J_{q}^{(2)}(p,\mu)+J_{q,\mathrm{sec}}^{(2),\mathrm{uf}}(p^{2},m^{2},\mu)\right]+\mathcal{O}(\alpha_{s}^{3})\,,$
(105)

with

$\displaystyle
J_{q,\mathrm{sec}}^{(2),\mathrm{uf}}(p^{2},m^{2},\mu)=C_{F}T_{F}\Bigg{\\{}\Bigg{[}\frac{8}{9}L_{m}^{3}+\frac{58}{9}L_{m}^{2}+\left(\frac{718}{27}-\frac{8\pi^{2}}{9}\right)L_{m}$
$\displaystyle\qquad+\frac{1531}{54}+\frac{10\pi^{2}}{27}-\frac{80\zeta_{3}}{9}\Bigg{]}\delta(p^{2})+\left[-\frac{8}{3}L_{m}^{2}-\frac{116}{9}L_{m}-\frac{224}{27}\right]\mathcal{L}_{0}(p^{2})+\frac{16}{3}L_{m}\mathcal{L}_{1}(p^{2})$
$\displaystyle\quad+G_{\mathrm{sec}}\Theta\left(p^{2}-(2m)^{2}\right)\Bigg{\\}}\,.$
(106)

Both jet functions satisfy eq. (12) in the same way as the corresponding
primary massive jet functions.

The renormalized primary massless jet-function ME including contributions
related to 1-loop corrections and secondary massive quark production at 2-loop
in the computational scheme with an infinitesimal gluon mass $\Lambda$ is
given by

$\displaystyle\left.\\!\mathscr{J}\\!\\!\left(p^{2},0,\Lambda^{2},\mu,\frac{Q}{\nu}\right)\right|_{C_{F}T_{F}}=\delta(p^{2})$
(107)
$\displaystyle\quad+a_{s}^{(n_{\ell}+1)}C_{F}\Bigg{\\{}\left[2\log^{2}\left(\frac{\Lambda^{2}}{\mu^{2}}\right)+4\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)\log\left(\frac{Q}{\nu}\right)-\frac{2\pi^{2}}{3}+7\right]\delta(p^{2})$
$\displaystyle\qquad+\left[-4\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)-3\right]\mathcal{L}_{0}(p^{2})+4\mathcal{L}_{1}(p^{2})\Bigg{\\}}$
$\displaystyle\quad+\left(a_{s}^{(n_{\ell}+1)}\right)^{2}C_{F}T_{F}\Bigg{\\{}\Bigg{[}\left(\frac{8}{3}\log^{2}\left(\frac{\Lambda^{2}}{\mu^{2}}\right)+10\right)L_{m}+2L_{m}^{2}$
$\displaystyle\qquad\quad+\log\left(\frac{Q}{\nu}\right)\left(\frac{16}{3}\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)L_{m}-\frac{80}{9}L_{m}-\frac{8}{3}L_{m}^{2}-\frac{224}{27}\right)-\frac{8\zeta_{3}}{3}+\frac{20\pi^{2}}{27}+\frac{73}{18}\Bigg{]}\delta(p^{2})$
$\displaystyle\qquad+\left(-\frac{16}{3}\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)-4\right)L_{m}\mathcal{L}_{0}(p^{2})+\frac{16}{3}L_{m}\mathcal{L}_{1}(p^{2})+G_{\mathrm{sec}}\Theta(p^{2}-(2m)^{2})\Bigg{\\}}+{\mathcal{O}}(\alpha_{s}^{3})\,.$

The ${\cal O}(\alpha_{s})$ and ${\cal
O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$
corrections in the computational scheme with a massless gluon in dimensional
regularization analytically agree with those of the massless quark SCET jet
function, see Becher:2006qw ; Pietrulewicz:2014qza .

##  Appendix B $G$-Functions in the Limit $m\to 0$

In the limit that the massive primary quark $Q$ becomes massless the singular
limits of the non-distributional $G$-functions have the form

$\displaystyle G_{1}={}$
$\displaystyle\left(-2L_{m}^{2}-L_{m}-1-\frac{2\pi^{2}}{3}\right)\delta(p^{2})+(1+4L_{m})\mathcal{L}_{0}(p^{2})-4\mathcal{L}_{1}(p^{2})+\mathcal{O}\\!\left(\frac{m^{2}}{p^{4}}\right),$
(108) $\displaystyle G_{1}^{(\varepsilon)}={}$
$\displaystyle\left(-\frac{3}{2}L_{m}^{2}+\frac{2\pi^{2}}{3}L_{m}-2-\frac{\pi^{2}}{3}-4\zeta_{3}\right)\delta(p^{2})$
$\displaystyle+(1+4L_{m}+2L_{m}^{2})\mathcal{L}_{0}(p^{2})+(-5-8L_{m})\mathcal{L}_{1}(p^{2})+6\mathcal{L}_{2}(p^{2})+\mathcal{O}\\!\left(\frac{m^{2}}{p^{4}}\right),$
(109) $\displaystyle G_{1}^{(\varepsilon^{2})}={}$
$\displaystyle\Bigg{[}-\frac{1}{6}L_{m}^{4}-\frac{1}{6}L_{m}^{3}+\left(-4+\frac{\pi^{2}}{6}\right)L_{m}^{2}+\left(\frac{7\pi^{2}}{12}+4\zeta_{3}\right)L_{m}$
(110)
$\displaystyle\quad-4-\frac{13\pi^{2}}{12}-\frac{\pi^{4}}{30}-2\zeta_{3}\Bigg{]}\delta(p^{2})$
$\displaystyle+\left[\frac{2}{3}L_{m}^{3}+2L_{m}^{2}+(8-\pi^{2})L_{m}+2-\frac{\pi^{2}}{4}\right]\mathcal{L}_{0}(p^{2})$
$\displaystyle+[-4L_{m}^{2}-8L_{m}-9+\pi^{2}]\mathcal{L}_{1}(p^{2})+\Big{[}8L_{m}+\frac{13}{2}\Big{]}\mathcal{L}_{2}(p^{2})-\frac{14}{3}\mathcal{L}_{3}(p^{2})+\mathcal{O}\\!\left(\frac{m^{2}}{p^{4}}\right),$
$\displaystyle
G_{C_{F}}^{(1Q)}={}\Biggl{[}-2L_{m}^{4}-2L_{m}^{3}+\left(2\pi^{2}-\frac{35}{2}\right)L_{m}^{2}+L_{m}\left(56\zeta_{3}-\frac{27}{2}+\pi^{2}\right)$
(111)
$\displaystyle\qquad-68\zeta_{3}+\frac{43\pi^{4}}{30}+\frac{125}{2}+\pi^{2}\left(24\log
2-\frac{131}{6}\right)\Biggr{]}\delta(p^{2})$
$\displaystyle\quad+\left[8L_{m}^{3}+12L_{m}^{2}+\left(38-\frac{20\pi^{2}}{3}\right)L_{m}-72\zeta_{3}+\frac{\pi^{2}}{3}+\frac{33}{2}\right]\mathcal{L}_{0}(p^{2})$
$\displaystyle\quad+\left[-32L_{m}^{2}-40L_{m}+\frac{20\pi^{2}}{3}-45\right]\mathcal{L}_{1}(p^{2})+6[8L_{m}+5]\mathcal{L}_{2}(p^{2})-24\mathcal{L}_{3}(p^{2})+\mathcal{O}\\!\left(\frac{m^{2}}{p^{4}}\right),$
$\displaystyle
G_{C_{A}}^{(1Q)}={}\Bigg{[}\frac{1}{9}\left(6\pi^{2}-179\right)L_{m}^{2}+\frac{1}{9}L_{m}\left(-180\zeta_{3}-139+56\pi^{2}\right)+42\zeta_{3}$
$\displaystyle\qquad+\frac{7\pi^{4}}{180}-\frac{508}{9}+\pi^{2}\left(\frac{25}{27}-12\log
2\right)\Bigg{]}\delta(p^{2})$
$\displaystyle\quad+\frac{1}{9}\left[66L_{m}^{2}+\left(391-12\pi^{2}\right)L_{m}+180\zeta_{3}-34\pi^{2}+172\right]\mathcal{L}_{0}(p^{2})$
$\displaystyle\quad-\frac{4}{9}\left[66L_{m}-3\pi^{2}+106\right]\mathcal{L}_{1}(p^{2})+22\mathcal{L}_{2}(p^{2})+\mathcal{O}\\!\left(\frac{m^{2}}{p^{4}}\right),$
(112) $\displaystyle
G_{T_{F}}={}\left[\frac{58}{9}L_{m}^{2}+\left(\frac{74}{9}-\frac{16\pi^{2}}{9}\right)L_{m}+\frac{110}{9}+\frac{46\pi^{2}}{27}\right]\delta(p^{2})$
$\displaystyle\quad+\left[-\frac{8}{3}L_{m}^{2}-\frac{128}{9}L_{m}-\frac{86}{9}+\frac{8\pi^{2}}{9}\right]\mathcal{L}_{0}(p^{2})+\left[\frac{32}{3}L_{m}+\frac{140}{9}\right]\mathcal{L}_{1}(p^{2})$
$\displaystyle\quad-8\mathcal{L}_{2}(p^{2})+\mathcal{O}\\!\left(\frac{m^{2}}{p^{4}}\right),$
(113) $\displaystyle
G^{(3Q)}={}\left[L_{m}^{2}+7L_{m}+88\zeta_{3}-91+\pi^{2}(11-16\log
2)\right]\delta(p^{2})+[-2L_{m}-7]\mathcal{L}_{0}(p^{2})$
$\displaystyle\quad+2\mathcal{L}_{1}(p^{2})+\mathcal{O}\\!\left(\frac{m^{2}}{p^{4}}\right),$
(114) $\displaystyle
G_{\mathrm{sec}}^{(3Q)}={}\frac{2}{81}\left[-36L_{m}^{3}-261L_{m}^{2}+3\left(12\pi^{2}-359\right)L_{m}+216\zeta_{3}+171\pi^{2}-2783\right]\delta(p^{2})$
$\displaystyle\quad+\frac{2}{27}\left[36L_{m}^{2}+174L_{m}-12\pi^{2}+359\right]\mathcal{L}_{0}(p^{2})-\frac{4}{9}[12L_{m}+29]\mathcal{L}_{1}(p^{2})$
$\displaystyle\quad+\frac{8}{3}\mathcal{L}_{2}(p^{2})+\mathcal{O}\\!\left(\frac{m^{2}}{p^{4}}\right).$
(115)

##  Appendix C Anomalous Dimensions

The $\mu$-anomalous dimensions as defined in eq. (4) are given by

$\displaystyle\gamma_{c}^{(n_{f})}(Q,\mu)={}-\gamma_{Q}^{(n_{f})}(Q,\mu)\,,$
(116)
$\displaystyle\gamma_{Q}^{(n_{f})}(Q,\mu)={}\Gamma_{\mathrm{cusp}}^{(n_{f})}[\alpha_{s}]\log\Bigl{(}\frac{-Q^{2}-i0}{\mu^{2}}\Bigr{)}+\gamma_{Q}^{(n_{f})}[\alpha_{s}]\,,$
(117)
$\displaystyle\gamma_{Q}^{(n_{f})}[\alpha_{s}]={}-6\,C_{F}a_{s}^{(n_{f})}+\left(a_{s}^{(n_{f})}\right)^{2}\Bigg{[}C_{F}^{2}\left(-3+4\pi^{2}-48\zeta_{3}\right)$
$\displaystyle\qquad+C_{F}C_{A}\left(-\frac{961}{27}-\frac{11\pi^{2}}{3}+52\zeta_{3}\right)+C_{F}T_{F}n_{f}\left(\frac{260}{27}+\frac{4\pi^{2}}{3}\right)\Bigg{]}+\mathcal{O}(\alpha_{s}^{3})\,,$
(118)
$\displaystyle\gamma_{c_{B}}^{(n_{f})}(Q,m,\mu)={}-\Gamma_{\mathrm{cusp}}^{(n_{f})}[\alpha_{s}]\log\Bigl{(}\frac{-Q^{2}-i0}{m^{2}}\Bigr{)}+\gamma_{c_{B}}^{(n_{f})}[\alpha_{s}]\,,$
(119)
$\displaystyle\gamma_{c_{B}}^{(n_{f})}[\alpha_{s}]={}4\,C_{F}a_{s}^{(n_{f})}+\left(a_{s}^{(n_{f})}\right)^{2}\Bigg{[}C_{F}C_{A}\left(\frac{196}{9}-\frac{4\pi^{2}}{3}+8\zeta_{3}\right)+C_{F}T_{F}n_{f}\left(-\frac{80}{9}\right)\Bigg{]}$
$\displaystyle\quad+\mathcal{O}(\alpha_{s}^{3})\,,$ (120)
$\displaystyle\gamma_{J}^{(n_{f})}(p^{2},\mu)={}-2\Gamma_{\mathrm{cusp}}^{(n_{f})}[\alpha_{s}]\mathcal{L}_{0}(p^{2})+\gamma_{J}^{(n_{f})}[\alpha_{s}]\delta(p^{2})\,,$
(121)
$\displaystyle\gamma_{J}^{(n_{f})}[\alpha_{s}]={}6\,C_{F}a_{s}^{(n_{f})}+\left(a_{s}^{(n_{f})}\right)^{2}\Bigg{[}C_{F}^{2}\left(3-4\pi^{2}+48\zeta_{3}\right)$
$\displaystyle\qquad+C_{F}C_{A}\left(\frac{1769}{27}+\frac{22\pi^{2}}{9}-80\zeta_{3}\right)+C_{F}T_{F}n_{f}\left(-\frac{484}{27}-\frac{8\pi^{2}}{9}\right)\Bigg{]}+\mathcal{O}(\alpha_{s}^{3})\,,$
(122)
$\displaystyle\gamma_{J_{B}}^{(n_{f})}(\hat{s},\mu)={}-2\Gamma_{\mathrm{cusp}}^{(n_{f})}[\alpha_{s}]\mathcal{L}_{0}(\hat{s})+\gamma_{J_{B}}^{(n_{f})}[\alpha_{s}]\delta(\hat{s})\,,$
(123)
$\displaystyle\gamma_{J_{B}}^{(n_{f})}[\alpha_{s}]={}4\,C_{F}a_{s}^{(n_{f})}+\left(a_{s}^{(n_{f})}\right)^{2}\Bigg{[}C_{F}C_{A}\left(\frac{1396}{27}-\frac{23\pi^{2}}{9}-20\zeta_{3}\right)$
$\displaystyle\qquad+C_{F}T_{F}n_{f}\left(\frac{4\pi^{2}}{9}-\frac{464}{27}\right)\Bigg{]}+\mathcal{O}(\alpha_{s}^{3})\,,$
(124)
$\displaystyle\gamma_{S}^{(n_{f})}(\ell^{+},\ell^{-},\mu)={}\delta(\ell^{+})\gamma_{s}^{(n_{f})}(\ell^{-},\mu)+\delta(\ell^{-})\gamma_{s}^{(n_{f})}(\ell^{+},\mu)\,,$
(125)
$\displaystyle\gamma_{s}^{(n_{f})}(\ell,\mu)={}2\Gamma_{\mathrm{cusp}}^{(n_{f})}[\alpha_{s}]\mathcal{L}_{0}(\ell)+\gamma_{s}^{(n_{f})}[\alpha_{s}]\delta(\ell)\,,$
(126)
$\displaystyle\gamma_{s}^{(n_{f})}[\alpha_{s}]={}\left(a_{s}^{(n_{f})}\right)^{2}\Bigg{[}C_{F}C_{A}\left(-\frac{808}{27}+\frac{11\pi^{2}}{9}+28\zeta_{3}\right)+C_{F}T_{F}n_{f}\left(\frac{224}{27}-\frac{4\pi^{2}}{9}\right)\Bigg{]}$
$\displaystyle\quad+\mathcal{O}(\alpha_{s}^{3})\,,$ (127)
$\displaystyle\Gamma_{\mathrm{cusp}}^{(n_{f})}={}4\,C_{F}a_{s}^{(n_{f})}+\left(a_{s}^{(n_{f})}\right)^{2}\left[C_{F}C_{A}\left(\frac{268}{9}-\frac{4\pi^{2}}{3}\right)+C_{F}T_{F}n_{f}\left(-\frac{80}{9}\right)\right]+\mathcal{O}(\alpha_{s}^{3}).$
(128)

The $\mu$\- and $\nu$-anomalous dimensions of the jet-function ME
$\\!\mathscr{J}_{f}^{(n_{\ell}+1)}(p^{2},m^{2})$ and the collinear-soft MEs
$\mathscr{S}^{(n_{\ell}+1)}(\ell,m)$ and $\mathscr{S}^{(n_{\ell})}(\ell)$,
defined in analogy to eqs. (4) and (57), read

$\displaystyle\left.\gamma_{\\!\mathscr{J}}(p^{2},\mu,Q/\nu)\right|_{C_{F}T_{F}}=a_{s}^{(n_{\ell}+1)}C_{F}\left[6-8\log\left(\frac{Q}{\nu}\right)\right]\delta(p^{2})$
$\displaystyle\quad+\left(a_{s}^{(n_{\ell}+1)}\right)^{2}C_{F}T_{F}\left[\frac{160}{9}\log\left(\frac{Q}{\nu}\right)-\frac{16\pi^{2}}{9}-\frac{4}{3}\right]\delta(p^{2})+{\mathcal{O}}(\alpha_{s}^{3})\,,$
(129)
$\displaystyle\gamma^{(n_{\ell})}_{\mathscr{S}}(\ell,\mu,\nu)=a_{s}^{(n_{\ell})}C_{F}\left[8\log\left(\frac{\nu}{\mu}\right)\delta(\ell)+8\mathcal{L}_{0}(\ell)\right]+{\mathcal{O}}(\alpha_{s}^{2})\,,$
(130)
$\displaystyle\left.\gamma^{(n_{\ell}+1)}_{\mathscr{S}}(\ell,\mu,\nu)\right|_{C_{F}T_{F}}=a_{s}^{(n_{\ell}+1)}C_{F}\left[8\log\left(\frac{\nu}{\mu}\right)\delta(\ell)+8\mathcal{L}_{0}(\ell)\right]$
$\displaystyle\quad+\left(a_{s}^{(n_{\ell}+1)}\right)^{2}C_{F}T_{F}\Bigg{[}\Bigg{(}-\frac{160}{9}\log\left(\frac{\nu}{\mu}\right)-\frac{8\pi^{2}}{9}+\frac{448}{27}\Bigg{)}\delta(\ell)-\frac{160}{9}\mathcal{L}_{0}(\ell)\Bigg{]}+{\mathcal{O}}(\alpha_{s}^{3})\,,$
(131)
$\displaystyle\left.\gamma_{\nu,\\!\mathscr{J}}(m,\Lambda,\mu)\right|_{C_{F}T_{F}}=a_{s}^{(n_{\ell}+1)}C_{F}\left[-4\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)\right]$
$\displaystyle\quad+\left(a_{s}^{(n_{\ell}+1)}\right)^{2}C_{F}T_{F}\Bigg{[}-\frac{16}{3}\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)L_{m}+\frac{8}{3}L_{m}^{2}+\frac{80}{9}L_{m}+\frac{224}{27}\Bigg{]}+{\mathcal{O}}(\alpha_{s}^{3})\,,$
(132)
$\displaystyle\gamma^{(n_{\ell})}_{\nu,\mathscr{S}}(\Lambda,\mu)=a_{s}^{(n_{\ell})}C_{F}\left[-4\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)\right]+{\mathcal{O}}(\alpha_{s}^{2})\,,$
(133)
$\displaystyle\left.\gamma^{(n_{\ell}+1)}_{\nu,\mathscr{S}}(m,\Lambda,\mu)\right|_{C_{F}T_{F}}=a_{s}^{(n_{\ell}+1)}C_{F}\left[-4\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)\right]$
$\displaystyle\quad+\left(a_{s}^{(n_{\ell}+1)}\right)^{2}C_{F}T_{F}\Bigg{[}-\frac{16}{3}\log\left(\frac{\Lambda^{2}}{\mu^{2}}\right)L_{m}+\frac{8}{3}L_{m}^{2}+\frac{80}{9}L_{m}+\frac{224}{27}\Bigg{]}+{\mathcal{O}}(\alpha_{s}^{3})\,.$
(134)

The ${\cal
O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{A}C_{F},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$
contributions to the $\mu$\- and $\nu$-anomalous dimensions of the jet-
function and collinear-soft MEs are unknown. We therefore only displayed the
${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ corrections coming from the massive quark
$Q$, which is indicated by subscript “$C_{F}T_{F}$”.

The virtuality $\mu$-anomalous dimensions for the collinear-soft function
$S_{c}$ and the mass mode jet function $J_{f}^{\mathrm{mf},(n_{\ell}+1)}$ are

$\displaystyle\gamma_{S_{c}}(\ell,m,\mu,\nu)={}a_{s}^{2}C_{F}T_{F}\Bigg{[}\left(-\frac{32}{3}L_{m}-\frac{160}{9}\right)\mathcal{L}_{0}(\ell)$
$\displaystyle\qquad+\left(-\frac{32}{3}\log\left(\frac{\nu}{\mu}\right)L_{m}-\frac{160}{9}\log\left(\frac{\nu}{\mu}\right)-\frac{8\pi^{2}}{9}+\frac{448}{27}\right)\delta(\ell)\Bigg{]}+{\mathcal{O}}(\alpha_{s}^{3}),$
(135)
$\displaystyle\left.\gamma_{J^{\mathrm{mf}}}(p^{2},m^{2},\mu,Q/\nu)\right|_{C_{F}T_{F}}={}a_{s}^{(n_{\ell}+1)}C_{F}\left[-8\mathcal{L}_{0}(p^{2})+6\delta(p^{2})\right]$
$\displaystyle\quad+\left(a_{s}^{(n_{\ell}+1)}\right)^{2}C_{F}T_{F}\Bigg{[}\left(\frac{32}{3}L_{m}\log\left(\frac{Q}{\nu}\right)+\frac{160}{9}\log\left(\frac{Q}{\nu}\right)-\frac{16\pi^{2}}{9}-\frac{4}{3}\right)\delta(p^{2})$
$\displaystyle\qquad-\frac{32}{3}L_{m}\mathcal{L}_{0}(p^{2})\Bigg{]}+{\mathcal{O}}(\alpha_{s}^{3})\,,$
(136)

where the remaining color factor contributions to $\gamma_{J^{\mathrm{mf}}}$
at ${\cal O}(\alpha_{s}^{2})$ (which do not arise from diagrams containing the
massive quark $Q$ loop) coincide with the ${\cal O}(\alpha_{s}^{2})$
contributions of $\gamma_{J}^{(n_{\ell})}$ given in eq. (121).

##  Appendix D Master Integrals

The family A master integrals depicted in figure 2 evaluate to ($y\equiv
m^{2}/p^{2}$, $\bar{y}\equiv 1-y$)

$\displaystyle
M_{1}={}\mathcal{I}(-y;0,0,0,1,1,0,0,0)=(-y/\bar{y})^{2-2\varepsilon}\,\Gamma(\varepsilon-1)^{2}\,,$
(137) $\displaystyle
M_{2}={}\mathcal{I}(-y;1,0,0,1,1,0,0,0)=-(-y/\bar{y})^{2-2\varepsilon}\,\Gamma(\varepsilon-1)^{2}\,{}_{2}F_{1}(1,2-2\varepsilon;2-\varepsilon;1/\bar{y})\,,$
(138) $\displaystyle
M_{3}={}\mathcal{I}(-y;1,1,0,1,1,0,0,0)=(-y/\bar{y})^{2-2\varepsilon}\,\Gamma(\varepsilon-1)^{2}\,{}_{2}F_{1}(1,2-2\varepsilon;2-\varepsilon;1/\bar{y})^{2}\,,$
(139) $\displaystyle
M_{4}={}\mathcal{I}(-y;1,0,1,0,1,0,0,0)=-2\Gamma(1-\varepsilon)\Gamma(\varepsilon)\Gamma(2\varepsilon-2)\,_{2}F_{1}\left(2-2\varepsilon,2\varepsilon-1;2-\varepsilon;1/\bar{y}\right),$
(140) $\displaystyle
M_{5}={}\mathcal{I}(-y;1,0,1,0,2,0,0,0)=-\Gamma(1-\varepsilon)\Gamma(\varepsilon-1)\Gamma(2\varepsilon)\,{}_{2}F_{1}\left(2-2\varepsilon,2\varepsilon;2-\varepsilon;1/\bar{y}\right),$
(141) $\displaystyle
M_{6}={}\mathcal{I}(-y;1,0,1,1,1,0,0,1)=\mathrm{e}^{-2\varepsilon\gamma_{E}}\bigg{[}\frac{\pi^{2}}{12\varepsilon^{2}}$
$\displaystyle\quad+\frac{1}{\varepsilon}\left(\frac{1}{6}\pi^{2}H(-1;-y)-H(-2,-1;-y)+H(-2,0;-y)+\frac{7\zeta_{3}}{2}\right)$
$\displaystyle\quad+7\zeta_{3}H(-1;-y)-\frac{2}{3}\pi^{2}H(-2;-y)+H(-3,-1;-y)-H(-3,0;-y)+3H(-2,-2;-y)$
$\displaystyle\quad+\frac{1}{3}\pi^{2}H(-1,-1;-y)-3H(-2,0,0;-y)-2H(-1,-2,-1;-y)+2H(-1,-2,0;-y)$
$\displaystyle\quad+\frac{11\pi^{4}}{72}+\mathcal{O}(\varepsilon)\bigg{]}\,,$
(142) $\displaystyle
M_{7}={}\mathcal{I}(-y;1,0,1,1,1,0,0,2)=\mathrm{e}^{-2\varepsilon\gamma_{E}}\bigg{[}-\frac{1}{4\varepsilon^{3}}$
$\displaystyle\quad+\frac{1}{\varepsilon}\left(-H(-1,-1;-y)+H(-1,0;-y)+\frac{5\pi^{2}}{24}\right)-\frac{1}{6}\pi^{2}H(-1;-y)-2H(-2,-1;-y)$
$\displaystyle\quad+2H(-2,0;-y)+3H(-1,-2;-y)-3H(-1,0,0;-y)+\frac{29\zeta_{3}}{3}$
$\displaystyle\quad+\varepsilon\bigg{(}-2\zeta_{3}H(-1;-y)-\pi^{2}H(-2;-y)+2H(-3,-1;-y)-2H(-3,0;-y)$
$\displaystyle\qquad+4H(-2,-2;-y)-7H(-1,-3;-y)+\frac{5\pi^{2}}{6}H(-1,-1;-y)+\frac{\pi^{2}}{3}H(-1,0;-y)$
$\displaystyle\qquad-4H(-2,0,0;-y)+4H(-1,-2,-1;-y)-4H(-1,-2,0;-y)$
$\displaystyle\qquad-4H(-1,-1,-1,-1;-y)+4H(-1,-1,-1,0;-y)+7H(-1,0,0,0;-y)+\frac{497\pi^{4}}{1440}\bigg{)}$
$\displaystyle\quad+\mathcal{O}(\varepsilon^{2})\bigg{]}\,,$ (143)

where the $H(\vec{n};-y)$ denote HPLs according to the conventions of ref.
Maitre:2005uu .

The relevant imaginary parts and $\bar{y}$-expansions related to the family B
master integrals depicted in figure 3 were determined or directly taken from
refs. Laporta:2004rb ; Bauberger:1994by ; Remiddi:2016gno  and read

$\displaystyle
M_{8}={}\mathcal{I}(-y;0,0,0,1,1,1,0,0)=\mathrm{e}^{-2\varepsilon\gamma_{E}}\left(\frac{-y}{\bar{y}}\right)^{1-2\varepsilon}\biggl{[}-\frac{3}{2\varepsilon^{2}}-\frac{\frac{1}{-y}+18}{4\varepsilon}$
$\displaystyle\quad-\frac{1}{8}\left(59+2\pi^{2}\right)-\frac{3\bar{y}}{8}-\frac{1}{8}\left(\pi^{2}-7\right)\bar{y}^{2}-\frac{1}{96}\left(9\pi^{2}-56\right)\bar{y}^{3}+\mathcal{O}\\!\left(\varepsilon,\bar{y}^{4}\right)\biggr{]},$
(144)
$\displaystyle\mathrm{Im}\big{[}(-s-i0)^{-2\varepsilon}M_{8}\big{]}\overset{s>0}{=}{}\Theta(1/9-y)\frac{\pi}{2\bar{y}}\sqrt{1-\sqrt{y}}\sqrt{3\sqrt{y}+1}$
$\displaystyle\quad\times\Biggl{[}8y\,\mathrm{K}\left(\frac{\left(\sqrt{y}+1\right)^{3}\left(3\sqrt{y}-1\right)}{\left(\sqrt{y}-1\right)^{3}\left(3\sqrt{y}+1\right)}\right)+\left(\sqrt{y}-1\right)(1+3y)\mathrm{\,}\mathrm{E}\left(\frac{\left(\sqrt{y}+1\right)^{3}\left(3\sqrt{y}-1\right)}{\left(\sqrt{y}-1\right)^{3}\left(3\sqrt{y}+1\right)}\right)\Biggr{]}$
$\displaystyle\quad+\mathcal{O}(\varepsilon),$ (145) $\displaystyle
M_{9}={}\mathcal{I}(-y;0,0,0,2,1,1,0,0)=\mathrm{e}^{-2\varepsilon\gamma_{E}}\left(\frac{-y}{\bar{y}}\right)^{-2\varepsilon}\biggl{[}\frac{1}{2\varepsilon^{2}}+\frac{1}{2\varepsilon}-\frac{1}{12}\left(6-\pi^{2}\right)$
$\displaystyle\quad-\left(\frac{\pi^{2}}{12}-1\right)\bar{y}-\frac{1}{32}\left(4-\pi^{2}\right)\bar{y}^{2}+\mathcal{O}\\!\left(\varepsilon,\bar{y}^{3}\right)\biggr{]},$
(146)
$\displaystyle\mathrm{Im}\big{[}(-s-i0)^{-2\varepsilon}M_{9}\big{]}\overset{s>0}{=}{}-\Theta(1/9-y)\frac{\pi\sqrt{3\sqrt{y}+1}}{3\sqrt{1-\sqrt{y}}}$
$\displaystyle\quad\times\Biggl{[}4\left(2\sqrt{y}-1\right)\,\mathrm{K}\left(\frac{\left(\sqrt{y}+1\right)^{3}\left(3\sqrt{y}-1\right)}{\left(\sqrt{y}-1\right)^{3}\left(3\sqrt{y}+1\right)}\right)+3\left(\sqrt{y}-1\right)^{2}\,\mathrm{E}\left(\frac{\left(\sqrt{y}+1\right)^{3}\left(3\sqrt{y}-1\right)}{\left(\sqrt{y}-1\right)^{3}\left(3\sqrt{y}+1\right)}\right)\Biggr{]}$
$\displaystyle\quad+\mathcal{O}(\varepsilon)\,,$ (147) $\displaystyle
M_{10}={}\mathcal{I}(-y;1,1,0,1,1,1,0,0)=\bar{y}\Bigg{[}-\frac{1}{6}\pi^{2}H(1;1/y)-2H(0,1,1;1/y)$
$\displaystyle\quad+H(1,0,1;1/y)+\frac{2}{3}\int_{9}^{\infty}\\!\mathrm{d}t\,\frac{t-9}{t-1}\log\left(1-\frac{1}{yt}\right)\frac{2\,\mathrm{K}\left(\frac{\left(\sqrt{t}-3\right)\left(\sqrt{t}+1\right)^{3}}{\left(\sqrt{t}-1\right)^{3}\left(\sqrt{t}+3\right)}\right)}{\sqrt{\left(\sqrt{t}-1\right)^{3}\left(\sqrt{t}+3\right)}}\,\Bigg{]}+\mathcal{O}(\varepsilon)$
(148)
$\displaystyle\,={}\bar{y}\Bigg{[}\bigg{(}\frac{17\zeta_{3}}{4}-\frac{3}{2}\pi^{2}\log
2\bigg{)}+\bar{y}\left(3\log\left(-\bar{y}\right)-\log^{2}\left(-\bar{y}\right)-\frac{3\pi^{2}}{16}-\frac{5}{2}\right)\Bigg{]}+\mathcal{O}(\varepsilon,\bar{y}^{3})\,,$
(149)
$\displaystyle\mathrm{Im}\big{[}(-s-i0)^{-2\varepsilon}M_{10}\big{]}\overset{s>0}{=}{}\bar{y}\,\pi\Bigg{[}\Theta(1-y)\left[2\log(y)\log(\bar{y}/y)-3\text{Li}_{2}(-\bar{y}/y)\right]$
$\displaystyle\quad-\Theta(1/9-y)\frac{4}{3}\int_{9}^{1/y}\mathrm{d}t\,\frac{t-9}{t-1}\frac{\mathrm{K}\left(\frac{\left(\sqrt{t}-3\right)\left(\sqrt{t}+1\right)^{3}}{\left(\sqrt{t}-1\right)^{3}\left(\sqrt{t}+3\right)}\right)}{\sqrt{\left(\sqrt{t}-1\right)^{3}\left(\sqrt{t}+3\right)}}\,\Bigg{]}+\mathcal{O}(\varepsilon)\,.$
(150)

The remaining integral in eq. (150) can be evaluated numerically very
efficiently. It originates from the branch cut of the logarithm in the
explicit integral in the full $\varepsilon$ expanded result for $M_{10}$ given
in eq. (148) when taking the imaginary part. The expansion of $M_{10}$ around
$s\sim 0$, i.e. $y\sim 1$, in eq. (149) was computed by expanding the
integrand and solving the resulting integrals using a PSLQ-type algorithm
implemented in Mathematica.

All results stated in this appendix were checked numerically using sector
decomposition.

## References

  * (1) _Les Houches 2017: Physics at TeV Colliders Standard Model Working Group Report_ , 2018.
  * (2) CLICdp collaboration, H. Abramowicz et al., _Top-Quark Physics at the CLIC Electron-Positron Linear Collider_ , 1807.02441.
  * (3) HL-LHC, HE-LHC Working Group collaboration, P. Azzi et al., _Standard Model Physics at the HL-LHC and HE-LHC_ , 1902.04070.
  * (4) C. W. Bauer, S. Fleming, D. Pirjol and I. W. Stewart, _An Effective field theory for collinear and soft gluons: Heavy to light decays_ , _Phys. Rev._ D63 (2001) 114020, [hep-ph/0011336].
  * (5) C. W. Bauer, D. Pirjol and I. W. Stewart, _Soft collinear factorization in effective field theory_ , _Phys. Rev._ D65 (2002) 054022, [hep-ph/0109045].
  * (6) C. W. Bauer and I. W. Stewart, _Invariant operators in collinear effective theory_ , _Phys. Lett._ B516 (2001) 134–142, [hep-ph/0107001].
  * (7) E. Lunghi, D. Pirjol and D. Wyler, _Factorization in leptonic radiative $B\to\gamma e\nu$ decays_, _Nucl. Phys._ B649 (2003) 349–364, [hep-ph/0210091].
  * (8) T. Becher and M. Neubert, _Toward a NNLO calculation of the $\bar{B}\to X_{s}\gamma$ decay rate with a cut on photon energy. II. Two-loop result for the jet function_, _Phys. Lett._ B637 (2006) 251–259, [hep-ph/0603140].
  * (9) R. Brüser, Z. L. Liu and M. Stahlhofen, _Three-Loop Quark Jet Function_ , _Phys. Rev. Lett._ 121 (2018) 072003, [1804.09722].
  * (10) A. K. Leibovich, Z. Ligeti and M. B. Wise, _Comment on quark masses in SCET_ , _Phys. Lett._ B564 (2003) 231–234, [hep-ph/0303099].
  * (11) I. Z. Rothstein, _Factorization, power corrections, and the pion form-factor_ , _Phys. Rev._ D70 (2004) 054024, [hep-ph/0301240].
  * (12) S. Fleming, A. H. Hoang, S. Mantry and I. W. Stewart, _Jets from massive unstable particles: Top-mass determination_ , _Phys. Rev._ D77 (2008) 074010, [hep-ph/0703207].
  * (13) B. Dehnadi, _Heavy quark mass determinations with sum rules and jets_.  PhD thesis, Vienna U., 2016.
  * (14) S. Fleming, A. H. Hoang, S. Mantry and I. W. Stewart, _Top Jets in the Peak Region: Factorization Analysis with NLL Resummation_ , _Phys. Rev._ D77 (2008) 114003, [0711.2079].
  * (15) B. Dehnadi, A. H. Hoang, V. Mateu, M. Preisser and I. W. Stewart _, in preparation_ .
  * (16) A. H. Hoang, S. Mantry, A. Pathak and I. W. Stewart, _Extracting a Short Distance Top Mass with Light Grooming_ , 1708.02586.
  * (17) Y. Makris and V. Vaidya, _Transverse Momentum Spectra at Threshold for Groomed Heavy Quark Jets_ , _JHEP_ 10 (2018) 019, [1807.09805].
  * (18) C. Lee, P. Shrivastava and V. Vaidya, _Predictions for energy correlators probing substructure of groomed heavy quark jets_ , 1901.09095.
  * (19) M. Butenschoen, B. Dehnadi, A. H. Hoang, V. Mateu, M. Preisser and I. W. Stewart, _Top Quark Mass Calibration for Monte Carlo Event Generators_ , _Phys. Rev. Lett._ 117 (2016) 232001, [1608.01318].
  * (20) A. H. Hoang, S. Plätzer and D. Samitz, _On the Cutoff Dependence of the Quark Mass Parameter in Angular Ordered Parton Showers_ , _JHEP_ 10 (2018) 200, [1807.06617].
  * (21) A. Jain, I. Scimemi and I. W. Stewart, _Two-loop Jet-Function and Jet-Mass for Top Quarks_ , _Phys. Rev._ D77 (2008) 094008, [0801.0743].
  * (22) S. Gritschacher, A. H. Hoang, I. Jemos and P. Pietrulewicz, _Secondary Heavy Quark Production in Jets through Mass Modes_ , _Phys. Rev._ D88 (2013) 034021, [1302.4743].
  * (23) J.-Y. Chiu, A. Jain, D. Neill and I. Z. Rothstein, _A Formalism for the Systematic Treatment of Rapidity Logarithms in Quantum Field Theory_ , _JHEP_ 05 (2012) 084, [1202.0814].
  * (24) A. V. Manohar and I. W. Stewart, _The Zero-Bin and Mode Factorization in Quantum Field Theory_ , _Phys. Rev._ D76 (2007) 074002, [hep-ph/0605001].
  * (25) J.-y. Chiu, A. Fuhrer, A. H. Hoang, R. Kelley and A. V. Manohar, _Soft-Collinear Factorization and Zero-Bin Subtractions_ , _Phys. Rev._ D79 (2009) 053007, [0901.1332].
  * (26) P. Pietrulewicz, S. Gritschacher, A. H. Hoang, I. Jemos and V. Mateu, _Variable Flavor Number Scheme for Final State Jets in Thrust_ , _Phys. Rev._ D90 (2014) 114001, [1405.4860].
  * (27) C. W. Bauer, F. J. Tackmann, J. R. Walsh and S. Zuberi, _Factorization and Resummation for Dijet Invariant Mass Spectra_ , _Phys. Rev._ D85 (2012) 074006, [1106.6047].
  * (28) M. Procura, W. J. Waalewijn and L. Zeune, _Resummation of Double-Differential Cross Sections and Fully-Unintegrated Parton Distribution Functions_ , _JHEP_ 02 (2015) 117, [1410.6483].
  * (29) C. Lee and G. F. Sterman, _Momentum Flow Correlations from Event Shapes: Factorized Soft Gluons and Soft-Collinear Effective Theory_ , _Phys. Rev._ D75 (2007) 014022, [hep-ph/0611061].
  * (30) A. Idilbi and T. Mehen, _On the equivalence of soft and zero-bin subtractions_ , _Phys. Rev._ D75 (2007) 114017, [hep-ph/0702022].
  * (31) A. H. Hoang, P. Pietrulewicz and D. Samitz, _Variable Flavor Number Scheme for Final State Jets in DIS_ , _Phys. Rev._ D93 (2016) 034034, [1508.04323].
  * (32) P. Pietrulewicz, D. Samitz, A. Spiering and F. J. Tackmann, _Factorization and Resummation for Massive Quark Effects in Exclusive Drell-Yan_ , _JHEP_ 08 (2017) 114, [1703.09702].
  * (33) D. Samitz, _in preparation_.  PhD thesis, Vienna U.
  * (34) S. Gritschacher, A. Hoang, I. Jemos and P. Pietrulewicz, _Two loop soft function for secondary massive quarks_ , _Phys. Rev._ D89 (2014) 014035, [1309.6251].
  * (35) T. Becher and M. D. Schwartz, _Direct photon production with effective field theory_ , _JHEP_ 02 (2010) 040, [0911.0681].
  * (36) T. Becher and G. Bell, _The gluon jet function at two-loop order_ , _Phys. Lett._ B695 (2011) 252–258, [1008.1936].
  * (37) P. Banerjee, P. K. Dhani and V. Ravindran, _Gluon jet function at three loops in QCD_ , _Phys. Rev._ D98 (2018) 094016, [1805.02637].
  * (38) I. W. Stewart, F. J. Tackmann and W. J. Waalewijn, _Factorization at the LHC: From PDFs to Initial State Jets_ , _Phys. Rev._ D81 (2010) 094035, [0910.0467].
  * (39) J.-y. Chiu, A. Jain, D. Neill and I. Z. Rothstein, _The Rapidity Renormalization Group_ , _Phys. Rev. Lett._ 108 (2012) 151601, [1104.0881].
  * (40) J. Ablinger, J. Blümlein, A. De Freitas, M. van Hoeij, E. Imamoglu, C. G. Raab et al., _Iterated Elliptic and Hypergeometric Integrals for Feynman Diagrams_ , _J. Math. Phys._ 59 (2018) 062305, [1706.01299].
  * (41) A. H. Hoang, A. Pathak, P. Pietrulewicz and I. W. Stewart, _Hard Matching for Boosted Tops at Two Loops_ , _JHEP_ 12 (2015) 059, [1508.04137].
  * (42) T. Becher and G. Bell, _Analytic Regularization in Soft-Collinear Effective Theory_ , _Phys. Lett._ B713 (2012) 41–46, [1112.3907].
  * (43) S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner and J. Schlenk, _Numerical evaluation of two-loop integrals with pySecDec_ , _Acta Phys. Polon. Supp._ 11 (2018) 375, [1712.05755].
  * (44) K. Melnikov and T. van Ritbergen, _The Three loop on-shell renormalization of QCD and QED_ , _Nucl. Phys._ B591 (2000) 515–546, [hep-ph/0005131].
  * (45) A. V. Smirnov, _FIRE5: a C++ implementation of Feynman Integral REduction_ , _Comput. Phys. Commun._ 189 (2015) 182–191, [1408.2372].
  * (46) T. Gehrmann and E. Remiddi, _Two loop master integrals for $\gamma^{*}\to 3\text{ jets}$: The Planar topologies_, _Nucl. Phys._ B601 (2001) 248–286, [hep-ph/0008287].
  * (47) E. Remiddi and J. A. M. Vermaseren, _Harmonic polylogarithms_ , _Int. J. Mod. Phys._ A15 (2000) 725–754, [hep-ph/9905237].
  * (48) T. Huber and D. Maitre, _HypExp: A Mathematica package for expanding hypergeometric functions around integer-valued parameters_ , _Comput. Phys. Commun._ 175 (2006) 122–144, [hep-ph/0507094].
  * (49) D. Maitre, _HPL, a mathematica implementation of the harmonic polylogarithms_ , _Comput. Phys. Commun._ 174 (2006) 222–240, [hep-ph/0507152].
  * (50) A. V. Kotikov, _Differential equations method: New technique for massive Feynman diagrams calculation_ , _Phys. Lett._ B254 (1991) 158–164.
  * (51) E. Remiddi, _Differential equations for Feynman graph amplitudes_ , _Nuovo Cim._ A110 (1997) 1435–1452, [hep-th/9711188].
  * (52) T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys._ B580 (2000) 485–518, [hep-ph/9912329].
  * (53) A. V. Smirnov, _FIESTA4: Optimized Feynman integral calculations with GPU support_ , _Comput. Phys. Commun._ 204 (2016) 189–199, [1511.03614].
  * (54) S. Laporta and E. Remiddi, _Analytic treatment of the two loop equal mass sunrise graph_ , _Nucl. Phys._ B704 (2005) 349–386, [hep-ph/0406160].
  * (55) S. Bauberger, F. A. Berends, M. Bohm and M. Buza, _Analytical and numerical methods for massive two loop selfenergy diagrams_ , _Nucl. Phys._ B434 (1995) 383–407, [hep-ph/9409388].
  * (56) E. Remiddi and L. Tancredi, _Differential equations and dispersion relations for Feynman amplitudes. The two-loop massive sunrise and the kite integral_ , _Nucl. Phys._ B907 (2016) 400–444, [1602.01481].
  * (57) E. C. G. Stueckelberg, _Theory of the radiation of photons of small arbitrary mass_ , _Helv. Phys. Acta_ 30 (1957) 209–215.
  * (58) G. Feldman and P. T. Matthews, _Massive Electrodynamics_ , _Phys. Rev._ 130 (1963) 1633–1638.