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

TL;DR

This paper develops an effective-theory framework using soft-collinear effective theory to incorporate massive quark effects into resummed differential distributions for exclusive processes like Drell-Yan, covering various mass and scale hierarchies.

Contribution

It introduces a comprehensive calculation of quark-mass effects at NNLL' order for $q_T$ and beam thrust, including primary and secondary contributions across different hierarchies.

Findings

Calculated all ingredients for primary and secondary mass effects at NNLL' order.

Analyzed rapidity divergences and evolution in the massive quark case.

Enabled detailed studies of quark-mass effects in $W$ and $Z$ boson spectra.

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.

Equations

$$\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}\}\,.$$

$$n_a\text{-collinear:}$$

$$n_b\text{-collinear:}$$

$$p_s^{\mu }\sim (q_T,q_T,q_T),$$

$$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 }),$$

$$\frac{\text{d}\sigma }{\text{d}{q}_T^2\text{d}{Q}^2\text{d}Y}$$

$$\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{]}\,,$$

$${\omega }_a=Q{e}^Y,{\omega }_b=Q{e}^{-Y},x_{a,b}=\frac{{\omega }_{a,b}}{E_{\mathrm{cm}}},$$

$$\displaystyle B^{(n_{f})}_{i}\Bigl{(}\vec{p}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}$$

$$\equiv \sum _k{\mathcal{I}}_{ik}^{(n_f)}({\vec{p}}_T,x,\mu ,\frac{\nu }{\omega }){\otimes }_x{f}_k^{(n_f)}(x,\mu ),$$

$$\mu \frac{\text{d}}{\text{d}\mu }{H}_{ij}^{(n_f)}(Q,\mu )$$

$$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}B^{(n_{f})}_{i}\Bigl{(}\vec{p}_{T},x,\mu,\frac{\nu}{\omega}\Bigr{)}$$

$$\mu \frac{\text{d}}{\text{d}\mu }{S}^{(n_f)}({\vec{p}}_T,\mu ,\nu )$$

$$\mu \frac{\text{d}}{\text{d}\mu }{f}_i^{(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{)}$$

$$\nu \frac{\text{d}}{\text{d}\nu }{S}^{(n_f)}({\vec{p}}_T,\mu ,\nu )$$

$$\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)$$

$$2{\gamma }_{\nu ,B}^{(n_f)}({\vec{p}}_T,\mu )+{\gamma }_{\nu ,S}^{(n_f)}({\vec{p}}_T,\mu )$$

$$\mu \frac{\text{d}}{\text{d}\mu }{\gamma }_{\nu ,S}^{(n_f)}({\vec{p}}_T,\mu )=\nu \frac{\text{d}}{\text{d}\nu }{\gamma }_S^{(n_f)}(\mu ,\nu )\delta ({\vec{p}}_T)$$

$$\frac{\text{d}\sigma }{\text{d}{q}_T^2\text{d}{Q}^2\text{d}Y}$$

$$\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{]}\,,$$

$$n_a\text{-collinear MM:}$$

$$n_b\text{-collinear MM:}$$

$$p_{m,s}^{\mu }\sim (m,m,m),$$

$$\frac{\text{d}\sigma }{\text{d}{q}_T^2\text{d}{Q}^2\text{d}Y}$$

$$\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{]}\,.$$

$$\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)\,,$$

$${\gamma }_{\nu ,H_s}(m,\mu )=-2{\gamma }_{\nu ,H_c}(m,\mu )=\frac{\text{d}}{\text{d}\ln \nu }\ln H_s(m,\mu ,\nu ).$$

$$\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{]}\,.$$

$$\frac{\text{d}\sigma }{\text{d}{q}_T^2\text{d}{Q}^2\text{d}Y}$$

$$\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{]}\,,$$

$$\displaystyle\mu\frac{\textrm{d}}{\textrm{d}\mu}\,\mathcal{I}_{ik}\Bigl{(}\vec{p}_{T},m,z,\mu,\frac{\nu}{\omega}\Bigr{)}$$

$$\displaystyle\nu\frac{\textrm{d}}{\textrm{d}\nu}B_{i}^{(n_{l}+1)}\Bigl{(}\vec{p}_{T},m,\mu,\frac{\nu}{\omega}\Bigr{)}$$

$$\nu \frac{\text{d}}{\text{d}\nu }{S}^{(n_l+1)}({\vec{p}}_T,m,\mu ,\nu )$$

$$\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{]}\,.$$

$$\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{]}\,.$$

$$\frac{\text{d}\sigma }{\text{d}{q}_T^2\text{d}{Q}^2\text{d}Y}$$

Full text

††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 ,

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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)$.

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$,

[TABLE]

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}$,

[TABLE]

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.

[TABLE]

Written out explicitly, the factorization theorem reads

[TABLE]

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

