Two-Loop Massive Quark Jet Functions in SCET

TL;DR

This paper computes second-order corrections to massive quark jet functions in SCET, clarifies their definitions, and explores their role in factorization approaches for accurate jet mass predictions involving massive quarks.

Contribution

It introduces a new way to define SCET jet functions using collinear-soft matrix elements, avoiding zero-bin subtractions, and compares two approaches for secondary massive quark effects.

Findings

Calculated $ ext{O}( ext{α}_s^2)$ corrections to massive quark jet functions.

Established the equivalence of two factorization approaches for secondary quark effects.

Addressed technical issues like rapidity divergences and infrared regularization.

Abstract

We calculate the $\mathcal O(\alpha_s^2)$ corrections to the primary massive quark jet functions in Soft-Collinear Effective Theory (SCET). They are an important ingredient in factorized predictions for inclusive jet mass cross sections initiated by massive quarks emerging from a hard interaction with smooth quark mass dependence. Due to the effects coming from the secondary production of massive quark-antiquark pairs there are two options to define the SCET jet function, which we call universal and mass mode jet functions. They are related to whether or not a soft mass mode (zero) bin subtraction is applied for the secondary massive quark contributions and differ in particular concerning the infrared behavior for vanishing quark mass. We advocate that a useful alternative to the common zero-bin subtraction concept is to define the SCET jet functions through subtractions related to collinear-soft matrix elements. This avoids the need to impose additional power counting arguments as required for zero-bin subtractions. We demonstrate how the two SCET jet function definitions may be used in the context of two recently developed factorization approaches to treat secondary massive quark effects. We clarify the relation between these approaches and in which way they are equivalent. Our two-loop calculation involves interesting technical subtleties related to spurious rapidity divergences and infrared regularization in the presence of massive quarks.

Equations

$${\mathcal{J}}_f^{(n_{\ell }+1)}(p^2,m^2)$$

$$=\delta (p^2-m^2)+\mathcal{O}({\alpha }_s).$$

$${\chi }_f^n(x)=W_n^{†}(x){\xi }_f^n(x),$$

$$W_{n}(x)=\sum_{\mathrm{perms}}\exp\bigg{(}\frac{-g}{\bar{n}\cdot\mathcal{P}}\,\bar{n}\cdot A_{n}(x)\bigg{)}\,.$$

$${\mathcal{J}}_f^{(n_{\ell }+1)}$$

$$\displaystyle=\mathrm{Im}\bigg{[}\frac{i}{2\pi N_{c}\,(\bar{n}\cdot p)}\int\!\mathrm{d}^{4}x\;e^{-ip\cdot x}\,\mathrm{Tr}\,\Big{\langle}0\Big{|}\mathrm{T}\bigg{\{}\frac{\bar{n}\!\!\!/}{2}W^{\dagger}(0)\,\psi_{f}(0)\,\overline{\psi}_{f}(x)\,W(x)\bigg{\}}\Big{|}0\Big{\rangle}\bigg{]}\,,$$

$$W(x)$$

$$J_f^{(n_{\ell }+1)}(p^2,0)\equiv$$

$$\mathscr{S}^{(n_{\ell}+1)}(\ell,m)=\frac{1}{N_{c}}\mathrm{Tr}\left\langle 0\left|\overline{\mathrm{T}}\bigl{[}V_{n}^{(0)\dagger}(0)X_{n}^{(0)}(0)\bigr{]}\,\delta(\ell-\hat{p}^{+})\,\mathrm{T}\bigl{[}X_{n}^{(0)\dagger}(0)V_{n}^{(0)}(0)\bigr{]}\right|0\right\rangle\,,$$

$$X_n(x)=$$

$$V_n(x)=$$

$$\int\mathrm{d}\ell^{\prime}\,\mathscr{S}^{(n_{\ell}+1)}(\ell-\ell^{\prime},m)\Bigl{(}\mathscr{S}^{(n_{\ell}+1)}(\ell^{\prime},m)\Bigr{)}^{-1}=\delta(\ell)\,,$$

$$J_{f}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},m^{2})\equiv\int\mathrm{d}\ell\,\!\mathscr{J}_{f}^{(n_{\ell}+1)}(p^{2}-(\bar{n}\cdot p)\ell,m^{2})\Bigl{(}\mathscr{S}^{(n_{\ell}+1)}(\ell,m^{2})\Bigr{)}^{-1}\,.$$

$$J_{f}^{\mathrm{mf},(n_{\ell}+1)}(p^{2},m^{2})=\int\mathrm{d}\ell\,\!\mathscr{J}_{f}^{(n_{\ell}+1)}(p^{2}-(\bar{n}\cdot p)\ell,m^{2})\Bigl{(}\mathscr{S}^{(n_{\ell})}(\ell)\Bigr{)}^{-1},$$

$$S_c(\ell ,m,\mu ,\nu )=$$

$$=$$

$$a_s^{(n_f)}\equiv \frac{{\alpha }_s^{(n_f)}(\mu )}{4\pi }\text{and}s\equiv p^2-m^2,$$

$$L_m\equiv \log \frac{m^2}{{\mu }^2},$$

$${\mathcal{L}}_n(x)\equiv \frac{1}{{\mu }^{\lambda }}{\left[\frac{\Theta (x){\log }^n(x/{\mu }^{\lambda })}{x/{\mu }^{\lambda }}\right]}_{+},$$

$$y\equiv \frac{m^2}{p^2},\bar{y}\equiv \frac{s}{p^2}=1-y,L_y\equiv \log (y),L_{\bar{y}}\equiv \log (\bar{y}).$$

$$Q\equiv (\bar{n}\cdot p)=p^{-},$$

$$J_Q^{\mathrm{mf},(n_{\ell }+1)}(p^2,m^2,\mu ,Q/\nu )$$

$$\displaystyle+\Bigl{(}a_{s}^{(n_{\ell})}\Bigr{)}^{2}\left(J_{Q}^{(2)}(p^{2},m^{2},\mu)+J_{Q,\mathrm{sec}}^{(2),\mathrm{mf}}(p^{2},m^{2},\mu,Q/\nu)\right)+\mathcal{O}\left(\alpha_{s}^{3}\right)\,,$$

$$\displaystyle J_{Q}^{(1)}(p^{2},m^{2},\mu)=C_{F}\Bigg{\{}\left(2L_{m}^{2}+L_{m}-\frac{\pi^{2}}{3}+8\right)\delta(s)+(-4L_{m}-4)\mathcal{L}_{0}(s)+8\mathcal{L}_{1}(s)$$

$$+G_1\Theta (p^2-m^2)$$

$$\displaystyle+\varepsilon\Bigg{[}\left(\frac{3}{2}L_{m}^{2}-\frac{2\pi^{2}}{3}L_{m}+16-\frac{5\pi^{2}}{12}-\frac{16\zeta_{3}}{3}\right)\delta(s)+\left(-2L_{m}^{2}-4L_{m}-8+\pi^{2}\right)\mathcal{L}_{0}(s)$$

$$\displaystyle\quad+\left(8+8L_{m}\right)\mathcal{L}_{1}(s)-8\mathcal{L}_{2}(s)+G_{1}^{(\varepsilon)}\Theta(p^{2}-m^{2})\Bigg{]}$$

$$\displaystyle+\varepsilon^{2}\Bigg{[}\Bigg{(}\frac{1}{6}L_{m}^{4}+\frac{1}{6}L_{m}^{3}+\left(4-\frac{\pi^{2}}{6}\right)L_{m}^{2}+\left(-\frac{7\pi^{2}}{12}-4\zeta_{3}\right)L_{m}+32-\frac{2\pi^{2}}{3}-\frac{\pi^{4}}{120}-5\zeta_{3}\Bigg{)}\delta(s)$$

$$+\left(-\frac{2}{3}{L}_m^3-2L_m^2+(-8+{\pi }^2)L_m-16+{\pi }^2+\frac{28{\zeta }_3}{3}\right){\mathcal{L}}_0(s)$$

$$\displaystyle\;+\left(4L_{m}^{2}+8L_{m}+16-2\pi^{2}\right)\mathcal{L}_{1}(s)+(-8-8L_{m})\mathcal{L}_{2}(s)+\frac{16}{3}\mathcal{L}_{3}(s)+G_{1}^{(\varepsilon^{2})}\Theta(p^{2}-m^{2})\Bigg{]}\Bigg{\}}.$$

$$G_1=$$

$$G_1^{(\epsilon )}=$$

$$=$$

$$G_1^{({\epsilon }^2)}=$$

$$+L_y^2\left(\frac{1}{2}\left(y^2-2y+13\right)-12L_{\bar{y}}\right)+2{\bar{y}}^2{L}_{\bar{y}}^2-2{\bar{y}}^2{L}_{\bar{y}}-\frac{1}{4}\left({\pi }^2-8\right){\bar{y}}^2+\frac{14}{3}{L}_y^3$$

$$\displaystyle+L_{y}\Big{(}-2\left(y^{2}-2y+5\right)L_{\bar{y}}+8L_{\bar{y}}^{2}+y^{2}-2y-\pi^{2}+9\Big{)}\bigg{]}$$

$$=$$

