[7]G^ #1,#2_#3,#4(#5 #6— #7) aainstitutetext: Theory Group, Deutsches
Elektronen-Synchrotron (DESY), D-22607 Hamburg, Germany bbinstitutetext:
University of Vienna, Faculty of Physics, Boltzmanngasse 5, A-1090 Wien,
Austria ccinstitutetext: Institut für Physik, Humboldt-Universität zu Berlin,
Zum Großen Windkanal 6, D-12489 Berlin, Germany

# Factorization and Resummation for Massive Quark Effects in Exclusive Drell-
Yan

Piotr Pietrulewicz b Daniel Samitz c Anne Spiering a and Frank J. Tackmann
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>berlin.de<EMAIL_ADDRESS>
###### Abstract

Exclusive differential spectra in color-singlet processes at hadron colliders
are benchmark observables that have been studied to high precision in theory
and experiment. We present an effective-theory framework utilizing soft-
collinear effective theory to incorporate massive (bottom) quark effects into
resummed differential distributions, accounting for both heavy-quark initiated
primary contributions to the hard scattering process as well as secondary
effects from gluons splitting into heavy-quark pairs. To be specific, we focus
on the Drell-Yan process and consider the vector-boson transverse momentum,
$q_{T}$, and beam thrust, ${\mathcal{T}}$, as examples of exclusive
observables. The theoretical description depends on the hierarchy between the
hard, mass, and the $q_{T}$ (or ${\mathcal{T}}$) scales, ranging from the
decoupling limit $q_{T}\ll m$ to the massless limit $m\ll q_{T}$. The
phenomenologically relevant intermediate regime $m\sim q_{T}$ requires in
particular quark-mass dependent beam and soft functions. We calculate all
ingredients for the description of primary and secondary mass effects required
at NNLL′ resummation order (combining NNLL evolution with NNLO boundary
conditions) for $q_{T}$ and ${\mathcal{T}}$ in all relevant hierarchies. For
the $q_{T}$ distribution the rapidity divergences are different from the
massless case and we discuss features of the resulting rapidity evolution. Our
results will allow for a detailed investigation of quark-mass effects in the
ratio of $W$ and $Z$ boson spectra at small $q_{T}$, which is important for
the precision measurement of the $W$-boson mass at the LHC.

††preprint:  DESY 16-199 UWTHPH 2017-3 HU-EP-16/42

##  1 Introduction

Differential cross sections for the production of color-singlet states (e.g.
electroweak vector bosons or the Higgs boson) represent benchmark observables
at the LHC. For the Drell-Yan process, the measurements of the transverse
momentum ($q_{T}$) spectrum of the vector boson (and related variables) have
reached uncertainties below the percent level Aad:2011fp ; Chatrchyan:2011wt ;
Aad:2014xaa ; Khachatryan:2015oaa ; Aad:2015auj ; Khachatryan:2016nbe ,
allowing for stringent tests of theoretical predictions from both analytic
resummed calculations and parton-shower Monte-Carlo programs. An accurate
description of the $q_{T}$ spectrum is also a key ingredient for a precise
measurement of the $W$-boson mass at the LHC, which requires a thorough
understanding of the $W$-boson and $Z$-boson spectra and in particular their
ratio Buge:2006dv ; Besson:2008zs ; CMS:2016nnd ; Atlas:2014 . The associated
uncertainties are one of the dominant theoretical uncertainties in the recent
$m_{W}$ determination by the ATLAS collaboration Aaboud:2017svj .

So far, mass effects from charm and bottom quarks in the initial state have
been discussed extensively for inclusive heavy-quark induced cross sections,
leading to the development of several variable-flavor number schemes in deep
inelastic scattering and $pp$ collisions (see e.g. refs. Aivazis:1993kh ;
Aivazis:1993pi ; Thorne:1997ga ; Kretzer:1998ju ; Collins:1998rz ;
Forte:2010ta ; Bonvini:2015pxa ). On the other hand, analogous heavy-quark
mass effects from initial-state radiation have received little attention so
far in the context of resummed exclusive (differential) cross sections, i.e.
where the measurement of an additional (differential) observable restricts the
QCD radiation into the soft-collinear regime requiring the resummation of the
associated logarithms. While e.g. for $m\ll q_{T}$ the mass effects in the
resummed $q_{T}$ distribution are simply encoded by the matching between the
parton distribution functions across a flavor threshold (e.g. matching four-
flavor PDFs onto five-flavor PDFs including a $b$-quark PDF at the scale
$m_{b}$, which happens much below the scale $q_{T}$), this description breaks
down for $q_{T}\sim m$ or $q_{T}\ll m$. A comprehensive treatment of these
regimes in resummed predictions has been missing so far. This concerns in
particular also parton-shower Monte-Carlo generators, which include massive
quark effects primarily as kinematic effects and by using massive splitting
functions. Since heavy-quark initiated corrections are one of the main
differences between the $W$ and $Z$ boson spectra, this issue can play
therefore an important role for $m_{W}$ measurements at the LHC.

In general, one can distinguish two types of mass effects as illustrated in
fig. 1, which have different characteristics: Contributions where the heavy-
quark enters the hard interaction process are called primary mass effects.
Contributions from a gluon splitting into a massive quark-antiquark pair with
light quarks entering the hard interaction are called secondary. For the
$q_{T}$ spectrum, earlier treatments of the heavy-quark initiated primary
contributions for $m\lesssim q_{T}$ have been given in refs. Nadolsky:2002jr ;
Berge:2005rv ; Belyaev:2005bs , essentially combining the ACOT scheme with the
standard CSS $q_{T}$ resummation. A complete setup also requires to account
for secondary mass effects. Their systematic description for differential
spectra in the various relevant hierarchies between mass and other physical
scales has been established in the context of event shapes in $e^{+}e^{-}$
collisions Gritschacher:2013pha ; Pietrulewicz:2014qza  and for threshold
resummation in DIS Hoang:2015iva , see also refs. Butenschoen:2016lpz ;
Dehnadi:2016snl  for a recent utilization in the context of boosted heavy
quark initiated jets. The application to differential spectra in $pp$
collisions will be part of the present paper.

We present a systematic effective-theory treatment of quark mass effects
including both types of mass effects and all possible scale hierarchies using
soft-collinear effective theory (SCET) Bauer:2000ew ; Bauer:2000yr ;
Bauer:2001ct ; Bauer:2001yt . We focus on the Drell-Yan process, $pp\to
Z/\gamma^{*}\to\ell^{+}\ell^{-}$, and consider two types of observables that
resolve additional QCD radiation and are used to constrain the process to the
exclusive region, namely the transverse momentum $q_{T}$ of the gauge boson
and beam thrust Stewart:2009yx ,

$\displaystyle
q_{T}=|\vec{q}_{T}|=\lvert\vec{p}_{T\ell}+\vec{p}_{T\bar{\ell}}\rvert=\Bigl{\lvert}\sum_{i}\vec{p}_{Ti}\Bigr{\rvert}\,,\qquad{\mathcal{T}}=\sum_{i}\min\\{n_{a}\\!\cdot\\!p_{i},n_{b}\\!\cdot\\!p_{i}\\}\,.$
(1)

Here, $p_{i}$ are all hadronic final-state momenta (i.e. excluding the color-
singlet final state), and $n_{a,b}^{\mu}=(1,\pm\hat{z})$ are lightlike vectors
along the beam axes. Due to transverse momentum conservation $q_{T}$ measures
the total transverse momentum of the final state hadronic radiation, while
beam thrust measures the momentum projections of all hadronic particles onto
the beam axis. The exclusive regime we are interested in corresponds to
$q_{T}\ll Q$ or ${\mathcal{T}}\ll Q$, where $Q=\sqrt{q^{2}}$ is the dilepton
invariant mass. These two observables restrict the allowed QCD radiation into
the collinear and soft regime in different ways, leading to different
effective-theory setups with distinct factorization and resummation
properties, which are well-known in the massless limit up to high orders in
the logarithmic counting (see e.g. refs. Collins:1981uk ; Collins:1981va ;
Collins:1984kg ; Balazs:1997xd ; Balazs:2000wv ; Bozzi:2010xn ;
deFlorian:2011xf ; Becher:2010tm ; Chiu:2011qc ; Chiu:2012ir ;
GarciaEchevarria:2011rb  and refs. Stewart:2009yx ; Stewart:2010pd ;
Berger:2010xi ). These two cases provide simple prototypical examples, which
cover the essential features of the factorization with massive quarks that
will also be relevant for including massive quark effects for other more
complicated jet resolution variables whose factorization is known in the
massless limit. Throughout the paper we always consider the limit
$\Lambda_{\mathrm{QCD}}\ll q_{T},{\mathcal{T}}$ allowing for a perturbative
description of the physics at these kinematic scales. We then consider all
relevant relative hierarchies between the heavy-quark mass $m$ and the
kinematic scales set by the measurement of $q_{T}$ or ${\mathcal{T}}$,
respectively.

Figure 1: Primary (a) and secondary (b) heavy-quark mass effects for $Z$-boson
production.

In the second part of the paper, we explicitly compute all required
ingredients for incorporating $m_{b}$ effects at NNLL′ order, which combines
NNLL evolution with the full NNLO singular boundary conditions (hard, beam,
and soft functions). For $Z$-boson production at NNLL′, primary effects
contribute via $\mathcal{O}(\alpha_{s})\times\mathcal{O}(\alpha_{s})$ heavy-
quark initiated contributions, illustrated in fig. 1. Secondary effects
contribute as $\mathcal{O}(\alpha_{s}^{2})$ corrections to light-quark
initiated hard interactions, illustrated in fig. 1. Due to the strong CKM
suppression primary $m_{b}$-effects do not play any significant role for
$W$-production, which represents a key difference to $Z$-boson production.
Primary $m_{c}$-effects enter $W$-production in the (sizeable) $cs$-channel,
where they start already at NLL′ via
$\mathcal{O}(\alpha_{s})\times\mathcal{O}(1)$ corrections. For this case, our
explicit results for the regime $q_{T}\sim m_{c}$ allows for up to NNLL
resummation. (Here, the resummation at NNLL′ would require the
$\mathcal{O}(\alpha_{s}^{2})$ primary massive contributions.)

The paper is organized as follows: We first discuss in detail the effective
field theory setup for the different parametric regimes for the case of
$q_{T}$ in sec. 2 and for ${\mathcal{T}}$ in sec. 3. Here, we elaborate on the
relevant mode setup in SCET, the resulting factorization formulae, and all-
order relations between the factorization ingredients in the different
regimes. In sec. 4, we give the $\mathcal{O}(\alpha_{s})$ and
$\mathcal{O}(\alpha_{s}^{2})$ results for the various ingredients for NNLL′
resummation. We also verify the consistency of our results with the associated
results in the massless limit. Further details on all calculations are given
in the appendices, where we also give the analytic results at fixed-order for
the massive quark effects in the $q_{T}$ and ${\mathcal{T}}$ distributions in
the singular limit $q_{T},{\mathcal{T}}\ll Q$. In sec. 5, we discuss the
consequences of the secondary mass effects on the rapidity evolution, in
particular for the $q_{T}$ distribution in the regime $q_{T}\sim m_{b}$. As an
outlook we provide in sec. 6 an estimate of the potential size of the bottom
quark effects for low-$q_{T}$ Drell-Yan measurements. In sec. 7 we conclude.

##  2 Factorization of quark mass effects for the $q_{T}$ spectrum

###  2.1 Factorization for massless quarks

Before discussing the massive quark corrections, we first briefly summarize
the EFT setup and factorization for massless quarks. The relevant modes for
the measurement of $q_{T}$ in the limit $q_{T}\ll Q$ are $n_{a}$-collinear,
$n_{b}$-collinear, and soft modes with the scaling

$\displaystyle n_{a}\text{-collinear:}$ $\displaystyle\quad
p_{n_{a}}^{\mu}\sim\Bigl{(}\frac{q_{T}^{2}}{Q},Q,q_{T}\Bigr{)}\,,$
$\displaystyle n_{b}\text{-collinear:}$ $\displaystyle\quad
p_{n_{b}}^{\mu}\sim\Bigl{(}Q,\frac{q_{T}^{2}}{Q},q_{T}\Bigr{)}\,,$ soft:
$\displaystyle\quad p_{s}^{\mu}\sim(q_{T},q_{T},q_{T})\,,$ (2)

which we have written in terms of light-cone coordinates along the beam axis,

$\displaystyle
p^{\mu}=n_{a}\\!\cdot\\!p\,\frac{n_{b}^{\mu}}{2}+n_{b}\\!\cdot\\!p\,\frac{n_{a}^{\mu}}{2}+p_{\perp}^{\mu}\equiv(n_{a}\\!\cdot\\!p,n_{b}\\!\cdot\\!p,p_{\perp})\equiv(p^{+},p^{-},p_{\perp})\,,$
(3)

with $\bar{n}_{a}\equiv n_{b}$. Besides these perturbative modes there are
also nonperturbative collinear modes with the scaling
$(\Lambda_{\mathrm{QCD}}^{2}/Q,Q,\Lambda_{\mathrm{QCD}})$ and
$(Q,\Lambda_{\mathrm{QCD}}^{2}/Q,\Lambda_{\mathrm{QCD}})$, which describe the
initial-state protons at the scale $\mu\sim\Lambda_{\mathrm{QCD}}$, and which
are unrelated to the specific jet resolution measurement. The typical
invariant mass of the soft modes is parametrically the same as for the
collinear modes, $p_{n_{a}}^{2}\sim p_{n_{b}}^{2}\sim p_{s}^{2}\sim
q_{T}^{2}$, which is the characteristic feature of a ${\rm SCET}_{\rm II}$
theory. The soft and collinear modes are only separated in rapidity leading to
the emergence of rapidity divergences and associated rapidity logarithms. The
traditional approach for their resummation in QCD relies on the work by
Collins, Soper, and Sterman Collins:1981uk ; Collins:1981va ; Collins:1984kg .
In SCET the factorization and resummation were devised in refs. Becher:2010tm
; Chiu:2011qc ; Chiu:2012ir ; GarciaEchevarria:2011rb .

Here we will use the rapidity renormalization approach of refs. Chiu:2011qc ;
Chiu:2012ir , where the rapidity divergences are regularized by a symmetric
regulator and are renormalized by appropriate counterterms (by a
$\overline{\rm MS}$-type subtraction). The rapidity logarithms are then
resummed by solving the associated rapidity renormalization group equations.
Within this framework the factorized differential cross section with $n_{f}$
massless quarks reads111In principle there is also a corresponding
contribution for a gluon initiated hard interaction. However, taking into
account the decay of the electroweak boson into massless leptons this
correction vanishes for onshell gluons and only contributes to the power
suppressed terms of $\mathcal{O}(q_{T}/Q)$.

$\displaystyle\frac{\textrm{d}\sigma}{\textrm{d}q^{2}_{T}\,\textrm{d}Q^{2}\,\textrm{d}Y}$
$\displaystyle=\sum_{i,j\in\\{q,\bar{q}\\}}H^{(n_{f})}_{ij}(Q,\mu)\int\\!\textrm{d}^{2}p_{Ta}\,\textrm{d}^{2}p_{Tb}\,\textrm{d}^{2}p_{Ts}\,\delta(q_{T}^{2}-|\vec{p}_{Ta}+\vec{p}_{Tb}+\vec{p}_{Ts}|^{2})$
(4) $\displaystyle\quad\times
B^{(n_{f})}_{i}\Bigl{(}\vec{p}_{Ta},x_{a},\mu,\frac{\nu}{\omega_{a}}\Bigr{)}\,B^{(n_{f})}_{j}\Bigl{(}\vec{p}_{Tb},x_{b},\mu,\frac{\nu}{\omega_{b}}\Bigr{)}\,S^{(n_{f})}(\vec{p}_{Ts},\mu,\nu)\,\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}}{Q}\Bigr{)}\Bigr{]}\,,$

where

$\displaystyle\omega_{a}=Qe^{Y}\,,\quad\omega_{b}=Qe^{-Y}\,,\quad
x_{a,b}=\frac{\omega_{a,b}}{E_{\mathrm{cm}}}\,,$ (5)

with $Y$ denoting the rapidity of the color-singlet state.

In eq. (4), the superscript $(n_{f})$ on all functions indicates that the
associated EFT operators and the strong coupling constant in these functions
are renormalized with $n_{f}$ active quark flavors, which matters for the
evolution already at LL. $H_{ij}$ denotes the process-dependent (but
measurement-independent) hard function. It encodes the tree-level result and
hard virtual corrections of the partonic process $ij\to Z/W/\gamma^{*}$ at the
scale $\mu\sim Q$. Following refs. Collins:1981uw ; Fleming:2006cd ;
Stewart:2009yx , the renormalized transverse-momentum dependent (TMD) beam
functions $B_{i}$, which are essentially equivalent to TMDPDFs, can be matched
onto PDFs as

$\displaystyle
B^{(n_{f})}_{i}\Bigl{(}\vec{p}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=\sum_{k}\int_{x}^{1}\frac{\textrm{d}z}{z}\,\mathcal{I}^{(n_{f})}_{ik}\Bigl{(}\vec{p}_{T},z,\mu,\frac{\nu}{\omega}\Bigr{)}\,f^{(n_{f})}_{k}\Bigl{(}\frac{x}{z},\mu\Bigr{)}\biggl{[}1+\mathcal{O}\biggl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{|\vec{p}_{T}|^{2}}\biggr{)}\biggr{]}$
$\displaystyle\equiv\sum_{k}\mathcal{I}^{(n_{f})}_{ik}(\vec{p}_{T},x,\mu,\frac{\nu}{\omega})\otimes_{x}f^{(n_{f})}_{k}(x,\mu)\,,$
(6)

where the perturbative matching coefficients $\mathcal{I}_{ik}$ describe the
collinear initial-state radiation at the invariant mass scale $\mu\sim q_{T}$
and rapidity scale $\nu\sim Q$, and the nonperturbative parton distribution
functions (PDFs) are denoted by $f_{k}$. In the following, we abbreviate the
Mellin-type convolution in $x$ as in the second line above. Finally, the soft
function $S$ describes the wide-angle soft radiation at the invariant mass and
rapidity scale $\mu\sim\nu\sim q_{T}$. The matching coefficients
$\mathcal{I}_{ik}$ and the soft function are process-independent and have been
computed to $\mathcal{O}(\alpha_{s}^{2})$ in refs. Catani:2012qa ;
Gehrmann:2012ze ; Gehrmann:2014yya ; Luebbert:2016itl  allowing for a full
NNLL′ analysis of Drell-Yan for massless quarks. The three-loop noncusp
rapidity anomalous dimension required for the resummation at N3LL has recently
become available Li:2016axz ; Li:2016ctv ; Vladimirov:2016dll .

In eq. (4), the logarithms of $q_{T}/Q$ are resummed by evaluating all
functions at their characteristic renormalization scales and evolving them to
common final scales $\mu$ and $\nu$ by solving the set of coupled evolution
equations

$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}H^{(n_{f})}_{ij}(Q,\mu)$
$\displaystyle=\gamma_{H}^{(n_{f})}(Q,\mu)\,H^{(n_{f})}_{ij}(Q,\mu)\,,$
$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}B^{(n_{f})}_{i}\Bigl{(}\vec{p}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=\gamma_{B}^{(n_{f})}\Bigl{(}\mu,\frac{\nu}{\omega}\Bigr{)}\,B^{(n_{f})}_{i}\Bigl{(}\vec{p}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}\,,$
$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}S^{(n_{f})}(\vec{p}_{T},\mu,\nu)$
$\displaystyle=\gamma_{S}^{(n_{f})}(\mu,\nu)\,S^{(n_{f})}(\vec{p}_{T},\mu,\nu)\,,$
$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}f_{i}^{(n_{f})}(x,\mu)$
$\displaystyle=\sum_{k}\gamma_{f,ik}^{(n_{f})}(x,\mu)\otimes_{x}f_{k}^{(n_{f})}(x,\mu)\,,$
$\displaystyle\nu\frac{\textrm{d}}{\textrm{d}\nu}B^{(n_{f})}_{i}\Bigl{(}\vec{p}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=\int\\!\textrm{d}^{2}k_{T}\,\gamma_{\nu,B}^{(n_{f})}(\vec{p}_{T}-\vec{k}_{T},\mu)\,B^{(n_{f})}_{i}\Bigl{(}\vec{k}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}\,,$
$\displaystyle\nu\frac{\textrm{d}}{\textrm{d}\nu}S^{(n_{f})}(\vec{p}_{T},\mu,\nu)$
$\displaystyle=\int\\!\textrm{d}^{2}k_{T}\,\gamma_{\nu,S}^{(n_{f})}(\vec{p}_{T}-\vec{k}_{T},\mu)\,S^{(n_{f})}(\vec{k}_{T},\mu,\nu)\,.$
(7)

Only the evolution of the PDF leads to flavor mixing. Consistency of RG
running implies that

$\displaystyle\gamma_{H}^{(n_{f})}(Q,\mu)+\gamma_{B}^{(n_{f})}\Bigl{(}\mu,\frac{\nu}{\omega_{a}}\Bigr{)}+\gamma_{B}^{(n_{f})}\Bigl{(}\mu,\frac{\nu}{\omega_{b}}\Bigr{)}+\gamma_{S}^{(n_{f})}(\mu,\nu)$
$\displaystyle=0\,,$ $\displaystyle
2\gamma_{\nu,B}^{(n_{f})}(\vec{p}_{T},\mu)+\gamma_{\nu,S}^{(n_{f})}(\vec{p}_{T},\mu)$
$\displaystyle=0\,,$
$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}\gamma_{\nu,S}^{(n_{f})}(\vec{p}_{T},\mu)=\nu\frac{\textrm{d}}{\textrm{d}\nu}\gamma_{S}^{(n_{f})}(\mu,\nu)\,\delta(\vec{p}_{T})$
$\displaystyle=-4\Gamma^{(n_{f})}_{\mathrm{cusp}}[\alpha_{s}(\mu)]\,\delta(\vec{p}_{T})\,.$
(8)