[TABLE]

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 .,

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

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}$,

[TABLE]

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.

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.

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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}$,

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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).

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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,

[TABLE]

where

[TABLE]

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,

[TABLE]

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

[TABLE]

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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)}$,

[TABLE]

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

[TABLE]

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}$,

[TABLE]

where here and in the following we abbreviate

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Its anomalous dimensions are

[TABLE]

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,

[TABLE]

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

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}$)

[TABLE]

with the splitting function

[TABLE]

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 .

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$

[TABLE]

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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,

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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}$.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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),

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Its solution is given by

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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$,

[TABLE]

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

[TABLE]

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

[TABLE]

The mass-dependent result has the limiting behavior

[TABLE]

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

[TABLE]

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.

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))

[TABLE]

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

[TABLE]

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 percent for $q_{T}\sim m_{b}$, which roughly corresponds to the peak region of the distribution where the cross section is largest.

The same can also be seen in fig. 18, where we show the mass nonsingular corrections to the massless limit for primary and secondary contributions as well as their sum. They are shown for $\mu=m_{b}$ on the left and for $\mu=q_{T}$ on the right. We see that these corrections are (at fixed order) indeed only weakly dependent on the value of $\mu$ (for $q_{T}\gtrsim 2$ GeV). All in all, the bottom quark mass can have a relevant effect for high precision predictions of the $q_{T}$-spectrum at the order of percent around the peak of the distribution ($\sim 5$ GeV). Below the peak of the distribution the fixed-order result is of course not expected to give a reliable quantitative result, and furthermore nonperturbative corrections become important in this regime. Nevertheless, we expect the qualitative features like the sign and order of magnitude of the mass effects to provide an indication for the behaviour of the full resummed result.

For $W$ production sizable corrections from bottom quark effects arise only through secondary contributions (due to the strong CKM suppression of the primary contributions), which have a similar impact as for $Z$-production. On the other hand, charm-initiated production plays an important role and enters already at $\mathcal{O}(\alpha_{s})$. Estimating the nonsingular mass corrections for $q_{T}\sim m_{c}$ is more subtle, since higher-order corrections in the strong coupling and nonperturbative effects are likely to dominate the effect from the known beam function at $\mathcal{O}(\alpha_{s})$ at these low scales. Thus, we do not attempt to determine their characteristic size here and leave this to future work. An analysis based on the leading-order matrix element and its potential impact on the determination of $m_{W}$ can be found in ref. Berge:2005rv .

7 Conclusions

Massive quark effects provide a challenge for high-precision predictions at colliders. Using a SCET-based factorization framework, we have discussed how to systematically incorporate massive quark corrections into exclusive differential cross sections at the LHC, using the measurement of the transverse momentum $q_{T}$ and beam thrust for Drell-Yan production as prototypical examples. We have discussed the relevant factorization setup for the different hierarchies between the mass scale and the other relevant kinematic scales. We find that the presence of (secondary) massive quarks can lead to the emergence or alteration of rapidity logarithms thus changing the resummation structure in a nontrivial way.

The generic framework for the description of mass effects generalizes to other exclusive cross sections with different jet-resolution measurements and final-state kinematic cuts, which will require additional calculations of the relevant factorization ingredients. Our results for the beam thrust spectrum allow for a systematic inclusion of massive quark effects into the Geneva Monte-Carlo program Alioli:2012fc ; Alioli:2015toa at NNLL${}^{\prime}+$NNLO in its underlying jet resolution variable. Several of our results are also immediately relevant for other processes besides Drell-Yan. The massive quark beam functions are relevant for any heavy-quark initiated process, for example exclusive $b\bar{b}H$-production. The mass-dependent soft function and rapidity anomalous dimension at $\mathcal{O}(\alpha_{s}^{2})$ satisfy Casimir scaling and can be therefore also utilized for the description of gluon-fusion processes, e.g. the Higgs $q_{T}$-spectrum.

An important application of our framework is to the precise theoretical description of the Drell-Yan $q_{T}$ spectrum. To this end, we have computed all required mass-dependent beam and soft functions up to $\mathcal{O}(\alpha_{s}^{2})$ allowing for the description of massive quark effects in the Drell-Yan $q_{T}$ spectrum at NNLL*′*. In particular, our results provide an important ingredient 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.

Acknowledgements.

We thank Markus Ebert and Maximilian Stahlhofen for useful discussions. This work was supported by the German Science Foundation (DFG) through the Emmy-Noether Grant No. TA 867/1-1 and the Collaborative Research Center (SFB) 676 Particles, Strings and the Early Universe. D.S. is supported by the Austrian Science Fund (FWF) under the Doctoral Program No. W1252-N27 Particles and Interactions. We also thank the Erwin-Schrödinger Institute (ESI) for partial support in the framework of the ESI program Challenges and Concepts for Field Theory and Applications in the Era of LHC Run-2.