$$\displaystyle J_{Q}^{(2)}(p^{2},m^{2},\mu)=C_{F}^{2}\Bigg{\{}\Bigg{[}2L_{m}^{4}+2L_{m}^{3}+\left(\frac{33}{2}-2\pi^{2}\right)L_{m}^{2}+\left(-56\zeta_{3}+\frac{13}{2}-\pi^{2}\right)L_{m}$$

$$\displaystyle\qquad-38\zeta_{3}-\frac{\pi^{4}}{2}+\frac{433}{8}+\pi^{2}\left(-\frac{1}{3}-8\log 2\right)\Bigg{]}\delta(s)$$

Full text

††preprint:

UWThPh-2019-13

MITP/19-025

1 Introduction

Gaining higher precision in theoretical predictions for the production of massive quarks in hard particle collisions represents an important field of research in the context of the LHC as well as future colliders (see e.g. refs. Bendavid:2018nar ; Abramowicz:2018rjq ; Azzi:2019yne ). Factorized predictions are of special relevance since they provide a separation of physical effects from different momentum scales for cases where the scale hierarchies are large, such as for kinematic edges or endpoint regions or when there are hierarchies between particle masses and dynamical scales. Once factorization is established for a given observable, it can be written as a product or a convolution of functions encoding the dynamics of particular phase space regions. The resulting factorization formulae (or theorems) provide an approximation in the limit where one can expand in the ratios of the hierarchical scales. This allows to sum large logarithms related to these scales and in addition may provide the basis for a field theoretic treatment of low-energy hadronization effects for processes dominated by the strong interaction. As such, factorization theorems provide important pieces of information that are not accessible directly in corresponding calculations obtained in fixed-order perturbation theory within full QCD. Factorization also entails that the various functions occurring in the factorization theorems frequently have a universal character and can be applied for different processes.

For the description of the dynamics related to energetic radiation that is emitted collinearly to a quark produced at high energy, the so-called quark jet function is an essential ingredient in factorization theorems for a variety of processes where the inclusive invariant mass of all the collinear radiation (including the energetic quark) enters. The quark jet function was first introduced in the context of the factorization framework of the Soft-Collinear Effective Theory (SCET) Bauer:2000yr ; Bauer:2001yt ; Bauer:2001ct for the theoretical description of inclusive semileptonic or radiative $B$ meson decays in the kinematic regions where the produced hadrons form a jet. For the case of massless quarks the jet function is defined as a non-local correlation function of two SCET massless quark jet fields. It was calculated at ${\cal O}(\alpha_{s})$ and ${\cal O}(\alpha_{s}^{2})$ in refs. Lunghi:2002ju and Becher:2006qw , respectively, and recently also the ${\mathcal{O}}(\alpha_{s}^{3})$ corrections became available in ref. Bruser:2018rad . These results are applicable for the treatment of light quarks but can also be applied for massive quarks as long as the jet invariant mass $\sqrt{p^{2}}$ is much larger than the quark mass $m$, i.e. $p^{2}\gg m^{2}$, and the associated mass corrections are negligible.

When quark mass effects are considered, it is useful to distinguish two types of mass effects. One is related to the quark produced by the hard interaction and which initiates the jet, called primary, and the other is related to quark-antiquark vacuum polarization effects, called secondary. Interestingly, primary as well as secondary quark mass effects introduce new kinds of subtleties.

For an inclusive jet initiated by a primary heavy quark and with invariant mass $p^{2}$, two kinematic regions where quark mass effects are important emerge: $p^{2}\sim p^{2}-m^{2}\sim m^{2}$ and $p^{2}-m^{2}\ll m^{2}\approx p^{2}$. The region $p^{2}\sim p^{2}-m^{2}\sim m^{2}$, called SCET region, is relevant for processes where the jet invariant mass is close to the primary quark mass $m$ but is also allowed to fluctuate at the same level. The corresponding jet function(s) can be formulated with SCET massive collinear quark jet fields and require that the massless quark SCET formalism is extended to account for mass effects modifying the propagation and interactions of collinear quarks Leibovich:2003jd ; Rothstein:2003wh ; Fleming:2007qr . Jet functions of this kind are therefore called SCET jet functions. In the context of bottom and charm quark production in hard collisions the SCET regime $p^{2}\sim p^{2}-m^{2}\sim m^{2}$ is essentially the only one that can ever arise in practical applications where quark mass effects are important due to the sizable momentum smearing coming from non-perturbative effects Dehnadi:2016snl . The region $p^{2}-m^{2}\ll m^{2}\approx p^{2}$, also called bHQET limit or heavy quark limit, is relevant for processes where the jet invariant mass is very close to the primary quark mass $m$ and only allowed to fluctuate at a scale much smaller than $m$. Here, the SCET description does not provide a full separation of all dynamical effects because $\hat{s}\equiv(p^{2}-m^{2})/m\ll m\approx\sqrt{p^{2}}$ emerges as a new relevant scale. This separation is achieved in the context of boosted Heavy Quark Effective Theory (bHQET) where the corresponding collinear SCET sector is matched on a theory containing a super-heavy quark with its small offshell field component being integrated out Fleming:2007qr ; Fleming:2007xt . Jet functions of this kind are therefore called bHQET jet functions. In practice, the bHQET limit $p^{2}-m^{2}\ll m^{2}\approx p^{2}$ can only arise in the context of top quark production due to the large size of the top quark mass. Here, scales below the top quark mass may be resolved and not smeared out by hadronization effects Dehnadi:2016snl ; Dehnadi:unpub .111For similar discussions in the context of groomed jets see e.g. refs. Hoang:2017kmk ; Makris:2018npl ; Lee:2019lge . The bHQET regime is particularly important for the theoretical understanding of top quark mass determinations from observables tied to kinematic regions with top invariant masses close to the reconstructed top quark resonance Butenschoen:2016lpz ; Hoang:2017kmk ; Hoang:2018zrp . The ${\cal O}(\alpha_{s})$ primary massive quark SCET and the bHQET jet functions were calculated in ref. Fleming:2007xt . The ${\cal O}(\alpha_{s}^{2})$ corrections to the bHQET jet function were computed in ref. Jain:2008gb .

Secondary massive quark effects in a jet function occur at ${\cal O}(\alpha_{s}^{2})$ (NNLO) in the fixed-order expansion and thus only become relevant in high precision predictions. Their treatment in a factorization approach that separates collinear and soft quantum fluctuations, however, inevitably leads to $\mathrm{SCET_{II}}$-type rapidity singularities and associated large logarithms in factorization theorems for physical cross sections Gritschacher:2013pha . These large logarithms contribute at the next-to leading logarithmic (NLL) order. They arise from modes of equal virtuality (i.e. invariant mass $k^{2}\sim k_{+}k_{-}\sim m^{2}$) but widely separated rapidity (i.e. the size of the ratio $k_{+}/k_{-}$), see e.g. ref. Chiu:2012ir . The treatment of these rapidity singularities in the context of jet functions is tied to the definition of jet functions being infrared finite (matching) functions describing collinear fluctuations with virtuality of order $p^{2}$. Traditionally the required infrared (IR) subtractions are associated to so-called zero-bin subtractions Manohar:2006nz (that can be understood from the point of view of the SCET label formalism Bauer:2001ct , where the collinear label momentum zero is removed as it may not be power-counted as containing collinear modes) or are simply ignored Becher:2006qw (since frequently they are related to scaleless integrals that vanish identically in dimensional regularization). However, when secondary massive quark effects (or, conceptually equivalent, the exchange of massive gauge bosons Gritschacher:2013pha ) are considered, these subtractions are associated to non-vanishing contributions that need to be computed and specified in a well-defined and systematic fashion. Here, the concept of a zero-bin subtraction gets tied up with the concept of a soft mass mode bin subtraction Chiu:2009yx ; Gritschacher:2013pha ; Pietrulewicz:2014qza . The ${\cal O}(\alpha_{s}^{2})$ secondary massive quark corrections for the SCET jet function for primary massless quarks with imposing zero-bin as well as soft mass mode bin subtractions were calculated in ref. Pietrulewicz:2014qza . However, the difference between both types of subtractions makes the treatment of secondary quark mass effects subtle since the mass mode subtractions may or may not be associated to lower virtuality modes that are, depending on the application, located in separated phase space regions (i.e. having different power counting).

We propose that an alternative and in fact more transparent view is to associate the subtractions needed to define the jet functions to so-called collinear-soft matrix elements Bauer:2011uc ; Procura:2014cba , an idea that was already explored in different contexts in refs. Lee:2006nr ; Idilbi:2007ff . This leads to jet function results that agree with the results based on the zero-bin and mass mode bin subtraction concepts, but avoids the need to impose power counting arguments to determine whether a bin subtraction for a given diagram is relevant. Interestingly, there are two alternative ways to define the subtractions which entail two options to define jet functions, both of which differ analytically once secondary quark mass effects at ${\cal O}(\alpha_{s}^{2})$ are accounted for. Both options differ in whether or not the collinear-soft matrix element used for subtractions accounts for secondary massive quark corrections. The former option leads to what we call “universal” jet functions which merge into the well known massless quark jet functions in the limit $m\to 0$ and which were already introduced in refs. Pietrulewicz:2014qza ; Hoang:2015iva . The latter option leads to what we call “mass mode” jet functions, which are divergent in the limit $m\to 0$ and contain (universal) rapidity singularities. The mass mode jet functions are analogous to a corresponding definition for invariant mass dependent beam functions introduced in ref. Pietrulewicz:2017gxc . Both types of jet functions are useful to gain a transparent view on the treatment of secondary massive quark effects and can be employed in practical applications. They are also useful to better understand the relation between the secondary massive quark factorization approaches provided in refs. Gritschacher:2013pha ; Pietrulewicz:2014qza ; Hoang:2015iva and in ref. Pietrulewicz:2017gxc which lead to equivalent results. We stress that the concepts of the “universal” and the “mass mode” jet functions also applies to the analogous situation when the effects of massive (e.g. electroweak) gauge bosons are accounted for. This is because the structure of the relevant dynamical modes is equivalent Gritschacher:2013pha . For the exchange of massive gauge bosons, however, the mass effects arise already at ${\cal O}(\alpha_{em})$, $\alpha_{em}$ being the electromagnetic coupling, which will be explored elsewhere. Furthermore, the same concepts can also be applied for invariant mass dependent beam functions Samitz:thesis .

The main aim of this paper is to present the ${\cal O}(\alpha_{s}^{2})$ results (and some details of their computation) for the primary massive quark SCET jet functions. They represent important ingredients in factorization theorems for massive quark production within an inclusive jet where the massive quark is produced by the hard interaction and where the jet masses are similar in size to the mass of the quark. To illustrate the application of the universal and the mass mode jet functions we also discuss their role in the equivalent factorization approaches of refs. Gritschacher:2013pha ; Pietrulewicz:2014qza and Pietrulewicz:2017gxc using as an example the double differential hemisphere mass distribution in $e^{+}e^{-}$ annihilation (for the hard production of a boosted massive quark-antiquark pair). This discussion also clarifies the relation between both approaches and furthermore emphasizes that a smooth dependence on the value of the quark mass $m$ (in the sense that assumptions on the hierarchy of $m$ with respect to the other physical scales are not mandatory) can be achieved with them. The latter issue was already fully accounted for in refs. Gritschacher:2013pha ; Pietrulewicz:2014qza ; Hoang:2015iva but not explored in ref. Pietrulewicz:2017gxc . On the other hand, in ref. Pietrulewicz:2017gxc a more systematic representation of the building blocks to sum virtuality and rapidity logarithms on an equal footing was provided.

The content of the paper is organized as follows. In section 2 we set up our notation and provide the definition of the universal and the mass mode SCET jet functions in the context of subtractions based on collinear-soft matrix elements. In section 3 we present our results for both SCET jet functions for massive primary quarks at ${\cal O}(\alpha_{s}^{2})$ and show that (and in which way) the results are consistent with respect to available ${\cal O}(\alpha_{s}^{2})$ results in the limit of massless quarks Becher:2006qw and the bHQET limit Jain:2008gb . The practical use of the jet functions in the context of the factorization approaches of refs. Gritschacher:2013pha ; Pietrulewicz:2014qza ; Hoang:2015iva and of ref. Pietrulewicz:2017gxc is discussed in section 4 for the example of the double differential hemisphere mass distribution in $e^{+}e^{-}$ annihilation. The discussion should be sufficient to illustrate how the jet functions can be used for other processes in the context of QCD as well as the electroweak theory. Here, we also show how both approaches, considered together, are useful to gain a clearer view on the summation of virtuality and rapidity logarithms, the structure of collinear and soft mass mode corrections and on how a smooth dependence on the quark mass $m$ can be achieved in practical applications. Section 5 contains details of the calculations for the ${\cal O}(\alpha_{s}^{2})$ corrections to the massive primary quark SCET jet functions. This section is self-contained and focuses on computational details because the calculations are involved and require specific methods for different types of diagrams that may be useful for similar future work. The readers not interested in these technical details may skip this section except for section 5.1 which summarizes the results for the jet-function and collinear-soft matrix elements with different infrared regulators. In section 6 we provide a brief numerical analysis of the universal primary massive quark SCET jet function at ${\cal O}(\alpha_{s}^{2})$ in comparison to the corresponding known result in the massless quark and the bHQET limits. Finally, in section 7 we conclude. The paper also contains a number of appendices. In appendix A we provide the expressions for the universal and the mass mode jet functions at ${\cal O}(\alpha_{s}^{2})$ for massless primary quarks for the sake of future use. The latter can be extracted from ref. Pietrulewicz:2014qza but was not given there. The massless limit of the non-distributional corrections of the primary massive SCET jet function is provided in appendix B, and a number of virtuality and rapidity anomalous dimensions used in the discussion of section 4 are collected in appendix C. Finally, in appendix D we provide the results for the two-loop master integrals needed for our calculations.

2 Jet Function Definitions and Notation

Quark jet functions in SCET are based on the vacuum correlator of two SCET jet fields, encoding the inclusive collinear dynamics of quark fields coherently accounting for the collinear gluon radiation from all other color sources of a process. We work in QCD with $n_{\ell}$ massless quark flavors and one quark flavor $Q$ with mass $m$. The SCET jet-function matrix element (ME) for a jet initiated by the primary quark $f$ in light-like direction $n^{\mu}$ ($\vec{n}^{2}=1$) is given by Bauer:2001yt ; Fleming:2007qr 222 We note that in this paper we write jet functions and jet-function MEs as functions of the invariant mass $p^{2}\geq m^{2}$ in order to make all dependence on the quark mass $m$ explicit.

[TABLE]

We decompose four-vectors in lightcone components according to $p^{\mu}=p^{-}n^{\mu}/2+p^{+}\bar{n}^{\mu}/2+p_{\perp}^{\mu}$ with $n^{2}=\bar{n}^{2}=0$, $\bar{n}\cdot n=2$, and $n\cdot p_{\perp}=\bar{n}\cdot p_{\perp}=0$. The SCET label momentum operators $\overline{\mathcal{P}}_{n}$ and ${\mathcal{P}}_{n\perp}$ yield the sum of the large minus and perpendicular light-cone momentum labels of the fields on their right, respectively Bauer:2001ct , while $\hat{p}^{+}$ is the momentum operator of the small plus momentum. The trace is over color ($N_{c}=3$) and spinor indices. The SCET jet fields $\chi_{f}^{n}(x)$ are defined by

[TABLE]

where $\xi_{f}^{n}(x)$ is the SCET primary quark field and the $n$-collinear Wilson line can be written as

[TABLE]

The Wilson lines in eq. (1) encode the universal coherent dynamics of collinear gluons coming from other color sources, which are boosted in direction $\vec{\bar{n}}=-\vec{n}$ with respect to the direction of motion of the primary quark. They entail that the jet-function ME is gauge-invariant with respect to collinear gauge transformations Bauer:2000yr .

After BPS field redefinition the collinear sectors of the (leading-order) SCET Lagrangian are equivalent to boosted versions of the QCD Lagrangian Bauer:2000yr ; Bauer:2001yt . We can thus express the $n$-collinear jet-function ME also in terms of QCD quark fields as Becher:2006qw

[TABLE]

with the corresponding QCD definition for the $n$-collinear Wilson line

[TABLE]

In eq. (4) we have reexpressed the jet-function ME as the imaginary part of a jet field correlator. Taking the imaginary part is equivalent to half the discontinuity w.r.t. $s=p^{2}-m^{2}>0$. In practice, for the computation of the two-loop corrections it is more convenient to work with eq. (4) rather than with the corresponding SCET expression of eq. (1), because the QCD Feynman rules are simpler than the SCET Feynman rules and lead to less diagrams.

In eqs. (1) and (4) the quark flavor index $f$ stands for the flavor of the incoming primary quark and can either be one of the $n_{\ell}$ massless quarks, generically referred to as $q$, or the massive quark $Q$. The argument $m$ represents the dependence on the mass of quark $Q$ and arises either from secondary $Q$ effects (due to gluon splitting) if the primary quark is a massless flavor $q$ or from primary and secondary $Q$ effects. The superscript $(n_{\ell}+1)$ in eqs. (1) and (4) (and for ME definitions and factorization functions used throughout this paper) indicates that we consider the MEs in the context of QCD with $n_{\ell}$ massless quarks and the massive quark $Q$. We stress, however, that here and throughout this work the superscript $(n_{\ell}+1)$ (or $n_{\ell}$) on MEs or factorization functions does not automatically imply that also the corresponding flavor scheme for the $\overline{\rm MS}$ strong coupling is used. The latter can in principle be chosen independently, particularly for renormalization scales $\mu$ close to the quark mass $m$. We therefore specify explicitly the flavor scheme of $\alpha_{s}$ when we quote results parametrized by either $\alpha_{s}^{(n_{\ell})}$ or $\alpha_{s}^{(n_{\ell}+1)}$.

It is well known Fleming:2007xt ; Fleming:2007qr that the jet-function ME $\!\mathscr{J}_{f}^{(n_{\ell}+1)}(p^{2},m^{2})$ defined in eqs. (1) and (4) is per se not infrared-finite and commonly said to implicitly require as part of its definition zero-bin subtractions Manohar:2006nz to avoid double counting concerning the collinear overlap with softer or lower virtuality regions described by other functions in factorization theorems where the jet function appears. Alternatively, one can consider the jet function as an IR-finite matching function that describes fluctuations with virtualities $p^{2}-m^{2}\sim p^{2}\gtrsim m^{2}$ and that is defined by the jet-function ME $\!\mathscr{J}_{f}^{(n_{\ell}+1)}(p^{2},m^{2})$ with additional (i.e. independent) ME subtractions related to (soft and possibly collinear) modes at smaller virtualities. So far in the literature the required subtractions, when they contributed non-vanishing terms, were implemented mostly via the zero-bin Fleming:2007xt ; Fleming:2007qr or soft mass mode bin Gritschacher:2013pha ; Pietrulewicz:2014qza ; Gritschacher:2013tza prescriptions that are imposed by hand. However, from the point of view of being matching functions, the jet functions can also be defined in a conceptually more systematic way by explicitly dividing out fluctuations with momenta $n\cdot k=k^{+}\sim p^{2}/Q$ and virtualities $k^{2}\ll p^{2}-m^{2}$ in terms of a so-called collinear-soft ME. This also has the advantage that no additional power counting arguments are mandatory to identify which contributions have to be accounted for (and which not), as it is necessary when implementing a bin subtraction.

In case that all $(n_{\ell}+1)$ quarks are massless this results in the jet function definition

[TABLE]

Here $\mathscr{S}^{(n_{\ell}+1)}(\ell,m)$ is the collinear-soft ME of our $(n_{\ell}+1)$ flavor theory (with $n_{\ell}$ massless quarks $q$ and one quark $Q$ with mass $m$) and defined as333Note that the collinear-soft matrix element is universal in the sense that it only depends on the color representations of the partons involved as noted in ref. Bauer:2011uc . Therefore the result is equivalent to the one given in eq. (B.51) of ref. Pietrulewicz:2017gxc where massive quark effects in exclusive Drell-Yan were considered.

[TABLE]

with the collinear-soft Wilson lines Bauer:2011uc ; Procura:2014cba

[TABLE]

including the symmetric analytic rapidity regularization of ref. Chiu:2012ir . The inverse of the collinear-soft ME is defined by the convolution

[TABLE]

and we use the analogue relations for defining the inverse of all (factorization) functions depending on dynamical variables. The combination of MEs in eq. (6) yields the well known IR finite massless quark jet function at one Lunghi:2002ju , two Becher:2006qw and three loops Bruser:2018rad . The analogous definition (with the appropriate color representation for the Wilson lines) can be applied also for the gluon jet function Becher:2009th ; Becher:2010pd ; Banerjee:2018ozf . For $(n_{\ell}+1)$ massless quarks the subtraction of the low virtuality modes is from the purely computational point of view, however, not strictly mandatory for computations of the jet function when dimensional regularization for UV and IR divergences is used because then $\mathscr{S}^{(n_{\ell}+1)}(\ell,0)$ yields vanishing scaleless integrals beyond tree-level.

On the other hand, when considering the quark $Q$ having a finite mass $m$, $\mathscr{S}^{(n_{\ell}+1)}(\ell,m)$ is not anymore scaleless and yields non-vanishing ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ contributions in dimensional regularization due to mass mode fluctuations with virtualities $k^{2}\sim m^{2}$. We note that the mass mode contributions in the collinear-soft ME $\mathscr{S}^{(n_{\ell}+1)}(\ell,m)$ are only related to secondary effects of the massive quark $Q$ due to closed $Q$ loops, so that it is the same for the jet functions for a massless and a massive primary quark.

Interestingly, secondary massive quark effects do not lead to IR divergences (as long as one does not take the massless limit) that must always be subtracted. It is therefore not mandatory to use the massive quark analogue of eq. (6) to define the jet function. This is also related to the emergence of rapidity singularities, which one may not associate to be of either UV or IR type and leads to two options to define the SCET jet function in the context of QCD with $n_{\ell}$ massless flavors and one quark $Q$ with mass $m$.

The first option is just the analogue of eq. (6) but with the $(n_{\ell}+1)$st quark $Q$ having mass $m$,

[TABLE]

The primary jet-initiating quark flavor $f$ can either be one of the $n_{\ell}$ massless quarks $q$ or the massive quark $Q$. Here the collinear-soft ME is calculated in the full $(n_{\ell}+1)$ flavor theory and thus encodes in particular the secondary effects coming from the massive quark $Q$. This jet function yields by construction the massless quark jet function $J_{f}^{(n_{\ell}+1)}(p^{2},0)$ of eq. (6) in the limit $m\to 0$. $J_{f}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},m^{2})$ and $J_{f}^{(n_{\ell}+1)}(p^{2},0)$ have the same anomalous dimension and are furthermore free of rapidity divergences as well as of any large rapidity logarithms. For $J_{f}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},m^{2})$ both cancel between $\!\mathscr{J}_{f}^{(n_{\ell}+1)}$ and $\mathscr{S}^{(n_{\ell}+1)}$. Upon renormalizing $J_{f}^{\mathrm{uf},(n_{\ell}+1)}$ the only dependence on a renormalization scale is that on the virtuality scale $\mu$. Following the philosophy of the SCET jet function being a matching function describing fluctuations with virtuality $p^{2}-m^{2}\sim p^{2}$, the jet function of eq. (10) is the mandatory definition for the case where $m^{2}\ll p^{2}$ where the quark mass $m$ is an IR scale and mass effects represent power corrections. In that regime, however, one may as well use the massless quark jet function $J_{f}^{(n_{\ell}+1)}(p^{2},0)$. The jet function $J_{f}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},m^{2})$ can, however, also be employed for $p^{2}\sim p^{2}-m^{2}\sim m^{2}$ following the approach of refs. Gritschacher:2013pha ; Pietrulewicz:2014qza (originally developed for massless primary quark jet functions) where universal hard, jet and soft (and in principle also invariant mass dependent beam) functions can be defined valid for any value of $m^{2}$, with respect to the other dynamic scales and the renormalization scale $\mu$. This approach makes the $\mu$-evolution in $n_{\ell}$ or $(n_{\ell}+1)$ flavor schemes to sum virtuality logarithms particularly transparent and allows to formulate factorization theorems that smoothly interpolate all possible hierarchies as far as the value of $m$ with respect to the other kinematic scales is concerned. The approach of refs. Gritschacher:2013pha ; Pietrulewicz:2014qza also entails that rapidity singularities (and the respective resummation of logarithms) only arise within the mass threshold matching conditions for the individual hard, jet (as well as beam) and soft factorization functions. Since these threshold matching conditions consist of a few universal collinear and soft mass mode factors, there can be cancellations among them in the context of a factorization theorem related to consistency relations. So mass mode factors can arise in the matching conditions that are physically not relevant for a particular observable and eventually cancel in the corresponding factorization theorem. We call this the “universal factorization function” approach and therefore give the jet function definition in eq. (10) the superscript “uf”. In the following we call it the “universal jet function”.