Note that in practice, the evolution is usually performed in Fourier space,
such that one actually resums the conjugate logarithms $\ln(b\mu)$ where
$b=\lvert\vec{b}_{T}\rvert\sim 1/q_{T}$ is the Fourier-conjugate variable to
$q_{T}$. The $q_{T}$ spectrum is then obtained as the inverse Fourier
transform of the resummed $b$-spectrum. The exact solution and evolution
directly in $q_{T}$ space, which directly resums the (distributional)
logarithms in $q_{T}$, has been recently discussed in Ebert:2016gcn  (see also
ref. Monni:2016ktx ), and turns out to be significantly more involved due to
the intrinsic two-dimensional nature of $\vec{q}_{T}$.

In the following subsections, we discuss how the mode and factorization setup
changes when massive quark flavors are involved. These lead to the appearance
of additional modes related to fluctuations around the mass shell as discussed
extensively in refs. Gritschacher:2013pha ; Pietrulewicz:2014qza . For the
different hierarchies between the mass scale $m$ and the scales $Q$ and
$q_{T}$ the relevant modes are illustrated in fig. 2. In the first case,
$q_{T}\ll m\sim Q$, the massive flavor is integrated out at the hard scale,
which leads to the above massless case with $n_{l}$ massless flavors, as
discussed in sec. 2.2. The second case, $q_{T}\ll m\ll Q$, where the quark
mass is larger than the jet resolution variable, is analogous to the
corresponding case for thrust in $e^{+}e^{-}\to$ dijets in refs.
Gritschacher:2013pha ; Pietrulewicz:2014qza  and DIS in the $x\to 1$ limit
Hoang:2015iva . We refer to these papers for details and only summarize
briefly the main features for this regime in sec. 2.3. Our main focus is on
the hierarchies $q_{T}\sim m\ll Q$ and $m\ll q_{T}\ll Q$, which are important
for bottom and charm quark mass effects at the LHC, and which are discussed in
secs. 2.4 and 2.5.

(a) ${q_{T}\ll m\sim Q}$

(b) ${q_{T}\ll m\ll Q}$

(c) ${q_{T}\sim m\ll Q}$

(d) ${m\ll q_{T}\ll Q}$

Figure 2: Effective theory modes for the $q_{T}$ spectrum with massive quarks
for $q_{T}\ll Q$ and $m\gg\Lambda_{\rm QCD}$.

###  2.2 Quark mass effects for $m\sim Q$

If the quark mass represents a large scale $\sim Q$ (which concerns the top
quark at the LHC), this quark flavor does not play a dynamic role in the low-
energy effective theory and is integrated out at the hard scale in the
matching from QCD to SCET. The relevant modes are shown in fig. 2(a). The
massive quark only contributes via mass-dependent contributions to the hard
function. This yields the factorization theorem

$\displaystyle\frac{\textrm{d}\sigma}{\textrm{d}q^{2}_{T}\,\textrm{d}Q^{2}\,\textrm{d}Y}$
$\displaystyle=\sum_{i,j\in\\{q,\bar{q}\\}}H_{ij}(Q,m,\mu)\,\int\textrm{d}^{2}p_{Ta}\,\textrm{d}^{2}p_{Tb}\,\textrm{d}^{2}p_{Ts}\,\delta(q_{T}^{2}-|\vec{p}_{Ta}+\vec{p}_{Tb}+\vec{p}_{Ts}|^{2})$
(9) $\displaystyle\quad\times
B^{(n_{l})}_{i}\Bigl{(}\vec{p}_{Ta},x_{a},\mu,\frac{\nu}{\omega_{a}}\Bigr{)}\,B^{(n_{l})}_{j}\Bigl{(}\vec{p}_{Tb},x_{b},\mu,\frac{\nu}{\omega_{b}}\Bigr{)}\,S^{(n_{l})}(\vec{p}_{Ts},\mu,\nu)\biggl{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}^{2}}{m^{2}},\frac{q_{T}}{Q}\Bigr{)}\biggr{]}\,,$

which is essentially equivalent to the massless case in the previous
subsection with $n_{l}$ massless flavors. The hard function $H_{ij}(Q,m,\mu)$
can be evaluated either in the $(n_{f}=n_{l})$ or the $(n_{f}=n_{l}+1)$ flavor
scheme for $\alpha_{s}$, where $n_{l}$ is the number of light (massless) quark
flavors. The associated massive quark corrections are directly related to the
virtual contributions to the quark form factors, e.g. given at
$\mathcal{O}(\alpha_{s}^{2})$ by the virtual diagrams in fig. 1. In general
both primary and secondary corrections contribute for initial (massless)
quarks. Using the $(n_{l})$ flavor scheme for $\alpha_{s}$ these vanish as
$\mathcal{O}(Q^{2}/m^{2})$ in the decoupling limit $m\gg Q$ for the conserved
vector current. For the axial-vector current, contributing to $Z$-boson
production, there are in addition also anomaly corrections starting at
$\mathcal{O}(\alpha_{s}^{2})$ from the massive quark triangle in fig. 1 that
do not decouple.222Instead, for $m\gg Q$ the heavy quark can be integrated out
around its mass scale and the axial current can be evolved between $m$ and $Q$
to resum the associated logarithms $\ln(m/Q)$. Since the massive quark does
not appear as a dynamic flavor in the EFT below the hard scale $Q$, the entire
RG evolution to sum the logarithms of $q_{T}$ is performed with $n_{l}$
massless flavors as in eq. (4).

###  2.3 Quark mass effects for $q_{T}\ll m\ll Q$

Next, we consider the hierarchies where the quark mass is parametrically
smaller than the hard scale, $m\ll Q$. These require a different factorization
setup than $m\sim Q$ since fluctuations around the mass-shell are now
parametrically separated from hard fluctuations, which would lead to large
unresummed logarithms inside the hard function $H_{ij}(Q,m,\mu)$. In this
subsection, we start with the case where the transverse momentum is much
smaller than the mass, $q_{T}\ll m\ll Q$, while $q_{T}\sim m\ll Q$ and $m\ll
q_{T}\ll Q$ are considered in the following subsections.

In a first step the QCD current is matched onto the SCET current with
$n_{l}+1$ dynamic quark flavors at the scale $\mu\sim Q$. Since $m\ll Q$ this
matching can be performed (at leading order in the expansion parameter $m/Q$)
only with massless quarks, leading to the hard function with $n_{l}+1$
massless flavors, $H_{ij}^{(n_{l}+1)}$, with the strong coupling inside it
renormalized with $n_{l}+1$ flavors. The matching is performed onto SCET
containing $n_{a}$-collinear, $n_{b}$-collinear, and soft mass modes with the
scaling

$\displaystyle n_{a}\text{-collinear MM:}$ $\displaystyle\quad
p_{m,n_{a}}^{\mu}\sim\Bigl{(}\frac{m^{2}}{Q},Q,m\Bigr{)}\,,$ $\displaystyle
n_{b}\text{-collinear MM:}$ $\displaystyle\quad
p_{m,n_{b}}^{\mu}\sim\Bigl{(}Q,\frac{m^{2}}{Q},m\Bigr{)}\,,$ soft MM:
$\displaystyle\quad p_{m,s}^{\mu}\sim(m,m,m)\,,$ (10)

as illustrated in fig. 2(b). These mass-shell fluctuations arise here purely
from secondary virtual contributions.

In a second step at the scale $\mu\sim m$, the mass modes are integrated out
and the SCET with $n_{l}$ massless and one massive flavor is matched onto SCET
with $n_{l}$ massless flavors with the usual scaling as in the massless case
in eq. (2.1). Since the soft and collinear mass modes have the same invariant
mass set by the quark mass and are only separated in rapidity, there are
rapidity divergences in their (unrenormalized) collinear and soft
contributions. Their renormalization and the resummation of the associated
logarithms can be again handled using the rapidity RG approach in refs.
Chiu:2011qc ; Chiu:2012ir , which has been explicitly carried out in ref.
Hoang:2015vua .333The matching in ref. Hoang:2015vua  was performed with
massive primary quarks yielding the matching functions denoted as $H_{m,n}$,
$H_{m,\bar{n}}$ and $H_{m,s}$ there. However, this does not affect the
structure of the rapidity logarithms arising from the secondary mass effects,
which are independent of the primary quarks being massive or massless. In
addition, all renormalized parameters like the strong coupling constant are
matched at the mass scale from $n_{l}+1$ to $n_{l}$ flavors taking into
account that the massive flavor is removed as a dynamic degree of freedom.

After these steps, the factorization at the low scale $\sim q_{T}$ proceeds as
in the massless case with all operator matrix elements depending on the
$n_{l}$ massless flavors, which yields the factorization theorem

$\displaystyle\frac{\textrm{d}\sigma}{\textrm{d}q^{2}_{T}\,\textrm{d}Q^{2}\,\textrm{d}Y}$
$\displaystyle=\sum_{i,j\in\\{q,\bar{q}\\}}H^{(n_{l}+1)}_{ij}(Q,\mu)\,H_{c}\Bigl{(}m,\mu,\frac{\nu}{\omega_{a}}\Bigr{)}H_{c}\Bigl{(}m,\mu,\frac{\nu}{\omega_{b}}\Bigr{)}H_{s}(m,\mu,\nu)$
$\displaystyle\quad\times\int\textrm{d}^{2}p_{Ta}\,\textrm{d}^{2}p_{Tb}\,\textrm{d}^{2}p_{Ts}\,\delta(q_{T}^{2}-|\vec{p}_{Ta}+\vec{p}_{Tb}+\vec{p}_{Ts}|^{2})\,B^{(n_{l})}_{i}\Bigl{(}\vec{p}_{Ta},x_{a},\mu,\frac{\nu}{\omega_{a}}\Bigr{)}$
$\displaystyle\quad\times
B^{(n_{l})}_{j}\Bigl{(}\vec{p}_{Tb},x_{b},\mu,\frac{\nu}{\omega_{b}}\Bigr{)}\,S^{(n_{l})}(\vec{p}_{Ts},\mu,\nu)\,\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}}{Q},\frac{q_{T}^{2}}{m^{2}},\frac{m^{2}}{Q^{2}}\Bigr{)}\Bigr{]}\,.$
(11)

Here $H_{c}$ and $H_{s}$ denote the hard functions that arise from the
matching at the mass scale $\mu\sim m$. Their natural rapidity scales are
$\nu\sim Q$ for the collinear contributions and $\nu\sim m$ for the soft ones.
They can be evaluated in either the $(n_{l})$ or $(n_{l}+1)$ scheme for
$\alpha_{s}$. We will give their expressions at $\mathcal{O}(\alpha_{s}^{2})$
in sec. 4.1. The resummation of all logarithms of ratios of $q_{T}$, $m$, and
$Q$ is achieved by performing the evolution in $\mu$ and $\nu$ of all
functions appearing in eq. (2.3) from their natural scales.

In principle, the $\mu$ evolution can be performed by evolving all functions
with their respective number of quark flavors without switching the flavor
scheme, i.e. with $n_{l}+1$ flavors for $H$, $n_{l}$ flavors for $B$ and $S$
and an additional evolution for the collinear and soft matching functions
$H_{c}$ and $H_{s}$. The consistency of RG running for the factorization
theorems in eqs. (2.3) and (9), and eq. (2.1) with $n_{l}$ massless flavors,
implies that the $\mu$-dependence of the mass-dependent hard functions $H_{c}$
and $H_{s}$ is precisely given by the difference between $n_{l}$ and $n_{l}+1$
active quark flavors in the evolution of the hard function $H_{ij}$,

$\displaystyle\gamma_{H_{c}}\Bigl{(}m,\mu,\frac{\nu}{\omega_{a}}\Bigr{)}+\gamma_{H_{c}}\Bigl{(}m,\mu,\frac{\nu}{\omega_{b}}\Bigr{)}+\gamma_{H_{s}}(m,\mu,\nu)=\gamma^{(n_{l})}_{H}(Q,\mu)-\gamma^{(n_{l}+1)}_{H}(Q,\mu)\,,$
(12)

where $\gamma^{(n_{f})}_{H}$ is defined in eq. (2.1), and $\gamma_{H_{c}}$ and
$\gamma_{H_{s}}$ are defined analogously. At two loops this relation can be
checked explicitly using the results in eqs. (4.1.2), (4.1.2) and (A.1). As a
result, the $\mu$ evolution for the hard functions can be conveniently
implemented as illustrated in fig. 3(a), by carrying out the $\mu$ evolution
with $n_{l}$ active quark flavors below the matching scale $\mu_{m}\sim m$ and
with $n_{l}+1$ flavors above $\mu_{m}$, providing in this sense a “variable-
flavor number scheme” Pietrulewicz:2014qza ; Hoang:2015iva . (This effectively
corresponds to using operator running for the hard scattering current, which
is renormalized with $n_{l}+1$ flavors above the mass scale and with $n_{l}$
flavors below the mass scale.) In addition there is also a rapidity evolution,
which is carried out at $\mu_{m}=m$, i.e. at the border between the
$(n_{l}+1)$ and $(n_{l})$-flavor theories (see ref. Hoang:2015vua ), which is
governed by the mass-dependent rapidity anomalous dimensions for $H_{s}$ and
$H_{c}$,

$\displaystyle\gamma_{\nu,H_{s}}(m,\mu)=-2\gamma_{\nu,H_{c}}(m,\mu)=\frac{\textrm{d}}{\textrm{d}\ln\nu}\ln
H_{s}(m,\mu,\nu)\,.$ (13)

(a) ${q_{T}\ll m\ll Q}$

(b) ${q_{T}\sim m\ll Q}$

(c) ${m\ll q_{T}\ll Q}$

Figure 3: Illustration of the renormalization group evolution for $q_{T}$ of
the hard, beam, soft, and parton distribution functions in invariant mass and
rapidity. The anomalous dimensions for each evolution step involve the
displayed number of active quark flavors. The label $m$ indicates that the
corresponding evolution is mass dependent.

###  2.4 Quark mass effects for $q_{T}\sim m\ll Q$

If the $q_{T}$ scale is of the order of the quark mass, $q_{T}\sim m$, the
massive quark becomes a dynamic degree of freedom, which contributes to the
$q_{T}$ spectrum via real radiation effects. The mass modes in eq. (2.3) are
now the same as the usual massless ${\rm SCET}_{\rm II}$ modes for the $q_{T}$
measurement in eq. (2.1), since their parametrically scaling coincides for
$q_{T}\sim m$, as illustrated in fig. 2(c). In this case, there is only a
single matching at the hard scale $\mu\sim Q$ from QCD onto SCET with these
common soft and collinear modes. This hard matching gives again rise to the
(mass-independent) hard function $H_{ij}^{(n_{l}+1)}$ for $n_{l}+1$ massless
flavors. The SCET operator matrix elements at the scale $\mu\sim q_{T}$, i.e.
the beam and soft functions, now encode the effects of the massive quark. They
are now renormalized with $n_{l}+1$ quark flavors and contain an explicit
dependence on the quark mass. When integrating out the modes with the
virtuality $q_{T}$ also the massive quark is integrated out and the collinear
matching functions $\mathcal{I}_{ik}$ between the beam functions and the PDFs
thus also contain the effect from changing from $n_{l}+1$ to $n_{l}$ flavors,
i.e.

$\displaystyle
B_{i}^{(n_{l}+1)}\Bigl{(}\vec{p}_{T},m,x,\mu,\frac{\nu}{\omega}\Bigr{)}=\sum_{k\in\\{q,\bar{q},g\\}}\mathcal{I}_{ik}\Bigl{(}\vec{p}_{T},m,x,\mu,\frac{\nu}{\omega}\Bigr{)}\otimes_{x}f^{(n_{l})}_{k}(x,\mu)\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{m^{2}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{q_{T}^{2}}\Bigr{)}\Bigr{]}\,.$
(14)

Written out explicitly, the factorization theorem reads

$\displaystyle\frac{\textrm{d}\sigma}{\textrm{d}q^{2}_{T}\,\textrm{d}Q^{2}\,\textrm{d}Y}$
$\displaystyle=\sum_{i,j\in\\{q,\bar{q},Q,\bar{Q}\\}}H^{(n_{l}+1)}_{ij}(Q,\mu)\,\int\textrm{d}^{2}p_{Ta}\,\textrm{d}^{2}p_{Tb}\,\textrm{d}^{2}p_{Ts}\,\delta(q_{T}^{2}-|\vec{p}_{Ta}+\vec{p}_{Tb}+\vec{p}_{Ts}|^{2})$
$\displaystyle\quad\times\biggl{[}\sum_{k\in\\{q,\bar{q},g\\}}\mathcal{I}_{ik}\Bigl{(}\vec{p}_{Ta},m,x_{a},\mu,\frac{\nu}{\omega_{a}}\Bigr{)}\otimes_{x}f^{(n_{l})}_{k}(x_{a},\mu)\biggr{]}$
$\displaystyle\quad\times\biggl{[}\sum_{k\in\\{q,\bar{q},g\\}}\mathcal{I}_{jk}\Bigl{(}\vec{p}_{Tb},m,x_{b},\mu,\frac{\nu}{\omega_{b}}\Bigr{)}\otimes_{x}f^{(n_{l})}_{k}(x_{b},\mu)\biggr{]}$
$\displaystyle\quad\times
S(\vec{p}_{Ts},m,\mu,\nu)\,\biggl{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}}{Q},\frac{m^{2}}{Q^{2}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{m^{2}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{q_{T}^{2}}\Bigr{)}\biggr{]}\,,$
(15)

where $i,j=Q,\bar{Q}$ denotes the massive quark flavor in the sum over
flavors. We stress that the renormalization of the bare soft and beam function
with $n_{l}$ massless and one massive flavor is carried out in the $n_{l}+1$
flavor scheme for $\alpha_{s}$, while the strong coupling in the PDFs (which
are defined in the lower theory with $n_{l}$ massless flavors) is renormalized
with $n_{l}$ flavors. The _renormalized_ soft function and beam function
coefficients $\mathcal{I}_{ik}$ can then be expressed in terms of either the
$(n_{l}+1)$ or the $(n_{l})$ flavor scheme for $\alpha_{s}$ without
introducing large logarithms.

In this hierarchy quark mass effects enter in eq. (2.4) at
$\mathcal{O}(\alpha_{s}^{2})$ in two ways: There are secondary radiation
effects appearing in the two-loop soft function $S^{(2)}$ and the flavor-
diagonal beam function matching coefficients $\mathcal{I}^{(2)}_{qq}$. In
addition, there are primary mass effects arising from a massive-quark
initiated hard process. For $Z/\gamma^{*}$ production, this requires the
production of the massive quarks via gluon splitting in both collinear
sectors, which manifests itself in two one-loop collinear matching
coefficients $\mathcal{I}^{(1)}_{Qg}\times\mathcal{I}^{(1)}_{\bar{Q}g}$. For
$W$-boson production, primary charm quark effects enter already at
$\mathcal{O}(\alpha_{s})$ from a single $\mathcal{I}^{(1)}_{Qg}$ with $Q=c$.

The resummation of logarithms $\ln(q_{T}/Q)$ and $\ln(m/Q)$ is again obtained
by performing the RG evolution for eq. (2.4), which is illustrated in fig.
3(b). While the evolution of the PDFs proceeds in $n_{l}$ flavors, the
$\mu$-evolution for the hard, beam, and soft functions above the scale $m$ is
now carried out purely with $n_{l}+1$ flavors. Consistency of RG running for
eq. (14) implies that the matching coefficients $\mathcal{I}_{ik}$ satisfy