Appendix A Results for massless quarks

Here we summarize the relevant results with massless quarks for the hard, beam, and soft functions.

A.1 Hard function

The massless quark hard function is directly related to the QCD form factor and has been computed at $\mathcal{O}(\alpha_{s}^{2})$ in ref. Matsuura . The $\mathcal{O}(\alpha_{s})$ and $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ corrections read in an expansion in terms of $\alpha_{s}=\alpha_{s}^{(n_{f})}(\mu)$ in analogy to eq. (49) (with $L_{Q}=\ln(Q^{2}/\mu^{2})$)

[TABLE]

where $H^{(0)}$ is the tree-level contribution. Note that for a single quark flavor there is in addition a nonvanishing correction to the axial current contribution relevant for $Z$-boson production, but cancels within an isospin doublet for massless quarks.

The anomalous dimensions are

[TABLE]

A.2 Beam functions

A.2.1 TMD beam function

The matching coefficients entering the TMD beam function have been computed at $\mathcal{O}(\alpha_{s}^{2})$ in various schemes Catani:2012qa ; Gehrmann:2012ze ; Gehrmann:2014yya ; Echevarria:2016scs and are obtained for the symmetric $\eta$-regulator in ref. Luebbert:2016itl . The results at $\mathcal{O}(\alpha_{s})$ are

[TABLE]

The splitting functions are

[TABLE]

At $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ the massless matching coefficient is given by

[TABLE]

The anomalous dimensions of the massless quark TMD beam function, as defined in eq. (2.1), are given at $\mathcal{O}(\alpha_{s})$ and $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ by

[TABLE]

A.2.2 Virtuality-dependent beam function

The virtuality-dependent beam functions for massless quarks are known to two loop order Gaunt:2014xga ; Gaunt:2014cfa . The matching coefficients at $\mathcal{O}(\alpha_{s})$ read

[TABLE]

The massless matching coefficient at order $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ for one quark flavor reads

[TABLE]

The anomalous dimension of the massless quark beam function at order $\mathcal{O}(\alpha_{s})$ and $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ are given by

[TABLE]

A.3 Soft functions

A.3.1 TMD soft function

The TMD soft function for massless quarks with the symmetric $\eta$-regulator has been computed at two loops in ref. Luebbert:2016itl . At $\mathcal{O}(\alpha_{s})$ and $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ it is given by

[TABLE]

The corresponding anomalous dimensions are

[TABLE]

A.3.2 Thrust soft function

The thrust soft function is known to two loops Kelley:2011ng ; Monni:2011gb . At $\mathcal{O}(\alpha_{s})$ and $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ it is given by

[TABLE]

The corresponding $\mu$ anomalous dimension is given by

[TABLE]

Appendix B Calculations of massive quark corrections

We calculate the quark mass dependent beam and soft functions for primary and secondary contributions at one and two loops, respectively. The final renormalized results are given and discussed in sec. 4. For the computation of the collinear massive quark corrections we use the Feynman rules determined from the collinear massive quark Lagrangian Rothstein:2003wh ; Leibovich:2003jd ; Fleming:2007qr . For the secondary corrections we use in practice regular QCD Feynman rules, since the collinear sector is essentially just a boosted version of QCD. (The interactions of the massive quarks in the soft sector are anyway given by the usual QCD Feynman rules.) First, we calculate the massive quark beam function in sec. B.1, before discussing the computation of the secondary corrections for the massless quark beam and soft functions in secs. B.2 – B.5. All computations are carried out in Feynman gauge.

B.1 Massive quark beam function at $\mathcal{O}(\alpha_{s})$

The massive quark beam function operator for a measurement function $\mathcal{M}$ is defined as (see e.g. refs. Stewart:2009yx ; Stewart:2010qs ; Jain:2011iu ; Chiu:2012ir )

[TABLE]

where $\chi_{n,m}$ indicates a massive collinear quark field, $\mathcal{P}^{\mu}$ is the label momentum operator, and $\hat{p}^{+}$ extracts the residual momentum component $n\cdot k$. For the transverse momentum dependent (TMD), virtuality dependent, and fully differential case the measurement functions are

[TABLE]

For convenience we discuss also the fully differential case here, from which the other two cases can be obtained by an integration over the respective other variable. The beam functions are proton matrix elements of the operators $\mathcal{O}_{Q}$. To compute the (perturbative) matching coefficients onto the PDFs, we take matrix elements with partonic states, denoting e.g.

[TABLE]

for an initial collinear gluon state with momentum $p^{\mu}=p^{-}n^{\mu}/2$.

At $\mathcal{O}(\alpha_{s})$ the only contribution to the massive quark beam function originates from an initial collinear gluon splitting into a heavy quark-antiquark pair. The corresponding diagram is given in fig. 19. The kinematics of the on-shell final state is fully constrained at one loop, so that the diagram can be evaluated without performing any integration. For the fully differential case we obtain

