Heavy Quark Fragmenting Jet Functions

Heavy quark fragmenting jet functions describe the fragmentation of a parton into a jet containing a heavy quark, carrying a fraction of the jet momentum. They are two-scale objects, sensitive to the heavy quark mass, $m_Q$, and to a jet resolution variable, $\tau_N$. We discuss how cross sections for heavy flavor production at high transverse momentum can be expressed in terms of heavy quark fragmenting jet functions, and how the properties of these functions can be used to achieve a simultaneous resummation of logarithms of the jet resolution variable, and logarithms of the quark mass. We calculate the heavy quark fragmenting jet function $\mathcal G_Q^Q$ at $\mathcal O(\alpha_s)$, and the gluon and light quark fragmenting jet functions into a heavy quark, $\mathcal G_g^Q$ and $\mathcal G_l^Q$, at $\mathcal O(\alpha_s^2)$. We verify that, in the limit in which the jet invariant mass is much larger than $m_Q$, the logarithmic dependence of the fragmenting jet functions on the quark mass is reproduced by the heavy quark fragmentation functions. The fragmenting jet functions can thus be written as convolutions of the fragmentation functions with the matching coefficients $\mathcal J_{i j}$, which depend only on dynamics at the jet scale. We reproduce the known matching coefficients $\mathcal J_{i j}$ at $\mathcal O(\alpha_s)$, and we obtain the expressions of the coefficients $\mathcal J_{g Q}$ and $\mathcal J_{l Q}$ at $\mathcal O(\alpha_s^2)$. Our calculation provides all the perturbative ingredients for the simultaneous resummation of logarithms of $m_Q$ and $\tau_N$.


Introduction
The production of heavy flavors and heavy flavored jets, where by heavy flavor we here mean charm or bottom, plays an important role in collider experiments. These processes are interesting in themselves, as a probe of QCD dynamics, since for heavy quarks one expects the closest correspondence between calculations at the parton level and experimentally measured hadrons. Furthermore, b jets are found in interesting electroweak processes; for instance, one of the most important channel to probe fermionic couplings of the Higgs boson is H → bb. Finally, a quick look at the ATLAS and CMS public results pages shows the almostomnipresence of b jets in searches of Beyond Standard Model physics. It is therefore important to have a good theoretical understanding of heavy flavor production (where a heavy flavored hadron is directly observed) and heavy flavored jets (where a jet is tagged by demanding that it contains at least one heavy flavored hadron) in collider experiments.
The current state of the art of fixed order calculations for heavy flavor hadroproduction is next-to-leading order (NLO) accuracy, and such NLO calculations have a long history [1][2][3][4]. By now, several processes, including heavy quark pair production, associated production of weak bosons and heavy quarks, and Higgs production with decay into bb, are implemented in the program MCFM [5] at NLO accuracy, and any distribution can be obtained for these processes.
In such fixed order calculations, the dependence on the heavy quark mass typically enters through the ratio were p T is the transverse momentum of the heavy quark. Many of the available NLO calculations include the full dependence of the heavy quark mass, which is important at small to moderate values of p T . For large values of p T , the ratio r Q becomes negligible, and one might want to perform calculations with massless heavy quarks, for which NLO calculations are significantly simpler. However, besides a dependence on powers of r Q , there is also a logarithmic dependence on that ratio, which arises from infrared divergences in the massless calculation which are regulated by the heavy quark mass. Thus, at higher orders in perturbation theory more powers of these logarithms appear, requiring resummation. This is accomplished by introducing a heavy quark fragmentation function [6]. As was discussed for hadroproduction in Ref. [7], the heavy quark fragmentation function can be calculated perturbatively at scales µ ∼ m Q without encountering any large logarithms. Running the fragmentation function from this low scale to µ ∼ p T using the familiar DGLAP evolution then resums all logarithms of m Q /p T . A combination of both approaches is needed to describe heavy flavor production for both large and small values of p T . Such a combination, named "Fixed Order plus Nextto-Leading-Log" (FONLL), has been proposed in Ref. [8], and applied to single inclusive production of heavy flavored hadrons. The general idea of FONLL is to add the massive fixed order calculation to the resummed calculation, and then subtract the overlap of the two. The overlap can be calculated either as the massless limit of the fixed order calculation (keeping the logarithms), or as the expansion of the resummed calculation to the appropriate order. The FONLL approach has been successfully compared to Tevatron and LHC data, for a recent discussion see Ref. [9].
Over the past decade we have learned how to combine NLO calculations with parton shower algorithms. This provides final states which are fully showered and hadronized, but which still provide NLO accuracy for predicted observables. Since such calculations can be compared much more directly to experimental data, this is used a great deal in analyses. The most popular available methods are MC@NLO [10] and POWHEG [11,12], with several other approaches being pursued. Both MC@NLO and POWHEG include heavy flavor production in their list of available processes [13,14]. Since parton showers resum leading logarithms in their evolution variable t, any calculation that is interfaced with such a shower needs to provide at least the same amount of resummation. In fact, as was discussed in detail in Ref. [15], any combination of a perturbative calculation with a parton shower algorithm requires at least LL resummation of the dependence on an infrared safe jet resolution variable, however this jet resolution variable does not necessarily have to be equal to the evolution variable of the parton shower. Thus, one can choose a resolution variable for which one has good theoretical control, such as N -jettiness. We will denote a general dimensionful N -jet resolution variable by τ N and define the dimensionless ratio with p T denoting the transverse momentum of the hadron, or of the jet in which the hadron is found, which we assume to be not too different. It follows from the above discussion that combining heavy quark production at large p T with parton shower algorithms requires the simultaneous resummation of logarithms of r Q and r τ . 1 Logarithmic dependence on a second ratio r τ can also arise from explicit experimental cuts restricting the size of the jet resolution variable τ N . For example, observables that explicitly restrict extra jet activity through jet vetoes will have logarithmic dependence on the jet veto scale τ cut N . Extending the discussion to heavy-flavor tagged jets, one might expect jet observables to be less sensitive to the heavy quark fragmentation function, and to logarithms of m Q [16]. This is because heavy-flavor tagged jets are essentially agnostic to the flavor of the heavy hadron and its energy fraction. Indeed heavy quark jets initiated by heavy quarks produced directly in the hard interaction do not have a logarithmic dependence on m Q . In this case it is not necessary to introduce a fragmentation function, and to resum log r Q . The resummation of log r τ can be achieved with methods similar to those used for light quark jets.
However, heavy-flavor tagging algorithms also tag jets initiated by gluons or light quarks, where heavy quarks are produced through g → QQ. In this case, Q-tagging introduces an infrared dependence on logarithms of m Q , and large uncertainties [17]. A possible way to deal with final state logarithms is not to label jets with gluon or light quark splittings into QQ as heavy-flavor jets. Banfi, Salam and Zanderighi in Ref. [17] explored the interesting possibility of using an IR safe jet flavor algorithm [18], which would label jets with no net heavy flavor as gluon or light quark jets. Alternatively, one can improve the theoretical description of Q-tagged jet cross sections by resumming log r Q , and, in the presence of another small ratio r τ , by simultaneously resumming log r τ . 1 It should be noted that the leading log resummation in the shower resums a subset of logarithms of the heavy quark mass, for example all terms originating from emissions from the heavy quark or antiquark in the final state. However, not all log rQ are included at leading logarithmic accuracy. In particular, gluon splitting in QQ, with almost collinear quark and antiquark, are included only at fixed order. For a full discussion, see [13,14].
In this paper we develop a formalism that allows to simultaneously resum the logarithmic dependence on the heavy quark mass as well as on the additional small ratio r τ . This opens the door to deal with vetoed heavy flavor production in a systematic way, and perhaps more importantly to interface FONLL-type calculations with a parton shower. Our formalism is based on the idea of fragmenting jet functions (FJFs), introduced in Refs. [19,20], and first applied to heavy quarks in Ref. [21]. The FJFs G j i describe the fragmentation of a parton j inside a jet initiated by the parton i, and contain information both on the jet dynamics, and on the parton fragmentation function. Therefore, FJFs encode the dependence on both m Q , as well as τ N . An important property of the heavy flavor FJFs is that the renormalization group evolution is independent of the heavy quark mass, with the anomalous dimension being identical to that of an inclusive jet function. Thus, the dependence on the jet resolution variable τ N can be resummed in the same way as for processes with only light jets.
To perform a simultaneous resummation of r Q and r τ requires to separate these scales in the factorization theorem, and therefore factorize the FJFs themselves. For τ N m 2 Q /Q this is accomplished by integrating out the degrees of freedom responsible for the τ N scale, with the remaining long-distance physics (and therefore the entire m Q dependence) determined by the heavy-quark fragmentation function.
The main part of this paper is an explicit calculation of the heavy flavor FJFs in fixed order perturbation theory. We calculate the heavy quark initiated FJF at O(α s ), and the gluon and light-quark initiated FJFs at O(α 2 s ). Besides being an important ingredient to obtain the resummed expressions, it also allows us to check explicitly various properties of the heavy flavor FJFs. In particular, we verify that the heavy quark fragmentation functions reproduce the logarithmic dependence on the heavy quark mass and that the anomalous dimension of the heavy-quark FJFs are independent of m Q . Our calculations are performed using Soft Collinear Effective Theory (SCET) [22][23][24][25][26].
The paper is organized as follows. In Section 2 we recall the main SCET ingredients needed in the rest of the paper. We define and state important properties of heavy quark fragmentation functions in Section 2.1, and of inclusive jet functions in Section 2.2. In Section 2.3 we introduce the fragmenting jet functions, extending the definition of Refs. [19,20] to heavy quarks. After reviewing the resummation of log r Q in single inclusive observables, and of log r τ in jet observables in Section 3, we describe how to achieve the simultaneous resummation of logarithms of the quark mass and the jet resolution variable τ N in Section 4. In Section 5 we calculate the FJFs G Q Q and G Q g at O(α s ) in the massless limit, we give the expressions with full mass dependence in Appendix B. In Sections 6.1 and 6.2 we carry out the calculation of G Q g and G Q l at O(α 2 s ). We draw our conclusions in Section 7. In Appendix A we discuss some additional details on how to take the massless limit m Q → 0. In Appendix C we give the analytic expression of the function g C F T R , defined in Section 6.1.