$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}\,\mathcal{I}_{ik}\Bigl{(}\vec{p}_{T},m,z,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=\Bigl{[}\gamma^{(n_{l}+1)}_{B_{i}}\times\mathcal{I}_{ik}\Bigr{]}\Bigl{(}\vec{p}_{T},m,z,\mu,\frac{\nu}{\omega}\Bigr{)}-\sum_{j\in
q,\bar{q},g}\Bigl{[}\mathcal{I}_{ij}\otimes\gamma^{(n_{l})}_{f,jk}\Bigr{]}\Bigl{(}\vec{p}_{T},m,z,\mu,\frac{\nu}{\omega}\Bigr{)}\,.$
(16)

Since the renormalization of the beam functions does not involve parton
mixing, the one-loop primary mass contributions to ${\mathcal{I}}_{Qg}^{(1)}$
cannot give rise to rapidity divergences and associated logarithms. On the
other hand, the secondary mass effects change the rapidity evolution. In
particular, the beam and soft $\nu$-anomalous dimensions become mass
dependent444The fact that quark masses can affect the evolution was already
pointed out in ref. Collins:1984kg .,

$\displaystyle\nu\frac{\textrm{d}}{\textrm{d}\nu}B_{i}^{(n_{l}+1)}\Bigl{(}\vec{p}_{T},m,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=\int\\!\textrm{d}^{2}k_{T}\,\gamma^{(n_{l}+1)}_{\nu,B}(\vec{p}_{T}-\vec{k}_{T},m,\mu)\,B_{i}^{(n_{l}+1)}\Bigl{(}\vec{k}_{T},m,\mu,\frac{\nu}{\omega}\Bigr{)}\,,$
$\displaystyle\nu\frac{\textrm{d}}{\textrm{d}\nu}S^{(n_{l}+1)}(\vec{p}_{T},m,\mu,\nu)$
$\displaystyle=\int\\!\textrm{d}^{2}k_{T}\,\gamma^{(n_{l}+1)}_{\nu,S}(\vec{p}_{T}-\vec{k}_{T},m,\mu)\,S^{(n_{l}+1)}(\vec{k}_{T},m,\mu,\nu)\,.$
(17)

We discuss the implications of the mass dependence for the rapidity evolution
in sec. 5.

###  2.5 Quark mass effects for $m\ll q_{T}\ll Q$

If $q_{T}$ is much larger than the mass, the fluctuations around the mass-
shell take place at a different scale than the jet resolution measurement.
There are no relevant soft fluctuations scaling like
$p_{m,s}^{\mu}\sim(m,m,m)$, since the measurement of $q_{T}$ is IR safe and is
thus insensitive to the lower mass scale. (In other words, if we were to
explicitly distinguish such soft mass modes their contribution would cancel as
for an inclusive observable since they are not constrained by the
$q_{T}$-measurement.) This means that the soft modes are described by a soft
function with $n_{l}+1$ massless flavors at the scale $\mu\sim q_{T}$. Due to
the collinear sensitivity of the initial-state radiation there are still
relevant collinear mass modes scaling like
$p_{m,n_{a}}^{\mu}\sim(m^{2}/Q,Q,m)$ and $p_{m,n_{b}}^{\mu}\sim(Q,m^{2}/Q,m)$,
as illustrated in fig. 2(d). Thus there are collinear modes in SCET at
different invariant mass scales, which can be disentangled by a multistage
matching. First, the beam functions are matched onto the PDFs with $n_{l}$
massless and one massive flavor. Since this matching takes place at the scale
$\mu_{B}\sim q_{T}\gg m$ this gives just the matching coefficients
$\mathcal{I}_{ik}$ for $n_{l}+1$ massless flavors,

$\displaystyle
B_{i}^{(n_{l}+1)}\Bigl{(}\vec{p}_{T},m,x,\mu,\frac{\nu}{\omega}\Bigr{)}=\sum_{k\in\\{q,\bar{q},Q,\bar{Q},g\\}}\mathcal{I}^{(n_{l}+1)}_{ik}\Bigl{(}\vec{p}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}\otimes_{x}f^{(n_{l}+1)}_{k}(x,m,\mu)\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{m^{2}}{q_{T}^{2}}\Bigr{)}\Bigr{]}\,.$
(18)

In a second step, at the mass scale $\mu_{m}\sim m$, the PDFs including the
massive quark effects are matched onto PDFs with $n_{l}$ massless quarks, and
with $\alpha_{s}$ in the $(n_{l})$ flavor scheme,

$\displaystyle
f^{(n_{l}+1)}_{i}(x,m,\mu)=\sum_{k\in\\{q,\bar{q},g\\}}\mathcal{M}_{ik}(x,m,\mu)\otimes_{x}f^{(n_{l})}_{k}(x,\mu)\biggl{[}1+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{m^{2}}\Bigr{)}\biggr{]}\,.$
(19)

The PDF matching functions $\mathcal{M}_{ik}$ can be expressed in either the
$(n_{l})$ or the $(n_{l}+1)$ flavor scheme for $\alpha_{s}$.

In total, the factorization theorem reads

$\displaystyle\frac{\textrm{d}\sigma}{\textrm{d}q^{2}_{T}\,\textrm{d}Q^{2}\,\textrm{d}Y}$
$\displaystyle=\sum_{i,j\in\\{q,\bar{q},Q,\bar{Q}\\}}H^{(n_{l}+1)}_{ij}(Q,\mu)\,\int\textrm{d}^{2}p_{Ta}\,\textrm{d}^{2}p_{Tb}\,\textrm{d}^{2}p_{Ts}\,\delta(q_{T}^{2}-|\vec{p}_{Ta}+\vec{p}_{Tb}+\vec{p}_{Ts}|^{2})$
$\displaystyle\quad\times\biggl{[}\sum_{k\in\\{q,\bar{q},Q,\bar{Q},g\\}}\sum_{l\in\\{q,\bar{q},g\\}}\mathcal{I}^{(n_{l}+1)}_{ik}\Bigl{(}\vec{p}_{Ta},x_{a},\mu,\frac{\nu}{\omega_{a}}\Bigr{)}\otimes_{x}\mathcal{M}_{kl}(x_{a},m,\mu)\otimes_{x}f^{(n_{l})}_{l}(x_{a},\mu)\biggr{]}$
$\displaystyle\quad\times\biggl{[}\sum_{k\in\\{q,\bar{q},Q,\bar{Q},g\\}}\sum_{l\in\\{q,\bar{q},g\\}}\mathcal{I}^{(n_{l}+1)}_{jk}\Bigl{(}\vec{p}_{Tb},x_{b},\mu,\frac{\nu}{\omega_{b}}\Bigr{)}\otimes_{x}\mathcal{M}_{kl}(x_{b},m,\mu)\otimes_{x}f^{(n_{l})}_{l}(x_{b},\mu)\biggr{]}$
$\displaystyle\quad\times
S^{(n_{l}+1)}(\vec{p}_{Ts},\mu,\nu)\biggl{[}1+\mathcal{O}\Bigl{(}\frac{q_{T}}{Q},\frac{m^{2}}{q_{T}^{2}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{m^{2}}\Bigr{)}\biggr{]}\,.$
(20)

As in sec. 2.4, massive quark corrections can arise at
$\mathcal{O}(\alpha_{s}^{2})$ either via primary mass effects involving the
product of two one-loop PDF matching corrections $\mathcal{M}^{(1)}_{Qg}$ (for
$Z/\gamma^{*}$) generating a massive quark-antiquark pair that initiates the
hard interaction, or via secondary mass effects involving one two-loop
contribution $\mathcal{M}^{(2)}_{qq}$. Note that also the running of the light
quark and gluon PDFs above $\mu_{m}$ generates an effective massive quark PDF
via evolution factors $U^{(n_{l}+1)}_{f,Qk}(\mu_{B},\mu_{m})$, which for large
hierarchies $m\ll q_{T}$ can give $\mathcal{O}(1)$ contributions.

The evolution of the hard, beam, and soft functions can be performed purely
with $n_{l}+1$ massless flavors, while the PDFs are evolved with $n_{l}$
flavor below $\mu_{m}\sim m$ and with $n_{l}+1$ flavors above $\mu_{m}\sim m$,
see fig. 3(c), as usual in any variable-flavor number scheme. The $\mu_{m}$
dependence is canceled order-by-order by the matching factors
$\mathcal{M}_{ij}$,

$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}\,\mathcal{M}_{ik}(m,z,\mu)$
$\displaystyle=\sum_{j\in\\{q,\bar{q},Q,\bar{Q},g\\}}\Bigl{[}\gamma^{(n_{l}+1)}_{f,ij}\otimes_{z}\mathcal{M}_{jk}\Bigr{]}(m,z,\mu)-\sum_{j\in\\{q,\bar{q},g\\}}\Bigl{[}\mathcal{M}_{ij}\otimes_{z}\gamma^{(n_{l})}_{f,jk}\Bigr{]}(m,z,\mu)\,.$
(21)

The absence of soft mass modes in this regime implies there is no rapidity
evolution at the mass scale, while the rapidity evolution between beam and
soft functions is the same as for $n_{l}+1$ massless flavors.

In this regime, the mass dependence is thus fully contained in the collinear
sectors. Within each collinear sector, the EFT setup is completely analogous
to that of the heavy-quark induced inclusive cross section discussed in detail
in ref. Bonvini:2015pxa , with the beam functions here playing the role of the
inclusive cross section there and the $q_{T}$ scale here playing the role of
the hard scale there.

###  2.6 Relations between hierarchies

After discussing all hierarchies separately, we now show how the ingredients
in each of the associated factorization theorems are related to each other.
This will also make it obvious how the mass-dependent fixed-order corrections
that are kept in one hierarchy but are dropped in another can be combined with
the resummation of logarithms to obtain a systematic inclusion of the mass
effects over the whole $q_{T}$ spectrum. The relations between the modes and
their contributions between the different regimes are summarized in fig. 4.

Figure 4: Relevant modes for the $q_{T}$ spectrum with $q_{T}\ll Q$ for
different hierarchies between the quark mass $m$ and the scales $q_{T}$ and
$Q$. The arrows indicate the relations between the modes and their associated
contributions.

The hard functions appearing in the hierarchy $q_{T}\ll m\ll Q$ in eq. (2.3)
are related to the hard function for $q_{T}\ll m\sim Q$ in sec. 2.2 as
follows555Here and in the following it is implied that at a specific fixed
order the functions on both sides have to be expanded in the same
renormalization scheme for $\alpha_{s}$. Since $H^{(n_{l}+1)}_{ij}$ is
generically written with $n_{l}+1$ (massless) flavors, the $(n_{l}+1)$-flavor
scheme is most convenient here and also to extract the
$\mathcal{O}(m^{2}/Q^{2})$ power corrections on the right-hand side of eq.
(22) from its fixed-order expansion.

$\displaystyle
H_{ij}(Q,m,\mu)=H^{(n_{l}+1)}_{ij}(Q,\mu)\,H_{c}\Bigl{(}m,\mu,\frac{\nu}{\omega_{a}}\Bigr{)}H_{c}\Bigl{(}m,\mu,\frac{\nu}{\omega_{b}}\Bigr{)}H_{s}(m,\mu,\nu)\biggl{[}1+\mathcal{O}\Bigl{(}\frac{m^{2}}{Q^{2}}\Bigr{)}\biggr{]}\,.$
(22)

In the product of functions on the right-hand side, which appear in eq. (2.3),
the logarithms $\ln(m/Q)$ can be resummed to all orders. One can construct a
smooth description of the cross section for $q_{T}\ll m$ that resums these
logarithms and also includes the associated mass-dependent
$\mathcal{O}(m^{2}/Q^{2})$ power corrections by simply adding the latter to
the hard function $H^{(n_{l}+1)}(Q,\mu)$ at the scale $\mu\sim Q$.

The fixed-order contributions to the operator matrix elements appearing in the
hierarchy $q_{T}\ll m$ are encoded in the ones for $q_{T}\sim m$. The mass-
dependent beam function matching coefficients for $q_{T}\sim m$ are related to
those for $q_{T}\ll m$ and the collinear mass-mode function $H_{c}$ by

$\displaystyle\mathcal{I}_{ik}\Bigl{(}\vec{p}_{T},x,m,\mu,\frac{\nu}{\omega}\Bigr{)}=H_{c}\Bigl{(}m,\mu,\frac{\nu}{\omega}\Bigr{)}\,\mathcal{I}_{ik}^{(n_{l})}\Bigl{(}\vec{p}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}\biggl{[}1+\mathcal{O}\Bigl{(}\frac{p_{T}^{2}}{m^{2}}\Bigr{)}\biggr{]}\,.$
(23)

Similarly, the mass-dependent soft function for $q_{T}\sim m$ is related to
the one for $q_{T}\ll m$ and the soft mass-mode function $H_{s}$ by

$\displaystyle
S(\vec{p}_{T},m,\mu,\nu)=H_{s}(m,\mu,\nu)\,S^{(n_{l})}(\vec{p}_{T},\mu,\nu)\biggl{[}1+\mathcal{O}\Bigl{(}\frac{p_{T}^{2}}{m^{2}}\Bigr{)}\biggr{]}\,.$
(24)

In the products on the right-hand sides, which appear in eq. (2.3), logarithms
$\ln(q_{T}/m)$ are resummed to all orders in the limit $q_{T}\ll m$. One can
include the associated $\mathcal{O}(q_{T}^{2}/m^{2})$ power corrections that
are important for $q_{T}\sim m$, by obtaining them from the fixed-order
expansions of eqs. (23) and (24) and adding them to the $(n_{l})$-flavor beam
function coefficients and soft function at the scale $\mu\sim q_{T}$.

Finally, the fixed-order contributions for the operator matrix elements
appearing in the hierarchy $m\ll q_{T}$ are also encoded in the corresponding
ones for $q_{T}\sim m$. Hence, the mass-dependent beam function matching
coefficients are related to those for $m\ll q_{T}$ and the PDF matching
functions by

$\displaystyle\mathcal{I}_{ik}\Bigl{(}\vec{p}_{T},m,x,\mu,\frac{\nu}{\omega}\Bigr{)}=\sum_{j=q,\bar{q},g}\mathcal{I}_{ij}^{(n_{l}+1)}\Bigl{(}\vec{p}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}\otimes_{x}\mathcal{M}_{jk}(m,x,\mu)\biggl{[}1+\mathcal{O}\Bigl{(}\frac{m^{2}}{p_{T}^{2}}\Bigr{)}\biggr{]}\,.$
(25)

Similarly, the mass-dependent and massless soft function are related by

$\displaystyle
S(\vec{p}_{T},m,\mu,\nu)=S^{(n_{l}+1)}(\vec{p}_{T},\mu,\nu)\biggl{[}1+\mathcal{O}\Bigl{(}\frac{m^{2}}{p_{T}^{2}}\Bigr{)}\biggr{]}\,,$
(26)

since there are no relevant soft IR fluctuations below the mass scale. In the
functions on the right-hand sides, which appear in eq. (2.5), logarithms
$\ln(m/q_{T})$ can be resummed to all orders in the limit $m\ll q_{T}$. This
can be combined with the associated $\mathcal{O}(m^{2}/q_{T}^{2})$ power
corrections relevant for $q_{T}\sim m$, by obtaining them from the fixed-order
expansions of eqs. (25) and (26) and adding them to the $(n_{l}+1)$-flavor
beam function matching coefficients and soft function at the scale $\mu\sim
q_{T}$.

By including the various power corrections, one combines the factorization
theorems in the different hierarchies and obtains a theoretical description
that is valid across the whole $q_{T}$ spectrum and includes the resummation
of logarithms in all relevant limits. This can be considered a variable-flavor
scheme for the resummed $q_{T}$ spectrum. (In addition one should of course
also include the usual $q_{T}/Q$ nonsingular corrections to reproduce the full
fixed-order result for $q_{T}\sim Q$.)

We stress that different specific ways of how to incorporate the various power
corrections are formally equivalent as long as the correct fixed-order
expansion and the correct resummation is reproduced in each limit. Any
differences then amount to resummation effects at power-suppressed level and
are thus beyond the formal (leading-power) resummation accuracy.

A particular scheme (“S-ACOT”) to merge the $m\ll q_{T}$ and $q_{T}\sim m$
regimes was discussed in ref. Nadolsky:2002jr  for the primary massive quark
corrections. In practice, for the numerical study of $b$-quark mass effects at
low $q_{T}\ll m\ll Q$ the off-diagonal evolution factor $U_{f,bg}$ and thus
the effective $b$-quark PDF at the scale $q_{T}$ are still quite small, so
that one may effectively count $f_{b}(\mu_{B})\sim\mathcal{O}(\alpha_{s})$. In
particular, this counting facilitates the seamless combination with the
nonsingular corrections for $m\sim q_{T}$ encoded in the beam function
matching coefficients in eq. (2.4). This was discussed in ref. Bonvini:2015pxa
in the context of the inclusive $b\bar{b}H$ production cross section, and the
analogous discussion applies here as well. In refs. Pietrulewicz:2014qza ;
Hoang:2015iva , the power corrections were included implicitly in the
construction of the variable-flavor number schemes for thrust in $e^{+}e^{-}$
and DIS in the endpoint region by applying different renormalization schemes
for the massive quark contributions to the EFT operators above and below the
mass scale.

##  3 Factorization of mass effects for beam thrust

###  3.1 Factorization for massless quarks

For the measurement of beam thrust with ${\mathcal{T}}\ll Q$ the relevant EFT
modes are $n_{a}$-collinear, $n_{b}$–collinear and usoft modes with the
scaling

$\displaystyle n_{a}\text{-collinear:}$ $\displaystyle\quad
p_{n_{a}}^{\mu}\sim({\mathcal{T}},Q,\sqrt{Q{\mathcal{T}}})\,,$ $\displaystyle
n_{b}\text{-collinear:}$ $\displaystyle\quad
p_{n_{b}}^{\mu}\sim(Q,{\mathcal{T}},\sqrt{Q{\mathcal{T}}})\,,$ usoft:
$\displaystyle\quad
p_{us}^{\mu}\sim({\mathcal{T}},{\mathcal{T}},{\mathcal{T}})\,.$ (27)

The usoft and collinear modes are now separated in invariant mass,
$p_{us}^{2}\sim{\mathcal{T}}^{2}\ll p_{n_{a}}^{2}\sim p_{n_{b}}^{2}\sim
Q{\mathcal{T}}$, which is the characteristic feature of a ${\rm SCET}_{\rm I}$
theory. In this case, there are no rapidity logarithms and the renormalization
and evolution is solely in invariant mass. The resulting factorization formula
reads Stewart:2009yx

$\displaystyle\frac{\textrm{d}\sigma}{\textrm{d}Q^{2}\,\textrm{d}Y\,\textrm{d}{\mathcal{T}}}$
$\displaystyle=\sum_{i,j\in\\{q,\bar{q}\\}}H^{(n_{f})}_{ij}(Q,\mu)\int\textrm{d}t_{a}\,\textrm{d}t_{b}\,B^{(n_{f})}_{i}(t_{a},x_{a},\mu)\,B^{(n_{f})}_{j}(t_{b},x_{b},\mu)$
$\displaystyle\quad\times
S^{(n_{f})}\Bigl{(}{\mathcal{T}}-\frac{t_{a}}{\omega_{a}}-\frac{t_{b}}{\omega_{b}},\mu\Bigr{)}\,\biggl{[}1+\mathcal{O}\Bigl{(}\frac{{\mathcal{T}}}{Q}\Bigr{)}\biggr{]}\,.$
(28)

This as well as the expressions including mass effects in the subsequent
subsections are valid for the primary hard scattering, and do not account for
spectator forward (multiparton) scattering effects, since the Glauber
Lagrangian of ref. Rothstein:2016bsq  has been neglected. (There are also
corrections from perturbative Glauber effects starting at ${\cal
O}(\alpha_{s}^{4})$ Gaunt:2014ska ; Zeng:2015iba , which are well beyond the
order we are interested in, but can be calculated and included using the
Glauber operator framework of ref. Rothstein:2016bsq .) This is sufficient for
our purposes of discussing the mass effects in a prototypical ${\rm SCET}_{\rm
I}$ scenario. Our results are also directly relevant to include massive quark
effects in the Geneva Monte-Carlo program Alioli:2012fc ; Alioli:2015toa ,
which employs ${\mathcal{T}}$ as the jet resolution variable for the primary
interaction and where multiparton effects are included Alioli:2016wqt  via the
combination with Pythia8 and its MPI model Sjostrand:1987su ; Sjostrand:2004pf
; Sjostrand:2004ef .

The hard function $H_{ij}$ in eq. (3.1) is measurement independent and the
same as in eq. (4). The beam and soft functions depend on the measurement and
are different from those in eq. (4). The virtuality-dependent beam functions
$B_{i}$ can be factorized into perturbative matching coefficients
$\mathcal{I}_{ik}$ at the scale $\mu\sim t\sim\sqrt{Q{\mathcal{T}}}$ and the
standard nonperturbative PDFs Stewart:2009yx ; Stewart:2010qs