[TABLE]

where $P_{qg}(z)=z^{2}+(1-z)^{2}$ is the leading-order gluon-quark splitting function. The correction $B_{Qg}$ at $\mathcal{O}(\alpha_{s})$ is UV and IR finite. It corresponds directly to the matching coefficient $\mathcal{I}^{(1)}_{Qg}$, given as the one-loop coefficient in an expansion in terms of $\alpha_{s}$ as in eq. (49). The matching coefficients for the TMD and virtuality-dependent beam functions can be obtained here by a trivial integration of this result,

[TABLE]

which yields the results in eqs. (63) and (69). Note that in general, this integration has to be performed for the bare result with the full dependence on the UV and rapidity regulator. However, in this case all matrix elements are finite and do not require any renormalization at this order.

B.2 Dispersive technique for secondary massive quark corrections

For observables where only the sum over the final-state hadronic momenta enters the measurement, one can use dispersion relations to obtain the results for secondary massive quark radiation at $\mathcal{O}(\alpha_{s}^{2})$ from the corresponding results for “massive gluon” radiation at $\mathcal{O}(\alpha_{s})$. This has been discussed in detail in ref. Pietrulewicz:2014qza . The key relation is that the insertion of a vacuum polarization function for massive quarks $\Pi_{\mu\nu}(m^{2},p^{2})$ between two gluon propagators can be written as

[TABLE]

The first term contains a gluon propagator with effective mass $M$ and the absorptive part of the vacuum polarization function, which reads in $d=4-2\epsilon$ dimensions

[TABLE]

To obtain the first term on the right-hand side in eq. (B.2) the vacuum polarization function (and thus the strong coupling) was renormalized in the on-shell scheme, i.e., with $n_{l}$ active quark flavors. The second term in eq. (B.2) translates back to an unrenormalized strong coupling and consists of a massless gluon propagator and the $\mathcal{O}(\alpha_{s})$ vacuum polarization function at zero momentum transfer, which is given by

[TABLE]

In the following we will first carry out the computation of the beam and soft functions at $\mathcal{O}(\alpha_{s})$ for the radiation of a “massive gluon” and in a second step use the relation in eq. (B.2) to obtain the associated results for massive quarks at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$. In our calculations we drop the contributions from the terms proportional to $p^{\mu}p^{\nu}$, which vanish in total due to gauge invariance.

B.3 Secondary mass effects in light-quark beam functions

We compute the massive quark corrections to the TMD and virtuality-dependent light-quark beam function at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ starting with the massive gluon case at $\mathcal{O}(\alpha_{s})$. Only the contributions to the matching coefficient $\mathcal{I}_{qq}$ are nontrivial, so we consider only diagrams with a quark in the initial state.

B.3.1 Quark beam function with a massive gluon at $\mathcal{O}(\alpha_{s})$

Contributions to the fully-differential beam function

As in sec. B.1 we start also here with the computation of the corrections for the fully-differential beam function. The contributing one-loop diagrams to the matrix element $B_{qq}$ with massless quarks in the initial state, defined in analogy to eq. (134), are displayed in fig. 20. They consist of a purely virtual and a real-radiation part,

[TABLE]

The virtual massive gluon contributions in fig. 20 are the same as for other collinear quark operators like the current or the PDF and have been computed e.g. in ref. Chiu:2012ir . Including the wave function renormalization diagrams the $d$-dimensional result reads Hoang:2015iva

[TABLE]

where $H_{\alpha}=\psi(1+\alpha)+\gamma_{E}$ is the Harmonic number. Here the rapidity divergences have been regulated using the symmetric $\eta$ regulator acting on the Wilson lines Chiu:2012ir ; Chiu:2011qc , while UV divergences are regulated with dimensional regularization as usual. Furthermore, the gluon mass provides an IR cutoff.

The real radiation contributions in figs. 20 and 20 can be easily evaluated, since all momentum components are fully determined by the measurement. For the first diagram we get

[TABLE]

Since UV divergences do not appear for the real radiation corrections and the gluon mass regulates all IR divergences we do not need to employ dimensional regularization here. The second diagram in fig. 20 yields

[TABLE]

While the fully-differential quark beam function itself does not contain any rapidity divergences, we have included here the $\eta$ regulator, since we will use this result to obtain the TMD beam function by integrating over the virtuality, which results in rapidity divergences for this real radiation correction. The full real radiation contributions at one loop yield

[TABLE]

For both virtual and real radiation corrections all soft-bin subtractions are parametrically power suppressed or scaleless and therefore do not contribute.

Contributions to the TMD beam function

The corrections for the TMD beam function with a massive gluon can be obtained by integrating the fully-differential beam function in eq. (140) over the virtuality $t$. We write them again as