Soft Collinear Effective Theory
In this paper we use the formalism of Soft Collinear Effective Theory (SCET) [22][23][24][25], generalized to massive quarks [26]. SCET is an effective theory for fast moving, almost light-like, quarks and gluons, and their interactions with soft degrees of freedom. It has been successfully applied to a variety of processes, from B physics to quarkonium, and there is a growing body of application to the study of jet physics and collider observables.
We are interested in processes sensitive to three scales, Q 2 , Qτ N and m 2 Q . Q is the hard scattering scale, represented by the p T of the hardest jet in the event. τ N defines the jet scale, so the typical size of Qτ N is the jet invariant mass, while m Q is the heavy quark mass. We are interested in the situation Q 2 Qτ N , m 2 Q . In this case, degrees of freedom with virtuality of order Q 2 can be integrated out by matching QCD onto SCET. The degrees of freedom of SCET are collinear quarks and gluons, with virtuality p 2 ∼ Q 2 λ 2 , and ultrasoft (usoft) quarks and gluons, with even smaller virtuality p 2 ∼ Q 2 λ 4 . λ is the SCET expansion parameter, λ ∼ m/Q 1, with m the next relevant scale in the problem, e.g. m 2 = Qτ N . In SCET different collinear sectors can only interact by exchanging usoft degrees of freedom. An important property of SCET is that usoft-collinear interactions can be moved from the SCET Lagrangian to matrix elements of external operators [24], greatly simplifying the proof of factorization theorems. Since the dynamics of different collinear sectors and of usoft degrees of freedom factorize, we can focus in this paper on jets in a single collinear sector.
If there is a large hierarchy between the remaining two scales, Qτ N m 2 Q , we can further lower the virtuality of the degrees of freedom in the effective theory by integrating out particles with virtuality Qτ N at the jet scale. This second version of SCET has collinear fields with p 2 ∼ m 2 Q . The additional matching step allows to factorize the dynamics of the two scales m Q and τ N , and to resum large logarithms of their ratio m 2 Q /(Qτ N ). We now summarize some SCET ingredients needed in the rest of the paper. For more details, we refer to the original papers [22][23][24][25][26]. We introduce two lightcone vectors n µ and n µ , satisfying n 2 =n 2 = 0, andn · n = 2. The momentum of a particle can be decomposed in lightcone coordinates according to Particles collinear to the jet axis have (p + , p − , p ⊥ ) ∼ Q(λ 2 , 1, λ), where λ 1 is the SCET expansion parameter. Usoft quarks and gluons have all components of the momentum roughly of the same size (p + , p − , p ⊥ ) ∼ Q(λ 2 , λ 2 , λ 2 ).
The SCET Lagrangian can be written as Each collinear sector is described by a copy of the collinear Lagrangian L n . For massless quarks, L n is L n =ξ n in · D n + gn · A us + ( / ξ n and A n are collinear quark and gluon fields, labeled by the lightcone direction n and by the large components of their momentum (p − , p ⊥ ). We leave the momentum label mostly implicit, unless explicitly needed. The label momentum operator P µ acting on collinear fields returns the value of the label, for example The collinear covariant derivative D n is defined as W n are Wilson lines, constructed with collinear gluon fields, A µ us is an usoft gluon field. At leading order in λ, usoft gluons couple to collinear quarks only through n · A us . This coupling can be eliminated from the Lagrangian via the BPS field redefinition [24]: with P denoting path ordering. The effect of the field redefinition is to eliminate the usoft gluon in Eq. (2.3), and to replace the collinear quark and gluon fields ξ n and A n with their non-interacting counterparts. The same field redefinition also decouples usoft from collinear gluons [24]. From here on we always use decoupled collinear fields, and drop the superscript (0). For fast moving massive particles there are additional mass terms [26], We work with one massive quark with mass m Q , n f − 1 massless quarks and assume that quarks heavier than m Q have been integrated out. In this paper we use q to denote both heavy and light quarks, when it is not necessary to specify the quark mass. Q (Q) is used exclusively for heavy quarks (antiquarks), while l (l) denotes the n l = n f − 1 light quarks (antiquarks). Using the Wilson line W n it is possible to construct gauge invariant combinations of collinear fields Collinear gauge invariant operators are expressed in terms of matrix elements of these building blocks [25]. In the next subsections, we discuss three such operators, heavy quark fragmentation functions, inclusive quark and gluon jet functions, and heavy quark fragmenting jet functions.

Heavy Quark Fragmentation Functions
Fragmentation functions describe the fragmentation of a parton into a hadron H, which carries a fraction z of the parton momentum. In SCET, the operator definitions of the fragmentation function of a quark or a gluon into a hadron H are given by [19] The trace in Eq. (2.12) is over Dirac and color indices. The sum over X denotes the integration over the phase space of all possible collinear final states. In Eqs. (2.12) and (2.13), ω denotes the large component of the momentum of the fragmenting parton, and the frame in which the fragmenting parton has zero p ⊥ has been chosen. The hadron H has momentum p h . The perpendicular component p ⊥,h is integrated over, while p − h , or, equivalently, the momentum fraction z = p − h /ω is measured. N c is the number of colors, N c = 3. The definitions in Eqs. (2.12) and (2.13) are equivalent to the classical definition of fragmentation function in QCD, in Ref. [27].
The evolution of the fragmentation functions is governed by the DGLAP equation [28] d d log where the P ji (ξ) are the time-like splitting functions. The splitting functions are computed in perturbation theory with, at one loop, [28] The color factors in Eqs. (2.16)-(2.19) are C F = 4/3, C A = 3, T R = 1/2, while β 0 is the leading order coefficient of the beta function, The space-like and time-like splitting functions at O(α 2 s ) are given in Refs. [29,30], and nicely summarized in Ref. [31]. Space-like splitting functions are known to O(α 3 s ) [32,33]. The nonsinglet component of the time-like splitting functions is also known to three-loops [34], while the singlet is, at the moment, unknown.
The fragmentation functions of light hadrons are non-perturbative matrix elements, which need to be extracted from data. In the case of heavy flavored hadrons, the heavy quark mass m Q is large compared to the hadronization scale Λ QCD . Neglecting corrections of order Λ QCD /m Q , one can identify the heavy hadron with a heavy quark or antiquark and the fragmentation function can be computed in perturbation theory at the scale m Q [6]. Expanding in α s , the fragmentation function for a heavy quark into a quark or a gluon, and for a gluon into a heavy quark at O(α s ) are The fragmentation functions of a heavy quark, heavy antiquark or light quark into a heavy quark were computed at O(α 2 s ) in Ref. [35], while the fragmentation of a gluon into a heavy quark in Ref. [36].
The fixed order expressions for the heavy quark fragmentation functions are reliable at scales µ ∼ m Q , where logarithms are small. The fragmentation functions at an arbitrary scale µ are obtained by taking the fixed order expressions as initial condition for the DGLAP evolution. The evolution of the one-loop initial condition (2.22)-(2.25) with O(α 2 s ) splitting functions resums all leading and next-to-leading logarithms (NLL), that is all terms of the form α n s log n (m 2 Q /µ 2 ) and α n s log n−1 (m 2 Q /µ 2 ). The knowledge of the initial condition at O(α 2 s ), and of the non-singlet splitting functions at O(α 3 s ), allows to achieve NNLL accuracy for non-singlet combinations of the quark fragmentation function, for example D Q Q − DQ Q . The evolution of the gluon distribution D Q g and of the singlet distribution require the time-like singlet splitting function to O(α 3 s ), which, at the moment, is not known. The picture obtained with the partonic initial conditions (2.22)-(2.25) and the DGLAP evolution is a valid description of the fragmentation functions of heavy hadrons, except in the endpoint region, 1 − z ∼ Λ QCD /m Q , corresponding to the peak of the quark distribution. In this region soft gluon resummation and non-perturbative effects become important and a model describing hadronization must be included, and fitted to data [37][38][39].
We conclude this section by mentioning two important sum rules obeyed by the fragmentation functions. The first is the momentum conservation sum rule, The sum is extended over a complete set of states. Eq. (2.26) is the statement that the total energy carried off by all the fragmentation products sums to that of the original parton. At the perturbative level, H ∈ {Q,Q, g, l,l}. Eq. (2.26) can be readily verified using the one loop results for the quark distributions in Eqs. (2.22)-(2.24). Using the one loop expression for the splitting functions, and the DGLAP equation (2.14), one can also check that the momentum conservation sum rule is not spoiled by renormalization. This is true at all orders [40]. In addition, there are flavor conservation sum rules. For heavy quarks, dz(D Q Q (z, µ 2 ) − DQ Q (z, µ 2 )) = 1 .

Inclusive Jet Functions
The gauge invariant quark and gluon fields, Eq. (2.11), are a natural ingredient for the description of jets in SCET. Quark and gluon inclusive jet functions are defined as matrix elements of χ n and B n⊥ [41][42][43], We work in a frame where the momentum is aligned with the jet direction, p µ J = (ω, r + , 0). ω is the large component of the momentum, of the size of the jet p T , and the jet invariant mass is ωr + p 2 T . To simplify the notation, in the rest of the paper we drop the superscript on the plus component of the jet momentum. The quark and gluon inclusive jet functions are infrared finite quantities, insensitive to the scale Λ QCD , and can be computed in perturbation theory. In the case the quark field χ n is massive, the quark jet function (2.28) depends on the quark mass, but the dependence is not singular [44,45].
Beyond leading order, the quark and gluon jet functions are UV divergent and require renormalization. The dependence on the renormalization scale µ is governed by the renor- The anomalous dimension is where the plus distribution of the dimensionful variable ωr is defined as and it is independent of the arbitrary cut-off ωκ. Γ q cusp and Γ g cusp are the quark and gluon cusp anomalous dimensions [46,47], which are known to three loops [32]. Up to this order, they are related by Γ g cusp /Γ q cusp = C A /C F . γ i is the non-cusp component of the anomalous dimension, known to O(α 2 s ) [48,49]. The form of the anomalous dimension (2.31), in particular its dependence on log µ 2 , allows to resum Sudakov double logarithms. The RGE (2.30) can be solved analytically, and, given an initial condition at the scale µ I , the jet function at the scale µ F is where U J is an evolution function, given, for example in Ref. [45].
In hadronic collisions, in addition to collinear radiation from final state particles, one has to account for initial state radiation from the incoming beams. Initial state radiation is described by beam functions B i [50,51]. The beam functions depend on the invariant mass t and also on the momentum fraction x of the incoming parton. They satisfy the same RGE as final state jets, Eq. (2.30). Large Sudakov logarithms induced by collinear radiation from the incoming beams are resummed in the same way as logarithms in the jet functions, Notice that the beam function evolution does not change the distribution in the momentum fraction x. The beam functions are perturbatively related to the parton distributions [50,51], according to In this case, the initial condition for the evolution (2.34) cannot be computed purely in perturbation theory, but it is obtained convoluting the perturbative matching coefficients I ij with the parton distributions evaluated at the scale µ 2 I . The last ingredient in factorization theorems for jet cross sections is a soft function, describing soft interactions between jets, and between jets and the beams. The precise definition of the soft function depends on the observable in consideration, but in general its RGE is of similar form as Eq. (2.30), and resums Sudakov double logarithms.

Heavy Quark Fragmenting Jet Functions
Fragmenting jet functions were introduced in Refs. [19,20] to describe the fragmentation of a hadron inside a quark or gluon jet. A first application to heavy quarks was discussed in Ref. [21]. FJFs combine the fragmentation function, given in Eqs. (2.12) and (2.13), with the inclusive jet function, given in Eqs. (2.28) and (2.29). This can be explicitly seen in their definition [19,20]: χ n and B n⊥ are the gauge-invariant fields defined in Eq. (2.11). Antiquark FJFs are defined in a similar way, by exchanging the fields χ n andχ n in Eq. (2.36). As in Eqs. (2.28) and (2.29), the large component of the jet momentum is ω, and ωr is the jet invariant mass. However, differently from inclusive jets, in the definition of FJF a heavy hadron H in the final state is singled out, and its momentum p − h = ωz is measured. The FJFs thus depend on the jet invariant mass ωr, on the momentum fraction z and on the heavy quark mass m Q .
The FJFs G H i have several important properties, which were proven for light partons in Ref. [19,20,52] and which we now discuss briefly.
The first relationship states that after integrating over z and summing over all the possible emitted particles, one should recover the inclusive jet function. This is guaranteed by the momentum [19,20] and flavor [52] sum rules obeyed by the FJFs. The momentum conservation sum rule states that 1 2(2π) 3 where the sum is over a complete set of states. The flavor sum rule for a quark is [52] These relations are valid both in the approximation ωr m 2 Q , and in the regime ωr ∼ m 2 Q . In the former case, J i are the inclusive quark and gluon jet functions, computed with massless quarks [41,43,49,53]. If ωr ∼ m 2 Q , J Q is the massive jet function of Ref. [44,45]. In both cases, the mass dependence of the jet function on the r.h.s. of Eqs. (2.38) and (2.39) is not singular.
The second property arises again due to the similarity of FJFs with inclusive jet functions. In the UV, the FJFs G H i look like inclusive jet functions, initiated by the parton i. In particular, the restriction on the final state, requiring the identification of the hadron H, does not affect the UV poles of the FJF, so that G H i have the same renormalization group equation as quark or gluon inclusive jet functions [19,20] The anomalous dimension γ J i is identical to the inclusive case, given in Eq. (2.31), and in particular is independent of the momentum fraction z and the mass of the heavy quark m Q . The resummation of log r τ proceeds as in the inclusive case discussed in Section 2.2, albeit with a different initial condition. The final relation is due to the fact that the IR sensitivity of the FJFs is completely captured by the unpolarized fragmentation functions. Therefore, if the jet scale and the heavy quark mass are well separated, ωr m 2 Q , at leading power in m 2 Q /ωr one can factorize the dynamics at the two scales by matching the FJFs onto fragmentation functions The coefficients J ij depend on the jet invariant mass, and on the momentum fraction, but are independent of the heavy quark mass, up to power corrections of the size m 2 Q /ωr. Eq. (2.41) is very similar to the relation between the beam functions and the parton distributions in Eq. (2.35).
We will further discuss the properties ( 3 Review of how to resum logarithms of r Q and r τ Consider single inclusive production of one (light) hadron h in pp collisions, pp → h + X. Beyond leading order, the partonic cross section contains collinear divergences, when additional emissions become collinear to initial or final state partons. The divergences are physically cut off by non-perturbative physics, and they need to be absorbed into non-perturbative matrix elements, parton distribution functions for the partons in the initial state, and fragmentation functions for the final state.
In the case of single inclusive hadroproduction of heavy hadrons, pp → H + X, the final state collinear divergences in the partonic cross section are cut off by the heavy quark mass m Q , a perturbative scale. Identifying the heavy hadron with a heavy quark, the cross section for the production of a heavy hadron, differential in the hadron p T and rapidity y, can be expressed as a convolution of the partonic cross section for the production of a heavy quark and two parton distribution functions for the incoming partons [2] dσ dp 2 T dy Here the functions f i (x a , µ 2 ) denote the standard parton distributions functions, while dσ ij denotes the partonic cross section for a parton i and parton j to scatter into a heavy quark with transverse momentum p T and rapidity y. We have omitted the dependence of the shortrange cross section on the momentum fractions x a,b , and on the heavy quark rapidity. The final state collinear divergences present in the massless case manifest as logarithms of the heavy quark mass in Eq. (3.1). As the energy increases, logarithms of r Q in the partonic cross section become large, threatening the validity of the perturbative expansion. In order to resum them, one needs to factorize the partonic cross section into two separate pieces, each of which depends on only one of the two scales m Q and p T . This is achieved by introducing a fragmentation function [7] dσ dp 2 T dy In Eq. (3.2), the short-range cross section is the cross section for the production of the parton k with transverse momentump T and rapidity y in the collision of partons i and j, computed with all partons considered massless. The parton k then fragments into a heavy hadron H, carrying a transverse momentum p T = zp T , and the same rapidity as the original parton y. If the short-range cross section and the fragmentation function are evaluated at their characteristic scale, respectively µ ∼ p T and µ ∼ m Q , no large logarithms arise in the perturbative expressions. Of course, in the end all functions have to be evaluated at a common scale µ, and one therefore has to use the RGE to evolve each function to this scale. The RG evolution of the fragmentation function is determined by the DGLAP equation, Eq. (2.14). By evolving the fragmentation function from µ ∼ m Q to µ ∼ p T one can sum the logarithms of m Q /p T .
As already noted, the short-range cross section in (3.2) is calculated in the limit m Q = 0. Thus, while this approach correctly resums the logarithms of m Q /p T , it does not contain any dependence on powers of the same ratio. Since the power dependence is correctly reproduced using (3.1), one can obtain an expression that correctly reproduces both the logarithmic and the power dependence on m Q /p T by combining the two ways of calculating. This is the approach taken in FONLL [8]. Now consider jet cross sections. As in the previous case, the starting point for a resummation of the large logarithms that arise in cross sections that are differential in a jet resolution parameter τ is the separation of the dimensionful variables whose ratio gives the value of τ . This is achieved by a factorization of the cross section. There are many jet resolution variables one can choose, and a large body of literature how to obtain the relevant factorization theorems. Since all approaches in the end contain the same physics, and the final factorization theorems look very similar, we simply state the result here for one specific resolution variable, namely N -jettiness τ N [54]. For the purposes of this discussion, the only relevant part of the definition of N -jettiness is that τ N has dimension one, τ N → 0 as we approach N pencil-like jets, and that τ N is linear in the contributions from each jet (both from initial and final state radiation) and soft physics τ N = τ (a) N . The factorization theorem can be written schematically as [54] dσ dp and we have only included the dependence on terms that are relevant for our discussion. H is the hard function for the production of the partons with flavor k 1 , . . . k N , and it depends on the p T of the N signal jets. Collinear radiation in the final state is described by the inclusive jet functions J k i , while two beam functions B describe initial state radiation from the incoming beams, initiated by the partons of flavor a and b. The jet and beam functions are function of the jet invariant mass Qτ N , where Q is of the size of the hard scattering scale. In addition, the beam functions depend on the momentum fraction x of the incoming partons. The soft function describes soft interactions between jets, and between jets and the beams. It depends on soft momenta ∼ τ (s) N . After evolving the parton distribution functions to the jet scale, as discussed in Section 2, each function in the factorization theorem (3.3) only depends on a single scale, thus one can again calculate each term at its characteristic scale without encountering any large logarithms, and then evolve them to a common scale using the RGEs discussed in Section 2.2.
The factorization formula in (3.3) is derived in the limit τ N → 0, such that no power corrections of the ratio τ N /p T can be included. In order to derive an expression that is valid in both the limits of small and large τ N /p T , one needs to combine the resummed result with the known fixed order expression, which includes this power dependence.

A combined resummation of r Q and r τ
In this section we give the factorization theorems that are required to combine both types of resummation, such that one can study the production of heavy flavor at high energy in the presence of jet vetoes, or perhaps more importantly, such that one can combine calculations which resum the dependence on the heavy quark mass with parton shower algorithms. The later sections in this paper are then devoted to calculating the new ingredients in the resulting factorization theorems perturbatively.
We will consider two separate cases. The first is the production of identified heavy flavored hadrons in hadronic collisions, with measured momentum of the heavy hadron. Examples are pp → H + X or the associated production of a heavy flavored hadron and a weak boson, pp → W + H + X. In particular, we consider the case in which the heavy hadron is part of an identified jet, and a jet veto limits the total number of jets in the event. The momentum of the heavy hadron is characterized by its fraction z of the total jet momentum it is part of. A second interesting application is the production of jets (identified by a regular jet algorithm), which are tagged as b jets, and again extra jet activity is vetoed. Since b-tagging algorithms rely on the presence of at least one weakly decaying b-flavored hadron, the situation is related to the previous case. The main difference is that the momentum of the B hadron is not measured in this case, and the momentum fraction z is therefore integrated over.
The factorization theorem is based on the FJFs, defined in Refs. [19,20] for light hadrons, and reviewed in Section 2.3. The main new ingredient in this work is to extend the idea of a FJF to the case of heavy quarks, in which case infrared singularities that were present in the light FJFs manifest themselves as a logarithmic dependence on the heavy quark mass.
In addition to cases we discuss, logarithmic dependence on m Q appears in the flavor excitation channel. In this channel, one heavy quark is present in the initial state and enters the hard collision. This is similar to the cases discussed above, with the difference that the production of the heavy flavor happens in the initial rather than the final state. In this case log r Q are resummed by introducing a perturbative b-quark parton distribution at the scale m Q , and running it with the DGLAP equation up to the hard scattering scale. Initial state radiation at a scale t m 2 Q can be studied using the same techniques developed in this paper, by introducing a heavy quark beam function. We leave a detailed discussion for future work, and will not discuss initial state splitting any further.

The production of an identified heavy hadron
We consider first the case of production of an identified heavy flavored hadron in the presence of a veto on extra jet activity. As already discussed, the momentum of the heavy hadron is measured to have a fraction z of the momentum of the jet it is part of. The extra jet activity is vetoed using a jet resolution variable τ N , where τ N is defined such that it goes to zero when there are at most N pencil-like jets present. Phenomenologically interesting applications are the two-jettiness cross section in pp → QQ + X, or the one-and two-jettiness cross sections for pp → W + Q + X.
In the limit of small τ N , the factorization theorem for the cross section differential in τ N and in the p T and rapidity of the observed hadron can schematically be written as dσ dp 2 This factorization theorem is almost identical to the one given in Eq. (3.3), and, as in that case, it holds up to power corrections in τ N /p T . The only difference is that a hadron H is observed inside the jet initiated by the parton k N , and its p T and rapidity are measured. To be able to describe this extra information, the inclusive jet function J k N needs to be replaced by the FJF G H k N . In addition to the argument τ N , describing the contribution of the inclusive jet to the jet resolution variable, the fragmenting jet function depends on the mass of the heavy quark m Q as well as the momentum fraction of the heavy hadron z.
As discussed in Section 2.3 the RGE of the FJF, Eq. (2.40), is identical to that of an inclusive quark or gluon jet function, so that the resummation of log τ N proceeds as in the inclusive case.
The FJFs are two-scale objects, sensitive to the jet invariant mass and to the heavy quark mass m Q . Differently from the light parton FJFs discussed in Refs. [19,20], the heavy quark FJFs can be computed purely in perturbation theory. If the jet scale Qτ N in Eq. (4.1) is close to m 2 Q , the fixed order expression for the FJFs at the scale µ 2 I ∼ Qτ N does not contain large logarithms, and the evolution (2.33) resums logarithms of τ N . On the other hand, if Qτ N m 2 Q , there is no choice of initial scale the minimizes the logarithms in the FJFs, and the initial condition for the jet evolution is still plagued by large logarithms. However, the IR sensitivity of the FJFs is completely captured by the unpolarized fragmentation functions. Therefore, if the jet scale and the heavy quark mass are well separated one can factorize the dynamics at the two scales by matching the FJFs onto heavy quark fragmentation functions D H i , as in Eq. (2.41). Since each term on the right-hand side of Eq. (2.41) depends only on a single scale, the logarithms of log m 2 Q /(Qτ N ) of the FJFs are reproduced through logarithms of Qτ N /µ 2 F and m 2 Q /µ 2 F on the right hand side, such that they can be resummed through RG evolution. Evolving the fragmentation function from the mass scale to a scale of order Qτ N , no large logarithms are left in the initial condition for the FJF evolution. Then, running the hard, beam, soft and jet functions to a common scale, all large logarithms in Eq. (4.1) are correctly resummed.
By matching the FJFs onto fragmentation functions and the beam functions onto parton distributions, we can recast Eq. (4.1) in a form that stresses the relation to single inclusive production discussed in Section 3. dσ dp 2 This form is very similar to Eq. Our goal is to define a heavy quark tagged (Q-tagged) jet function J Q i (ωr, m 2 Q , µ 2 ) with these features, such that one can use the factorization formula given in Eq. (3.3) and simply replace the standard jet function by its heavy quark tagged version. A heavy quark tagged jet function is agnostic as to which type of hadron gives rise to the long decay time and therefore the secondary vertex. Furthermore, it is insensitive to the momentum fraction of the heavy hadron, as long as the transverse momentum is above the minimum transverse momentum imposed in the b-tagging algorithm. Therefore, the b-tagged jet function can be obtained from the heavy flavor FJF by summing over all heavy flavored hadrons as well as integrating over the momentum fraction of the heavy hadron (down to a cutoff z 0 related to the minimum p T ).
As is the case for a fragmentation function, the hard interaction does not necessarily need to involve the production of a heavy quark Q, since this can be produced from the splitting g → QQ in the radiation happening within the jet. Thus, heavy quark tagged jets can be initiated by any possible flavor. In the case of heavy quark initiated jets we define where the subtraction of the antiquark contribution avoids the double counting of configurations in which the heavy quark splits in an additional QQ pair, Q → QQQ. When z 0 approaches 0, a heavy quark initiated jet should always be tagged. In virtue of the flavor sum rules obeyed by the quark fragmentation function and FJF, Eq. (4.3) does indeed guarantee that for z 0 → 0 one recovers an inclusive quark jet. In particular, any dependence on the fragmentation function, and thus on logarithms of the mass, disappears. For J Q Q there is no need to resum DGLAP logarithms, while Sudakov double logarithms of r τ are resummed by the evolution of the inclusive quark jet function. Powers of m 2 Q /(Qτ N ) can be retained by using inclusive massive jet functions [44,45].
For gluon and light quark initiated jets, i ∈ {g, l}, heavy quarks are always produced in pairs. In this case we define where, in the last step, we used charge conjugation invariance. With this definition, J Q g,l count the multiplicity of QQ pairs in a gluon or light quark jet of invariant mass ωr. 2 Since the anomalous dimension of the FJF is z-independent, the RGE of the Q-tagged jet is still identical to that of the inclusive jet function. For gluon or light quark initiated jets the dependence on the fragmentation function does not drop out. For small values of Qτ N ∼ 4m 2 Q , this does not cause problems. The only large logarithms in this case are log τ N , which are resummed by using the fixed order expression of J Q g,l as initial condition for the jet evolution (2.33). For larger values of 4m 2 Q Qτ N Q 2 , resummation of log m 2 Q /(Qτ N ) does become necessary and is achieved by running the fragmentation function to the scale Qτ N .

Heavy Quark Fragmenting Jet Functions at O(α s )
We now discuss in more detail the FJFs of massive quarks, illustrating the general properties discussed in Section 2.3 with one loop examples. We work in perturbation theory, identifying the hadron H with one of the parton species, {Q,Q, g}. For heavy quark production, the most interesting functions are those with an identified Q orQ. For completeness, and to  verify the cancellation of mass dependent terms in the momentum sum rule (2.38) for the quark and gluon FJFs, we also consider the effects of the heavy quark mass on the FJF of a heavy quark into a gluon, and of a gluon into a gluon, even though these FJFs are of less practical interest.
We expand the FJFs G and the matching coefficients J in powers of α s /(2π) , At tree level the heavy quark FJF, G Q Q , and the gluon FJF, G g g , are the product of a delta function on the jet invariant mass, and a delta function on the observed momentum fraction, 2) with the factor of 2(2π) 3 due to the choice of normalization of Refs. [19,20]. All other FJFs vanish at tree level.
At order O(α s ), G Q Q and G g g receive corrections from virtual one loop diagrams, and real diagrams, with the emission of an additional parton. We show the diagrams contributing to G Q Q in Fig. 1, and to G g g in Fig. 2. At this order, one finds the first contributions to G g Q and G Q g . They originate purely from real emissions, the real diagrams in Fig. 1 for G g Q , and the diagram in Fig. 3 for G Q g . We compute the diagrams in Figs. 1, 2 and 3 with a finite quark mass m Q , and take the limit m 2 Q ωr at the end of the calculation. We present here the results in this limit, which we refer to as "massless limit", and relegate the one loop expressions for finite m Q to Appendix B.
The individual diagrams contributing to G Q Q contain UV and IR divergences. We regulate UV divergences in dimensional regularization, and IR divergences by introducing a ∆ regulator, as defined in Ref. [57]. Double counting of the collinear and usoft regions is eliminated via zero-bin subtractions [58]. We choose different regulators in the IR and UV to explicitly check that, after zero-bin subtraction, all infrared divergences cancel between real and virtual emission diagrams, and the remaining 1/ε poles are UV in nature.
The UV divergences of the diagrams in Fig. 1 are canceled by introducing the counterterm Z Jq relating the renormalized and unrenormalized FJF with, at one loop, As expected Z Jq does not depend on m Q , which is an IR scale in SCET. Z Jq is also independent of the momentum fraction z, implying that the evolution of the quark FJF from the jet scale to the hard scale does not change the shape of the momentum fraction distribution.
Z 2 is the quark field renormalization. At one loop, in the MS scheme and in Feynman gauge, Eq. (5.5) has the same form as Eq. (2.31). The coefficient of the plus distribution is the one-loop value of the quark cusp anomalous dimension Γ q = α s /(4π)(4C F ) [46,47]. The coefficient of the delta function reproduces the one-loop value of the non-cusp component of the anomalous dimension of the inclusive quark jet function, γ q = α s /(4π)(6C F ). The anomalous dimension γ Jq is thus identical to the anomalous dimension that governs the evolution of inclusive quark jets.
In the massless limit, the quark FJF depends on the jet invariant mass ωr through plus distributions of the form d(ωr) θ(ωr) log n (ωr) ωr Since the integration range extends to infinity, one has to introduce an arbitrary cut-off ωκ in the subtraction term, ϕ(0). The cut-off dependence is canceled by the second line of Eq. (5.7), and the distributions do not depend on ωκ. It is also convenient to introduce the notation θ(ωr) log n (ωr/µ 2 ) ωr (µ 2 ) + ≡ θ(ωr) log n (ωr/µ 2 ) ωr + + (−1) n+1 1 n + 1 log(µ 2 )δ(ωr). (5.8) In terms of the distribution in Eq. (5.8), the renormalized heavy quark FJF at one loop is (5.9) The gluon FJF G g g is affected by the quark mass only via quark loop corrections to the gluon propagator, the last virtual diagram in Fig. 2. The remaining diagrams in Fig. 2 are unchanged with respect to the massless case discussed in Ref. [20], and we did not evaluate them. The correction to G g g is obtained by multiplying the tree level result by the contribution of massive quarks to the residue of the gluon propagator, and we obtain where G g g (ωr, z, µ 2 ) light is the gluon FJF computed with n l massless quarks. The same correction to the gluon propagator affects the fragmentation function for a gluon into a gluon, D g g , so that it cancels in the matching and the matching coefficients are mass independent. The heavy quark mass does not affect the anomalous dimension of G g g , which, as showed in Ref. [20], is the same as that of an inclusive gluon jet.
At order O(α s ), the first contributions to heavy quark fragmentation into a gluon, G g Q , and gluon fragmentation into heavy quark, G Q g , arise. G g Q receives contributions from the real emission diagrams in Fig. 1, when one integrates over the heavy quark phase space and fixes the gluon momentum fraction to z. The lowest order diagram contributing to G Q g is showed in Fig. 3. The diagrams are UV and IR finite, and, in the massless limit, they give Notice that Eqs. (5.11) and (5.12) are independent of µ. The evaluation of the anomalous dimension of G g Q and G Q g requires the calculation of UV poles at O(α 2 s ). We explicitly verify in Section 6.1 that G Q g has, as expected, the same RGE as an inclusive gluon jet. FJFs for heavy antiquarks, GQ Q , G ḡ Q and GQ g have the same expressions as the quark FJFs in Eqs. (5.9), (5.11) and (5.12). FJFs of a heavy antiquarkQ, or of a light quark (or antiquark) l into a heavy quark Q vanish at O(α s ).
The only dependence on the quark mass in Eqs.
At one loop, the quark matching coefficients J QQ are Eq. (5.14) reproduces the results in Ref. [20,21], as one expects, since all the IR dependence should cancel in the matching. Similarly, we obtain which agree with Ref. [20]. The gluon matching coefficient J gg is also unaffected by m Q , since the correction to gluon propagator cancels between the gluon FJF and fragmentation function. For completeness, and since it is needed in the O(α 2 s ) calculation, we report the result of Ref. [20].
with p gg (z) given by The explicit expressions for the FJFs, Eqs. (5.9)-(5.12), allow to check the sum rules (2.38) and (2.39). At the perturbative level, the sum over a complete set of states in Eq. (2.38) is a sum over partons, H ∈ {Q,Q, g, l,l}. Integrating the fixed order results (5.9), (5.11) over z one finds that that agrees with the massless quark inclusive jet function at one loop, given in Ref. [42]. In Appendix B we prove the analogous relations in the regime ωr ∼ m 2 Q . The momentum and flavor conservation sum rules are not affected by the evolution of the fragmentation functions. Using the momentum and flavor sum rules for the fragmentation function, Eqs. The dependence on the heavy quark mass cancels in the inclusive gluon jet function. The combination dzz G Q g (ωr, z, m 2 Q , µ 2 ) + GQ g (ωr, z, m 2 Q , µ 2 ) + G g g (ωr, z, m 2 Q , µ 2 ) + 2n l G l g (ωr, z, µ 2 ) (5.21) is indeed equal to the inclusive gluon jet function, and mass independent, up to power corrections. However, if one insists on tagging the heavy quark, she is left with some mass dependence. Retaining terms of O(α s ) in the matching coefficients J ij , the Q-tagged jet function is (5.22) Logarithms of the ratio m 2 Q /(Qτ N ) are resummed at NLL accuracy by solving the DGLAP equation, with two-loop splitting functions and one loop initial conditions for D Q g and D Q Q at µ 0 ∼ m Q , given in Eqs. (2.25) and (2.23).
Setting z 0 = 0 in Eq. (5.22) would express J Q g in terms of the first Mellin moment of the quark and gluon fragmentation function. However, the DGLAP equation for the first moment of D Q g is not well defined, since the Mellin transform of P gg has a pole for N = 1. In common approaches for the study of heavy quark multiplicity in gluon jets [59,60], the DGLAP equation is modified at small z to regulate the singularity of P gg , by including coherence effects. In this work, we will assume that z 0 is large enough that small z effects can be neglected.

Gluon and light quark fragmentation into heavy quarks at O(α 2 s )
We have seen in Section 5 that gluon initiated Q-tagged jets are sensitive to the scale of the quark mass, and that large logs of the ratio m 2 Q /(Qτ N ) can be resummed by solving the DGLAP equation for the fragmentation function. In this section we calculate the gluon and light quark FJFs into a heavy quark at O(α 2 s ), both of which involve gluons splitting into QQ pairs. In the case of the gluon FJF G Q g , the leading order is O(α s ), and the calculation of this section amounts to the NLO contribution. The light quark FJF starts at O(α 2 s ), and we give here the leading order term.
There are several reason to go beyond the lowest order in processes involving gluon splitting into heavy quark pairs. First of all, one can study the renormalization group properties of G Q g . We explicitly show that the RGE for G Q g is the same as for the gluon inclusive jet. Furthermore, knowing the fixed order expression of G Q g at NLO, together with two loop cusp anomalous dimension and one loop non-cusp anomalous dimension, allows to resum Sudakov double logarithms of r τ at NNLL accuracy. Finally is interesting to explicitly check that at O(α 2 s ) the infrared sensitivity is exactly reproduced by the heavy quark fragmentation functions, computed at O(α 2 s ) in Ref. [36]. But perhaps the most important reason for this calculation is that for many interesting SM processes involving b jets fixed order calculations are available at NLO accuracy [61,62]. In order to match the resummed result to these fixed order results, the knowledge of the gluon and the light quark FJFs into heavy quarks at O(α 2 s ) is required.  6.1 Gluon fragmentation into heavy quarks at O(α 2 s ) The virtual and real diagrams contributing to G Q g at O(α 2 s ) are shown, respectively, in Figs. 4 and 5. We decompose G Q g in terms of color factors where n l = n f − 1 is the number of light quarks. At O(α 2 s ), the matching condition Eq. (2.41) reads where we used that at tree level D , has been computed in Ref. [36]. Therefore the calculation of G Q g allows for the extraction of the matching coefficient J gQ at O(α 2 s ). We stress that the matching is meaningful if the coefficients are independent of m Q , which is an infrared scale in the problem, and, thus, all the singular mass dependence of G Q g needs to be reproduced by the fragmentation functions. We will see that this is indeed the case.
We use dimensional regularization to regulate UV and IR divergences and work in Feynman gauge. We compute the diagrams in Figs. 4 and 5 analytically, for finite value of m Q , and take the massless limit, m 2 Q ωr, at the end of the calculation. This limit has to be taken carefully, and we refer the reader to Appendix A for details. We find that G Q g depends on ωr through plus distributions [log n (ωr)/(ωr)] + , defined in Eq. (5.7), with n ≤ 2. As discussed in more detail in Appendix A, the coefficients of the plus distributions and the logarithms of the mass in the δ(ωr) piece are determined by the "naive" massless limit of G Q g . The mass independent component of δ(ωr), on the other hand, requires to integrate the result for G Q g , obtained at fixed m Q , from the minimum invariant mass required for the production of a QQ pair, ωr = m 2 Q /(z(1 − z)), to ∞. Most of the integrals can be performed analytically, and expressed in terms of polylogarithms up to rank three. A few integrals from the real emission diagrams with color structure C A T R had to be solved numerically.
G C F T R g receives contributions from the virtual diagrams (a), (b) and (c) in Fig. 4, and from the square of the real diagrams (a), (b) and (c) in Fig. 5. The ultraviolet divergences in these diagrams are canceled by charge and mass renormalization, while infrared divergences cancel between the virtual and real emission diagrams. The color structures T 2 R and T 2 R n l receive contributions from heavy and light quark loop corrections to the gluon propagator, diagram (d) in Fig. 4. The diagrams are IR finite, and the UV divergences are renormalized by charge renormalization.
The situation is more interesting for G C A T R g . Even after charge renormalization, the result is still divergent and are rendered finite only by an operator renormalization. Defining the jet renormalization Z Jg , in the same way as in Eq. (5.3), one finds This leads to the RGE for G Q g d d log µ G Q g (ωr, z, m 2 Q , µ 2 ) = d (ωs) γ Jg ωr − ωs, µ 2 G Q g (ωs, z, m 2 Q , µ 2 ), (6.4) with anomalous dimension Z 3 is the gluon field strength renormalization in the MS scheme, which at one loop and in Feynman gauge is given by Eq. (6.5) is of the form (2.31), and it is the same anomalous dimension that governs the running of an inclusive gluon jet function. As in the case of the quark FJF, the anomalous dimension of G Q g is mass and momentum fraction independent. We now give the expression of the renormalized gluon FJF G Q(2) g , in the massless limit. We work in the pole mass scheme, and use as mass counterterm which subtracts the entire one loop correction to the quark mass. The color structure G C F T R g is given by The massless limit of the color structure C A T R is given by The functions g C F T R (z) and g C A T R (z) are functions of z only, independent on the mass. We plot them in Fig. 6. We give the analytic expression of g C F T R in Eq. (C.1). A few contributions to g C A T R from the interference of the real diagrams (d) and (e) with (a), (b), and (c) had to be computed numerically, so that we do not have the full analytic expression of this function. g C F T R (z) and g C A T R (z) are not smooth at z = 1/2, and they exhibit the same behavior as the gluon fragmentation function, D Q(2) g , described in Ref. [36]. The non-smooth terms can be traced back to the virtual diagram 4(a), with color structure C F − C A /2 and are related to the production of the heavy quark pair at threshold. In matching the FJFs onto fragmentation functions, the non-smooth terms cancel and the matching coefficients are smooth functions of z.
The color structures T 2 R n l and T 2 R arise from light and heavy quark loop corrections to the gluon propagator, respectively. In the massless limit, the plus distribution and the mass dependence of G T 2 R n l g and G T 2 R g are the same, the only difference between the two functions is in the mass-independent part of the coefficient of the delta function. Thus, we can write  qualitatively similar to G T 2 R g . To capture the contribution of terms proportional to δ(ωr), we integrate massive and massless FJFs against a constant test function ϕ(ωr), with cut-off ωκ ϕ(ωr) = θ(ωκ − ωr). (6.13) ωκ is representative of the jet scale, and we set the renormalization scale µ 2 = ωκ. The integration of FJFs in the massless limit can be carried out very easily, while we integrate the massive FJFs numerically. The massive jet has an additional theta function, that constrains the jet invariant mass to be larger than m 2 Q /(z(1−z)), the minimum invariant mass to produce a QQ pair at fixed z. We show results for two choices of jet invariant masses, ωκ = (50 GeV By comparing our results with the expressions for the fragmentation functions given in Ref. [36], one finds that the mass dependence of G Q g is exactly reproduced by the one and two loop fragmentation functions, as it should. It therefore cancels in the matching, leaving mass independent matching coefficients. For the color structure C F T R , the matching coefficient (6.14) with (3)). (6.15) For the color structure C A T R , we find +δ(ωr) I C A T R (z). The matching condition (6.2) implies that the function I C A T R (z) is obtained subtracting from g C A T R the C A T R components of the fragmentation function D Q(2) g , evaluated at µ = m Q . In the notation of Ref. [36] where F (C A T R ) g is given in Eq. (21) of Ref. [36]. Setting µ 0 = m Q eliminates the logarithmic terms in the fragmentation function, leaving only the finite pieces. The functions I C F T R and I C A T R are plotted in Figs. 6, and are smooth functions of z.
The difference between G T 2 R g and G T 2 R n l g is accounted for by the difference between the contributions of light and heavy quark loops to D Q(2) g , and thus it cancels in the matching. The matching coefficient is the same for heavy and light flavors and we find In Fig. 10, we show the effects of the DGLAP evolution on the Q-tagged jet function J Q g . We fix ω = 100 GeV, and integrate the FJF G Q g from a minimum z 0 = 0.05 to 1. To show the effects of terms proportional to δ(ωr), we integrate the jet function with a constant test function (6.13). We set the renormalization scale µ 2 = ωκ. The solid blue line denotes the fixed order result, obtained integrating Eq. (5.12) between z 0 and 1. The dashed-magenta line uses Eq. (5.22), with the fragmentation functions D Q g and D Q Q evolved from the scale µ 2 0 ∼ m 2 Q , to µ 2 = ωκ. We work at NLL accuracy, using two-loop time-like splitting functions, in the conventions of Ref. [31]. 3 The dot-dashed yellow and dotted green lines show the fixed-order O(α 2 s ) expression, obtained integrating Eq. (5.12), and Eqs. (6.8), (6.9), (6.10), from z 0 to 1. The dotted green line is the complete O(α 2 s ) result, while the yellow dot-dashed line includes only the logarithmic enhanced terms. From Fig. 10 we see that the logarithmically enhanced terms in the NLO result are, for the moderate value of the jet invariant mass showed here, reproduced by the DGLAP evolution of the quark and gluon fragmentation functions with O(α s ) initial condition. The non-logarithmic terms in G Q g also provide an important correction to the Q-tagged jet function.
6.2 Light quark fragmentation into heavy quarks at O(α 2 s ) We compute in this section the FJF for light quark fragmenting into a heavy quark at O(α 2 s ). In the case of light quark fragmentation, the only color structure at this order is C F T R , and we write G The matching condition (2.41), specified to H = Q, reads ll are delta functions. The one loop result for the matching coefficient J lg is the same as for heavy quark, and it is given in Eq. (5.15), while the O(α s ) fragmentation function D Q(2) l is given in Ref. [35]. Our calculation, then, allows to extract J lQ at O(α 2 s ). Figure 11. Diagrams contributing to the light quark FJF G Q l . The plain line denotes the light quark that initiates the jet, while the dashed lines heavy collinear Q andQ.
The diagrams in Fig. 11 are UV and IR finite, and we find (1 + z) − 7 9z + 1 9 −60 + 42z + 25z 2 log(1 − z) The logarithmic dependence on m Q is canceled by the fragmentation function, and the match-ing coefficient is (6.23)

Conclusion
In this paper we outlined a framework for the simultaneous resummation of logarithms of the heavy quark mass m Q and of a jet resolution variable τ N . These logarithms arise in heavy quark production, when, for experimental or theoretical reasons, the final state is not fully inclusive. Examples are hadroproduction of heavy flavored hadrons, when one hadron is detected and its momentum measured, and the final state is restricted to contain a maximum number of jets, or b jets cross sections, when the tagged b orb is found in a jet initiated by a light parton. Furthermore, a resolution variable is always needed when combining NLO calculations with parton shower MonteCarlo. The possibility to carry out the two resummations simultaneously is an important step towards the combination of FONLL-like calculations, which resum the logarithmic dependence on the heavy quark mass, with parton shower algorithms, which require a resummation on a jet resolution variable. Simultaneously resumming the logarithms of m Q and τ N requires a factorization theorem that separates these two scales. We discussed the generic structure of factorization theorems, using the inclusive event shape N -jettiness as the jet resolution variable. The crucial ingredient in this factorization theorem is the fragmenting jet function for a heavy quark, which describes a jet of particles with fixed invariant mass ωr containing a heavy-flavored hadron of given momentum fraction z. The heavy quark FJFs are generalizations of the massless parton FJFs introduced in Ref. [19,20]. Up to power corrections of size m 2 Q /(Qτ N ), the heavy quark FJFs can then be further factorized into a heavy quark fragmentation function that contains only the dependence on m Q and a matching coefficient that contains only the dependence on the jet resolution scale τ N . We computed the heavy quark FJFs at O(α s ), giving their expressions in the massless limit m 2 Q Qτ N in Eqs. (5.9), (5.11), and (5.12), and their full mass dependence in Eqs. (B.1), (B.2) and (B.3). Our calculation explicitly verified that to the calculated order the heavy quark fragmentation functions reproduce the entire dependence on the heavy quark mass, such that the matching coefficients are independent of m Q and reproduce the results of Ref. [20,21]. We have also shown that the anomalous dimension of the FJFs are independent of the mass m Q , and are in fact equal to the anomalous dimension of the inclusive jet functions. Thus, the resummation of the jet resolution variables is identical to the more familiar case of massless jet functions.
Final state splitting of gluons into heavy quarks is a particularly important source of uncertainty in heavy quark production. They are included at NLL accuracy in FONLL [8], which, however, can be applied only to fully inclusive final states. For more exclusive observables, one might rely on NLO plus parton shower Monte Carlo programs. While the shower evolution resums all the leading logarithms of m Q originating by emission of light partons from massive legs, MC@NLO and POWHEG only include gluon splitting at fixed order [13,14]. Fixed-order NLO calculations of b-jet cross sections are also affected by log m Q originating in final state splittings of gluons, or light quarks, into heavy quarks [16,61]. While the fixed-order approach is sufficient for b jets with moderate p T , at high p T we expect the resummation to become important.
To compare with fixed order NLO calculations, and to improve on Monte Carlo parton shower, we computed the gluon and light quark FJFs into heavy quarks to O(α 2 s ). We give their expressions in the massless limit in Eqs. (6.8), (6.9), (6.10) and (6.22). The calculation of the gluon FJF at O(α 2 s ), that is at NLO, is particularly interesting because it allows to explicitly check that the RGE of G Q g is that of an inclusive gluon jet. We verified that the singular mass dependence of G Q g and G Q l is reproduced by the heavy quark fragmentation functions at O(α 2 s ), and that the matching coefficients J (2) gQ and J (2) lQ are independent on the quark mass, up to power corrections. The FJFs G Q(2) g in Eqs. (6.8), (6.9) and (6.10), and G Q(2) l in Eq. (6.22), and the matching coefficients in Eqs. (6.14), (6.16), (6.19) and (6.23) are the main results of this work.
Our calculation of the heavy quark FJFs provides a key ingredient for using the factorization formulae (4.1) and (4.2) to describe phenomenologically interesting heavy quark production cross sections. In most cases, the remaining functions in Eqs. (4.1) and (4.2) can be computed with massless quarks without encountering divergences, and exist in the literature. If some of the inclusive quark jets in Eq. (4.1) and (4.2) are massive, one should use massive quark jet functions [44,45]. The effects of the heavy quark mass on inclusive light quark jet functions, and on the thrust hemisphere soft function, which start at O(α 2 s ) and are relevant in the case the jet or soft scale are close to m Q , have been considered in Refs. [64,65] The exception is the flavor excitation channel. In this case logarithms of m Q originate in initial state splittings of gluons and light quarks into QQ pairs, and one of the two heavy quarks enters the hard collision. These logarithms are resummed by introducing a perturbative b-quark parton distribution, and evolving it to the hard scale. In the presence of a resolution variable m 2 Q Qτ N p 2 T , one needs to introduce a heavy quark beam function, using techniques very similar to those discussed in this paper.
While we have given all perturbative ingredients required to perform the simultaneous resummation of logarithms of m Q and τ N , we have left this resummation for future work.
large and need not to be resummed. In this regime, it might be more important to retain the full dependence on m Q . In this section we provide the expressions of the quark and gluon FJF at one loop, keeping all powers of m 2 Q /(ωr). A first important observation is that the renormalization of the FJFs is not affected by assumptions on the relative size of the jet invariant mass and m Q . The UV divergences of diagrams Fig. 1, as well as those of the O(α 2 s ) diagrams, Figs. 4 and 5, are not changed. Consequently, even in the limit ωr ∼ m 2 thus will not be completely canceled by the fragmentation function. However this is just an artifact of the form of Eq. (B.1), that still contain expressions, as 1/(s + m 2 Q ) 2 , which need to be regulated in the m Q → 0 limit.
We now define the distributions in Eq. (B.1), applied to a test function factorized as ψ(s)ϕ(z) Eq. (B.11) reproduces the result for an inclusive quark massive jet [44,45]. Also in this case, one can take the limit m Q → 0 without encountering divergences. We also computed the expressions of G Q g and G Q l at O(α 2 s ), at fixed m Q . Their dependence on the jet invariant mass is not as simple as in Eqs. (5.11), (6.8), (6.9), (6.10), and (6.22). They are expressed analytically in terms of logarithms and polylogarithms up to rank two. Their expression is too lengthy to be reproduced here.