The second option is to define the jet function in the presence of a massive quark by

[TABLE]

where the collinear-soft ME is calculated in the $n_{\ell}$ flavor theory, i.e. containing the $n_{\ell}$ massless quarks but not accounting for the massive quark $Q$. Like for the jet function $J_{f}^{\mathrm{uf},(n_{\ell}+1)}$, the primary jet-initiating quark $f$ can either be one of the $n_{\ell}$ massless quarks $q$ or the massive quark $Q$. The jet function $J_{f}^{\mathrm{mf},(n_{\ell}+1)}$ is IR divergent for $m\to 0$ and furthermore has rapidity singularities that are treated by rapidity renormalization. So upon renormalization, $J_{f}^{\mathrm{mf},(n_{\ell}+1)}$ depends on the (virtuality) renormalization scale $\mu$ and in addition on the (rapidity) renormalization scale $\nu$. Following the philosophy of the SCET jet function being a matching function describing the fluctuations with virtuality $p^{2}-m^{2}\sim p^{2}$, the jet function of eq. (11) is the natural option for the case where $m^{2}\sim p^{2}\sim p^{2}-m^{2}$, when the mass is not an IR scale. It has also been advocated for this scale hierarchy in ref. Pietrulewicz:2017gxc (originally developed for the massless primary quark and invariant mass dependent beam function Stewart:2009yx ). In the factorization approach of ref. Pietrulewicz:2017gxc separate factorization theorems are formulated for each possible hierarchy of the mass scale $m$ with respect to the other kinematic scales. The approach is not designed to provide smooth dependence on $m$, but is more economical (in the sense of being minimalistic) concerning the appearance of collinear and soft mass mode contributions so that no cancellations arise between them in the context of a factorization theorem. Furthermore, the collinear and soft mass mode contributions adopt the status of genuine factorization functions. The resulting structure of factorization theorems makes the structure of $\nu$-evolution to sum rapidity logarithms particularly transparent as it appears at the same level as the virtuality $\mu$-evolution. We call this the “mass mode factorization” approach and therefore give the jet function definition in eq. (11) the superscript “mf”. In the following we call it the “mass mode jet function”.