[TABLE]

where $B^{(1,\rm{bare})}_{qq,\rm{virt}}$ is given in eq. (141) and

[TABLE]

Here it is necessary to keep a nonvanishing value for $\eta$ in the second term to regularize the rapidity divergence for $z\to 1$. Expanding for $\eta\to 0$ we get

[TABLE]

Contributions to the virtuality-dependent beam function

The virtuality-dependent beam function with a massive gluon can be obtained by integrating the results for the fully-differential beam function over $\vec{p}_{T}$. We decompose the corrections again into a virtual and real radiation part,

[TABLE]

where $B^{(1,\rm{bare})}_{qq,\rm{virt}}$ is given in eq. (141) and

[TABLE]

with the fully-differential real radiation contributions in eq. (144). Here the $\eta$ regulator has already been dropped, since for the virtuality-dependent beam function no rapidity divergences arise from the real radiation contributions.

B.3.2 Secondary massive quark effects in the TMD beam function

To obtain the secondary massive quark corrections from the one-loop results with a massive gluon, we first convolve the one-loop results with the imaginary part of the vacuum polarization function according to eq. (B.2) and define

[TABLE]

The results from these dispersion integrations are

[TABLE]

with

[TABLE]

and $\hat{m}$, $c$, $d$ defined in eq. (66). Using eq. (B.3.2) entails that the massive quark corrections to the strong coupling are renormalized in the on-shell scheme, i.e., the expansion is in terms of $\alpha_{s}=\alpha_{s}^{(n_{l})}$. Since the beam function matrix element has to be renormalized entirely in the $n_{l}+1$ flavor theory, we need to account for the second term in eq. (B.2) (which switches back to an unrenormalized $\alpha_{s}$) and renormalize the massive quark corrections to the strong coupling in the $\overline{\rm MS}$ scheme, such that the expansion is in terms of $\alpha_{s}=\alpha_{s}^{(n_{l}+1)}$. The beam function operator is renormalized according to

[TABLE]

where the counterterm encodes also the rapidity divergences. This yields for the renormalized matrix element with initial state quarks at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ in terms of $\alpha_{s}=\alpha_{s}^{(n_{l}+1)}$

[TABLE]

where the (bare) vacuum polarization function $\Pi^{(1)}(m^{2},0)$ is given in eq. (139). The one-loop counterterm reads

[TABLE]

The two-loop counterterm $Z_{B}^{(2)}$ absorbs all remaining UV and rapidity divergences in eq. (155) and is given by

[TABLE]

This yields the anomalous dimensions in eq. (4.2.1). The renormalized one-loop partonic beam function $B_{qq}^{(1)}$ still contains IR divergences, so its exact form depends on the choice of the IR regulator.

The beam function matching coefficient $\mathcal{I}_{qq}$ as defined in (14) can be now easily obtained. Note that the PDFs are renormalized in an $n_{l}$-flavor theory with $\alpha_{s}=\alpha_{s}^{(n_{l})}$ in contrast to the beam function. Thus, there is a contribution coming from the scheme change of $\alpha_{s}$ to $n_{l}+1$ flavors for the (renormalized) one-loop PDF correction, i.e.

[TABLE]

Here the IR divergences cancel between the one-loop beam function and the PDF to give the finite one-loop matching coefficient $\mathcal{I}_{qq}^{(1)}$, which is given in eq. (119). Using eqs. (151), (152), (156) and (B.3.2) we obtain the full result for the secondary massive quark corrections to the beam function matching coefficient given in eq. (4.2.1).

B.3.3 Secondary massive quark effects in the virtuality-dependent beam function

We proceed with the virtuality-dependent beam function. While the virtual contributions are the same as for the TMD beam function given in eq. (151), the dispersion integration for the real radiation terms yields

[TABLE]

with

[TABLE]

and $\hat{m}_{t}$, $u$, $v$ as in eq. (71).

To obtain the quark mass dependent matching coefficient $\mathcal{I}^{(2,h)}_{qq}$ we carry out our calculation using a gluon mass $\Lambda\ll\sqrt{Q{\mathcal{T}}}\sim m$ as IR regulator. Although the result is independent of the regulator, this is technically most convenient, since this allows us to match two ${\rm SCET}_{\rm II}$ theories with each other in a straightforward way.777Alternatively, one can also perform the matching between theories where the fluctuations related to the $n_{l}$ massless flavors are described within a ${\rm SCET}_{\rm I}$ theory. In this setup, there is no csoft function on the right-hand side of the matching relation in contrast to eq. (161). However, in this case the zero-bin subtractions for the collinear fields with respect to the ultrasoft modes in the ${\rm SCET}_{\rm I}$ $n_{l}$ flavor theory yield a nontrivial contribution to the beam-function matrix element on the left-hand side of the matching relation. Their contribution is equivalent to the inverse of the csoft function in eq. (161), such that the resulting matching coefficient $\mathcal{I}_{qq}$ is the same. While the ${\rm SCET}_{\rm II}$ theory with $n_{l}+1$ flavors (i.e. above the mass scale) contains collinear modes, the ${\rm SCET}_{\rm II}$ theory with $n_{l}$ flavors (i.e. below the mass scale) contains collinear and csoft modes like in the mode setup of sec. 3.3. The matching relation reads