$\displaystyle
B^{(n_{f})}_{i}(t,x,\mu)=\sum_{k}\mathcal{I}^{(n_{f})}_{ik}(t,x,\mu)\otimes_{x}f^{(n_{f})}_{k}(x,\mu)\biggl{[}1+\mathcal{O}\biggl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{t}\biggr{)}\biggr{]}\,.$
(29)

The matching coefficients $\mathcal{I}_{ik}$ have been calculated to
$\mathcal{O}(\alpha_{s}^{2})$ Gaunt:2014xga ; Gaunt:2014cfa . The soft
function at the scale $\mu\sim{\mathcal{T}}$ is equivalent to the thrust soft
function Kang:2015moa , which is known to $\mathcal{O}(\alpha_{s}^{2})$
Kelley:2011ng ; Monni:2011gb . The noncusp anomalous dimensions required at
N3LL are available from existing results Stewart:2010qs .

The resummation of logarithms $\ln({\mathcal{T}}/Q)$ is performed by
evaluating all functions at their characteristic scales and evolving them to a
common final scale $\mu$ using the solutions of the RGEs

$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}B^{(n_{f})}_{i}(t,x,\mu)$
$\displaystyle=\int\textrm{d}t^{\prime}\,\gamma_{B}^{(n_{f})}(t-t^{\prime},\mu)\,B^{(n_{f})}_{i}(t^{\prime},x,\mu)\,,$
$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}S^{(n_{f})}(\ell,\mu)$
$\displaystyle=\int\textrm{d}\ell^{\prime}\,\gamma_{S}^{(n_{f})}(\ell-\ell^{\prime},\mu)\,S^{(n_{f})}(\ell^{\prime},\mu)\,.$
(30)

In contrast to eq. (4), there is no rapidity evolution in ${\rm SCET}_{\rm I}$
for massless quarks. Consistency of the RG evolution implies that

$\displaystyle\omega_{a}\gamma_{B}^{(n_{f})}(\omega_{a}\ell,\mu)+\omega_{b}\gamma_{B}^{(n_{f})}(\omega_{b}\ell,\mu)+\gamma_{S}^{(n_{f})}(\ell,\mu)=\gamma_{H}^{(n_{f})}(Q,\mu)\,\delta(\ell)\,.$
(31)

For beam thrust the number of possible scale hierarchies with a massive quark
is larger due to the fact that the (massless) collinear and soft modes have
different invariant mass scales. The discussion for the hierarchies with
$\sqrt{Q{\mathcal{T}}}\ll m$ where the massive quark cannot be produced via
real emissions, is completely identical to $q_{T}\ll m$, since the quark mass
effects in these cases are independent of the low-energy measurement. For
$m\sim Q$, all mass effects are encoded by using the mass-dependent hard
function from sec. 2.2 in eq. (3.1) together with $n_{f}=n_{l}$ everywhere
else. Similarly, the case $\sqrt{Q{\mathcal{T}}}\ll m\ll Q$ is described by
using eq. (3.1) with $n_{f}=n_{l}$, and replacing the hard function by the
product of massless $(n_{l}+1)$-flavor hard function and the soft and
collinear mass-mode functions $H_{s}$ and $H_{c}$, as for the case $q_{T}\ll
m\ll Q$ in sec. 2.3. We therefore proceed directly to the hierarchies
$m\lesssim\sqrt{Q{\mathcal{T}}}$, where the massive quark can be produced in
collinear and/or soft real radiation. The four possible hierarchies and the
relevant EFT modes in the $p^{+}p^{-}$-plane are illustrated in fig. 5, and
are discussed in the following subsections.

(a) ${{\mathcal{T}}\ll m\sim\sqrt{Q{\mathcal{T}}}}$

(b) ${{\mathcal{T}}\ll m\ll\sqrt{Q{\mathcal{T}}}}$

(c) ${{\mathcal{T}}\sim m\ll\sqrt{Q{\mathcal{T}}}}$

(d) ${m\ll{\mathcal{T}}\ll\sqrt{Q{\mathcal{T}}}}$

Figure 5: Effective theory modes for the beam thrust spectrum with massive
quarks for $m^{2}/Q\lesssim{\mathcal{T}}\ll Q$ and $m\gg\Lambda_{\rm QCD}$.

###  3.2 Quark mass effects for $\sqrt{Q{\mathcal{T}}}\sim m\ll Q$

For $\sqrt{Q{\mathcal{T}}}\sim m\ll Q$ massive quarks can be produced via
collinear initial-state radiation, but not via soft real radiation. After the
hard matching, carried out with $n_{l}+1$ massless quark flavors as discussed
in sec. 2.3, the degrees of freedom in the EFT are collinear and soft modes
with the scaling

$\displaystyle n_{a}\text{-collinear + MM:}$ $\displaystyle\quad
p_{n_{a}}^{\mu}\sim({\mathcal{T}},Q,\sqrt{Q{\mathcal{T}}})\sim\Bigl{(}\frac{m^{2}}{Q},Q,m\Bigr{)}\,,$
$\displaystyle n_{b}\text{-collinear + MM:}$ $\displaystyle\quad
p_{n_{b}}^{\mu}\sim(Q,{\mathcal{T}},\sqrt{Q{\mathcal{T}}})\sim\Bigl{(}Q,\frac{m^{2}}{Q},m\Bigr{)}\,,$
soft MM: $\displaystyle\quad p_{s}^{\mu}\sim(m,m,m)\,,$ usoft:
$\displaystyle\quad
p_{us}^{\mu}\sim({\mathcal{T}},{\mathcal{T}},{\mathcal{T}})\,,$ (32)

as illustrated in fig. 5(a). While the usual usoft modes live at a lower
virtuality scale than the collinear modes, the soft mass-modes are separated
from the collinear modes only in rapidity, leading to a mix of ${\rm
SCET}_{\rm I}$ and ${\rm SCET}_{\rm II}$ features. In particular, there will
be mass-related rapidity divergences.

At the scale $\mu\sim\sqrt{Q{\mathcal{T}}}\sim m$ this theory with $n_{l}+1$
dynamical quark flavors is matched onto a theory with $n_{l}$ flavors
integrating out also fluctuations related to initial-state collinear radiation
of massless particles. The matching in the collinear sectors leads to mass-
dependent beam function coefficients $\mathcal{I}_{ik}$,

$\displaystyle
B_{i}^{(n_{l}+1)}\Bigl{(}t,m,x,\mu,\frac{\nu}{\omega}\Bigr{)}=\sum_{k\in\\{q,\bar{q},g\\}}\mathcal{I}_{ik}\Bigl{(}t,m,x,\mu,\frac{\nu}{\omega}\Bigr{)}\otimes_{x}f^{(n_{l})}_{k}(x,\mu)\biggl{[}1+\mathcal{O}\Bigl{(}\frac{\Lambda_{\mathrm{QCD}}^{2}}{m^{2}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{t}\Bigr{)}\biggr{]}\,,$
(33)

analogous to eq. (14). The dependence on the rapidity scale $\nu$ here arises
due to virtual secondary massive quark corrections and is the same as for the
collinear mass-mode function $H_{c}$ in eq. (2.3), i.e.,

$\displaystyle\nu\frac{\textrm{d}}{\textrm{d}\nu}B_{i}^{(n_{l}+1)}\Bigl{(}t,m,x,\mu,\frac{\nu}{\omega}\Bigr{)}=\gamma_{\nu,H_{c}}(m,\mu)\,B_{i}^{(n_{l}+1)}\Bigl{(}t,m,x,\mu,\frac{\nu}{\omega}\Bigr{)}\,.$
(34)

In the soft sector the soft mass modes are integrated out, leaving only the
usoft modes. This gives exactly the soft mass-mode function $H_{s}$ in eq.
(2.3), which encodes the effects of virtual secondary massive quark radiation.
As usual, also the strong coupling constant has to be matched from $n_{l}+1$
to $n_{l}$ flavors. The remaining contributions at the lower scales, the soft
function and the PDFs, are given in terms of $n_{l}$ massless flavors and in
the $(n_{l})$-scheme for $\alpha_{s}$. The resulting factorized cross section
reads

$\displaystyle\frac{\textrm{d}\sigma}{\textrm{d}Q^{2}\,\textrm{d}Y\,\textrm{d}{\mathcal{T}}}$
$\displaystyle=\sum_{i,j\in\\{q,\bar{q},Q,\bar{Q}\\}}H^{(n_{l}+1)}_{ij}(Q,\mu)\,H_{s}(m,\mu,\nu)\int\\!\textrm{d}t_{a}\,\textrm{d}t_{b}\,$
$\displaystyle\quad\times\biggl{[}\sum_{k\in\\{q,\bar{q},g\\}}\mathcal{I}_{ik}\Bigl{(}t_{a},m,x_{a},\mu,\frac{\nu}{\omega_{a}}\Bigr{)}\otimes_{x}f^{(n_{l})}_{k}(x_{a},\mu)\biggr{]}$
$\displaystyle\quad\times\biggl{[}\sum_{k\in\\{q,\bar{q},g\\}}\mathcal{I}_{jk}\Bigl{(}t_{b},m,x_{b},\mu,\frac{\nu}{\omega_{a}}\Bigr{)}\otimes_{x}f^{(n_{l})}_{k}(x_{b},\mu)\biggr{]}$
$\displaystyle\quad\times
S^{(n_{l})}\Bigl{(}{\mathcal{T}}-\frac{t_{a}}{\omega_{a}}-\frac{t_{b}}{\omega_{b}},\mu\Bigr{)}\,\biggl{[}1+\mathcal{O}\Bigl{(}\frac{{\mathcal{T}}}{Q},\frac{m^{2}}{Q^{2}},\frac{{\mathcal{T}}^{2}}{m^{2}},\frac{\Lambda_{\mathrm{QCD}}}{{\mathcal{T}}}\Bigr{)}\biggr{]}\,.$
(35)

The resummation of logarithms in eq. (3.2) is obtained by evolving all
functions from their natural scales, as illustrated in fig. 6(a). The mass-
dependent $\nu$ evolution, which resums the rapidity logarithms $\ln(Q/m)$, is
identical to the one for the hard functions $H_{c}$ and $H_{s}$ in sec. 2.3.
The $\mu$ evolution can be conveniently carried out by evolving the hard,
beam, and soft functions with $n_{l}+1$ active flavors above the mass scale
and with $n_{l}$ active flavors below the mass scale, which automatically
takes into account the $\mu$ dependence of $H_{S}$. To see this, the
consistency of RG running for eq. (3.2) together with the consistency relation
for $n_{l}+1$ massless quarks in eq. (31) implies

$\displaystyle\omega_{a}\gamma_{B,m}^{(n_{l}+1)}\Bigl{(}\omega_{a}\ell,m,\mu,\frac{\nu}{\omega_{a}}\Bigr{)}+\omega_{b}\gamma_{B,m}^{(n_{l}+1)}\Bigl{(}\omega_{b}\ell,m,\mu,\frac{\nu}{\omega_{b}}\Bigr{)}+\gamma_{S}^{(n_{l})}(\ell,\mu)+\gamma_{H_{s},\mu}(m,\mu,\nu)\,\delta(\ell)$
$\displaystyle\qquad=\omega_{a}\gamma_{B}^{(n_{l}+1)}(\omega_{a}\ell,\mu)+\omega_{b}\gamma_{B}^{(n_{l}+1)}(\omega_{b}\ell,\mu)+\gamma_{S}^{(n_{l}+1)}(\ell,\mu)\,,$
(36)

where $\gamma_{S}^{(n_{l})}$, $\gamma_{S}^{(n_{l}+1)}$,
$\gamma_{B}^{(n_{l}+1)}$ are the anomalous dimensions for the soft and beam
functions with $n_{l}$ and $n_{l}+1$ massless flavors as defined in eq. (3.1),
and $\gamma_{B,m}^{(n_{l}+1)}(t,m,\mu,\nu/\omega)$ is the anomalous dimension
of the mass-dependent beam function,

$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}B_{i}^{(n_{l}+1)}\Bigl{(}t,m,x,\mu,\frac{\nu}{\omega}\Bigr{)}=\int\textrm{d}t^{\prime}\,\gamma_{B,m}^{(n_{l}+1)}\Bigl{(}t-t^{\prime},m,\mu,\frac{\nu}{\omega}\Bigr{)}\,B_{i}^{(n_{l}+1)}\Bigl{(}t^{\prime},m,x,\mu,\frac{\nu}{\omega}\Bigr{)}\,.$
(37)

The consistency relation in eq. (3.2) can be confirmed explicitly at two loops
with the expressions in eqs. (A.2.2), (A.3.2), (4.1.2), and (73). Note that
this relation does not imply that
$\gamma_{B,m}^{(n_{l}+1)}(t,m,\mu,\nu/\omega)$ and
$\gamma_{B}^{(n_{l}+1)}(t,\mu)$ are the same, which is indeed not the case for
the massive quark corrections as we will see explicitly in sec. 4.2. The
reason is that the presence of the quark mass leads to a ${\rm SCET}_{\rm
II}$-type theory, in which the required rapidity regularization redistributes
the $\mu$ anomalous dimension between soft and collinear corrections with
individually regularization scheme dependent pieces. Only their sum, as given
on the left-hand side of eq. (3.2), is independent of the regularization
scheme and yields the combined running for beam and soft functions with
$n_{l}+1$ massless flavors above $\mu_{m}\sim m$, as on the right-hand side of
eq. (3.2).

(a) ${{\mathcal{T}}\ll m\sim\sqrt{Q{\mathcal{T}}}}$

(b) ${{\mathcal{T}}\ll m\ll\sqrt{Q{\mathcal{T}}}}$

(c) ${{\mathcal{T}}\sim m\ll\sqrt{Q{\mathcal{T}}}}$

Figure 6: Illustration of the renormalization group evolution for beam thrust
of the hard, beam, soft, and parton distribution function in invariant mass
and rapidity. The anomalous dimensions for each evolution step involve the
displayed number of active quark flavors. The label $m$ indicates that the
corresponding evolution is mass dependent.

###  3.3 Quark mass effects for ${\mathcal{T}}\ll m\ll\sqrt{Q{\mathcal{T}}}$

When the beam scale becomes larger than the mass scale, but the soft scale is
still larger than the mass, which happens for $m^{2}/Q\ll{\mathcal{T}}\ll m$,
the beam function matching coefficients $\mathcal{I}_{ik}$ encode only
fluctuations related to initial-state collinear radiation with $n_{l}+1$
massless quarks. The EFT below $\sqrt{Q{\mathcal{T}}}$ contains the usual
collinear and soft mass modes scaling as $p_{m,n_{a}}^{\mu}\sim(m^{2}/Q,Q,m)$,
$p_{m,n_{b}}^{\mu}\sim(Q,m^{2}/Q,m)$, and $p_{m,s}^{\mu}\sim(m,m,m)$, which do
not contribute to the beam thrust measurement. However, besides these there
are also additional modes with fluctuations around the mass scale which can
have a dynamic impact on the ${\mathcal{T}}$ spectrum in this hierarchy, as
illustrated in fig. 5(b). Their scaling is precisely determined by this
condition and the on-shell constraint, yielding the scaling

$\displaystyle n_{a}\text{-csoft MM:}$ $\displaystyle\quad
p_{cs,n_{a}}^{\mu}\sim\Bigl{(}{\mathcal{T}},\frac{m^{2}}{{\mathcal{T}}},m\Bigr{)}\,,$
$\displaystyle n_{b}\text{-csoft MM:}$ $\displaystyle\quad
p_{cs,n_{b}}^{\mu}\sim\Bigl{(}\frac{m^{2}}{{\mathcal{T}}},{\mathcal{T}},m\Bigr{)}\,.$
(38)

We refer to these intermediate modes as collinear-soft (csoft), since they are
simultaneously boosted (by a factor $m/{\mathcal{T}}$) but are softer than the
standard collinear modes, thus coupling to the latter via Wilson lines and
leading to a ${\rm SCET}_{+}$ theory Bauer:2011uc . This type of intermediate
${\rm SCET}_{+}$ modes have appeared in various contexts Bauer:2011uc ;
Procura:2014cba ; Larkoski:2015zka ; Pietrulewicz:2016nwo . The setup here is
similar to the case of double-differential distributions with a simultaneous
$q_{T}$ and beam thrust measurement discussed in ref. Procura:2014cba . Also
there, several hierarchies are possible ranging from a ${\rm SCET}_{\rm II}$
regime for $q_{T}\sim{\mathcal{T}}$ to a ${\rm SCET}_{\rm I}$ regime for
$q_{T}\sim\sqrt{Q{\mathcal{T}}}$ with a ${\rm SCET}_{+}$ regime in between.
The csoft modes in their ${\rm SCET}_{+}$ regime are separated from the
collinear modes only in rapidity. In our case here, the csoft mass modes are
separated in invariant mass from the standard ${\rm SCET}_{\rm I}$ soft and
collinear modes and in rapidity from their ${\rm SCET}_{\rm II}$-type soft
mass-mode cousins.

The matching in the collinear sector can be performed in two steps as in eqs.
(18) and (19). After integrating out all of the mass modes, the PDF and the
soft function are still given in a $(n_{l})$-flavor theory. Thus the
factorization formula reads

$\displaystyle\frac{\textrm{d}\sigma}{\textrm{d}Q^{2}\,\textrm{d}Y\,\textrm{d}{\mathcal{T}}}$
$\displaystyle=\sum_{i,j\in\\{q,\bar{q},Q,\bar{Q}\\}}H^{(n_{l}+1)}_{ij}(Q,\mu)\,H_{s}(m,\mu,\nu)\int\\!\textrm{d}k_{a}\,\textrm{d}k_{b}\,\mathcal{S}_{c}(k_{a},m,\mu,\nu)\,\mathcal{S}_{c}(k_{b},m,\mu,\nu)$
$\displaystyle\quad\times\int\\!\textrm{d}t_{a}\,\biggl{[}\sum_{k\in\\{q,\bar{q},Q,\bar{Q},g\\}}\sum_{l\in\\{q,\bar{q},g\\}}\mathcal{I}^{(n_{l}+1)}_{ik}(t_{a},x_{a},\mu)\otimes_{x}\mathcal{M}_{kl}(x_{a},m,\mu)\otimes_{x}f^{(n_{l})}_{k}(x_{a},\mu)\biggr{]}$
$\displaystyle\quad\times\int\\!\textrm{d}t_{b}\,\biggl{[}\sum_{k\in\\{q,\bar{q},Q,\bar{Q},g\\}}\sum_{l\in\\{q,\bar{q},g\\}}\mathcal{I}^{(n_{l}+1)}_{jk}(t_{b},x_{b},\mu)\otimes_{x}\mathcal{M}_{kl}(x_{b},m,\mu)\otimes_{x}f^{(n_{l})}_{n}(x_{b},\mu)\biggr{]}$
$\displaystyle\quad\times
S^{(n_{l})}\Bigl{(}{\mathcal{T}}-\frac{t_{a}}{\omega_{a}}-\frac{t_{b}}{\omega_{b}}-k_{a}-k_{b},\mu\Bigr{)}\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{{\mathcal{T}}}{Q},\frac{m^{2}}{Q{\mathcal{T}}},\frac{{\mathcal{T}}^{2}}{m^{2}},\frac{\Lambda_{\mathrm{QCD}}}{{\mathcal{T}}}\Bigr{)}\Bigr{]}\,.$
(39)

The functions $\mathcal{S}_{c}$ here are the csoft matching functions encoding
the interactions of the collinear-soft radiation at the invariant mass scale
$\mu\sim m$ and the rapidity scale $\nu\sim m^{2}/{\mathcal{T}}$. The
$\mathcal{M}_{ij}$ correspond to the well-known PDF matching correction
incorporating the effect of the collinear mass modes, as in eq. (2.5). The
virtual soft massive quark corrections are still described by the function
$H_{s}$ at the rapidity scale $\nu\sim m$ as in eq. (3.2).

The RG evolution for eq. (3.3) is illustrated in fig. 6(b). The csoft function
satisfies the same rapidity RGE as the collinear mass-mode function $H_{c}$ in
eq. (2.3) and the massive beam functions in eq. (34), i.e.,

$\displaystyle\nu\frac{\textrm{d}}{\textrm{d}\nu}\mathcal{S}_{c}(k,m,\mu,\nu)=\gamma_{\nu,H_{c}}(m,\mu)\,\mathcal{S}_{c}(k,m,\mu,\nu)\,.$
(40)