The difference between both $\overline{\rm MS}$-renormalized massive quark jet function definitions yields444The inverse jet function is defined in analogy to eq. (9).

[TABLE]

and is called the “collinear-soft function” in ref. Pietrulewicz:2017gxc . It is IR-finite for finite $m$ and has a dependence on the virtuality and the rapidity renormalization scales $\mu$ and $\nu$, respectively.

In section 4 we illustrate the practical use of both types of jet functions and the two factorization approaches of refs. Pietrulewicz:2014qza ; Gritschacher:2013pha and Pietrulewicz:2017gxc exemplarily for the double differential hemisphere mass distribution in $e^{+}e^{-}$ collisions for the hierarchies $m^{2}\sim p^{2}\sim p^{2}-m^{2}$ and $m^{2}\ll p^{2}$. We also stress that the form of the collinear-soft ME of eq. (7) and the two options of using either eq. (10) or eq. (11) for IR finite jet functions in the presence of massive quarks can also be applied in a fully equivalent way for invariant mass dependent beam functions Stewart:2009yx and the bHQET jet function Fleming:2007xt ; Fleming:2007qr .555In the context of bHQET Fleming:2007xt ; Fleming:2007qr we refer to the situation that some of the $n_{\ell}$ light flavors have finite masses smaller than the mass $m$ of the primary HQET quark, but larger than the hadronization scale $\Lambda_{\mathrm{QCD}}$. The equivalence of the collinear-soft MEs (or the zero-bin subtractions) for SCET and bHQET is a consequence of a rescaling property of the bHQET and SCET soft gluon eikonal couplings that can e.g. be trivially seen from the structure of the zero-bin diagrams in both effective theories. The approach also applies in an analogous way for the treatment of massive gauge bosons which is conceptually related to the issues of secondary quark mass effects as was discussed in Gritschacher:2013pha .

3 Analytic Results: Primary Massive Quark SCET Jet Functions

In this section we present the analytic results for the $\mathcal{O}(\alpha_{s}^{2})$ primary quark SCET jet functions. Details on the computation of the $\mathcal{O}(\alpha_{s}^{2})$ corrections, which are all new, are given in section 5. We use the abbreviations

[TABLE]

for the strong coupling in the $n_{f}$-flavor scheme and the squared jet function invariant mass $p^{2}$, respectively,

[TABLE]

as well as

[TABLE]

for plus-distributions for a variable $x$ having mass dimension $\lambda$, where $\mu$ is the common $\overline{\mathrm{MS}}$ renormalization scale. They are defined such that $\int_{0}^{1}\mathrm{d}x\,\mathcal{L}_{n}(x)=0$. Furthermore we define

[TABLE]

In the physical kinematic region we have $0\leq y,\bar{y}\leq 1$. Note that we also use the notation

[TABLE]

and that we carry out all computations in $d=4-2\varepsilon$ dimensions.

3.1 Mass Mode Jet Function

The result for the renormalized primary massive mass mode jet function $J_{Q}^{\mathrm{mf},(n_{\ell}+1)}$, defined in eq. (11), in the pole mass scheme at $\mathcal{O}(\alpha_{s}^{2})$ can be written in the form

[TABLE]

where we adopt the $n_{\ell}$ flavor scheme for the strong coupling. This can be considered the natural scheme choice for the mass mode jet function since it matches directly to the bHQET jet function, see section 3.3. It is of course straightforward to switch to the $(n_{\ell}+1)$ flavor scheme using eq. (42). At $\mathcal{O}(\alpha_{s}^{2})$ we distinguish between the corrections from diagrams containing only one single quark $Q$ line ($J_{Q}^{(2)}$) and those from diagrams containing the $Q\bar{Q}$ vacuum polarization correction ($J_{Q,\mathrm{sec}}^{(2)}$). The mass mode jet function contains rapidity singularities from the $Q\bar{Q}$ vacuum polarization diagrams which are renormalized following the approach of Chiu:2011qc ; Chiu:2012ir and which lead to the dependence on the rapidity renormalization scale $\nu$. For the prototypical application of the mass mode jet function within factorization theorems where all logarithms are resummed by explicit renormalization group evolution factors the natural choice for the rapidity renormalization scale is $\nu\sim Q$. This is indicated by the argument $Q/\nu$. Throughout this paper we adopt this convention for the rapidity scale $\nu$ in an analogous way. In case the natural choice for $\nu$ is process (or factorization theorem) dependent, only $\nu$ appears as the argument.

The $\mathcal{O}(\alpha_{s})$ contribution accounting for terms up to $\mathcal{O}(\varepsilon^{2})$ in order to maintain all information mandatory for $\mathcal{O}(\alpha_{s}^{2})$ calculations reads

[TABLE]

Here, all distributional contributions are displayed explicitly and the non-distributional contributions are contained in the functions

[TABLE]

These functions do not contribute to the singular behavior in the bHQET limit $s=p^{2}-m^{2}\ll m^{2}$ and thus represent power corrections in this kinematic region. They constitute the so-called $\mathcal{O}(\alpha_{s})$ bHQET non-singular corrections that are important for physical factorization predictions in the region $s\lesssim m^{2}$ to achieve a smooth transition between bHQET and SCET factorization theorems Fleming:2007qr ; Dehnadi:2016snl . For completeness we also quoted the $s\to 0$ limit of the $G_{1}$-functions. The $\varepsilon\to 0$ limit of eq. (19) was already given in ref. Fleming:2007xt , and we refer to this reference for pictures of the Feynman diagrams and details of the calculation. The determination of the ${\cal O}(\varepsilon)$ and ${\cal O}(\varepsilon^{2})$ contributions is straightforward.

The $\mathcal{O}(\alpha_{s}^{2})$ contributions arising from Feynman diagrams containing only one single quark $Q$ line (see all diagrams in figure 1 with the diagrams (o)-(r) containing only massless quark bubbles) read

[TABLE]

The distributional contributions are again displayed explicitly. There are two sets of non-distributional contributions. One set ($\propto\Theta(p^{2}-m^{2})$) is related to terms that are non-vanishing above the one-particle cut and the other ($\propto\Theta(p^{2}-(3m)^{2})$) is related to terms that are non-vanishing above the three particle $QQ\bar{Q}$ cut. Even though the term $J_{Q}^{(2)}$ exclusively arises from diagrams with one single quark $Q$ line, see figure 1, three particle $QQ\bar{Q}$ cuts can arise due to the antiparticle component of the quark $Q$ propagator. The non-distributional contributions $\propto\Theta(p^{2}-m^{2})$ read

[TABLE]

[TABLE]

The non-distributional contributions $\propto\Theta(p^{2}-(3m)^{2})$ read

[TABLE]

The non-distributional contributions $\propto\Theta(p^{2}-m^{2})$ coming from diagrams (b) and (c) in figure 1 could only be computed numerically for general values of $p^{2}$ and $m^{2}$, and are parametrized by the fit function $G_{\mathrm{fit}}$ which has the convenient form666We note that the function $G_{\mathrm{fit}}$ contains contributions arising both from $Q$ and $QQ\bar{Q}$ final states.

[TABLE]

with

[TABLE]

With the values of the fitting parameters $c_{ij}$ quoted in eq. (29) the fit function $G_{\mathrm{fit}}$ approximates the corresponding numerical results obtained in section 5.2.2 with a relative precision of better than $0.5\%$. Also accounting for the intrinsic uncertainty of the numerical results we reach an accuracy of at least $1\%$ for $G_{\mathrm{fit}}$. Besides that we note that the contributions parametrized in the fit function account for less then $0.5\%$ of the $\mathcal{O}(\alpha_{s}^{2})$ corrections for all physical values of $y=m^{2}/p^{2}$. So for all practical purposes the uncertainties associated to $G_{\mathrm{fit}}$ are negligible, and we did not quote uncertainties in eq. (29).

The $\mathcal{O}(\alpha_{s}^{2})$ contributions arising from Feynman diagrams containing the $Q\bar{Q}$ vacuum polarization subdiagrams (see diagrams (o)-(r) in figure 1 with only massive quark bubbles) read

[TABLE]

where again the distributional contributions are displayed explicitly. The dependence on the rapidity renormalization scale $\nu$ is displayed explicitly in terms of the logarithm $\log(Q/\nu)$. In the $(n_{\ell}+1)$ flavor scheme for the strong coupling the coefficient $J_{Q,\mathrm{sec}}^{(2),\mathrm{mf}}$ is the only one that is modified: $J_{Q,\mathrm{sec}}^{(2),\mathrm{mf}}\to J_{Q,\mathrm{sec}}^{(2),\mathrm{mf}}+4/3L_{m}J_{Q}^{(1)}$.

The non-distributional contributions arise entirely from the three particle $QQ\bar{Q}$ cut and its result reads

[TABLE]

The non-distributional three particle $QQ\bar{Q}$ cut contributions in $J_{Q}^{(2)}$ and $J_{Q,\mathrm{sec}}^{(2)}$ involve the elliptic functions

[TABLE]

and the integral functions

[TABLE]

which can be readily evaluated numerically. For related discussions see e.g. ref. Ablinger:2017bjx . In practical applications it is, however, useful to have an approximate representation of $G^{(3Q)}$ and $G_{\mathrm{sec}}^{(3Q)}$ involving only elementary functions. Such approximate representations are provided by

[TABLE]

[TABLE]

with $c_{1}^{C_{F}(3Q)}=-0.7298$, $c_{2}^{C_{F}(3Q)}=-0.2085$, $c^{T_{F}(3Q)}=-0.3183$. They both approximate the exact expressions to better than $0.5\%$ over the physical kinematic range and provide exact analytic results in the limit $p^{2}\to 9m^{2}$ and $m\to 0$. We also quoted a fully numerical version of the expressions useful for practical applications.

In analogy to the non-distributional terms at $\mathcal{O}(\alpha_{s})$ the non-distributional $G$-functions at $\mathcal{O}(\alpha_{s}^{2})$ do not contribute to the singular contributions in the bHQET limit $s=p^{2}-m^{2}\ll m^{2}$. They represent power corrections in this kinematic region and constitute the $\mathcal{O}(\alpha_{s}^{2})$ bHQET non-singular corrections. For completeness we also quoted the respective results in the limit $p^{2}\to m^{2}$ (eqs. (24) – (26)) and in the limit $p^{2}\to 9m^{2}$ (eqs. (27) and (31)).

For completeness we also provide the result for the renormalized mass mode jet function $J_{q}^{\mathrm{mf},(n_{\ell}+1)}(p^{2},m^{2},\mu,Q/\nu)$ for massless primary quarks in eq. (100), which can be extracted from the calculations in ref. Pietrulewicz:2014qza , but was not given in the literature before. The $\mu$ (virtuality) and $\nu$ (rapidity) anomalous dimensions of the mass mode jet functions for massive and massless primary quarks agree identically.