[TABLE]

where $\mathscr{B}_{qq}^{(n_{l}+1)}$ corresponds to the pure ${\rm SCET}_{\rm II}$ beam function matrix element and $\mathscr{S}^{(n_{l})}$ represents the csoft matrix element.

In close analogy to eq. (155) the renormalized ${\rm SCET}_{\rm II}$ matrix element $\mathscr{B}_{qq}^{(n_{l}+1)}$ is given at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ by

[TABLE]

To separate UV, rapidity, and IR divergences properly from each other, we also employ the ${\rm SCET}_{\rm II}$-type IR regulator (here a gluon mass $\Lambda$) for the one-loop expressions, and at this stage the renormalized matrix elements and the counterterms still depend on this IR regulator. The matching coefficient $\mathcal{I}_{qq}$ can now be calculated as (in an expansion in terms of $\alpha_{s}^{(n_{l}+1)}$)

[TABLE]

Here the IR divergences cancel between the one-loop beam function, the PDF, and the csoft matrix element and yield the finite one-loop matching coefficient $\mathcal{I}_{qq}^{(1)}$ given in eq. (123). The counterterm $Z_{\mathscr{B}}$ in ${\rm SCET}_{\rm II}$ is defined via

[TABLE]

Using the results in eqs. (148) and (141) for a massive gluon gives the associated expression for $Z_{\mathscr{B}}^{(1)}$ (expanded in $\eta$ and $\epsilon$)

[TABLE]

The two-loop counterterm $Z^{(2,h)}_{\mathscr{B}}$ cancels all divergences in eq. (B.3.3) and reads888While the $1/\eta$-divergences in the counterterm of the beam function matrix element still contain IR sensitivity, this also happens for the counterterm of the csoft matrix element in eq. (161), such that the resulting rapidity anomalous dimension for the running at the boundary between the $n_{l}+1$ and $n_{l}$ theory is IR finite.

[TABLE]

Using eqs. (141), (B.3.3), (165), (B.3.3) and (123) in eq. (B.3.3) we obtain the full two-loop result for the matching coefficient in eq. (4.2.2).

B.4 Secondary mass effects in the TMD soft function

The TMD soft function is defined as

[TABLE]

with the soft Wilson line $S_{n}$ given by Chiu:2012ir

[TABLE]

and in analogy for the others. Again we will first calculate the one-loop corrections to the soft function with a massive gluon, which is used in a second step to obtain the corrections from secondary massive quarks at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$.

B.4.1 TMD soft function with a massive gluon at $\mathcal{O}(\alpha_{s})$

We decompose the soft function with a massive gluon at one loop in terms of virtual and real radiation corrections,

[TABLE]

The virtual contributions from the diagram in fig. 21 (and its mirror diagram) are the same as for the Sudakov form factor computed in ref. Chiu:2012ir and yield

[TABLE]

The UV-finite and IR-finite real radiation diagram in Fig. 21 gives

[TABLE]

After expanding in $\eta$ and adding the mirror diagram, the real radiation contribution to the TMD soft function at one loop then reads

[TABLE]

B.4.2 Secondary corrections at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$

To obtain the secondary massive quark corrections from the one-loop results with a massive gluon, we first convolve the one-loop results with the imaginary part of the vacuum polarization function,

[TABLE]

The results from these dispersion integrations are

[TABLE]

with

[TABLE]

and $\hat{m}$ and $c$ as in eq. (66). Using eq. (B.4.2) entails that the massive quark corrections to the strong coupling are renormalized in the on-shell scheme, i.e., the expansion is in terms of $\alpha_{s}=\alpha_{s}^{(n_{l})}$. Since the soft function matrix element has to be renormalized entirely in the $n_{l}+1$ flavor theory, we need to account for the second term in eq. (B.2) (which switches back to an unrenormalized $\alpha_{s}$) and renormalize the massive quark corrections to the strong coupling in the $\overline{\rm MS}$ scheme, such that the expansion is in terms of $\alpha_{s}=\alpha_{s}^{(n_{l}+1)}$. The soft function is renormalized according to

[TABLE]

This yields for the renormalized matrix element with initial state quarks at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ in terms of $\alpha_{s}=\alpha_{s}^{(n_{l}+1)}$

[TABLE]