The only difference with respect to the rapidity evolution in eq. (3.2) is
that it now happens between $H_{s}$ and $\mathcal{S}_{c}$ with
$\nu_{\mathcal{S}_{c}}\sim m^{2}/{\mathcal{T}}$ rather than between $H_{s}$
and the beam functions with $\nu_{B}\sim Q$, such that now the (smaller)
rapidity logarithms $\ln(m/{\mathcal{T}})$ are resummed. The $\mu$ evolution
can be performed with $n_{l}+1$ flavors for the hard function $H_{ij}$, the
beam and soft function above the mass scale and with $n_{l}$ flavors below.
This automatically accounts for the $\mu$ dependence of $\mathcal{S}_{c}$ and
$H_{s}$ above $\mu_{m}\sim m$, which precisely gives the difference between
the evolution of the soft function with $n_{l}+1$ and $n_{l}$ flavors, as
implied by the consistency of RG running for eq. (3.3) and the relation in eq.
(31) with $n_{l}+1$ massless quarks,

$\displaystyle\gamma_{S}^{(n_{l})}(\ell,\mu)+2\gamma_{\mathcal{S}_{c}}(\ell,m,\mu,\nu)+\delta(\ell)\,\gamma_{H_{s}}(m,\mu,\nu)=\gamma_{S}^{(n_{l}+1)}(\ell,\mu)\,,$
(41)

where

$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}\mathcal{S}_{c}(k,m,\mu,\nu)=\int\textrm{d}k^{\prime}\,\gamma_{\mathcal{S}_{c}}(k-k^{\prime},m,\mu,\nu)\,\mathcal{S}_{c}(k,m,\mu,\nu)\,.$
(42)

At two loops, the consistency relation eq. (41) can be explicitly confirmed
with the expressions in eqs. (A.3.2), (80), and (4.1.2).

###  3.4 Quark mass effects for ${\mathcal{T}}\sim m$ and $m\ll{\mathcal{T}}$

For ${\mathcal{T}}\sim m$ the csoft and soft mass modes in the previous
section merge with the usual usoft modes,

usoft: $\displaystyle\quad
p_{s}^{\mu}\sim({\mathcal{T}},{\mathcal{T}},{\mathcal{T}})\sim\Bigl{(}{\mathcal{T}},\frac{m^{2}}{{\mathcal{T}}},m\Bigr{)}\sim(m,m,m)\,.$
(43)

In this hierarchy massive quarks can be also produced in soft real radiation
leading to a soft function at the scale $\mu\sim{\mathcal{T}}$ that depends on
the quark mass. In addition, there are the usual collinear modes as well as
the collinear mass modes, as illustrated in fig. 5(c). Since we still have
$m\ll\sqrt{Q{\mathcal{T}}}$, the matching in the collinear sectors is the same
as in the previous subsection. The factorization formula reads

$\displaystyle\frac{\textrm{d}\sigma}{\textrm{d}Q^{2}\,\textrm{d}Y\,\textrm{d}{\mathcal{T}}}$
$\displaystyle=\sum_{i,j\in\\{q,\bar{q},Q,\bar{Q}\\}}H^{(n_{l}+1)}_{ij}(Q,\mu)\int\\!\textrm{d}t_{a}\,\textrm{d}t_{b}$
$\displaystyle\quad\times\biggl{[}\sum_{k\in\\{q,\bar{q},Q,\bar{Q},g\\}}\sum_{l\in\\{q,\bar{q},g\\}}\mathcal{I}^{(n_{l}+1)}_{ik}(t_{a},x_{a},\mu)\otimes_{x}\mathcal{M}_{kl}(x_{a},m,\mu)\otimes_{x}f^{(n_{l})}_{k}(x_{a},\mu)\biggr{]}$
$\displaystyle\quad\times\biggl{[}\sum_{k\in\\{q,\bar{q},Q,\bar{Q},g\\}}\sum_{l\in\\{q,\bar{q},g\\}}\mathcal{I}^{(n_{l}+1)}_{jk}(t_{b},x_{b},\mu)\otimes_{x}\mathcal{M}_{kl}(x_{b},m,\mu)\otimes_{x}f^{(n_{l})}_{k}(x_{b},\mu)\biggr{]}$
$\displaystyle\quad\times
S\Bigl{(}{\mathcal{T}}-\frac{t_{a}}{\omega_{a}}-\frac{t_{b}}{\omega_{b}},m,\mu\Bigr{)}\biggl{[}1+\mathcal{O}\Bigl{(}\frac{{\mathcal{T}}}{Q},\frac{m^{2}}{Q{\mathcal{T}}},\frac{\Lambda_{\mathrm{QCD}}}{{\mathcal{T}}},\frac{\Lambda_{\mathrm{QCD}}^{2}}{m^{2}}\Bigr{)}\biggr{]}\,.$
(44)

Now all rapidity divergences cancel within the soft function and do not leave
behind any potentially large rapidity logarithms. The RG evolution for this
case is illustrated in fig. 6(c).

Finally, for $m\ll{\mathcal{T}}$ the mass dependence in the IR insensitive
soft function vanishes, if expressed in terms of the $(n_{l}+1)$-flavor scheme
for $\alpha_{s}$. Otherwise, eq. (3.4) remains unchanged, such that now the
only dependence on the mass scale arises in the PDF matching corrections
$\mathcal{M}_{ij}$. The hard, beam, and soft functions can now be always
evolved with $n_{l}+1$ massless flavors and only the evolution of the PDF
changes, when crossing the flavor threshold.

###  3.5 Relations between hierarchies

Figure 7: Relevant modes for the beam thrust spectrum with ${\mathcal{T}}\ll
Q$ for different hierarchies between the quark mass $m$ and the scales
${\mathcal{T}}$, $\sqrt{Q{\mathcal{T}}}$ and $Q$. The arrows indicate the
relations between the modes and their associated contributions.

We now discuss how the ingredients appearing in the different factorization
formulae are related to each other. The relations between the modes and their
contributions are illustrated in fig. 7 for the different possible
hierarchies. As in sec. 2.6, these relations show how one can combine the
resummation of logarithms relevant in one regime with the power-suppressed
fixed-order content that becomes important in the neighboring regimes,
enabling a systematic inclusion of mass corrections across the entire
${\mathcal{T}}$ spectrum.

Similar to eq. (23), the mass-dependent beam function coefficients appearing
for $\sqrt{Q{\mathcal{T}}}\sim m$ (incorporating massive quark fluctuation as
discussed in sec. 3.2) are related to those for $\sqrt{Q{\mathcal{T}}}\ll m$
with $n_{l}$ massless quarks and the collinear mass-mode function $H_{c}$ by

$\displaystyle\mathcal{I}_{ik}\Bigl{(}t,m,x,\mu,\frac{\nu}{\omega}\Bigr{)}=H_{c}\Bigl{(}m,\mu,\frac{\nu}{\omega}\Bigr{)}\,\mathcal{I}_{ik}^{(n_{l})}(t,x,\mu)\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{t}{m^{2}}\Bigr{)}\Bigr{]}\,.$
(45)

At the same time, the mass-dependent beam function also encodes information
about the fixed-order content for ${\mathcal{T}}\ll
m\ll\sqrt{Q{\mathcal{T}}}$. Comparing eqs. (3.2) and (3.3), they are related
to those with $n_{l}+1$ massless flavors, the PDF matching functions, and the
csoft function $\mathcal{S}_{c}$ by

$\displaystyle\mathcal{I}_{ik}\Bigl{(}t,m,x,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=\\!\\!\sum_{j=\\{q,\bar{q},Q,\bar{Q},g\\}}\int\\!\textrm{d}\ell\,\mathcal{I}_{ij}^{(n_{l}+1)}(t-Q\ell,x,\mu)\otimes_{x}\mathcal{M}_{jk}(x,m,\mu)\,\mathcal{S}_{c}(\ell,m,\mu,\nu)$
$\displaystyle\quad\times\biggl{[}1+\mathcal{O}\Bigl{(}\frac{m^{2}}{Q{\mathcal{T}}}\Bigr{)}\biggr{]}\,.$
(46)

The mass-dependent soft function for ${\mathcal{T}}\sim m$ in eq. (3.4)
contains massive quark fluctuations that for ${\mathcal{T}}\ll m$ get split
into the massless soft function with $n_{l}$ flavors, the soft mass mode
function $H_{s}$, and the csoft function $\mathcal{S}_{c}$ in eq. (3.3) as

$\displaystyle
S(\ell,m,\mu)=H_{s}(m,\mu,\nu)\int\\!\textrm{d}\ell^{\prime}\,S^{(n_{l})}(\ell-\ell^{\prime},\mu)\,\mathcal{S}_{c}(\ell^{\prime},m,\mu,\nu)\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{\ell^{2}}{m^{2}}\Bigr{)}\Bigr{]}\,.$
(47)

Finally, as already mentioned below eq. (3.4), the soft function approaches
its massless limit for $m\ll{\mathcal{T}}$,

$\displaystyle
S(\ell,m,\mu)=S^{(n_{l}+1)}(\ell,\mu)\Bigl{[}1+\mathcal{O}\Bigl{(}\frac{m^{2}}{\ell^{2}}\Bigr{)}\Bigr{]}\,.$
(48)

###  3.6 Relation to previous literature

Here, we briefly comment on the connection of the factorization setup
presented here for beam thrust to the closely related ${\rm SCET}_{\rm I}$
setup in refs. Gritschacher:2013pha ; Pietrulewicz:2014qza  for thrust in
$e^{+}e^{-}$-collisions (or similarly also for DIS with $x\to 1$ Hoang:2015iva
). Besides the fact that the jet functions appearing for thrust in
$e^{+}e^{-}$ are replaced by virtuality-dependent beam functions for beam
thrust in $pp$ collisions, there are also some differences in the description
of the different regimes. While we have discussed each possible hierarchy in a
strict EFT sense identifying a single operator matrix element or matching
function with each EFT mode, refs. Gritschacher:2013pha ; Pietrulewicz:2014qza
already set up their factorization theorems in a way that they apply for
neighboring hierarchies (e.g. ${\mathcal{T}}\ll m\sim\sqrt{Q{\mathcal{T}}}$
and ${\mathcal{T}}\ll m\ll\sqrt{Q{\mathcal{T}}}$). Using appropriate
renormalization conditions, the mass dependent corrections to the jet and soft
functions were assigned such that they directly give the massless results in
the small mass limit and decouple in the infinite mass limit. In addition, the
factorization theorems contained mass mode matching functions for hard, jet,
and soft function, whenever the evolution of one of the matrix elements
crossed the mass scale. In our setup this essentially amounts to a specific
practical choice how to incorporate the power corrections in eqs. (45)–(48).
Although the final outcome is thus essentially the same once the correct
rapidity scales are chosen in the mass mode matching functions, it is perhaps
more transparent conceptionally to first distinguish all hierarchies with the
associated modes as we do here, and separately discuss the possible ways to
add the nonsingular corrections later. In particular, for the hierarchy
${\mathcal{T}}\ll m\ll\sqrt{Q{\mathcal{T}}}$, this leads us to identify the
csoft modes as a relevant degree of freedom with a corresponding function
evaluated naturally at the rapidity scale $\nu\sim m^{2}/{\mathcal{T}}$. In
contrast, refs. Gritschacher:2013pha ; Pietrulewicz:2014qza  the corresponding
corrections appeared inside the mass mode matching functions as soft-bin
contributions that had to be evaluated at this rapidity scale to minimize
large rapidity logarithms.

##  4 Results for massive quark corrections

In this section we present our results for the contributions from primary
massive quarks at $\mathcal{O}(\alpha_{s})$ and from secondary massive quarks
at $\mathcal{O}(\alpha_{s}^{2})$ to all components of the various
factorization theorems discussed in secs. 2 and 3, providing all required
ingredients for the Drell-Yan spectrum at NNLL′. The results in this section
are only given for a single massive quark flavor and with the rapidity
divergences regularized by the symmetric Wilson line regulator introduced in
refs. Chiu:2011qc ; Chiu:2012ir . The actual computations of the primary and
secondary massive quark corrections to the beam and soft functions are carried
out in some detail in app. B. In sec. 4.4, we show explicitly that the results
satisfy the small and large mass limits, and illustrate the numerical size of
the mass-dependent corrections for the case of $b$ quarks.

The fixed-order results for the mass-dependent corrections can be expanded
either in terms of the $(n_{l})$-flavor or $(n_{l}+1)$-flavor scheme for
$\alpha_{s}$. For definiteness we expand in this section any function $F(m)$
using $\alpha_{s}^{(n_{l}+1)}$,

$\displaystyle
F(m)=\sum_{n=0}^{\infty}\biggl{(}\frac{\alpha_{s}^{(n_{l}+1)}(\mu)}{4\pi}\biggr{)}^{n}\,F^{(n)}(m)\;.$
(49)

The different two-loop contributions to $F^{(2)}(m)$ are written as

$\displaystyle
F^{(2)}(m)=T_{F}n_{l}\,F^{(2,l)}+T_{F}F^{(2,h)}(m)+\,.\,.\,.\;,$ (50)

where $F^{(2,h)}$ contains all mass dependent two-loop corrections and
$F^{(2,l)}$ the associated contributions for massless flavors. The expansion
of $F$ in terms of $\alpha_{s}^{(n_{l})}$ can be easily obtained by using the
matching relation for $\alpha_{s}$,

$\displaystyle\alpha_{s}^{(n_{l}+1)}(\mu)=\alpha_{s}^{(n_{l})}(\mu)\biggl{[}1-\frac{\alpha_{s}^{(n_{l})}(\mu)T_{F}}{4\pi}\,\frac{4}{3}L_{m}+\mathcal{O}(\alpha_{s}^{2})\biggr{]}\,,$
(51)

where here and in the following we abbreviate

$\displaystyle L_{m}\equiv\ln\frac{m^{2}}{\mu^{2}}\,.$ (52)

###  4.1 Hard matching functions

All hard matching functions, i.e. the hard function $H$ at the scale $Q$ and
the mass mode matching functions $H_{c}$ and $H_{s}$ at the scale $m\ll Q$,
are insensitive to the measurement performed at a lower scale and are
therefore the same for $q_{T}$ and beam thrust ${\mathcal{T}}$. Since the QCD
and SCET currents are the same as for $e^{+}e^{-}\to 2\;\text{jets}$, the
results can be read off from the corresponding ones in refs.
Pietrulewicz:2014qza ; Hoang:2015vua .

####  4.1.1 Massive quark corrections to the hard function

The secondary massive quark corrections to the hard function in eq. (9) read

$\displaystyle H^{(2,h)}(Q,m,\mu)$
$\displaystyle=H^{(0)}(Q)\biggl{\\{}C_{F}h_{\rm
virt}\Bigl{(}\frac{m^{2}}{Q^{2}}\Bigr{)}+\frac{4}{3}L_{m}\,H^{(1)}(Q,\mu)\biggr{\\}}\,,$
(53)

where $H^{(0)}$ denotes the tree-level normalization and $H^{(1)}$ the
massless one-loop contribution given in eq. (116). The function $h_{\rm virt}$
contains the $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ virtual massive quark
bubble correction in full QCD shown in fig. 1. It has been calculated in refs.
Kniehl1990 ; Hoang:1995fr  and is given by

$\displaystyle h_{\rm virt}(x)$
$\displaystyle=\Bigl{(}16x^{2}-\frac{8}{3}\Bigr{)}\biggl{[}-4\textrm{Li}_{3}\Bigl{(}\frac{r-1}{r+1}\Bigr{)}-\frac{1}{3}\ln^{3}\frac{r-1}{r+1}+\frac{2\pi^{2}}{3}\ln\frac{r-1}{r+1}+4\zeta_{3}\biggr{]}$
$\displaystyle\quad+r\Bigl{(}\frac{184}{9}x+\frac{76}{9}\Bigr{)}\biggl{[}4\textrm{Li}_{2}\Bigl{(}\frac{r-1}{r+1}\Bigr{)}+\ln^{2}\frac{r-1}{r+1}-\frac{2\pi^{2}}{3}\biggr{]}+\Bigl{(}\frac{880}{9}x+\frac{1060}{27}\Bigr{)}\ln
x$ $\displaystyle\quad+\frac{1904}{9}x+\frac{6710}{81}\,,$ (54)

with $r=\sqrt{1+4x}$. For $m\to\infty$ the massive quark decouples such that
$h_{\rm virt}(x)\to 0$ for $x\to\infty$.

For $Z$-boson production there is an additional primary massive quark
contribution to the axial vector current, namely the massive quark triangle
correction in fig. 1, which we denote by $\Delta h_{\rm axial}$ with the same
prefactor as for $h_{\rm virt}$ using the narrow width approximation for
notational simplicity. It has been computed in refs. Kniehl:1989qu ;
Gonsalves:1991qn ; Bernreuther:2005rw  and is given by

$\displaystyle\Delta h_{\rm axial}(Q,m,\mu)$
$\displaystyle=\frac{8a_{q}a_{Q}}{v_{q}^{2}+a_{q}^{2}}\biggl{[}3\ln\frac{Q^{2}}{\mu^{2}}-9+\frac{\pi^{2}}{3}+\theta(Q^{2}-4m^{2})\,G_{1}\Bigl{(}\frac{m^{2}}{Q^{2}}\Bigr{)}$
$\displaystyle\qquad\qquad+\theta(4m^{2}-Q^{2})\,G_{2}\Bigl{(}\frac{m^{2}}{Q^{2}}\Bigr{)}\biggr{]}\,,$
(55)

where the vector and axial vector couplings for up- and down-type quarks are
proportional to $v_{u}=1-8/3\sin^{2}\theta_{W}$,
$v_{d}=-1+4/3\sin^{2}\theta_{W}$, $a_{u}=1$, $a_{d}=-1$. The functions $G_{1}$
and $G_{2}$ are given in eqs. (2.8) and (2.9) of ref. Gonsalves:1991qn . In
the small mass limit $m\ll Q$ the function $G_{1}(m^{2}/Q^{2})$ vanishes, such
that $\Delta h_{\rm axial}$ gives the same result as for a massless flavor in
the loop,

$\displaystyle\Delta h_{\rm
axial}(Q,m,\mu)=\frac{8a_{q}a_{Q}}{v_{q}^{2}+a_{q}^{2}}\biggl{[}3\ln\frac{Q^{2}}{\mu^{2}}-9+\frac{\pi^{2}}{3}+\mathcal{O}\Bigl{(}\frac{m^{2}}{Q^{2}}\Bigr{)}\biggr{]}\,.$
(56)

For a massless isospin partner this correction is thus canceled within the
$SU(2)_{L}$ doublet, while for different masses (as for $m_{b}\ll m_{t}$)
there is a ($\mu$-independent) remainder. Note that for $Q\ll m$ the function
$\Delta h_{\rm axial}$ gives a nonvanishing contribution

$\displaystyle\Delta h_{\rm
axial}(Q,m,\mu)=\frac{8a_{q}a_{Q}}{v_{q}^{2}+a_{q}^{2}}\biggl{[}3\ln\frac{m^{2}}{\mu^{2}}+\frac{3}{2}+\mathcal{O}\Bigl{(}\frac{Q^{2}}{m^{2}}\Bigr{)}\biggr{]}\,.$
(57)

In this case one would integrate out the heavy quark at the scale $\mu_{m}\sim
m$ and evolve the axial current to $\mu_{H}\sim Q$ to resum logarithms
$\ln(m^{2}/Q^{2})$.

####  4.1.2 Soft and collinear mass-mode matching functions

The contributions to the mass-mode matching functions originate only from
secondary radiation. The soft mass-mode function $H_{s}$ appearing in eqs.
(2.3), (3.2), and (3.3) has been computed at two loops with the symmetric
$\eta$-regulator in ref. Hoang:2015vua . It is given by

$\displaystyle H_{s}(m,\mu,\nu)$
$\displaystyle=1+\frac{\alpha_{s}^{2}C_{F}T_{F}}{16\pi^{2}}\biggl{[}-\Bigl{(}\frac{16}{3}L_{m}^{2}+\frac{160}{9}L_{m}+\frac{448}{27}\Bigr{)}\ln\frac{\nu}{\mu}+\frac{8}{9}L_{m}^{3}+\frac{40}{9}L_{m}^{2}$
$\displaystyle\quad+\Big{(}\frac{448}{27}-\frac{4\pi^{2}}{9}\Big{)}L_{m}+\frac{656}{27}-\frac{10\pi^{2}}{27}-\frac{56\zeta_{3}}{9}\biggr{]}+\mathcal{O}(\alpha_{s}^{3})\,.$
(58)

Since there are no $\mathcal{O}(\alpha_{s})$ corrections, the flavor scheme
for $\alpha_{s}$ does not affect the results at $\mathcal{O}(\alpha_{s}^{2})$.
Its anomalous dimensions are