3.2 Universal Jet Function

The result for the renormalized primary massive universal jet function $J_{Q}^{\mathrm{uf},(n_{\ell}+1)}$, defined in eq. (10), in the pole mass scheme at $\mathcal{O}(\alpha_{s}^{2})$ can be written in the form

[TABLE]

where we adopt the $(n_{\ell}+1)$ flavor scheme for the strong coupling, which is the natural scheme choice for the universal jet function. The universal jet function does not contain rapidity divergences and thus only depends on the UV renormalization scale $\mu$ coming from dimensional regularization. The $\mathcal{O}(\alpha_{s})$ coefficient $J_{Q}^{(1)}$ and the $\mathcal{O}(\alpha_{s}^{2})$ contribution from diagrams containing only one single quark $Q$ line $J_{Q}^{(2)}$ (see all diagrams in figure 1 with the diagrams (o)-(r) containing only massless quark bubbles) are identical to the mass mode jet function and given in eqs. (19) and (23), respectively.

The $\mathcal{O}(\alpha_{s}^{2})$ contributions arising from Feynman diagrams containing the $Q\bar{Q}$ vacuum polarization subdiagrams (see diagrams (o)-(r) in figure 1 containing only massive quark bubbles) read

[TABLE]

where again the distributional contributions are displayed explicitly. The non-distributional contributions come from the single heavy quark $Q$ cut and from the three particle $QQ\bar{Q}$ cut. The $Q$ cut term is given in eq. (20) and arises from using the $(n_{\ell}+1)$ flavor scheme for the strong coupling. The $QQ\bar{Q}$ cut term $G_{\mathrm{sec}}^{(3Q)}$ is given in eq. (31).

The non-distributional $G$-functions represent power corrections in the bHQET limit $s=p^{2}-m^{2}\ll m^{2}$ and represent $\mathcal{O}(\alpha_{s}^{2})$ bHQET non-singular corrections. For completeness we also quoted the respective results in the limit $p^{2}\to m^{2}$ (eqs. (24) – (26)) and in the limit $p^{2}\to 9m^{2}$ (eqs. (27) and (31)).

For completeness we also quote the result for the renormalized universal jet function $J_{q}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},m^{2},\mu)$ for primary massless quarks in eq. (105), which was computed in ref. Pietrulewicz:2014qza . The renormalization $Z$-factor as well as the $\mu$ (virtuality) anomalous dimension of of the universal jet functions (to all orders as well as for primary massive and massless quarks) agree with those of the well known massless quark jet functions Becher:2006qw .

3.3 Consistency Checks and Kinematic Limits

There are a number of essential consistency properties the $\mathcal{O}(\alpha_{s}^{2})$ primary massive quark jet functions have to satisfy and which provide important consistency checks. These concern the relation between the mass mode and universal jet functions, the limit of massless quarks and the bHQET limit of a supermassive quark. We discuss them in the following.

As explained in section 2, the difference between the mass mode and the universal jet functions is related to the collinear-soft function $S_{c}$ Pietrulewicz:2017gxc via eq. (12). Using the results of eqs. (18) and (39) and accounting for the result for the renormalized collinear-soft function determined in ref. Pietrulewicz:2017gxc , which in our notation has the form

[TABLE]

it is straightforward to check the validity of eq. (12). Here one has to account for the decoupling relation

[TABLE]

for $\alpha_{s}$ in the $n_{\ell}$ and $(n_{\ell}+1)$ flavor schemes (where we also displayed the $\mathcal{O}(\varepsilon)$ and $\mathcal{O}(\varepsilon^{2})$ terms useful for divergent expressions). The sole dependence of the collinear-soft function on the rapidity scale $\nu$ indicates that the natural choice for $\nu$ in general depends on the process and the structure of the factorization theorem.

The definition of the universal jet function $J_{f}^{\mathrm{uf},(n_{\ell}+1)}$ given in eq. (10) entails that for $m\to 0$ it converges to the jet function $J_{f}^{(n_{\ell}+1)}$ for $(n_{\ell}+1)$ massless quarks, see eq. (6), which at $\mathcal{O}(\alpha_{s}^{2})$ was computed in ref. Becher:2006qw . In the limit $m\to 0$, the non-distributional $G$-functions also yield distributions and the corresponding limiting expressions are given in appendix B. Accounting for these results our expression for universal jet function correctly approaches the massless quark two-loop jet function for $m\to 0$.

Finally, we discuss the bHQET limit $s=p^{2}-m^{2}\ll m^{2}$, the jet mass threshold region. Taking the double differential hemisphere jet mass distribution in $e^{+}e^{-}$ annihilation (w.r.t. the thrust axis) as an example, one can consider the kinematic region where both jet masses are close to threshold. The corresponding factorization theorem for boosted top quarks ($Q\gg m$) has the form Fleming:2007qr ; Hoang:2015vua 777Note that throughout this paper we do not account for renormalization group functions to sum (virtuality or rapidity) logarithms in our discussions of factorization theorems in order to keep the notations brief.

[TABLE]

where $M_{Q,\bar{Q}}$ stand for the hemisphere masses, $\sigma_{0}$ denotes the tree level cross section for $e^{+}e^{-}\to Q\bar{Q}$. The term $H_{Q}^{(n_{\ell}+1)}$ is the hard function related to the matching of the QCD $Q\bar{Q}$ current to the SCET massive quark dijet current at the scale $Q$ and the term $H_{m}$ is the mass threshold hard function related to the matching of the SCET dijet current to the bHQET current at the scale $m$. The terms $J_{B}^{(n_{\ell})}$ and $S^{(n_{\ell})}$ denote the bHQET jet function Fleming:2007qr ; Jain:2008gb and dijet soft function, respectively. The soft function $S^{(n_{\ell})}$ describes ultrasoft cross talk between the two jets at the scale $s_{Q}/Q$. The bHQET jet function $J_{B}^{(n_{\ell})}$ Fleming:2007qr ; Jain:2008gb describes the effects of the ultra-collinear radiation inside the jets (i.e. soft radiation in the rest frame of the massive quarks) at the scales $s_{Q,\bar{Q}}/m$ and contain all the dynamical (ultra) collinear effects for $s_{Q,\bar{Q}}\ll m^{2}$. The mass threshold hard function $H_{m}(m,Q/m,\mu)$ consists of $n$-collinear, $\bar{n}$-collinear and soft mass mode contributions which can be written in factorized form as Hoang:2015vua

[TABLE]

Each of the collinear ($n,\bar{n}$) and soft ($s$) mass mode factors depends on the rapidity renormalization scale $\nu$. The natural scaling for the collinear factors is $\nu\sim Q$ and for the soft factor it is $\nu\sim m$. They were calculated at $\mathcal{O}(\alpha_{s}^{2})$ in ref. Hoang:2015vua (see eqs. (5.7) and (5.8) in ref. Hoang:2015vua and note that $H_{m,n}=H_{m,\bar{n}}$ for the symmetric regulator of ref. Chiu:2012ir , see eq. (8)).

Up to power-suppressed (for $s_{Q,\bar{Q}}\ll m^{2}$) and non-distributional contributions the factorization theorem in eq. (43) is also valid in the kinematic region $s_{Q,\bar{Q}}=p_{Q,\bar{Q}}^{2}-m^{2}\sim m^{2}\lesssim p_{Q,\bar{Q}}^{2}$. Interestingly, in ref. Pietrulewicz:2017gxc it was pointed out in the context of invariant mass dependent beam functions for exclusive Drell-Yan gauge boson production, that for this kinematic region the collinear-soft function $S_{c}$ describes dynamical fluctuations located in the same phase space region as those contained in the jet function. The most economical formulation (w.r.t. the number of factorization functions containing the dynamical effects of the collinear and soft mass mode contributions) for $s_{Q,\bar{Q}}=p_{Q,\bar{Q}}^{2}-m^{2}\sim m^{2}\sim p^{2}$ is therefore the one where the contributions of the collinear-soft function are included in the definition of the (invariant mass dependent) beam function, which corresponds to the mass mode definition of eq. (11). For the double differential hemisphere jet function this implies that for $s_{Q,\bar{Q}}=p_{Q,\bar{Q}}^{2}-m^{2}\sim m^{2}\sim p_{Q,\bar{Q}}^{2}$ one can formulate a factorization theorem of the form

[TABLE]

This in turn implies the relation

[TABLE]

between the primary massive quark mass mode jet function and the bHQET jet function up to non-singular terms integrable in $p^{2}$ at the threshold $m^{2}$. Using the $\mathcal{O}(\alpha_{s}^{2})$ results for the bHQET jet function (naturally given in the $n_{\ell}$-flavor scheme for $\alpha_{s}$) in ref. Jain:2008gb (see eq. (39)), the collinear mass mode hard functions $H_{m,n}$ in ref. Hoang:2015vua (see eq. (5.7)) and our result for the mass mode jet function in eq. (11), accounting for the fact that all $G$-functions are power-suppressed in the bHQET limit, it is straightforward to see that the relation in eq. (46) is indeed satisfied. As discussed in more detail in section 5.2.2, our calculations of the ${\cal O}(\alpha_{s}^{2})$ corrections for the mass mode jet function provide a high precision semi-analytic check of eq. (46) as some loop integrals could only be determined by numerical methods. The relation (46) was then used as an input to analytically parametrize the jet functions result of eqs. (11) and (10).

4 Alternative Factorization Approaches and Practical Use

The non-trivial conceptual aspect of the $\mathcal{O}(\alpha_{s}^{2})$ primary massive quark SCET jet function concerns the mass mode contribution due to secondary radiation. Since these are tied to corrections arising from collinear as well as soft phase space sectors with invariant masses of order the quark mass $m$, their treatment entails rapidity singularities (and the resummation of associated logarithms) within a $\mathrm{SCET_{II}}$ type approach regardless of whether the observable is $\mathrm{SCET_{I}}$ or $\mathrm{SCET_{II}}$. As already pointed out in section 2, in the literature two alternative factorization formulations to account for mass mode corrections have been advocated which are (or can be rendered) conceptually and phenomenologically equivalent due to consistency relations. We emphasize that none of these approaches may be considered superior to the respective other and, by emphasizing different and complementary aspects of the mass mode effects, together provide a thorough and comprehensive treatment and understanding from the conceptual as well as practical point of view.

In this section we discuss the relation between the approaches employing the mass mode and the universal jet functions for primary massive quarks and using the double differential hemisphere mass distribution already introduced in section 3.3 as an example. To keep the discussion brief we also restrict the discussion to the kinematic regions $s^{2}/Q^{2}\ll m^{2}\ll s\sim p^{2}\ll Q^{2}$ and $s=p^{2}-m^{2}\sim p^{2}\sim m^{2}$ where the quark mass effects calculated in this paper are most relevant. We emphasize, however, that the discussion can be extended in a straightforward way to all kinematic situations and also to other types of factorization theorems where mass mode effects shall be treated (and where factorization functions other than jet functions arise).

The universal factorization approach of refs. Gritschacher:2013pha ; Gritschacher:2013tza ; Pietrulewicz:2014qza ; Hoang:2015iva was devised with the motivation to obtain universal hard, jet and soft functions defined such that they can smoothly cover all possible choices for the quark mass $m$ in the context of the strict hierarchy $s^{2}/Q^{2}\ll s\ll Q^{2}$. The resulting factorization theorems therefore contain power corrections w.r.t. mass mode contributions in at least some of the factors in any kinematic region, which, in case they are known, may be expanded away accordingly if a smooth description of mass effects is not wanted. The approach was also constructed to allow for the formulation of universal and simple rules to implement the summation of logarithms through flavor number dependent renormalization evolution factors of the hard, jet and soft function for arbitrary choices of the global renormalization scale $\mu$. The approach also entails that only the hard, jet and soft functions and, depending on the choice of $\mu$, their respective threshold matching factors appear in the factorization theorems (in case they are evolved through the massive quark flavor threshold). By construction in the universal factorization approach, all rapidity logarithms (and their summations) are fully contained within these individual quark mass threshold matching factors at the scale $\mu_{m}\sim m$.