where the (bare) vacuum polarization function $\Pi^{(1)}(m^{2},0)$ is given in eq. (139) and the renormalized one-loop soft function $S^{(1)}$ is given in eq. (126). The one-loop counterterm reads

[TABLE]

The two-loop counterterm $Z_{S}^{(2)}$ absorbs all remaining UV and IR divergences in eq. (178) and is given by

[TABLE]

This yields the anomalous dimensions in eq. (4.3.1). Using eqs. (B.4.2), (175), (179), and (B.4.2) we obtain the full result for the secondary massive quark corrections to the TMD soft function in eq. (4.3.1).

B.5 Csoft function at two loops

We compute the csoft function $\mathcal{S}_{c}$ for beam thrust appearing in the hierarchy ${\mathcal{T}}\ll m\ll\sqrt{Q{\mathcal{T}}}$. As in the computation for the beam function matching coefficient in app. B.3.3 we carry out the calculation using a ${\rm SCET}_{\rm II}$ IR regulator (a gluon mass $\Lambda\ll m$). In this context the csoft function is the matching coefficient between the csoft matrix elements in the $n_{l}+1$ and $n_{l}$ flavor ${\rm SCET}_{\rm II}$ theories,

[TABLE]

The latter are defined for any direction $n$ as

[TABLE]

with the csoft Wilson lines given by (see e.g. refs. Bauer:2011uc ; Procura:2014cba )

[TABLE]

Besides replacing the soft fields by csoft fields we have also expanded the $\eta$ regulator according to the soft scaling as in ref. Hoang:2015iva .999If the regulator is not expanded, nonvanishing soft-bin subtraction appear which eliminate the overlap with soft mass mode momentum regions, see ref. Lustermans:2016nvk .

B.5.1 Csoft function with a massive gluon at $\mathcal{O}(\alpha_{s})$

We will first calculate the one-loop corrections to the csoft matrix elements $\mathscr{S}$ with a massive gluon, that can then be used to obtain the two-loop corrections with secondary massive quarks using the dispersion technique described in sec. B.2. The one-loop results for the csoft matrix elements can be written as

[TABLE]

The relevant contributions at one loop are displayed in the diagrams in fig. 21, with the soft Wilson lines $S_{n}$ and $S_{\bar{n}}$ replaced by the csoft Wilson lines $X_{n}$ and $V_{n}$. With the choice of regularization in eq. (183) the virtual diagram leads to a scaleless integral, such that $\mathscr{S}^{(1,\rm{bare})}_{\rm{virt}}=0$. The real radiation diagram corresponding to fig. 21 yields

[TABLE]

Including also the mirror diagram and expanding in $\eta$ the total real radiation contribution to the csoft matrix element with a massive gluon is

[TABLE]

B.5.2 Csoft function at $\mathcal{O}(\alpha_{s}^{2})$

We convolve the one-loop results with the imaginary part of the vacuum polarization function, which yields for the nonvanishing contributions

[TABLE]

The result of this dispersion integral is

[TABLE]

Using eq. (187) entails that the massive quark corrections to the strong coupling are renormalized in the on-shell scheme, i.e., the expansion is in terms of $\alpha_{s}=\alpha_{s}^{(n_{l})}$. To obtain the csoft function $\mathcal{S}_{c}$ we need to switch to $\alpha_{s}^{(n_{l}+1)}$ and furthermore subtract the correction $\mathscr{S}^{(2,n_{l})}$ (with a strong coupling in the $n_{l}$ flavor scheme) according to eq. (181). All purely massless contributions cancel each other and we obtain for the $\mathcal{O}(\alpha_{s}^{2})$ corrections in an expansion in terms of $\alpha_{s}^{(n_{l}+1)}$

[TABLE]

Here the ${\rm SCET}_{\rm II}$ counterterm is defined via

[TABLE]

Employing a gluon mass the associated expression for $Z_{\mathscr{S}}^{(1)}$ can be read off from eq. (186) and is given by (expanded in $\eta$ and $\epsilon$)

[TABLE]

The counterterm $Z^{(2,h)}_{\mathscr{S}}$ absorbs all divergences and is given by101010The anomalous dimension for the csoft function $\mathcal{S}_{c}$ can be obtained from the ratio of $Z^{(n_{l}+1)}_{\mathscr{S}}$ and $Z^{(n_{l})}_{\mathscr{S}}$, upon which the IR sensitivity cancels.

[TABLE]

Using eqs. (B.5.2), (191) and (B.5.2) in eq. (189) we obtain the full result for the renormalized csoft function at two loops in eq. (79).

Appendix C Massive quark effects at fixed order

The factorization formulae in the secs. 2 and 3 contain together all information about the singular massive quark corrections to the differential cross sections in QCD at fixed order (for any given hierarchy between the mass and $q_{T}$/${\mathcal{T}}$). Here we provide the results at $\mathcal{O}(\alpha_{s}^{2})$ for Drell-Yan for both primary and secondary corrections. We write for each of these contributions ($e=q_{T}^{2},{\mathcal{T}}$)