$\displaystyle\gamma_{H_{s}}(m,\mu,\nu)$
$\displaystyle=\frac{\alpha_{s}^{2}C_{F}T_{F}}{16\pi^{2}}\biggl{[}\Bigl{(}\frac{64}{3}L_{m}+\frac{320}{9}\Bigr{)}\ln\frac{\nu}{\mu}-\frac{448}{27}+\frac{8\pi^{2}}{9}\biggr{]}+\mathcal{O}(\alpha_{s}^{3})\,,$
$\displaystyle\gamma_{\nu,H_{s}}(m,\mu)$
$\displaystyle=\frac{\alpha_{s}^{2}C_{F}T_{F}}{16\pi^{2}}\biggl{[}-\frac{16}{3}L_{m}^{2}-\frac{160}{9}L_{m}-\frac{448}{27}\biggr{]}+\mathcal{O}(\alpha_{s}^{3})\,.$
(59)

The rapidity anomalous dimension is even known at
$\mathcal{O}(\alpha_{s}^{3})$, see ref. Hoang:2015iva .

The result for the collinear mass-mode function $H_{c}$ in eq. (2.3) can be
inferred at $\mathcal{O}(\alpha_{s}^{2})$ from the computations in refs.
Pietrulewicz:2014qza ; Hoang:2015vua  and reads