Let us consider the situation $s^{2}/Q^{2}\ll m^{2}\lesssim s\sim p^{2}\ll Q^{2}$. The corresponding leading power factorization theorem in the universal factorization function approach for $s^{2}/Q^{2}<\mu^{2}<m^{2}$ (when the SCET current and the jet functions are evolved below $m$) has the form

[TABLE]

where $H_{m}$ and $\mathcal{M}_{J_{Q}}$ are the quark mass threshold matching factors for the SCET current and the universal jet function $J_{Q}^{\mathrm{uf}}$ respectively. For $m^{2}<\mu^{2}<s$ (where the soft function is evolved above $m$) the same factorization theorem has the form

[TABLE]

where $\mathcal{M}_{S}$ is the mass threshold matching factor for the soft function $S$. It is straightforward to consider other choices of $\mu$ as well.

For a strong hierarchy $m^{2}\ll s$ one may replace the universal primary massive quark jet function $J_{Q}^{\mathrm{uf},(n_{\ell}+1)}$ by the jet function $J_{q}^{(n_{\ell}+1)}$ for $(n_{\ell}+1)$ massless quark flavors (calculated at $\mathcal{O}(\alpha_{s}^{2})$ in Becher:2006qw ). But using $J_{Q}^{\mathrm{uf},(n_{\ell}+1)}$ provides validity of eqs. (4) and (4) smoothly covering all cases between $m^{2}\ll s$ and $m^{2}\sim s$. The resummation of virtuality logarithms (with respect to $\mu$) is achieved by the flavor number dependent $\mu$-anomalous dimensions of the hard currents and jet functions known from SCET (with $(n_{\ell}+1)$ flavors) and bHQET (with $n_{\ell}$ flavors) and of the soft function (with $(n_{\ell}+1)$ flavors for scales above $m$ and $n_{\ell}$ flavors for scales below $m$). Furthermore one has to account for the natural scales of the hard, jet and soft function being $\mu_{Q}\sim Q$, $\mu_{J}\sim s^{1/2}$ and $\mu_{S}\sim s/Q$, respectively. In eq. (4) the mass threshold matching functions $H_{m}$ and $\mathcal{M}_{J}$ account for the change between $(n_{\ell}+1)$-flavor SCET evolution (above the mass threshold) and $n_{\ell}$-flavor bHQET evolution (below the mass threshold) for the hard current (see ref. Hoang:2015vua ) and universal jet function, respectively, at the mass threshold scale $\mu_{m}\sim m$. For the universal jet function for a massive primary quark this matching relation reads

[TABLE]

(The analogous relation for the massless primary quark case has been given in eq. (119) of ref. Pietrulewicz:2014qza .) In eq. (4) the mass threshold matching function $\mathcal{M}_{S}$ accounts for the corresponding change between the corresponding $n_{\ell}$- and $(n_{\ell}+1)$-flavor evolution of the soft function (see eqs. (136) and (139) of ref. Pietrulewicz:2014qza ). Rapidity logarithms and their summation take place only inside each of the mass threshold matching functions $H_{m}$, $\mathcal{M}_{J_{Q}}$ and $\mathcal{M}_{S}$ (see for example eq. (44)), and the jet and soft functions are free of any rapidity logarithms. The $\nu$-anomalous dimensions are the same for primary massive and primary massless quarks and were provided in ref. Hoang:2015iva ; Hoang:2015vua . Finally, consistency between eqs. (4) and (4) implies the consistency identity

[TABLE]

The mass mode factorization approach of ref. Pietrulewicz:2017gxc was devised with the motivation to formulate distinct and unique factorization theorems for the different possible scale hierarchies of the mass $m$ with respect to $Q$, $s^{1/2}$ and $s/Q$, where the physically relevant $n$- and $\bar{n}$-collinear and soft mass mode effects appear explicitly within three distinct factors to achieve a transparent and economic representation of their contributions. In contrast to the universal factorization approach, the mass mode factorization approach does not result in two options to formulate the same factorization theorem, such as eqs. (4) and (4), which are equivalent due to consistency relations such as eq. (50). The factorization theorems in the mass mode approach contain functions that have explicit $\mu$ and $\nu$ dependence, and the summation of virtuality and rapidity logarithms is carried out on equal footing. The summation of the associated logarithms, on the other hand (and in contrast to the universal factorization approach), does explicitly involve consistency relations among the anomalous dimensions to allow for a transparent treatment of renormalization group evolution with either $n_{\ell}$ or $(n_{\ell}+1)$ dynamical flavors with respect to the mass threshold scale $m$. In addition, in the formulation of the mass mode factorization approach of ref. Pietrulewicz:2017gxc the resulting factorization theorems have been free of any power corrections in the quark mass $m$ and did not provide a smooth dependence on the quark mass $m$. The latter is mandatory for smooth predictions of differential cross sections when the quark mass threshold is crossed.

Let us consider the situation $s^{2}/Q^{2}\ll m^{2}\ll s\sim p^{2}\ll Q^{2}$. The corresponding factorization theorem in the mass mode factorization approach has the form

[TABLE]

where $H_{m,s}$ is the soft mass mode contribution of the mass threshold hard function $H_{m}$ quoted in eq. (44), $S_{c}$ is the collinear-soft function of eqs. (12) and (41), and $J_{q}^{(n_{\ell}+1)}$ is the jet function for $(n_{\ell}+1)$ massless quarks which coincides with the universal jet function for vanishing quark mass, i.e. $J_{Q}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},0,\mu)$. Since the factorization theorems (4) and (4) as well as (4) provide identical descriptions for the situations $s^{2}/Q^{2}\ll m^{2}\ll s\sim p^{2}\ll Q^{2}$ (up to suppressed quark mass power corrections), they imply the following consistency relations for the mass threshold matching factors $\mathcal{M}_{J_{Q}}$ for the universal function $J_{Q}^{\mathrm{uf},(n_{\ell}+1)}$ and $\mathcal{M}_{S}$ for the soft function $S^{(n_{\ell})}$:

[TABLE]

where $H_{m,n}$ and $H_{m,s}$ are the collinear and soft mass mode contributions of the mass threshold hard function $H_{m}$, see eq. (44), which appears in the universal factorization theorem eq. (4).

In the situation $s^{2}/Q^{2}\ll m^{2}\sim s\sim p^{2}\ll Q^{2}$ the mass mode factorization approach states the factorization theorem

[TABLE]

Using the consistency relations (52) and (53) as well as the relation between the mass mode jet function $J_{Q}^{\mathrm{mf}}$, the universal jet function $J_{Q}^{\mathrm{uf}}$ and the collinear-soft function of eq. (12), it is easy to see that eq. (4) and the universal factorization theorems (4) and (4) provide identical descriptions for $s^{2}/Q^{2}\ll m^{2}\sim s\sim p^{2}\ll Q^{2}$. The relations and the analogous equivalence for eq. (4) show that the factorization theorem (4) – upon replacing the massless quark jet function $J_{q}^{(n_{\ell}+1)}$ by the universal jet function $J_{Q}^{\mathrm{uf},(n_{\ell}+1)}$ – fully accounts for the content of the factorization theorem (4) for $s^{2}/Q^{2}\ll m^{2}\ll s\sim p^{2}\ll Q^{2}$. Thus in this modified form, (4) is applicable for the entire region $s^{2}/Q^{2}\ll m^{2}\lesssim s\sim p^{2}\ll Q^{2}$ and exactly equivalent to the universal factorization theorems (4) and (4).

Thus from a practical application point of view the mass mode factorization theorems (4) and (4) can be considered as special cases of the universal factorization theorems (4) and (4). On the other hand, the mass mode factorization approach for the case $s^{2}/Q^{2}\ll m^{2}\ll s\sim p^{2}\ll Q^{2}$ makes explicit use of the collinear-soft function $S_{c}$ and provides a representation of the quark mass threshold matching factors for the jet and soft functions that disentangle the relevant process-independent collinear and soft mass mode contributions. So, eqs. (52) and (53) in connection with eq. (44) clarify the interplay of the different elementary (collinear and soft) mass mode functions $(H_{m,n},H_{m,\bar{n}},H_{m,s},S_{c})$. Interestingly, by simply replacing the collinear mass mode factors $H_{m,n}$ and $H_{m,\bar{n}}$ valid for primary massive quarks by the corresponding mass mode factors for primary massless quarks (see eqs. (4.10) and (4.12) in ref. Pietrulewicz:2017gxc for analytic expressions) the analogous versions of eqs. (44), (52) and (53) are also valid for the corresponding factorization theorems for primary massless quarks. In fact, upon making these modifications and replacing the primary massive quark jet functions $J_{Q}^{\mathrm{mf},(n_{\ell}+1)}$ and $J_{Q}^{\mathrm{uf},(n_{\ell}+1)}$ by the corresponding jet functions for primary massless quarks, eq. (4), (4) (and (4), (4)) agree identically with the results given in ref. Pietrulewicz:2014qza . This demonstrates the universal structure of the mass mode corrections and illustrates how to obtain them in an economic way in factorization theorems for other applications.

Our final comment concerns the summation of logarithms. Combining the information from the universal factorization and the mass mode factorization approaches provides a modular and transparent method to resum virtuality as well as rapidity logarithms (where we acknowledge, however, that complete treatments of correctly summing both types of logarithms is available in both approaches, see Pietrulewicz:2017gxc and Hoang:2015iva ).

From the point of view of summing virtuality logarithms, the universal factorization theorems (4) and (4) may be seen as more transparent and practical because all anomalous dimensions coincide with those known from factorization theorems where secondary massive quark effects (which entail rapidity singularities) are absent: one only has to keep track of the dependence on the number of dynamical flavors $n_{f}$ being either equal to $n_{\ell}$, when the evolution is below the quark mass $m$ or $(n_{\ell}+1)$ when it is above. If we write the virtuality $\mu$-anomalous dimensions of the SCET and bHQET quark-antiquark production currents ($\mathcal{J}_{\mathrm{SCET}}$ and $\mathcal{J}_{\mathrm{bHQET}}$, respectively Fleming:2007qr ; Fleming:2007xt ), the hard factor, the universal and bHQET jet function, and the soft function in the form ($\hat{s}=s/m$)

[TABLE]

the $\mu$-anomalous dimensions of the corresponding mass threshold matching factors appearing in the factorization theorems (4) and (4) read

[TABLE]

Their form expresses the fact that the mass threshold matching functions $H_{m}$, $\mathcal{M}_{J}$ and $\mathcal{M}_{S}$ simply account for the appropriate change between $(n_{\ell}+1)$ and $n_{\ell}$ flavor $\mu$-evolution. In this context one just has to take into account that for primary massive quarks the $n_{\ell}$-flavor evolution for the jet and the hard functions is carried out in bHQET. For sake of completeness we have displayed all $\mu$-anomalous dimensions up to $\mathcal{O}(\alpha_{s}^{2})$ in appendix C.

From the point of view of summing rapidity logarithms the structure provided by the mass mode factorization theorems (4) and (4) may be seen as most transparent because the universal structure of the $\nu$-anomalous dimensions of the collinear and soft mass mode factors is fully determined by quoting the results

[TABLE]