[TABLE]

and expand the partonic result in the $n_{l}$-flavor scheme for $\alpha_{s}$ as

[TABLE]

where $\textrm{d}\sigma^{(0)}_{q\bar{q}}/(\textrm{d}Q^{2}\textrm{d}Y)$ denotes the Born cross section for the corresponding Drell-Yan process $q\bar{q}\to Z/\gamma^{*}\to\ell\bar{\ell}$. In this context $\textrm{d}\sigma^{(0)}_{Q\bar{Q}}/(\textrm{d}Q^{2}\textrm{d}Y)$ indicates the Born cross section for a massless quark $q$ with the same charge and isospin as the heavy quark $Q$.

C.1 Fixed-order result for the $q_{T}$ spectrum

The singular fixed-order corrections for the $q_{T}$-spectrum (i.e. for $q_{T}\ll Q$) at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ consist of the virtual (full QCD) contributions encoded in eqs. (4.1.1) and (4.1.1) and the secondary collinear and soft real radiation corrections contained in eqs. (4.2.1) and (4.3.1). Setting common scales $\mu=\mu_{H}=\mu_{B}=\mu_{S}$ and $\nu=\nu_{B}=\nu_{S}$ yields for the corrections to virtual photon production

[TABLE]

where $h_{\rm virt}$, $b^{q_{T}}_{\rm{real}}$ and $s^{q_{T}}_{{\rm real}}$ are given in eqs. (4.1.1), (153) and (176). For $Z$-boson production one has to include in addition the anomalous axial current correction in eq. (4.1.1) as contribution to the $\delta(q_{T}^{2})$-term (which gives in conjunction with the isospin partner a $\mu$-independent result). Writing out the nontrivial terms in the spectrum explicitly we get

[TABLE]

where

[TABLE]

The singular fixed-order corrections for the $q_{T}$ spectrum at $\mathcal{O}(\alpha_{s}^{2}T_{F}^{2})$ consist of the primary collinear real radiation corrections in eq. (63) for both beam directions,

[TABLE]

Depending on the hierarchy between $m$ and $q_{T}$ and $Q$ some of the contributions in eqs. (C.1) and (C.1) are power-suppressed and therefore only appear via nonsingular corrections in the factorization formula for the associated parametric regime in sec. 2. Note also that virtual corrections are reshuffled among the components of the factorization theorem, which are in addition evaluated with $\alpha_{s}$ in different flavor number schemes. This essentially allows for a consistent factorization and the resummation of logarithms at higher orders.

C.2 Fixed-order result for the beam thrust spectrum

The singular fixed-order corrections for the ${\mathcal{T}}$ spectrum (i.e. for ${\mathcal{T}}\ll Q$) at $\mathcal{O}(\alpha_{s}^{2}C_{F}T_{F})$ consist of the virtual (full QCD) contributions encoded in eqs. (4.1.1) and (4.1.1) and the secondary collinear and soft real radiation corrections contained in eqs. (4.2.2) and (4.3.3). Setting common scales $\mu=\mu_{H}=\mu_{B}=\mu_{S}$ yields for the corrections to virtual photon production

[TABLE]

where $h_{\rm virt}$ and $b_{\rm{real}}^{{\mathcal{T}}}$ are given in eqs. (4.1.1) and (160), respectively, and $s_{{\rm real}}^{{\mathcal{T}}}$ is given implicitly by the nondistributive terms in eq. (4.3.3). Again, for $Z$-boson production the anomalous axial current correction in eq. (4.1.1) has to be included in the $\delta({\mathcal{T}})$ term. Writing out the nontrivial terms in the spectrum we get

[TABLE]

where

[TABLE]

The singular fixed-order corrections for the ${\mathcal{T}}$ spectrum at $\mathcal{O}(\alpha_{s}^{2}T_{F}^{2})$ consist of the collinear real radiation corrections in eq. (69) for both beam directions,

[TABLE]

Depending on the hierarchy between $m$ and ${\mathcal{T}}$ and $Q$ some of the contributions in eqs. (C.2) and (C.2) are power-suppressed and therefore only appear via nonsingular corrections in the factorization formula for the associated parametric regime in sec. 2.

Appendix D Plus distributions

The standard plus distribution for some dimensionless function $g(x)$ is defined as

[TABLE]

The special case used in this paper is

[TABLE]

The 2-dimensional plus distributions that appear in the TMD beam and soft functions are defined as

[TABLE]

For the Fourier transform we use the convention

[TABLE]

The Fourier transforms of the 2-dimensional distributions required here are

[TABLE]

with $L_{b}$ defined in eq. (111).