$\displaystyle H_{c}\Bigl{(}m,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=1+\frac{\alpha_{s}^{2}C_{F}T_{F}}{16\pi^{2}}\biggl{[}\left(\frac{8}{3}L_{m}^{2}+\frac{80}{9}L_{m}+\frac{224}{27}\right)\ln\frac{\nu}{\omega}$
$\displaystyle\quad+2L_{m}^{2}+\Bigl{(}\frac{2}{3}+\frac{8\pi^{2}}{9}\Bigr{)}L_{m}+\frac{73}{18}+\frac{20\pi^{2}}{27}-\frac{8\zeta_{3}}{3}\biggr{]}+\mathcal{O}(\alpha_{s}^{3})\,.$
(60)

Its anomalous dimensions are

$\displaystyle\gamma_{H_{c}}\Bigl{(}m,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=\frac{\alpha_{s}^{2}C_{F}T_{F}}{16\pi^{2}}\biggl{[}-\Bigl{(}\frac{32}{3}L_{m}+\frac{160}{9}\Bigr{)}\ln\frac{\nu}{\omega}-8L_{m}-\frac{4}{3}-\frac{16\pi^{2}}{9}\biggr{]}+\mathcal{O}(\alpha_{s}^{3})\,,$
$\displaystyle\gamma_{\nu,H_{c}}(m,\mu)$
$\displaystyle=\frac{\alpha_{s}^{2}C_{F}T_{F}}{16\pi^{2}}\biggl{(}\frac{8}{3}L_{m}^{2}+\frac{80}{9}L_{m}+\frac{224}{27}\biggr{)}+\mathcal{O}(\alpha_{s}^{3})\,.$
(61)

One can easily verify that the relation in eq. (22) between the massive hard
function in eq. (53), the hard function contribution for a massless flavor in
eq. (116), and the two mass-mode functions in eqs. (4.1.2) and (4.1.2) is
satisfied,

$\displaystyle H^{(2,h)}(Q,m,\mu)$
$\displaystyle=H^{(2,l)}(Q,\mu)+H_{c}^{(2)}\Bigl{(}m,\mu,\frac{\nu}{\omega_{a}}\Bigr{)}+H_{c}^{(2)}\Bigl{(}m,\mu,\frac{\nu}{\omega_{b}}\Bigr{)}+H_{s}^{(2)}(m,\mu,\nu)+\mathcal{O}\Bigl{(}\frac{m^{2}}{Q^{2}}\Bigr{)}.$
(62)

###  4.2 Beam functions

Here we give our results for the massive quark beam function coefficient
$\mathcal{I}_{Qg}$ at $\mathcal{O}(\alpha_{s})$ and the secondary massive
quark corrections to the light-quark coefficients $\mathcal{I}_{qq}$ at
$\mathcal{O}(\alpha_{s}^{2})$, which appear in eqs. (2.4) and (3.2) for the
$q_{T}$ and beam thrust measurement. We also give the massive quark
contributions to the beam function anomalous dimensions. We also give the
well-known results for the corresponding PDF matching coefficients
$\mathcal{\mathcal{M}}_{Qg}$ at $\mathcal{O}(\alpha_{s})$ and
$\mathcal{\mathcal{M}}_{qq}$ at $\mathcal{O}(\alpha_{s}^{2})$ appearing in
eqs. (2.5), (3.3) and (3.4).

####  4.2.1 TMD beam function coefficients

Figure 8: Feynman diagram for the massive quark beam function at one loop.

The matching coefficient $\mathcal{I}_{Qg}$ generating a massive beam function
from a gluon splitting is calculated at $\mathcal{O}(\alpha_{s})$ in sec. B.1
and corresponds to the diagram shown in fig. 8. The result reads
($p_{T}^{2}=|\vec{p}_{T}|^{2}$)

$\displaystyle\mathcal{I}_{Qg}(\vec{p}_{T},m,z)=\mathcal{I}_{\bar{Q}g}(\vec{p}_{T},m,z)=\frac{\alpha_{s}T_{F}}{4\pi^{2}}\,\theta(z)\,\theta(1-z)\,\frac{2}{p_{T}^{2}+m^{2}}\biggl{[}P_{qg}(z)+\frac{2m^{2}z(1-z)}{p_{T}^{2}+m^{2}}\biggr{]}+\mathcal{O}(\alpha_{s}^{2})\,,$
(63)

with the splitting function

$\displaystyle P_{qg}(z)=z^{2}+(1-z)^{2}\,.$ (64)

This result is equivalent to the Fourier transform of the mass-dependent
matching functions $\mathcal{C}_{h/G}$ in ref. Nadolsky:2002jr . After
performing an appropriate crossing it also agrees with the massive final-state
splitting functions Catani:2000ef ; Kang:2016ofv  or fragmenting jet function
Bauer:2013bza .

Figure 9: Secondary massive quark corrections for the light-quark beam
function at $\mathcal{O}(\alpha_{s}^{2})$. In addition, also the wave-function
renormalization correction and the mirror diagrams for (b) and (c) have to be
included.

The contributions from secondary massive quarks to the matching coefficient
$\mathcal{I}_{qq}$ are computed in sec. B.3 at $\mathcal{O}(\alpha_{s}^{2})$.
The corresponding diagrams are shown in fig. 9. The result is given by

$\displaystyle\mathcal{I}_{qq}^{(2,h)}\Bigl{(}\vec{p}_{T},m,z,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle\quad=\theta(z)C_{F}\biggl{\\{}\delta^{(2)}(\vec{p}_{T})\,\delta(1-z)\biggl{[}\Bigl{(}\frac{8}{3}L_{m}^{2}+\frac{80}{9}L_{m}+\frac{224}{27}\Bigr{)}\ln\frac{\nu}{\omega}+2L_{m}^{2}+\Bigl{(}\frac{2}{3}+\frac{8\pi^{2}}{9}\Bigr{)}L_{m}$
$\displaystyle\qquad\quad+\frac{73}{18}+\frac{20\pi^{2}}{27}-\frac{8\zeta_{3}}{3}\biggr{]}$
$\displaystyle\qquad+\frac{16}{9\pi
p_{T}^{2}}\Bigl{[}\mathcal{L}_{0}(1-z)-\delta(1-z)\ln\frac{\nu}{\omega}\Bigr{]}\Bigl{[}-5+12\hat{m}^{2}+3c(1-2\hat{m}^{2})\ln\frac{c+1}{c-1}\Bigr{]}$
$\displaystyle\qquad+\frac{16}{9\pi
p_{T}^{2}}\theta(1-z)\biggl{[}\frac{3}{2d(1-z)}\bigl{[}(1+z^{2})(1+2\hat{m}^{2}z)+4\hat{m}^{4}z^{2}(-5+6z-5z^{2})\bigr{]}\ln\frac{d+1}{d-1}$
$\displaystyle\qquad\quad-\frac{3c(1-2\hat{m}^{2})}{1-z}\ln\frac{c+1}{c-1}+1+4z+3\hat{m}^{2}(-4+z-5z^{2})\biggr{]}\biggr{\\}}$
$\displaystyle\qquad+\frac{4}{3}L_{m}\,\mathcal{I}_{qq}^{(1)}\Bigl{(}\vec{p}_{T},z,\mu,\frac{\nu}{\omega}\Bigr{)}\,,$
(65)

where

$\displaystyle\hat{m}\equiv\frac{m}{p_{T}}\,,\qquad
c=\sqrt{1+4\hat{m}^{2}}\,,\qquad d=\sqrt{1+4\hat{m}^{2}z}\,,$ (66)

and the one-loop term $\mathcal{I}_{qq}^{(1)}$ is given in eq. (119). Here
$\mathcal{L}_{n}(1-z)$ denotes the standard plus distribution as defined in
appendix D.

In the $(n_{l}+1)$-flavor scheme for $\alpha_{s}$ there is also a correction
from a virtual massive quark loop to the flavor-nondiagonal matching
coefficient $\mathcal{I}^{(2)}_{qg}$. This contribution is trivial, since it
factorizes into a vacuum polarization correction corresponding to the matching
of $\alpha_{s}$ between the $(n_{l})$ and $(n_{l}+1)$-flavor schemes, and the
one-loop contribution, such that

$\displaystyle\mathcal{I}_{qg}^{(2,h)}(\vec{p}_{T},m,z,\mu)=\frac{4}{3}L_{m}\,\mathcal{I}_{qg}^{(1)}(\vec{p}_{T},z,\mu)\,,$
(67)

with $\mathcal{I}_{qg}^{(1)}$ given in eq. (118). In the $(n_{l})$-flavor
scheme for $\alpha_{s}$ the $\mathcal{I}_{qg}^{(2,h)}$ contribution vanishes.

The contributions from a massive flavor to the beam function anomalous
dimensions are

$\displaystyle\gamma_{B}^{(2,h)}\Bigl{(}\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=C_{F}\,\biggl{(}-\frac{160}{9}\ln\frac{\nu}{\omega}-\frac{4}{3}-\frac{16\pi^{2}}{9}\biggr{)}\,,$
$\displaystyle\gamma_{\nu,B}^{(2,h)}(\vec{p}_{T},m,\mu)$
$\displaystyle=C_{F}\,\biggl{\\{}-\frac{16}{3}L_{m}\mathcal{L}_{0}(\vec{p}_{T},\mu)+\delta^{(2)}(\vec{p}_{T})\Bigl{(}\frac{8}{3}L_{m}^{2}+\frac{80}{9}L_{m}+\frac{224}{27}\Bigr{)}$
$\displaystyle\quad+\frac{16}{9\pi
p_{T}^{2}}\biggl{[}5-12\hat{m}^{2}-3c(1-2\hat{m}^{2})\ln\frac{c+1}{c-1}\biggr{]}\biggr{\\}}\,.$
(68)

The $\mathcal{L}_{0}(\vec{p}_{T},\mu)$ distribution is defined in appendix D.
The $\mu$ anomalous dimension here is the same as for a massless quark flavor,
$\gamma_{B}^{(2,h)}=\gamma_{B}^{(2,l)}$ [see eq. (A.2.1)]. The rapidity
anomalous dimension is explicitly mass dependent and only reproduces the
result for a massless flavor in the limit $m\ll p_{T}$.

####  4.2.2 Virtuality-dependent beam function coefficients

The massive quark-gluon virtuality beam function matching coefficient at
$\mathcal{O}(\alpha_{s})$ shown in fig. 8 is given by

$\displaystyle\mathcal{I}_{Qg}(t,m,z)=\frac{\alpha_{s}T_{F}}{4\pi}\,\theta(t)\,\theta(z)\,\theta\biggl{[}\frac{t(1-z)}{z}-m^{2}\biggr{]}\,\frac{2}{t}\biggl{[}P_{qg}(z)+\frac{2m^{2}z^{2}}{t}\biggr{]}+\mathcal{O}(\alpha_{s}^{2})\,.$
(69)

The contributions from secondary massive quarks to the light-quark coefficient
at $\mathcal{O}(\alpha_{s}^{2})$ as shown in fig. 9 are given by

$\displaystyle\mathcal{I}_{qq}^{(2,h)}\Bigl{(}t,m,z,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle\quad=\theta(z)C_{F}\biggl{\\{}\delta(t)\,\delta(1-z)\biggl{[}\Bigl{(}\frac{8}{3}L_{m}^{2}+\frac{80}{9}L_{m}+\frac{224}{27}\Bigr{)}\ln\frac{\nu}{\omega}+2L_{m}^{2}+\Bigl{(}\frac{2}{3}+\frac{8\pi^{2}}{9}\Bigr{)}L_{m}$
$\displaystyle\qquad\quad+\frac{73}{18}+\frac{20\pi^{2}}{27}-\frac{8\zeta_{3}}{3}\biggr{]}$
$\displaystyle\qquad+\theta\biggl{(}t-\frac{4m^{2}z}{1-z}\biggr{)}\frac{8}{9t(1-z)}\biggl{[}-\frac{3}{u}\Bigl{[}(1+z^{2})(1-2\hat{m}_{t}^{2}z)-4\hat{m}_{t}^{4}z^{2}(2-3z+5z^{2})\Bigr{]}\ln\frac{u-v}{u+v}$
$\displaystyle\qquad\quad-2v\Bigl{[}4-3z+4z^{2}+\frac{z(11-21z+29z^{2}-15z^{3})}{1-z}\,\hat{m}_{t}^{2}\Bigr{]}\biggr{]}\biggr{\\}}+\frac{4}{3}L_{m}\,\mathcal{I}_{qq}^{(1)}(t,z,\mu)\,,$
(70)

with

$\displaystyle\hat{m}_{t}=\frac{m}{\sqrt{t}}\,\qquad
u=\sqrt{1-4\hat{m}_{t}^{2}z}\,,\qquad
v=\sqrt{1-\frac{4\hat{m}_{t}^{2}z}{1-z}}\,,$ (71)

and the one-loop term $\mathcal{I}_{qq}^{(1)}$ is given in eq. (123).

In the $(n_{l}+1)$-flavor scheme for $\alpha_{s}$ there is also the analogous
contribution to eq. (67) to the flavor-nondiagonal coefficient

$\displaystyle\mathcal{I}_{qg}^{(2,h)}(t,m,z,\mu)=\frac{4}{3}L_{m}\,\mathcal{I}_{qg}^{(1)}(t,z,\mu)\,,$
(72)

with $\mathcal{I}_{qg}^{(1)}$ given in eq. (123). In the $(n_{l})$-flavor
scheme for $\alpha_{s}$ the $\mathcal{I}_{qg}^{(2,h)}$ contribution vanishes.

The contribution from the massive flavor to the $\mu$ anomalous dimension at
$\mathcal{O}(\alpha_{s}^{2})$ is given by

$\displaystyle\gamma^{(2,h)}_{B,m}\Bigl{(}t,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle=C_{F}\,\delta(t)\biggl{(}-\frac{160}{9}\ln\frac{\nu}{\omega}-\frac{4}{3}-\frac{16\pi^{2}}{9}\biggr{)}\,.$
(73)

We emphasize that the massive quark contribution to the $\mu$ anomalous
dimension is not the same as for a massless flavor, but is in fact the same as
for the TMD beam function in eq. (4.2.1). This is required by consistency with
the large mass limit $Q{\mathcal{T}},q_{T}\ll m$, where the massive flavor can
only contribute to the (local) running of the common current operators, which
are independent of the measurement. Only in combination with the soft mass-
mode function $H_{s}$ and the soft function, the combined $\mu$ evolution
above the mass scale is the same as for $n_{l}+1$ massless flavors as
discussed in eq. (3.2).

The secondary massive quarks introduce rapidity divergences and associated
logarithms also in the virtuality-dependent beam function. The $\nu$ anomalous
dimension induced by the secondary massive effects is the same as for the
collinear mass-mode function, see eq. (34), given in eq. (4.1.2).

####  4.2.3 PDF matching coefficients

The matching coefficients relating the PDFs in the $(n_{l}+1)$ and the
$(n_{l})$-flavor scheme are all known at two loops Buza:1996wv  and partially
beyond (see e.g. refs. Ablinger:2014nga ; Ablinger:2014vwa ; Ablinger:2014lka
and references therein). The matching coefficient for a primary massive quark
originating from an initial-state gluon at $\mathcal{O}(\alpha_{s})$ is

$\displaystyle\mathcal{M}_{Qg}(m,z,\mu)=-\frac{\alpha_{s}T_{F}}{4\pi}\,\theta(1-z)\theta(z)\,2P_{qg}(z)L_{m}+\mathcal{O}(\alpha_{s}^{2})\,.$
(74)

The matching coefficient coming from secondary massive quark corrections to
the light-quark PDFs reads up to $\mathcal{O}(\alpha_{s}^{2})$

$\displaystyle\mathcal{M}_{qq}(m,z,\mu)$
$\displaystyle=1+\frac{\alpha_{s}^{2}C_{F}T_{F}}{16\pi^{2}}\,\theta(z)\biggl{\\{}\mathcal{L}_{0}(1-z)\biggl{(}\frac{8}{3}L_{m}^{2}+\frac{80}{9}L_{m}+\frac{224}{27}\biggr{)}$
$\displaystyle\qquad+\delta(1-z)\biggl{[}2L_{m}^{2}+\left(\frac{2}{3}+\frac{8\pi^{2}}{9}\right)L_{m}+\frac{73}{18}+\frac{20\pi^{2}}{27}-\frac{8\zeta_{3}}{3}\biggr{]}$
$\displaystyle\qquad+\theta(1-z)\
\biggl{[}-\frac{4}{3}L_{m}^{2}(1+z)+L_{m}\biggl{(}\frac{8}{9}-\frac{88}{9}z+\frac{8}{3}\frac{1+z^{2}}{1-z}\,\ln
z\biggr{)}+\frac{2}{3}\frac{1+z^{2}}{1-z}\,\ln^{2}z$
$\displaystyle\qquad+\frac{\ln
z}{1-z}\biggl{(}\frac{44}{9}-\frac{16}{3}z+\frac{44}{9}z^{2}\biggr{)}+\frac{44}{27}-\frac{268}{27}z\biggr{]}\biggr{\\}}+\mathcal{O}(\alpha_{s}^{3})\,.$
(75)

The matching coefficient between the gluon PDF in the $(n_{l})$ and
$(n_{l}+1)$-flavor schemes at $\mathcal{O}(\alpha_{s})$, which is also
required for Drell-Yan at $\mathcal{O}(\alpha_{s}^{2})$, is equivalent to the
matching relation for $\alpha_{s}$

$\displaystyle\mathcal{M}_{gg}(m,z,\mu)=\delta(1-z)+\frac{\alpha_{s}T_{F}}{4\pi}\,\delta(1-z)\,\frac{4}{3}L_{m}+\mathcal{O}(\alpha_{s}^{2})\,.$
(76)

Note that taking into account the nondiagonal evolution of the PDFs the known
$\mathcal{O}(\alpha_{s}^{2})$ corrections for all matching factors
$\mathcal{M}_{ij}$ become relevant at NNLL′.

###  4.3 Soft and collinear-soft functions

Here we give all massive quark corrections at $\mathcal{O}(\alpha_{s}^{2})$ to
the soft and csoft functions. They arise exclusively from secondary radiation.
Note that the soft functions satisfy Casimir scaling at this order and can be
thus applied also to color-singlet production in gluon-fusion by replacing an
overall $C_{F}\to C_{A}$.

####  4.3.1 TMD soft function

Figure 10: Corrections from secondary massive quarks to the (c)soft function.
Also the mirror diagrams need to be included.

The contributions from secondary massive quarks to the TMD soft function,
which appears in eq. (2.4) for $q_{T}\sim m$, are calculated in app. B.4 at
$\mathcal{O}(\alpha_{s}^{2})$ and correspond to the diagrams shown in fig. 10.
The result reads

$\displaystyle S^{(2,h)}(\vec{p}_{T},m,\mu,\nu)$
$\displaystyle=C_{F}\,\biggl{\\{}\delta^{(2)}(\vec{p}_{T})\biggl{[}\Bigl{(}-\frac{16}{3}L_{m}^{2}-\frac{160}{9}L_{m}-\frac{448}{27}\Bigr{)}\ln\frac{\nu}{\mu}+\frac{8}{9}L_{m}^{3}+\frac{40}{9}L_{m}^{2}$
$\displaystyle\qquad+\Bigl{(}\frac{448}{27}-\frac{4\pi^{2}}{9}\Bigr{)}L_{m}+\frac{656}{27}-\frac{10\pi^{2}}{27}-\frac{56\zeta_{3}}{9}\biggr{]}$
$\displaystyle\quad+\frac{16}{9\pi
p_{T}^{2}}\biggl{[}2\Bigl{[}-5+12\hat{m}^{2}+3c(1-2\hat{m}^{2})\ln\frac{c+1}{c-1}\Bigr{]}\ln\frac{\nu}{m}$
$\displaystyle\qquad+3c(1-2\hat{m}^{2})\biggl{[}\textrm{Li}_{2}\biggl{(}\frac{(c-1)^{2}}{(c+1)^{2}}\biggr{)}+2\ln\frac{c+1}{c-1}\ln\frac{\hat{m}(c+1)}{2c}-\frac{\pi^{2}}{6}\biggr{]}$
$\displaystyle\qquad+c(5-16\hat{m}^{2})\ln\frac{c+1}{c-1}+8\hat{m}^{2}\biggr{]}\biggr{\\}}+\frac{4}{3}L_{m}\,S^{(1)}(\vec{p}_{T},\mu,\nu)\,,$
(77)

where $\hat{m}=m/p_{T}$ and $c=\sqrt{1+4\hat{m}^{2}}$ as in eq. (66) and the
one-loop soft function $S^{(1)}$ given in eq. (126).

The massive quark contributions to the anomalous dimensions of the soft
function are

$\displaystyle\gamma_{S}^{(2,h)}(\vec{p}_{T},\mu,\nu)$
$\displaystyle=C_{F}\,\biggl{(}\frac{320}{9}\ln\frac{\nu}{\mu}-\frac{448}{27}+\frac{8\pi^{2}}{9}\biggr{)}\,,$
$\displaystyle\gamma_{\nu,S}^{(2,h)}(\vec{p}_{T},m,\mu)$
$\displaystyle=C_{F}\,\biggl{\\{}\frac{32}{3}L_{m}\mathcal{L}_{0}(\vec{p}_{T},\mu)+\delta^{(2)}(\vec{p}_{T})\biggl{(}-\frac{16}{3}L_{m}^{2}-\frac{160}{9}L_{m}-\frac{448}{27}\biggr{)}$
$\displaystyle\quad+\frac{32}{9\pi
p_{T}^{2}}\biggl{[}-5+12\hat{m}^{2}+3c(1-2\hat{m}^{2})\ln\frac{c+1}{c-1}\biggr{]}\biggr{\\}}\,.$
(78)

The $\mu$ anomalous dimension here is the same as for an additional massless
flavor, $\gamma_{S}^{(2,h)}=\gamma_{S}^{(2,l)}$ [see eq. (A.3.1)]. The
rapidity anomalous dimension is explicitly mass dependent and only reduces to
the result for a massless flavor in the limit $m\ll p_{T}$.

####  4.3.2 Csoft function for beam thrust

The csoft function is a matching coefficient between an eikonal matrix element
in the $n_{l}+1$ and $n_{l}$ flavor theories appearing for the hierarchy
${\mathcal{T}}\ll m\ll\sqrt{Q{\mathcal{T}}}$ in eq. (3.3). The relevant
diagrams at $\mathcal{O}(\alpha_{s}^{2})$ are shown in fig. 10 and are
calculated in sec. B.5. The result is given by

$\displaystyle\mathcal{S}_{c}(\ell,m,\mu,\nu)$
$\displaystyle=\delta(\ell)+\frac{\alpha_{s}^{2}C_{F}T_{F}}{16\pi^{2}}\biggl{\\{}\frac{\nu}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{\ell\,\nu}{\mu^{2}}\Bigr{)}\Bigl{(}\frac{8}{3}L_{m}^{2}+\frac{80}{9}L_{m}+\frac{224}{27}\Bigr{)}$
(79)
$\displaystyle\quad+\delta(\ell)\biggl{[}-\frac{8}{9}L_{m}^{3}-\frac{40}{9}L_{m}^{2}+\Bigl{(}-\frac{448}{27}+\frac{4\pi^{2}}{9}\Bigr{)}L_{m}-\frac{656}{27}+\frac{10\pi^{2}}{27}+\frac{56\zeta_{3}}{9}\biggr{]}\biggr{\\}}+\mathcal{O}(\alpha_{s}^{3}).$

We can see that with the scale choices $\mu\sim m$ and
$\nu\sim\mu^{2}/\ell\sim m^{2}/{\mathcal{T}}$ all large logarithms (including
the implicit one inside the plus distribution) are minimized. The $\mu$
anomalous dimensions of the csoft matching function is given by

$\displaystyle\gamma_{\mathcal{S}_{c}}(\ell,m,\mu,\nu)$
$\displaystyle=\frac{\alpha_{s}^{2}C_{F}T_{F}}{16\pi^{2}}\biggl{[}-\frac{\nu}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{\ell\,\nu}{\mu^{2}}\Bigr{)}\Bigl{(}\frac{32}{3}L_{m}+\frac{160}{9}\Bigr{)}+\delta(\ell)\Bigl{(}\frac{448}{27}-\frac{8\pi^{2}}{9}\Bigr{)}\biggr{]}+\mathcal{O}(\alpha_{s}^{3})\,.$
(80)

The $\nu$ anomalous dimension is the same as for the collinear mass mode
function in eq. (4.1.2), $\gamma_{\nu,\mathcal{S}_{c}}=\gamma_{\nu,H_{c}}$.

####  4.3.3 (Beam) thrust soft function

The secondary massive quark corrections to the (beam) thrust soft function at
$\mathcal{O}(\alpha_{s}^{2})$ were calculated in ref. Gritschacher:2013tza
and are given by

$\displaystyle S^{(2,h)}(\ell,m,\mu)$
$\displaystyle=C_{F}\,\biggl{\\{}\frac{1}{\mu}\mathcal{L}_{0}\Bigl{(}\frac{\ell}{\mu}\Bigr{)}\biggl{(}\frac{16}{3}L_{m}^{2}+\frac{160}{9}L_{m}+\frac{448}{27}\biggr{)}$
$\displaystyle\quad+\delta(\ell)\biggl{[}-\frac{8}{9}L_{m}^{3}-\frac{40}{9}L_{m}^{2}+\Bigl{(}-\frac{448}{27}+\frac{4\pi^{2}}{9}\Bigr{)}L_{m}-\frac{656}{27}+\frac{10\pi^{2}}{27}+\frac{56}{9}\zeta_{3}\biggr{]}$
$\displaystyle\quad+\theta(\ell-2m)\frac{1}{\ell}\biggl{[}\frac{64}{3}\,\textrm{Li}_{2}\Bigl{(}\frac{w-1}{w+1}\Bigr{)}+\frac{16}{3}\ln^{2}\frac{1-w}{1+w}-\frac{64}{3}\ln\frac{1-w}{1+w}\ln\hat{m}_{\ell}$
$\displaystyle\qquad-\frac{160}{9}\ln\frac{1-w}{1+w}-w\Bigl{(}\frac{896}{27}+\frac{256}{27}\hat{m}_{\ell}^{2}\Bigr{)}+\frac{16\pi^{2}}{9}\biggr{]}\biggr{\\}}$
$\displaystyle\quad+\Delta
S_{\tau}(\ell,m)+\frac{4}{3}L_{m}\,S^{(1)}(\ell,\mu)\,,$ (81)

where

$\displaystyle\hat{m}_{\ell}\equiv\frac{m}{\ell}\,,\qquad
w=\sqrt{1-4\hat{m}_{\ell}^{2}}\,,$ (82)

and the one-loop soft function $S^{(1)}$ is given in eq. (129). The term
$\Delta S_{\tau}(\ell,m)$ contains the correction from two real final-state
emissions entering two opposite hemispheres, which vanishes both for $\ell\ll
m$ and $m\ll\ell$ and is currently only known numerically. The integral
expression for this numerically small contribution is given in eq. (61) of
ref. Gritschacher:2013tza , and a precise parametrization can be found in ref.
Pietrulewicz:2014qza .

The massive quark contribution to the anomalous dimension is the same as for a
massless flavor, $\gamma_{S}^{(2,h)}(\ell,\mu)=\gamma_{S}^{(2,l)}(\ell,\mu)$,
given in eq. (A.3.2).

###  4.4 Small and large mass limits

In secs. 2.6 and 3.5 we explained how the ingredients in the factorization
theorems for different hierarchies are related to each other. Here we verify
these relations for the beam and soft functions up to
$\mathcal{O}(\alpha_{s}^{2})$. We also scrutinize the numerical impact of the
power corrections for these functions. We focus in particular on the
$\mathcal{O}(m^{2}/q_{T}^{2})$ corrections the $q_{T}$ spectrum for $b$
quarks, which are contained in the factorization theorem eq. (2.4) for
$q_{T}\sim m$ but not in the massless limit for $m\ll q_{T}$ in eq. (2.5), as
these are phenomenologically important hierarchies for $b$-quark mass effects
at the LHC.

For the numerical results we use the MMHT2014 NNLO PDFs Harland-Lang:2014zoa
and evaluate the contributions for $\mu=m_{b}=4.8$ GeV, $\omega=m_{Z}$, and
$E_{\rm cm}=13$ TeV. The main qualitative features of the results do not
depend on these specific input parameters.

####  4.4.1 Limiting behavior for $q_{T}$

We first consider the primary mass effects at one loop, which are encoded in
the TMD beam function matching coefficient $\mathcal{I}_{Qg}^{(1)}$ in eq.
(63). In the limit $p_{T}\ll m$ the primary massive quarks decouple, which is
manifest in the result,

$\displaystyle\mathcal{I}_{Qg}^{(1)}(\vec{p}_{T},m,z)=\mathcal{O}\Bigl{(}\frac{p_{T}^{2}}{m^{2}}\Bigr{)}\,.$
(83)

On the other hand, in the opposite limit $m\ll p_{T}$ it becomes

$\displaystyle\mathcal{I}_{Qg}^{(1)}(\vec{p}_{T},m,z)$
$\displaystyle\quad=T_{F}\,\theta(1-z)\theta(z)\,\Bigl{\\{}2P_{qg}(z)\mathcal{L}_{0}(\vec{p}_{T},\mu)+\delta^{(2)}(\vec{p}_{T})\Bigl{[}-2P_{qg}(z)L_{m}+4z(1-z)\Bigr{]}\\!+\mathcal{O}\Bigl{(}\frac{m^{2}}{p_{T}^{2}}\Bigr{)}\biggr{\\}}$
$\displaystyle\quad=\mathcal{I}_{qg}^{(1)}(\vec{p}_{T},z,\mu)+\delta^{(2)}(\vec{p}_{T})\mathcal{M}_{Qg}^{(1)}(m,z,\mu)+\mathcal{O}\Bigl{(}\frac{m^{2}}{p_{T}^{2}}\Bigr{)}\,,$
(84)

confirming that the relation in eq. (25) is satisfied at
$\mathcal{O}(\alpha_{s})$. The massless one-loop matching coefficient
$\mathcal{I}_{qg}^{(1)}$ can be found in eq. (118) and the PDF matching
coefficient $\mathcal{M}_{Qg}^{(1)}$ in eq. (74). To account for the correct
distributive structure in $\vec{p}_{T}$ that emerges in the massless limit,
one can integrate the expressions with massive quarks and identify the
distributions at the cumulant level.

Figure 11: Massive $b$-quark beam function (left panel) and its cumulant
(right panel) together with its $m\to 0$ limit. The input parameters are
described in the text.

In fig. 11 we show the result for the massive quark beam function
$B^{(1)}_{Q}=\mathcal{I}^{(1)}_{Qg}\otimes_{x}f_{g}$ at
$\mathcal{O}(\alpha_{s})$ as function of $p_{T}$ using the full massive
matching coefficient $\mathcal{I}^{(1)}_{Qg}$ (solid orange) and its small
mass limit in eq. (4.4.1). Note that the results differential in $p_{T}$ are
not explicitly $\mu$-dependent at $\mathcal{O}(\alpha_{s})$. In the right
panel we show the corresponding results for the cumulant

$\displaystyle B_{Q}(p_{T}^{\rm cut},m)\equiv\int_{|\vec{p}_{T}|<p_{T}^{\rm
cut}}\textrm{d}^{2}p_{T}\,B_{Q}(\vec{p}_{T},m)\,,$ (85)

which also includes the $\delta^{(2)}(\vec{p}_{T})$ constant contribution. We
can see that in both cases the small mass limit is correctly approached for
$p_{T}^{(\rm cut)}\gg m_{b}$, while for $p_{T}^{(\rm cut)}\ll m_{b}$ the
primary mass effects decouple with the result going to zero. The corrections
to the small mass limit become sizeable for $p_{T}\sim m_{b}$ and vanish quite
fast for larger $p_{T}$.

Figure 12: Convolution of two massive $b$-quark beam functions together with
the result in the $m\to 0$ limit differential in the total
$p_{T}\equiv|\vec{p}_{T}|$ (left panel) and the corresponding cumulant (right
panel). This is proportional to the primary massive quark correction to the
$Z$-boson spectrum at $\mathcal{O}(\alpha_{s}^{2})$.

In fig. 12 we show the result for the convolution between two massive quark
beam functions,

$\displaystyle B^{(1)}_{Q}\otimes
B^{(1)}_{Q}(\vec{p}_{T},m)\equiv\int\\!\textrm{d}^{2}p_{T}^{\prime}\,B^{(1)}_{Q}(\vec{p}_{T}-\vec{p}^{\,\prime}_{T},m)\,B^{(1)}_{Q}(\vec{p}^{\,\prime}_{T},m)\,,$
(86)

which enters the result for $Z$-boson production at
$\mathcal{O}(\alpha_{s}^{2}T_{F}^{2})$ and NNLL′. The analytic expression for
the convolution between the two one-loop mass-dependent coefficients is given
in eq. (C.1). We see that now the corrections to small-mass limit remain
nonnegligible even for larger values of $p_{T}$. This is due to the fact that
the $\vec{p}_{T}$-convolution generates a logarithmic dependence in the
spectrum, such that the power corrections of
$\mathcal{O}(m_{b}^{2}/p_{T}^{2})$ become enhanced by logarithms
$\ln(p_{T}^{2}/m_{b}^{2})$.

Next, we consider the secondary massive quark corrections at
$\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$. The result for the mass-dependent TMD
beam function coefficient $\mathcal{I}^{(2,h)}_{qq}(\vec{p}_{T},m,z)$ is given
in eq. (4.2.1). In the decoupling limit $p_{T}\ll m$ all its terms without
distributions in $\vec{p}_{T}$ give $\mathcal{O}(p_{T}^{2}/m^{2})$ power-
suppressed contributions. Combining its remaining distributional terms with
the contributions arising from changing the $\alpha_{s}$ scheme from $n_{l}+1$
to $n_{l}$ flavors yields

$\displaystyle\mathcal{I}_{qq}^{(2,h)}\Bigl{(}\vec{p}_{T},m,z,\mu,\frac{\nu}{\omega}\Bigr{)}-\frac{4}{3}L_{m}\mathcal{I}_{qq}^{(1)}\Bigl{(}\vec{p}_{T},z,\mu,\frac{\nu}{\omega}\Bigr{)}=\delta^{(2)}(\vec{p}_{T})\,\delta(1-z)\,H_{c}^{(2)}\Bigl{(}m,\mu,\frac{\nu}{\omega}\Bigr{)}+\mathcal{O}\Bigl{(}\frac{p_{T}^{2}}{m^{2}}\Bigr{)}\,,$
(87)

confirming the relation in eq. (23) at this order. The massless one-loop
coefficient $\mathcal{I}_{qq}^{(1)}$ and the collinear mass-mode function
$H_{c}^{(2)}$ can be found in eqs. (119) and (4.1.2), respectively. On the
other hand, in the limit $m\ll p_{T}$ we get

$\displaystyle\mathcal{I}_{qq}^{(2,h)}\Bigl{(}\vec{p}_{T},m,z,\mu,\frac{\nu}{\omega}\Bigr{)}=\mathcal{I}^{(2,l)}_{qq}\Bigl{(}\vec{p}_{T},z,\mu,\frac{\nu}{\omega}\Bigr{)}+\delta^{2}(\vec{p}_{T})\,\mathcal{M}_{qq}^{(2)}(m,z,\mu)+\mathcal{O}\Bigl{(}\frac{m^{2}}{p_{T}^{2}}\Bigr{)}\;,$
(88)

such that all infrared mass dependence is given by the PDF matching, as
required by the relation in eq. (25). The results for the massless coefficient
and the PDF matching coefficient are given in eqs. (A.2.1) and (4.2.3),
respectively.

For the coefficient $\mathcal{I}_{qg}^{(2,h)}$ at
$\mathcal{O}(\alpha_{s}^{2}T_{F}^{2})$ the limiting behavior is trivial, since
it vanishes identically in the $(n_{l})$-flavor scheme for $\alpha_{s}$, and
in the $(n_{l}+1)$-flavor scheme for $\alpha_{s}$ it is exactly

$\displaystyle
T_{F}\,\mathcal{I}_{qg}^{(2,h)}(\vec{p}_{T},m,z,\mu)=\mathcal{I}_{qg}^{(1)}(\vec{p}_{T},z,\mu)\otimes_{z}\mathcal{M}^{(1)}_{gg}(m,z,\mu)\,.$
(89)

The mass-dependent TMD soft function is given in eq. (4.3.1). In the limit
$p_{T}\ll m$ all its terms without distributions in $\vec{p}_{T}$ become
$\mathcal{O}(p_{T}^{2}/m^{2})$ power suppressed, just as for the beam
function. Combining its remaining distributional terms with the contributions
arising from changing the scheme of the strong coupling from $n_{l}+1$ to
$n_{l}$ flavors yields

$\displaystyle
S^{(2,h)}(\vec{p}_{T},m,\mu,\nu)-\frac{4}{3}L_{m}S^{(1)}(\vec{p}_{T},\mu,\nu)=\delta^{(2)}(\vec{p}_{T})\,H_{s}^{(2)}(m,\mu,\nu)+\mathcal{O}\Bigl{(}\frac{p_{T}^{2}}{m^{2}}\Bigr{)}\,,$
(90)

confirming the relation in eq. (24). The massless one-loop TMD soft function
$S^{(1)}$ and the softmass-mode function $H_{s}^{(2)}$ are given in eqs. (126)
and (4.1.2), respectively. Since the soft function is free of IR
singularities, the limit $m\ll p_{T}$ just yields the massless soft function
in eq. (A.3.1),

$\displaystyle
S^{(2,h)}(\vec{p}_{T},m,\mu,\nu)=S^{(2,l)}(\vec{p}_{T},\mu,\nu)+\mathcal{O}\Bigl{(}\frac{m^{2}}{p_{T}^{2}}\Bigr{)}\,.$
(91)

Figure 13: Secondary massive bottom quark corrections to the $u$-quark beam
function (left panel) and its cumulant (right panel) at
$\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ (including the square-root of the soft
function here) for $\mu=m_{b}$ as a function of $q_{T}$.

We now discuss the numerical impact of the $\mathcal{O}(m^{2}/p_{T}^{2})$
terms from secondary mass effects. Since the individual results for the beam
and soft functions depend on the specific regularization scheme, we consider
their symmetrized combination

$\displaystyle\tilde{B}_{q}(\vec{p}_{T},m,\omega,x,\mu)=\int\textrm{d}^{2}p_{T}^{\prime}\,B_{q}\Bigl{(}\vec{p}_{T}-\vec{p}^{\,\prime}_{T},m,x,\mu,\frac{\nu}{\omega}\Bigr{)}\,\sqrt{S(\vec{p}^{\,\prime}_{T},m,\mu,\nu)}\,,$
(92)

which is independent of $\nu$.666This combination is sometimes used as
definition of a TMD PDF. $\tilde{B}_{q}$ contains large rapidity logarithms,
which are resummed once the soft and beam functions are evaluated at their
natural rapidity scales and evolved to a common scale $\nu$. For demonstrating
the size of the power corrections here, we evaluate it at fixed order. The
$\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ corrections explicitly depend on $\mu$
and the flavor-number scheme, but the difference between the full result and
the small mass limits given in eqs. (88) and (91) do not. In fig. 13 we show
the result for the $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ corrections (with
$\alpha_{s}=\alpha_{s}^{(n_{l}+1)}$) to the $u$-quark beam function, both
differential in $p_{T}$ and the corresponding cumulant. We see that the full
mass dependent results correctly reproduce the small and large mass limits.
The corrections to the massless are much larger than for the primary mass
effects. In particular, they are still of $\mathcal{O}(100\%)$ for
$p_{T}^{(\rm cut)}\sim 10$ GeV. This clearly indicates that for secondary
radiation involving two massive quarks in the final state the corrections are
rather of $\mathcal{O}(4m^{2}/p_{T}^{2})$, as one might expect.

####  4.4.2 Limiting behavior for ${\mathcal{T}}$

We carry out the discussion for beam thrust in close analogy. The virtuality-
dependent massive quark beam function coefficient at one loop is given in eq.
(69). In the limit $t\ll m^{2}$ the primary massive quarks correctly decouple,

$\displaystyle\mathcal{I}_{Qg}^{(1)}(t,m,z)=\mathcal{O}\Bigl{(}\frac{t}{m^{2}}\Bigr{)}\,.$
(93)

In the opposite limit $m^{2}\ll t$ we get

$\displaystyle\mathcal{I}_{Qg}^{(1)}(t,m,z)$
$\displaystyle=T_{F}\,\theta(1-z)\theta(z)\biggl{\\{}2P_{qg}(z)\,\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}+\delta(t)\biggl{[}2P_{qg}(z)\biggl{(}-L_{m}+\ln\frac{1-z}{z}\biggr{)}$
$\displaystyle\qquad+4z(1-z)\biggr{]}\biggr{\\}}+\mathcal{O}\left(\frac{m^{2}}{t}\right)$
$\displaystyle=\mathcal{I}_{qg}^{(1)}(t,z,\mu)+\delta(t)\,\mathcal{M}_{Qg}^{(1)}(m,z,\mu)+\mathcal{O}\Bigl{(}\frac{m^{2}}{t}\Bigr{)}\;,$
(94)

as required by the relation (3.5). The massless one-loop matching coefficient
$\mathcal{I}_{qg}^{(1)}$ and the PDF matching coefficient
$\mathcal{M}_{Qg}^{(1)}$ are given in eqs. (123) and (74), respectively.

The secondary massive quark corrections to the virtuality-dependent beam
function are given in eq. (4.2.2). In the decoupling limit $t\ll m^{2}$ all
its nondistributional terms become $\mathcal{O}(t/m^{2})$ power suppressed.
Combining the remaining distributional terms in $t$ with the contributions
arising from changing the scheme of the strong coupling from $n_{l}+1$ to
$n_{l}$ flavors yields

$\displaystyle\mathcal{I}_{qq}^{(2,h)}\Bigl{(}t,m,z,\mu,\frac{\nu}{\omega}\Bigr{)}-\frac{4}{3}L_{m}\mathcal{I}_{qq}^{(1)}(t,z,\mu)=\delta(t)\delta(1-z)\,H_{c}^{(2)}\Bigl{(}m,\mu,\frac{\nu}{\omega}\Bigr{)}+\mathcal{O}\Bigl{(}\frac{t}{m^{2}}\Bigr{)}\,,$
(95)

in agreement with eq. (45). The massless result for $\mathcal{I}_{qq}^{(1)}$
and the collinear mass-mode function $H_{c}^{(2)}$ are given in eqs. (123) and
(4.1.2), respectively. In the limit $m^{2}\ll t$ we get

$\displaystyle\mathcal{I}_{qq}^{(2,h)}\Bigl{(}t,m,z,\mu,\frac{\nu}{\omega}\Bigr{)}$
$\displaystyle\quad=\mathcal{I}^{(2,l)}_{qq}(t,z,\mu)+\delta(t)\,\mathcal{M}_{qq}^{(2)}(m,z,\mu)+\delta(1-z)\,\frac{1}{\omega}\,\mathcal{S}_{c}^{(2)}\Bigl{(}\frac{t}{\omega},m,\mu,\nu\Bigr{)}+\mathcal{O}\Bigl{(}\frac{m^{2}}{t}\Bigr{)}\,.$
(96)

All infrared mass dependence is contained in the PDF matching coefficient and
the csoft function, as required by eq. (3.5). The functions on the right-hand
side are given in eqs. (124), (4.2.3), and (79), respectively.

Figure 14: Massive $b$-quark beam function (left panel) and the convolution
between two of these (right panel) together with the $m\to 0$ limit as a
function of $\sqrt{t}\sim\sqrt{Q{\mathcal{T}}}$.

Figure 15: Secondary massive $b$-quark corrections to the $u$-quark beam
function for $Y=0$ (left panel) and the soft function (right panel) at
$\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ for $\mu=m_{b}$ as functions of
$\sqrt{t}\sim\sqrt{Q{\mathcal{T}}}$ and $\ell\sim{\mathcal{T}}$, respectively.

The mass-dependent corrections to the (beam) thrust soft function are given in
eq. (4.3.3). In the limit $\ell\ll m$ all its nondistributional terms become
$\mathcal{O}(\ell^{2}/m^{2})$ power suppressed. Combining the remaining
distributional terms with the contributions arising from changing the scheme
of the strong coupling from $n_{l}+1$ to $n_{l}$ flavors yields

$\displaystyle
S^{(2,h)}(\ell,m,\mu)-\frac{4}{3}L_{m}S^{(1)}(\ell,\mu)=\delta(\ell)\,H_{s}^{(2)}(m,\mu,\nu)+\mathcal{S}_{c}^{(2)}(\ell,m,\mu,\nu)+\mathcal{O}\Bigl{(}\frac{\ell^{2}}{m^{2}}\Bigr{)}\,,$
(97)

in agreement with eq. (47). The massless one-loop thrust soft function
$S^{(1)}$, the soft mass-mode function $H_{s}^{(2)}$, and the csoft function
$\mathcal{S}_{c}^{(2)}$ can be found in eqs. (129), (4.1.2), and (79),
respectively. For $m\ll\ell$ the correct massless result is recovered,

$\displaystyle
S^{(2,h)}(\ell,m,\mu)=S^{(2,l)}(\ell,\mu)+\mathcal{O}\Bigl{(}\frac{m^{2}}{\ell^{2}}\Bigr{)}\,,$
(98)

which was already checked in ref. Gritschacher:2013tza .

In fig. 14, we show the numerical results for the one-loop massive beam
function and the convolution between two of these (which is the leading order
correction from primary massive quarks for the $Z$-boson production) as a
function of $\sqrt{t}\sim\sqrt{Q{\mathcal{T}}}$. The mass effects become
relevant for $\sqrt{t}\sim m_{b}\sim 5$ GeV (corresponding to
${\mathcal{T}}\lesssim 1$ GeV for $Q=m_{Z}$). The corrections to the massless
limit for the convolution of two beam functions is nonnegligible also for
larger values. In fig. 15, we show the result for the secondary
$\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ corrections to the beam and soft
function. The corrections to the massless limit for the beam function remain
sizeable even for $\sqrt{t}\gtrsim 2m_{b}$. For the soft function, the mass
effects are important for ${\mathcal{T}}\sim\ell\sim m_{b}$ and become small
for $\ell>10\,{\rm GeV}\sim 2m_{b}$. Note that the small bump in the soft
function in fig. 15 originates from the correction term $\Delta S_{\tau}$ in
eq. (4.3.3). The associated correction in the massless limit is fully
contained in the $\delta(\ell)$ term.

##  5 Rapidity evolution

Here, we discuss the solutions of the rapidity RGEs in eq. (13), or
equivalently eqs. (34) and (40), and in particular the rapidity evolution for
the mass-dependent soft function in eq. (2.4) for $q_{T}\sim m$, where the
massive quark corrections give rise to a different running than for massless
flavors. Our primary aim here is to highlight the different features with
respect to the massless case, while leaving the practical implementation for
future work.

The rapidity evolution for the mass-mode matching functions $H_{s}$ and
$H_{c}$ according to eq. (13) has been discussed in ref. Hoang:2015vua . The
evolution for the beam thrust beam function and csoft function according to
eqs. (34) and (40) is completely analogous. For example, the $\nu$-evolved
soft matching function $H_{s}$ is given by

$\displaystyle H_{s}(m,\mu,\nu)$
$\displaystyle=V(m,\mu,\nu,\nu_{0})\,H_{s}(m,\mu,\nu_{0})\,,$ (99)
$\displaystyle V(m,\mu,\nu,\nu_{0})$
$\displaystyle=\exp\biggl{\\{}\Big{[}4\eta_{\Gamma}^{(n_{l})}(\mu_{0}(m),\mu)-4\eta_{\Gamma}^{(n_{l}+1)}(\mu_{0}(m),\mu)+\gamma_{\nu,H_{s}}(m,\mu_{0}(m))\Big{]}\ln\frac{\nu}{\nu_{0}}\biggr{\\}}\,.$

The evolution function $\eta_{\Gamma}$ is defined by

$\displaystyle\eta_{\Gamma}^{(n_{f})}(\mu_{0},\mu)$
$\displaystyle=\int_{\mu_{0}}^{\mu}\frac{\textrm{d}\mu^{\prime}}{\mu^{\prime}}\,\Gamma^{(n_{f})}_{\rm{cusp}}[\alpha_{s}^{(n_{f})}(\mu^{\prime})]\,,$
(100)

and resums the $\mu$-dependent logarithms inside the $\nu$ anomalous dimension
as required by consistency with the $\mu$ evolution to maintain the path
independence in $\mu$-$\nu$-space Chiu:2012ir . With the canonical scale
choice

$\displaystyle\mu_{0}(m)=m\,,$ (101)

all logarithmic terms in the boundary condition
$\gamma_{\nu,H_{s}}(m,\mu_{0}(m))$ are minimized.

The solution of the rapidity RGE for the soft function is substantially more
involved due to its two-dimensional convolution structure on $\vec{p}_{T}$.
The formal solution of the rapidity RGE for massless quarks in eq. (2.1) is
most conveniently found by Fourier transforming to impact parameter space with
$b=|\vec{b}|$, where the rapidity RGE becomes multiplicative

$\displaystyle\nu\frac{\textrm{d}}{\textrm{d}\nu}\,\tilde{S}^{(n_{f})}(b,\mu,\nu)=\tilde{\gamma}^{(n_{f})}_{\nu,S}(b,\mu)\,\tilde{S}^{(n_{f})}(b,\mu,\nu)\,.$
(102)

The consistency (path independence) between $\mu$ and $\nu$ evolution requires
the rapidity anomalous dimension in Fourier space to satisfy

$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}\,\tilde{\gamma}^{(n_{f})}_{\nu,S}(b,\mu)=-4\,\Gamma_{\rm{cusp}}[\alpha_{s}^{(n_{f})}(\mu)]\,.$
(103)

Its solution is given by

$\displaystyle\tilde{\gamma}^{(n_{f})}_{\nu,S}(b,\mu)=-4\eta_{\Gamma}^{(n_{f})}(\mu_{0}(b),\mu)+\tilde{\gamma}^{(n_{f})}_{\nu,S}(b,\mu_{0}(b))\,.$
(104)

The logarithms of $\ln(\mu\,b\,e^{\gamma_{E}}/2)$ in the second boundary term
are eliminated by the canonical scale choice

$\displaystyle\mu^{(l)}_{0}(b)=\frac{2\,e^{-\gamma_{E}}}{b}\,.$ (105)

With this choice, the $\nu$ evolution of the soft function in Fourier space at
any given scale $\mu$ is given by

$\displaystyle\tilde{S}(b,\mu,\nu)=\tilde{S}(b,\mu,\nu_{0})\exp\biggl{[}\tilde{\gamma}^{(n_{f})}_{\nu,S}(b,\mu)\,\ln\frac{\nu}{\nu_{0}}\biggr{]}\,.$
(106)

As is well known, the rapidity evolution kernel becomes intrinsically
nonperturbative at $1/b\ll\Lambda_{\rm QCD}$ Collins:1981uk ; Collins:1981va ;
Collins:1984kg . This nonperturbative sensitivity appears through the resummed
rapidity anomalous dimension, which with the canonical scale choice in eq.
(105) gets evaluated at $\alpha_{s}(1/b)$. It is important to note that this
is not an artefact of performing the evolution in Fourier space. Rather this
is a physical effect, which also happens when the $\nu$ evolution is
consistently performed in momentum space. As shown in ref. Ebert:2016gcn , in
this case the appropriate resummed result for
$\gamma_{\nu,S}(\vec{p}_{T},\mu)$ explicitly depends on $\alpha_{s}(p_{T})$,
which means it becomes nonperturbative for $p_{T}\lesssim\Lambda_{\rm QCD}$.

For the massive quark corrections in the regime $q_{T}\sim m$ the $\mu$
dependence of the rapidity anomalous dimension is the same as for the massless
quarks, i.e. eq. (103), such that

$\displaystyle\tilde{\gamma}^{(h)}_{\nu,S}(b,m,\mu)=4\eta_{\Gamma}^{(n_{l})}(\mu_{0}(b,m),\mu)-4\eta_{\Gamma}^{(n_{l}+1)}(\mu_{0}(b,m),\mu)+\tilde{\gamma}^{(h)}_{\nu,S}(b,m,\mu_{0}(b,m))\,.$
(107)

Here $\tilde{\gamma}^{(h)}_{\nu,S}$ denotes only the contributions of the
massive flavor to the full anomalous dimension. The explicit mass dependence
arises in the $\mu$-independent boundary contribution, which depends on both
$b$ and $m$. From the relations in eqs. (26) and (24) we can directly infer
the limiting behavior to the anomalous dimension,

$\displaystyle\tilde{\gamma}_{\nu,S}(b,m,\mu)$
$\displaystyle=\tilde{\gamma}_{\nu,S}^{(n_{l}+1)}(b,\mu)+\mathcal{O}(m^{2}b^{2})\,,$
$\displaystyle\tilde{\gamma}_{\nu,S}(b,m,\mu)$
$\displaystyle=\tilde{\gamma}_{\nu,S}^{(n_{l})}(b,\mu)+\gamma_{\nu,H_{s}}(m,\mu)+\mathcal{O}\Bigl{(}\frac{1}{m^{2}b^{2}}\Bigr{)}\,.$
(108)

This means that the massive quark corrections $\tilde{\gamma}^{(h)}_{\nu,S}$
are the same as for a massless flavor in the limit $m\ll 1/b$ and are the same
as the rapidity anomalous dimension of the soft mass mode function $H_{s}$ in
the limit $1/b\ll m$, provided one uses the $(n_{l}+1)$ and $(n_{l})$-flavor
scheme for $\alpha_{s}$, respectively. To eliminate the logarithms inside
$\tilde{\gamma}^{(h)}_{\nu,S}$, the canonical scale choice $\mu_{0}(b,m)$
should behave like the massless case for $m\ll 1/b$ and like the choice for
the mass-mode matching functions for $m\gg 1/b$,

$\displaystyle\mu^{(h)}_{0}(b,m)$
$\displaystyle\sim\mu^{(l)}_{0}(b)=\frac{2\,\mathrm{e}^{-\gamma_{E}}}{b}$
$\displaystyle\text{for}\quad 1/b\to\infty\,,$
$\displaystyle\mu^{(h)}_{0}(b,m)$ $\displaystyle\sim m$
$\displaystyle\text{for}\quad 1/b\to 0\,.$ (109)

Since $\mu^{(h)}_{0}(b,m)$ freezes out naturally at the perturbative mass
scale for $1/b\to 0$, the nonperturbative sensitivity in the $\nu$ evolution
gets regulated by the quark mass for the massive quark contributions.

We first illustrate this behavior in a simple one-loop toy example: We
consider the radiation of a massive gluon (with mass $M$) having the same
couplings as a (massless) gluon in QCD, which exhibits the main features of
the full results for secondary massive quarks. The associated corrections are
obtained in the calculations of app. B.4.1 as intermediate results for the
two-loop case. In $b$-space the one-loop rapidity anomalous dimensions for
massless and massive gluons are given by

$\displaystyle\tilde{\gamma}_{\nu,S}^{(1)}(b,\mu)$
$\displaystyle=-C_{F}\,8L_{b}\,,$
$\displaystyle\tilde{\gamma}_{\nu,S}^{(1)}(b,M,\mu)$
$\displaystyle=C_{F}\,\Bigl{[}8L_{M}+16K_{0}(bM)\Bigr{]}\,,$ (110)

where $K_{0}$ denotes the modified Bessel function of the second kind and

$\displaystyle
L_{b}\equiv\ln\frac{b^{2}\mu^{2}\mathrm{e}^{2\gamma_{E}}}{4}\,,\qquad
L_{M}\equiv\ln\frac{M^{2}}{\mu^{2}}\,.$ (111)

The mass-dependent result has the limiting behavior

$\displaystyle\tilde{\gamma}_{\nu,S}^{(1)}(b,M,\mu)$
$\displaystyle=-C_{F}\,8L_{b}+\mathcal{O}(M^{2}b^{2})\,,$
$\displaystyle\tilde{\gamma}_{\nu,S}^{(1)}(b,M,\mu)$
$\displaystyle=C_{F}\,8L_{M}+\mathcal{O}\Bigl{(}\frac{1}{M^{2}b^{2}}\Bigr{)}\;,$
(112)

in close analogy to eq. (5). A natural choice to eliminate any large terms in
eq. (5) in both limits is

$\displaystyle\mu^{(h)}_{0}(b,M)=M\,e^{K_{0}(bM)}\,.$ (113)

for which $\tilde{\gamma}^{(1)}_{\nu,S}(b,M,\mu^{(h)}_{0}(b,M))$ just
vanishes. The behavior of this choice as a function of $b$ compared to the
massless result is shown in fig. 16.

Figure 16: The canonical scales $\mu^{(h)}_{0}(b,M)$ for the massive case
(red, solid) and $\mu_{0}^{(l)}(b)=\mu_{0}(b,M=0)$ for the massless case
(blue, dashed) with $M=5$ GeV.

For the full secondary massive quark corrections at
$\mathcal{O}(\alpha_{s}^{2})$ the Fourier transform of eq. (4.3.1) reads
(expanded in terms of $\alpha_{s}^{(n_{l}+1)}$ as in eq. (49))

$\displaystyle\tilde{\gamma}_{\nu,S}^{(2,h)}(b,m,\mu)$
$\displaystyle=C_{F}\,\biggl{\\{}-\frac{32}{3}L_{b}L_{m}-\frac{16}{3}L_{m}^{2}-\frac{160}{9}L_{m}-\frac{448}{27}$
$\displaystyle\quad+\frac{8\sqrt{\pi}}{3}\Bigl{[}2\,G^{\,3,0}_{1,3}\Big{\lparen}\begin{smallmatrix}\frac{3}{2}\\\
0,0,0\end{smallmatrix}\Big{|}\,m^{2}b^{2}\Big{\rparen}+G^{\,3,0}_{1,3}\Big{\lparen}\begin{smallmatrix}\frac{5}{2}\\\
0,0,1\end{smallmatrix}\Big{|}\,m^{2}b^{2}\Big{\rparen}\Bigr{]}\biggr{\\}}\,,$
(114)

where $G$ denotes a Meijer G function. This result has the limiting behavior

$\displaystyle\tilde{\gamma}_{\nu,S}^{(2,h)}(b,m,\mu)$
$\displaystyle=C_{F}\,\biggl{(}\frac{16}{3}L_{b}^{2}+\frac{160}{9}L_{b}+\frac{448}{27}\biggr{)}+\mathcal{O}(m^{2}b^{2})\,,$
$\displaystyle\tilde{\gamma}_{\nu,S}^{(2,h)}(b,m,\mu)-\frac{4}{3}L_{m}\tilde{\gamma}_{\nu,S}^{(1)}(b,\mu)$
$\displaystyle=C_{F}\,\biggl{(}-\frac{16}{3}L_{m}^{2}-\frac{160}{9}L_{m}-\frac{448}{27}\biggr{)}+\mathcal{O}\Bigl{(}\frac{1}{m^{2}b^{2}}\Bigr{)}\,.$
(115)

Hence, the correct massless limit is recovered, while in the large-mass limit
one obtains the anomalous dimension in eq. (4.1.2). Note that one needs to
perform a change for the strong coupling between the $n_{l}+1$ and $n_{l}$
flavor schemes to obtain both limits correctly. To minimize the logarithms for
any regime one should thus adopt a canonical scale choice that satisfies eq.
(5), as for example in eq. (113).

##  6 Outlook: Phenomenological impact for Drell-Yan

Our results can be applied to properly take into account bottom quark mass
effects for the Drell-Yan $q_{T}$ spectrum at NNLL′. While a full resummation
analysis is beyond the scope of this paper, we can estimate the potential size
of the quark-mass effects by looking at the fixed-order $q_{T}$ spectrum.

In fig. 17, we show separately the contributions from primary and secondary
massive quarks to the cross section at $\mathcal{O}(\alpha_{s}^{2})$,
normalized to the $\mathcal{O}(\alpha_{s})$ spectrum $\textrm{d}\sigma^{(1)}$
including all flavors (treating the charm as a massless flavor). We utilize
the MMHT2014 NNLO PDFs Harland-Lang:2014zoa  and evaluate the contributions
for $\mu=m_{b}=4.8$ GeV, $Q=m_{Z}$, $Y=0$, and $E_{\rm cm}=13$ TeV. Note that
the secondary mass contributions at $\mathcal{O}(\alpha_{s}^{2})$ are
explicitly $\mu$-dependent and scheme-dependent, the nonsingular mass
correction, i.e. the difference between the full massive result for $\mu\sim
m_{b}$ and the massless limit (encoded partially in a massive PDF), is $\mu$
independent at this order. As can be seen, the relative contribution of the
$b\bar{b}$-initiated channel grows with larger $q_{T}$, while the impact of
the secondary contributions including the full mass dependence is at the sub-
percent level throughout the spectrum. As expected, the nonsingular mass
corrections are very small for $m_{b}\ll q_{T}$, but can reach the order of