Interestingly, owing to the consistency relations (44), (52) and (53), the $\nu$-anomalous dimensions of all collinear and soft mass mode factors are identical to all orders in $\alpha_{s}$ up to trivial factors and signs which are determined from the fact that the rapidity divergences within each of the threshold matching functions $H_{m}$, $\mathcal{M}_{J_{Q}}$ and $\mathcal{M}_{S}$ cancel. The consistency relations also entail that the $\nu$-anomalous dimensions of the collinear-soft function $S_{c}$ and the mass mode jet function $J_{Q}^{\mathrm{mf},(n_{\ell}+1)}$ are proportional to a $\delta$-function. The convolution in their anomalous dimension therefore reduces to a simple multiplication which we have already accounted for in the formulae above. The universal $\nu$-anomalous dimension up to $\mathcal{O}(\alpha_{s}^{2})$ reads

[TABLE]

To resum virtuality as well as rapidity logarithms in the universal factorization approach, one evaluates the factorization theorems of eqs. (4) or (4) and uses the $\mu$-anomalous dimensions in eq. (4) to sum all virtuality logarithms by evolving the hard ($H_{Q}^{(n_{\ell}+1)}$), jet ($J_{Q}^{\mathrm{uf},(n_{\ell}+1)}$) and soft ($S^{(n_{\ell})}$) functions from their natural scale to the global renormalization scale. Whenever one of their $\mu$-evolution crosses the mass threshold $m$, so that the anomalous dimension is modified, the corresponding matching factor $H_{m}$, $\mathcal{M}_{J_{Q}}$ or $\mathcal{M}_{S}$ has to be included. The rapidity logarithms are fully contained within each of the mass threshold matching factors $H_{m}$, $\mathcal{M}_{J_{Q}}$ and $\mathcal{M}_{S}$ evaluated at $\mu\sim m$ and can be resummed by using (44), (52) and (53) and the $\nu$-anomalous dimensions of eq. (57).

5 Calculation of Massive Quark Jet Functions at ${\cal O}(\alpha_{s}^{2})$

5.1 Comment on the Calculations and Summary of Matrix Element Results

In this section we provide details on the calculation of the SCET primary massive quark jet functions at ${\cal O}(\alpha_{s}^{2})$. The universal and mass modes SCET jet functions differ concerning the ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ corrections coming from diagrams with secondary massive quarks (which for the jet functions refers to diagrams with a closed massive quark-antiquark loop subdiagram). These corrections involve a number of conceptual and technical subtleties, and we briefly address these in the following before presenting details of the computations in the subsequent sections. We also summarize the individual results for the jet-function ME (calculated in the subsequent sections) and the two collinear-soft MEs (extracted from Pietrulewicz:2017gxc ) which are needed to determine the universal and mass mode SCET jet functions. We believe that this makes the conceptual subtleties of the calculations explicit.

The conceptual starting point is the calculation of the jet-function ME $\!\mathscr{J}_{f}(p^{2},m^{2})$, defined in eq. (1), which is infrared divergent and defined in the $(n_{\ell}+1)$ flavor theory. The infrared finite universal jet function $J_{f}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},m^{2})$ and the mass mode jet function $J_{f}^{\mathrm{mf},(n_{\ell}+1)}(p^{2},m^{2})$ are then given via subtractions related to the likewise infrared divergent collinear-soft MEs $\mathscr{S}^{(n_{\ell}+1)}(\ell,m)$ and $\mathscr{S}^{(n_{\ell})}(\ell,m)$, respectively, defined in the $(n_{\ell}+1)$ flavor theory (i.e. containing $n_{\ell}$ massless quarks and one quark $Q$ with mass $m$, see eq. (7)) and in the $n_{\ell}$ flavor theory (i.e. containing only $n_{\ell}$ massless quarks), respectively. The corresponding definitions are given in eq. (10) and eq. (11), respectively. We stress that the concept of zero-bin subtractions Manohar:2006nz does not play any role in our calculation, and we also believe that our approach can be generalized to other calculations in the context of SCET factorization theorems. As already explained in section 1, the universal jet function $J_{f}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},m^{2})$ by construction approaches the well-known massless quark jet function Becher:2006qw in the limit $m\to 0$. It can therefore be thought of as the mandatory jet function definition for $m^{2}\ll q^{2}$ that involves subtractions related to the mass singularity coming from the secondary quark mass effects in this limit. On the other hand, the mass mode jet function $J_{f}^{\mathrm{mf},(n_{\ell}+1)}(p^{2},m^{2})$ does not involve subtractions related to the secondary quark mass effects and can be thought of as the jet function definition suitable for $m^{2}\sim q^{2}$, where the quark mass $m$ does not represent an infrared scale. It is useful for the computations to keep these conceptual aspects of the two jet functions in mind, even though we emphasize that both jet functions can be applied in a broader kinematic context and are related through the collinear-soft function $S_{c}(\ell,m)$ via eq. (12) (see the discussion in section 4).888In previous literature on the factorization concerning secondary massive quark corrections these subtleties were hidden.

On the technical level an additional subtlety arises because it turns out that - at least at this time - the only feasible way to determine the full set of ${\cal O}(\alpha_{s}^{2})$ corrections is to use two completely separate calculational schemes: one for the purely gluonic and the massless quark corrections at ${\cal O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$ and another one for the secondary massive quark corrections at ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$. This separation arises because the ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ secondary massive quark corrections involve rapidity singularities. To quantify them analytically it is mandatory to introduce an explicit dimensionful infrared regulator (for which we adopt a gluon mass $\Lambda$), such that UV and rapidity divergences can be disentangled and identified unambiguously. The treatment of subdivergences at ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ then also requires a calculation of the ${\cal O}(\alpha_{s})$ corrections with a finite gluon mass for consistency. These calculations have to be carried out for the jet-function ME $\!\mathscr{J}_{f}(p^{2},m^{2})$ and the collinear-soft MEs $\mathscr{S}^{(n_{\ell}+1)}(\ell,m)$ and $\mathscr{S}^{(n_{\ell})}(\ell,m)$ which also entails infrared gluon mass dependent $\overline{\rm MS}$ renormalization $Z$-factors for each of them. For the renormalized universal and mass mode jet functions as well as for their virtuality and rapidity renormalization group equations contributing at ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ it can then be explicitly checked that the dependence on the infrared gluon mass cancels. In contrast, the purely gluonic as well as the massless quark corrections at ${\cal O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$ (as well as the ${\cal O}(\alpha_{s})$ corrections) do not involve rapidity singularities and just using dimensional regularization without an additional infrared regulator is in principle sufficient. In fact using this simpler regularization is the only feasible way to compute the gluonic ${\cal O}(\alpha_{s}^{2})$ corrections due to their complexity. For these corrections then all diagrams contributing to the collinear-soft MEs $\mathscr{S}^{(n_{\ell}+1)}(\ell,m)$ and $\mathscr{S}^{(n_{\ell})}(\ell,m)$ lead to vanishing scaleless integrals which just turn all $1/\varepsilon^{n}$ IR singularities in the infrared divergent jet-function ME $\!\mathscr{J}_{f}(p^{2},m^{2})$ into $1/\varepsilon^{n}$ UV singularities of the infrared finite jet functions $J_{f}^{\mathrm{uf},(n_{\ell}+1)}(p^{2},m^{2})$ and $J_{f}^{\mathrm{mf},(n_{\ell}+1)}(p^{2},m^{2})$. Both jet functions thus agree concerning their ${\cal O}(\alpha_{s})$ and ${\cal O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$ corrections. This justifies that for the ${\cal O}(\alpha_{s}^{2})$ corrections involving only diagrams with massless partons (in addition to the primary massive quark line), at least at the computational level, the subtractions associated to collinear-soft MEs or zero-bins can be ignored. The treatment of subdivergences in this calculation also requires the calculations of the ${\cal O}(\alpha_{s})$ corrections, which differ from the corresponding results with a gluon mass regulator. The lack of having an infrared regularization scheme that allows the computation of all ${\cal O}(\alpha_{s}^{2})$ corrections in a uniform and manageable way entails that - at least at the present time - the computation of the ${\cal O}(\alpha_{s}^{3})$ corrections to the SCET jet functions in the presence of massive quarks would represent an extremely challenging task.

We present the calculations of the purely gluonic and massless quark corrections at ${\cal O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$ for the jet-function ME $\!\mathscr{J}_{f}(p^{2},m^{2})$ in section 5.2, and those of the secondary massive quark corrections at ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ in section 5.3. The reader not interested in the computational technical details may skip these two subsections. In the following we summarize the analytic results for jet-function ME and the two collinear-soft MEs within the two calculational schemes. Results needed for the calculation of the ${\cal O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$ corrections to the jet function are computed in the massless gluon computational scheme (where the ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ corrections cannot be given) and are quoted with the subscript “$\cancel{C_{F}T_{F}}$”. In contrast all results needed for the calculation of the ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ terms from secondary mass effects are computed in the massive gluon computational scheme (where the ${\cal O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$ corrections cannot be given) and are quoted with the subscript “$C_{F}T_{F}$”. Note that we use this labeling also for the corresponding one-loop results that contribute (through the treatment of subdivergences) to the determination of the respective renormalized two-loop jet functions and their anomalous dimensions.999We note that all statements made before in section 5.1 apply for massless ($f=q$) and massive primary quarks ($f=Q$).

In the computational scheme with a massless gluon the renormalized jet-function ME can be written in the form

[TABLE]

where the ${\cal O}(\alpha_{s})$ coefficient $J_{Q}^{(1)}(p^{2},m^{2},\mu)$ agrees with the well-known one-loop SCET primary massive quark jet function result from ref. Fleming:2007xt and is given up to ${\cal O}(\varepsilon^{2})$ in eq. (19). The two-loop coefficient $J_{Q}^{(2)}$ containing the ${\cal O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$ corrections is given in eq. (23). Both are free of rapidity divergences and therefore do not carry an argument with respect to the rapidity renormalization scale $\nu$. For convenience, we also provide the corresponding renormalized result for the massless primary quark jet-function ME $\left.\!\mathscr{J}_{q}^{(n_{\ell}+1)}(p^{2},m^{2},\Lambda^{2}=0,\mu,Q/\nu)\right|_{\cancel{C_{F}T_{F}}}$ in appendix A.

The divergent $Z$-factor of the jet-function ME reads

[TABLE]

By definition it contains all $1/\varepsilon$-divergent terms and quantifies the difference between the renormalized and bare jet-function MEs,

[TABLE]

In the context of our jet-function ME computation it includes UV and IR divergent $1/\varepsilon$ terms. We stress that, upon (literally) replacing the color factor $C_{F}T_{F}n_{\ell}$ by $C_{F}T_{F}(n_{\ell}+1)$ in eq. (60), one obtains the full ${\cal O}(\alpha_{s}^{2})$ $\overline{\rm MS}$ $Z$-factor for the universal SCET jet functions $J_{f}^{\mathrm{uf},(n_{\ell}+1)}$ (for massless or massive primary quarks) as well as for the SCET jet function for $(n_{\ell}+1)$ massless quark flavors, where all $1/\varepsilon$ divergent terms are UV divergences.101010Note that that the $Z$-factor for the mass mode jet function $J_{f}^{\mathrm{mf},(n_{\ell}+1)}$ can be obtained from the results for the bare and renormalized jet-function ME $\!\mathscr{J}_{Q}^{(n_{\ell}+1)}$ and collinear-soft ME $\mathscr{S}^{(n_{\ell})}$. Because both are defined in theories with different number of quark flavors, the $Z$-factor for the mass mode jet function cannot be quoted in the $\overline{\rm MS}$ scheme and in general has finite contributions. This relation is required by consistency and can be used as a cross check for calculations.

In the massless gluon computational scheme the loop corrections to the bare collinear-soft ME $\mathscr{S}^{(n_{\ell})}$ are given by scaleless integrals and thus vanish to all orders in perturbation theory

[TABLE]

The corresponding bare collinear-soft ME $\mathscr{S}^{(n_{\ell}+1)}$, on the other hand, corresponds to scaleless integrals only at ${\cal O}(\alpha_{s})$ and ${\cal O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$:

[TABLE]

In the computational scheme with an infinitesimal gluon mass $\Lambda$ the terms of the renormalized jet-function ME relevant for the calculation of the ${\cal O}(\alpha_{s}^{2}C_{F}T_{F})$ secondary quark mass effects read

[TABLE]

where $G_{\mathrm{sec}}^{(3Q)}$ is given in eq. (31). For convenience, we also provide the analogous results for the renormalized massless primary quark jet-function ME $\left.\!\mathscr{J}_{q}^{(n_{\ell}+1)}(p^{2},m^{2},\Lambda^{2},\mu,Q/\nu)\right|_{C_{F}T_{F}}$ in appendix A. The ${\overline{\rm MS}}$ renormalization $Z$-factor associated with the jet-function ME contributions in eq. (64) (which is identical for massive and massless primary quarks) accounts for UV and rapidity divergences and reads

[TABLE]

The corresponding results for the two different renormalized collinear-soft MEs are

[TABLE]

and

[TABLE]

Their respective ${\overline{\rm MS}}$ renormalization $Z$-factors accounting for UV and rapidity divergences read

[TABLE]

and

[TABLE]

These results can be extracted from ref. Pietrulewicz:2017gxc , and we have explicitly cross checked their results.111111Note that eq. (B.55) of ref. Pietrulewicz:2017gxc has a typo due to a missing global factor of $4$.

From the results quoted above it is straightforward to obtain the results for the universal (see eqs. (10) and (39)) and the mass mode SCET jet functions (see eqs. (11) and (18)). The virtuality and rapidity anomalous dimensions for the universal and mass mode SCET jet functions (and the collinear-soft function) (see eqs. (4) and (57) and appendix C) are related to the corresponding anomalous dimensions for the jet-function and the collinear-soft MEs (see appendix C as well).

5.2 Gluonic and Massless Quark Corrections

In figure 1 all two-loop Feynman diagrams relevant for the calculation of the two-loop jet-function ME in eq. (4) using QCD Feynman rules and Feynman gauge are displayed. In this section we discuss the computation of the ${\cal O}(\alpha_{s}^{2}C_{F}^{2},\alpha_{s}^{2}C_{F}C_{A},\alpha_{s}^{2}C_{F}T_{F}n_{\ell})$ corrections using the computational scheme with a massless gluon, as explained in the previous subsection.

After evaluating the Dirac trace in eq. (4) the contribution of each diagram for the jet function forward scattering ME can be expressed as a linear combination of dimensionally regularized ($d=4-2\varepsilon$) scalar integrals of the generic form

[TABLE]

where all propagator denominators have the usual $-i0$ term and the prescription $s\to s+i0$ for the jet virtuality $s=p^{2}-m^{2}\geq 0$ is understood. In the second line we pulled out the conventional factor of $i\pi^{d/2}$ per loop integral. We also pulled out the factor $(\bar{n}\cdot p)^{-c_{1}-c_{2}}$ as this dependence is fixed by the behavior of the integral under the rescaling of $\bar{n}$ or Lorentz covariance. The power of $(-s)$ is then determined by choosing $\mathcal{I}$ to be dimensionless. This factor encodes the full $s$ dependence in the massless limit ($s\to p^{2}$), where $\mathcal{I}$ becomes a function of $\varepsilon$ only.

The set of propagator denominators in eq. (70) does not represent a linearly independent basis. Either $a_{3}$ or $b_{3}$ can be rendered zero by partial fraction decomposition. We therefore can arrange all integrals into two integral families: family A, where $b_{3}=0$, and family B, where $a_{3}=0$. Integrals where $a_{3}=b_{3}=0$ are for convenience assigned to family A. Some family B integrals with only two massive propagators can be mapped by loop momentum shifts to family A. Diagrams (g) and (c) contain the top sector integrals of family A and B, respectively (i.e. where all propagator powers are positive). We distinguish two types of diagrams: the “planar” diagrams (e)-(q), which contain the leading color contributions, and the “nonplanar” diagrams (a)-(d) which are proportional to $C_{F}^{2}-C_{F}C_{A}/2$ and therefore color-suppressed. The planar diagrams only involve family A integrals, while the “nonplanar” diagrams contain integrals of family A and B. Diagrams (r) are scaleless and vanish.

The classification in family A and B integrals is useful for two reasons. The first is that they are closed under integration by parts (IBP) reductions and require different sets of IBP identities. The second is related to the different mathematical properties of family A and B integrals. Family A integrals have at most two massive propagators and involve at most harmonic polylogarithms in the $\varepsilon$ expanded results. In contrast, family B integrals can contain elliptic functions in the $\varepsilon$ expanded expressions, which makes their evaluation more complicated. This is connected to the fact that they admit a triple massive-particle cut. Moreover, we observed that some family B integrals (including cases with only two massive propagators) are ill-defined in pure dimensional regularization due to rapidity singularities even though the Feynman diagrams are individually rapidity finite. Further, using IBP reductions on rapidity finite family B integrals can lead to rapidity singular integrals. Therefore, the computations of family B integrals in general require a rapidity regulator.

Let us briefly illustrate this issue for the following rapidity-finite family B integral

[TABLE]

where we introduced an analytic rapidity regulator $\eta$ for the $k$-integration, see e.g. refs. Becher:2011dz ; Chiu:2011qc ; Chiu:2012ir . Applying the IBP relation

[TABLE]

where the bold letters represent operators that increase ($+$) or decrease ($-$) the corresponding propagator powers in $\mathcal{I}$, we obtain

[TABLE]

The first three integrals on the right hand side turn out to be individually rapidity divergent, i.e. they have $1/\eta$ poles as can e.g. be checked numerically with the sector decomposition program pySecDec Borowka:2017esm . These $1/\eta$ poles cancel in the sum. Without a rapidity regulator, i.e. naively setting $\eta=0$ in eq. (73) literally, however, gives an incorrect relation because the term $\eta\,\mathcal{I}(-y;0,1,0,1,1,1,2+\eta,0)$ leads to finite contributions which are missed when the rapidity regularization is introduced after using the IBP relation. An analogous issue also arises when using a dimensionful rapidity regulator such as the $\Delta$ regulator proposed in ref. Chiu:2009yx . We stress that the IBP reduction with a rapidity regulator (no matter of what kind) in general leads to a substantially increased number of master integrals. Moreover, these are typically more complicated to compute than comparable integrals where this kind of problem does not arise and an additional rapidity regularization is not necessary.

Interestingly, the IBP reduction of family A integrals does not give rise to spurious rapidity divergences and works consistently without any rapidity regulator in the standard way. This is (currently) an empirical observation as we are lacking a simple systematic criterion that would allow us to identify when the problem of spurious rapidity divergences arises and to possibly avoid or treat them in an automatized way. We therefore explicitly verified all family A integral reductions used in our computations numerically in order to dispel any doubt about their validity. In this context the known (or easily computed) massless primary quark limit of each individual diagram provides a valuable cross check of our calculations.

The calculation of the planar diagrams (e)-(q) can thus be performed with standard modern multi-loop technology as explained in section 5.2.1. For the nonplanar diagrams (a)-(d), the issue of the spurious rapidity divergences and the elliptic nature of certain family B integrals forced us to use different and less uniform methods as explained in section 5.2.2. As a general strategy, we analytically compute the two-loop master integrals for the jet field correlator to the required order in $\varepsilon$ first, insert them into the different diagrams, and finally take the imaginary part of the total contribution to the jet field correlator according to eq. (4) for each color factor. The only exceptions are diagrams (b) and (c), for which we take a different semi-numerical approach, as explained in section 5.2.2. We renormalize the resulting jet-function MEs \!\mathscr{J}_{Q}^{(n_{\ell}+1)}\big{|}_{\cancel{C_{F}T_{F}}} according to eq. (61), i.e. we absorb all divergent terms into the $Z$-factor, and we use the pole mass scheme for the quark mass. The necessary counterterms for the pole mass $m$ (and the ${\overline{\mathrm{MS}}}$ coupling $\alpha^{(n_{\ell}+1)}_{s}$) can e.g. be found in ref. Melnikov:2000zc . We note that in the pole mass scheme derivatives of the distributions $\delta(s)$ and $\mathcal{L}_{n}(s)$ do not arise in the final result.

5.2.1 Planar Diagrams

The calculation of the planar diagrams (e)-(q) (including only massless partonic loops in the gluon self-energy) is rather straightforward. They involve ${\mathcal{O}}(100)$ family A integrals. We reduced this set of integrals to the seven master integrals depicted in figure 2 using the IBP program FIRE5 Smirnov:2014hma .

The master integrals $M_{1}$-$M_{5}$ can be computed using Feynman parameters for arbitrary $\varepsilon$ in terms of hypergeometric functions. They can be expanded in terms of harmonic polylogarithms (HPLs) Gehrmann:2000zt ; Remiddi:1999ew with the help of the Mathematica package HypExp Huber:2005yg . We also used the Mathematica package HPL Maitre:2005uu to exploit relations among the HPLs for the simplification of expressions or to make their singular behavior for $s\to 0$ explicit, see below.

For $M_{6}$, $M_{7}$ we used the method of differential equations Kotikov:1990kg ; Remiddi:1997ny ; Gehrmann:1999as . To this end we take their derivatives with respect to $y=m^{2}/p^{2}$. The result can be expressed via IBP reduction in terms of $M_{1}$-$M_{7}$. The coupled system of seven homogeneous first-order differential equations can also be rewritten in terms of a coupled system of two first-order differential equations for $M_{6}$, $M_{7}$ with linear combinations of $M_{1}$-$M_{5}$ as (known) inhomogeneous terms:

[TABLE]

These can be decoupled into two separate second-order inhomogeneous differential equations for $M_{6}$ and $M_{7}$. Upon expansion in $\varepsilon$ these can be iteratively rewritten in second-order inhomogeneous differential equations where only derivatives of expansion coefficients of $M_{6}$ or $M_{7}$ arise. These differential equations happen to have a particularly simple form because these derivatives appear in a combination, such that upon integration the resulting differential equations only involve first derivatives of the expansion coefficients. Another integration w.r.t. $y$ then yields the result depending on two integration constants (boundary values) that can be fixed by taking the (massless) limit $y\to 0$, where the integrals are known Becher:2006qw . We give the explicit solutions for $M_{1}$-$M_{7}$ in appendix D. We checked all of them numerically to the required order in $\varepsilon$ using the Sector Decomposition codes FIESTA4 Smirnov:2015mct and pySecDec Borowka:2017esm .

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

[TABLE]

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

[TABLE]

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

5.2.2 Nonplanar Diagrams

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

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

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

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

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

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

[TABLE]

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

[TABLE]

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

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

Furthermore, one can argue on physical grounds that

[TABLE]

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

[TABLE]

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

[TABLE]

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

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

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

5.3 Secondary Massive Quark Corrections

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

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

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

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

[TABLE]

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

[TABLE]

and the remainder term reads

[TABLE]

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

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

[TABLE]

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

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

[TABLE]

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

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

[TABLE]

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

6 Numerical Analysis

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

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

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

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

7 Conclusions

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

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

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

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

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

Acknowledgments

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

Appendix A Other Jet Functions

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

[TABLE]

with

[TABLE]

[TABLE]

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

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

[TABLE]

with

[TABLE]

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

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

[TABLE]

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

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

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

[TABLE]

[TABLE]

Appendix C Anomalous Dimensions

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

[TABLE]

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

[TABLE]

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

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

[TABLE]

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

Appendix D Master Integrals

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

[TABLE]

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

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

[TABLE]

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

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