B-hadron hadro-production in NNLO QCD: application to LHC t\bar{t} events with leptonic decays
Michal Czakon, Terry Generet, Alexander Mitov, Rene Poncelet
PPrepared for submission to JHEP
B-hadron hadro-production in NNLO QCD:application to LHC t ¯ t events with leptonic decays Micha(cid:32)l Czakon a Terry Generet a Alexander Mitov b and Rene Poncelet b a Institut f¨ur Theoretische Teilchenphysik und Kosmologie, RWTH Aachen University,D-52056 Aachen, Germany b Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, UK
E-mail: [email protected] , [email protected] , [email protected] , [email protected] Abstract:
We calculate, for the first time, the NNLO QCD corrections to identified heavyhadron production at hadron colliders. The calculation is based on a flexible numericframework which allows the calculation of any distribution of a single identified heavyhadron plus jets and non-QCD particles. As a first application we provide NNLO QCDpredictions for several differential distributions of B hadrons in t ¯ t events at the LHC.Among others, these predictions are needed for the precise determination of the top quarkmass. The extension of our results to other processes, like open or associated B and charmproduction is straightforward. We also explore the prospects for extracting heavy flavorfragmentation functions from LHC data. Keywords:
QCD, Top-quark physics, NNLO Calculations, Fragmentation
TTK-21-05, P3H-21-011, Cavendish-HEP-21/02 a r X i v : . [ h e p - ph ] F e b ontents b -fragmentation in top-quark decay 115.2 b -fragmentation in top-quark pair-production and decay at the LHC 115.3 Extraction of B -hadron FFs from t ¯ t events 16 t ¯ t production and top-quark decay in-cluding fragmentation 21B Collinear counterterms for processes involving fragmentation 22C Checks on our computational setup 23 C.1 b -fragmentation in e + e − collisions 23C.2 Sum rules in top decay 24 The production of heavy flavors, like bottom and charm, is a cornerstone high-energy col-lider process. It offers a wealth of information about the Standard Model and representsan excellent tool for probing the QCD dynamics. Heavy flavor production has been exten-sively studied at past and present high-energy lepton and/or hadron colliders as well asin nuclear collisions where heavy flavors are a prominent probe of the underlying nucleardynamics.Heavy flavors are copiously produced at the LHC. Indeed, the b ¯ b and c ¯ c cross sectionsare among the largest at this collider. Such large production rates enable detailed andvery precise measurements in wide kinematic ranges. The theoretical description of theseprocesses, currently at next-to leading order (NLO) in QCD, is lagging in precision behindthe experimental needs. For improving the precision of theory predictions the inclusion ofthe NNLO QCD corrections is mandatory.– 1 –hen discussing the production of a heavy flavor of mass m at a hadron collider, it isinstructive to distinguish two kinematic regimes: the low p T regime where p T ∼ m and thehigh p T one where p T (cid:29) m . The low p T production of a heavy flavor can be described infixed order perturbation theory as an expansion in powers of the strong coupling constantevaluated at the scale m , i.e. α s ( m ), and including the full dependence of the heavy quarkmass m . For bottom, and especially charm, this expansion converges slowly since α s ( m ) isnot much smaller than unity. This expectation was confirmed by the recent fully-differentialNNLO QCD calculation of b ¯ b production at the quark level [1]. Such behavior is to becontrasted with t ¯ t production which is very similar technically but the smallness of α s ( m t )leads to a well-converging perturbative expansion [2, 3] through NNLO in QCD.The description of heavy flavor production at high p T involves a different set of chal-lenges. Fixed order perturbation theory is no longer adequate there since large quasi-collinear logarithms log( p T /m ) appear to all orders in perturbation theory and need to beresummed. The resummation of these logs can be consistently carried out in the so-calledperturbative fragmentation function (PFF) formalism [4]. Unlike the low p T case, a calcu-lation of heavy flavor production at high p T is performed with a massless heavy quark sincein the high-energy limit all terms that are power suppressed with m are negligible whilethe mass-independent terms as well as the logarithmically enhanced ones are automaticallyaccounted for by the PFF formalism. The current state of the art is NLO with next toleading logarithmic (NLL) accuracy. The goal of the present paper is to extend, for thefirst time, this description at hadron colliders to NNLO in QCD.A generic application to heavy flavor production that is valid in all kinematic regimeswould require the merging of the low p T and high p T descriptions mentioned above. Thishas been achieved at NLO in QCD within the so called FONLL approach [5]. Some of therecent hadron collider applications include refs. [6, 7]. Its generalization to NNLO goesbeyond the scope of this paper.As a first application of this formalism at NNLO in QCD we compute several B -hadrondifferential distributions in top quark pair production and decay at the LHC pp → t ¯ t + X → B + X . The reason for choosing this process is twofold: first, B -production is central totop quark physics and B -hadron related observables are a great tool for precise top quarkmass determination at hadron colliders. Second, in t ¯ t events the top quark mass providesa natural large hard scale such that for almost all distributions of interest the powersuppressed effects ∼ ( m b ) n are negligible. This makes this process an ideal applicationfor the massless b quark PFF formalism used in this work. B -hadron production in otherprocesses, like open B production at high p T , would be a straightforward extension of thecurrent work and we hope to report on it in future publications.This work is organized as follows: in sec. 2 we discuss the general features of theformalism for calculations with an identified hadron. In sec. 3 we explain our calculationalframework. In sec. 4 we introduce the B fragmentation functions used in this work. Sec. 5 isdevoted to phenomenological LHC applications. We study in detail B -hadron distributionsin top quark decay and in t ¯ t production and decay. We also propose an observable whichwe find suitable for extracting B -hadron fragmentation functions from LHC data. Severalappendices contain additional results. In appendix A we give the structure of the NNLO– 2 –ross section for the process pp → t ¯ t + X → B + X . In appendix B we give in explicitform the general expressions for the collinear counterterms needed for any NNLO hadroncollider process with fragmentation. Appendices C.1 and C.2 present two highly non-trivialchecks of our calculational setup: the calculation of B production in e + e − collisions whichis compared to the exact analytic result and the fulfillment of sum rules in top quark decay. A typical calculation in perturbative QCD involves final states with QCD partons, whichare clustered into jets, and colourless particles such as leptons. By clustering particlesinto jets, information is lost about the properties of the individual particles. On theexperimental side, it also introduces jet energy scale uncertainties, which can dominate thetotal uncertainty on jet-based observables (see e.g. ref. [8]), but are largely absent wheninstead measuring a single hadron’s momentum (e.g. ref. [9]). As an alternative to thisusual approach of jet-based observables, it therefore seems appealing to instead considerobservables involving the momentum of a single hadron, h .Perturbation theory alone cannot describe non-perturbative phenomena like the tran-sition from partons to hadrons, called fragmentation. The solution is to factorise thenon-perturbative aspects into fragmentation functions [10] in analogy to how parton dis-tribution functions are introduced to describe transitions from hadrons to partons in theinitial state. The fragmentation functions depend on the hadron h but are otherwise uni-versal and can thus be extracted from experimental data.The theoretical description of the production of an identified hadron proceeds as fol-lows. Standard tools and techniques are used to describe the production of on-shell partons.The partonic calculation is then extended by fragmenting the final-state partons, one ata time, into the observed hadron h which has a well-defined momentum p h . In practice,fragmentation corresponds to multiplying the fragmenting parton’s momentum with a mo-mentum fraction between 0 and 1, and then integrating the partonic cross section over itwith a weight given by the corresponding fragmentation function. This procedure is equiv-alent to convolving the differential partonic cross sections with fragmentation functions: dσ h dE h ( E h ) = (cid:88) i (cid:18) D i → h ⊗ dσ i dE i (cid:19) ( E h ) ≡ (cid:88) i (cid:90) dxx D i → h ( x ) dσ i dE i (cid:18) E h x (cid:19) , (2.1)where the summation over i is over all partons in the final state. D i → h is the fragmentationfunction for the transition i → h . Although the hadron’s energy E h is used as an examplehere, any observable linear in the hadron’s momentum can be utilized.The kinematics of the collinear fragmentation process can be represented as follows i ( p i ) → h ( p h ) + X ( p i − p h ) , p µh = xp µi , x ∈ [0 , , (2.2)where the momenta of particles have been indicated in brackets and X represents theparticles produced in the fragmentation process which are not explicitly described by thefragmentation function, i.e. all particles in the jet initiated by i other than the observed– 3 –adron h . Essentially, this means that one relates the hadron’s momentum to that of asingle parton, the latter being an infrared-unsafe quantity.As the above discussion indicates, the partonic cross section for producing a parton i is infrared unsafe and therefore contains uncancelled divergences. These are collineardivergences which factorise into lower-order contributions to the cross section and process-independent splitting functions. Because of this general and process-independent structure,it is possible to absorb the uncancelled divergences into the fragmentation functions viacollinear renormalisation [11]: D bare i → h ( x ) = (cid:88) j (cid:0) ˆΓ ij ⊗ D i → h (cid:1) ( x ) , (2.3)where the sum is over all partons. The collinear counterterms ˆΓ ij are functions of x and canbe specified, not uniquely, within perturbation theory. In practice a choice is made aboutthe finite terms contained in these counterterms. Such a choice implies that the IR renor-malized coefficient and fragmentation functions, dσ i and D i → h , are individually schemedependent however their convolution dσ h is not, as one may expect from an observable.As for parton distribution functions, it is standard practice to define the counterterms ˆΓ ij in the MS scheme.The collinear renormalisation eq. (2.3) introduces scale dependence into the renor-malised fragmentation functions, which is described by the (time-like) Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [12–14]: µ F r dD i → h dµ F r ( x, µ F r ) = (cid:88) j (cid:0) P T ij ⊗ D j → h (cid:1) ( x, µ F r ) , (2.4)where P T ij are the time-like splitting functions, known through NNLO [15–17], and µ F r is the fragmentation factorisation scale, or simply the fragmentation scale. Because frag-mentation functions are extracted from experiment at a certain scale, it is necessary torelate fragmentation functions evaluated at two different scales. This is achieved by solv-ing the DGLAP equations eq. (2.4). The initial conditions necessary for fully specifying thesolution are discussed in sec. 4. The solution of the DGLAP equation has the additionalbenefit that any large logarithms of the ratio of two scales are resummed with a logarithmicaccuracy given by the order of the splitting functions used.
Fixed order calculations are typically performed using a subtraction scheme. The pur-pose of a subtraction scheme is to ensure that in any singular limit of the kinematics, thesingularities of physical cross sections are matched by those of the relevant subtractionterms and that in those limits, the corresponding final states are indistinguishable. Theseare the requirements for the numerical integrability of the cross section. If the singularbehaviour of the cross section is not matched by its subtraction terms, then a numeri-cally non-integrable singularity remains. If the singular behaviours of the contributions– 4 –atch, but the kinematics are distinct, then the fully inclusive cross section is numericallyintegrable, but differential and fiducial cross sections may not be. Schematically, a crosssection differential in some observable O can be written as dσdO = (cid:88) n (cid:90) n dσ n ( { y i } n ) δ ( O ( { y i } n ) − O )= (cid:88) n (cid:90) n (cid:18) dσ n δ ( O ( { y i } n ) − O ) − (cid:88) m dσ mn δ ( O m ( { y i } n ) − O )+ (cid:88) m dσ mn δ ( O m ( { y i } n ) − O ) (cid:19) , (3.1)where n denotes the number of final-state particles, dσ n is the fully differential n -particlecross section, { y i } n is the set of n -particle phase space parameters to be integrated over, e.g.the set of momentum components of the particles, dσ mn is a subtraction term, the integral (cid:82) n is over the full n -particle phase space and the dependence of dσ n and dσ mn on { y i } n hasbeen omitted on the second and third lines for brevity. As usual, the intent is to integratethe combination of terms on the second line fully numerically, while the integration overthe singular behaviour of the terms on the third line is performed analytically.Due to conceptual differences, a subtraction scheme has to be modified with respect tothe case without fragmentation in order to perform calculations involving fragmentation.Such modifications have been made to the sector-improved residue subtraction scheme [18–21] and its implementation in the Stripper library, enabling the calculations presentedhere. The additional complications due to fragmentation have been discussed in the pastin the context of NLO subtraction schemes, see e.g. ref. [22], and no further complicationsare introduced beyond NLO. Nonetheless, the required modifications will also be discussedhere for consistency and completeness. All of the necessary changes can be identified byconsidering which additional requirements fragmentation effectively puts on the calculation.Writing the fragmentation equivalent of eq. (3.1), these requirements become apparent: dσdO = (cid:88) n (cid:90) dx (cid:90) n dσ n ( { y i } n ) D ( x ) δ ( O ( { y i } n , x ) − O )= (cid:88) n (cid:90) dx (cid:90) n (cid:18) dσ n ( { y i } n ) D ( x ) δ ( O ( { y i } n , x ) − O ) − (cid:88) m dσ mn ( { y i } n ) D (˜ x m ( { y i } n , x )) δ ( O m ( { y i } n , x ) − O )+ (cid:88) m dσ mn ( { y i } n ) D (˜ x m ( { y i } n , x )) δ ( O m ( { y i } n , x ) − O ) (cid:19) , (3.2)where for simplicity a single fragmentation function contributes. A realistic cross sectionwould simply be a sum over such contributions. The functions ˜ x m will be discussed later.The two differences with respect to eq. (3.1) are multiplication by the fragmentation func-tion and the dependence of the observable on the momentum fraction x . For subtractionterms, the dependence of the observable on the phase space parameters changes as welland this point will be discussed first. – 5 – i p j p h p i +p j p h Figure 1 . Observable kinematics in a collinear limit. When the partons i and j become collinear,the jets they initiate become indistinguishable from a jet initiated by a single particle carrying thesum of their momenta (blue shaded regions). If a single hadron h is identified in the final state(red), then the fraction of the fragmenting particle’s momentum carried by h is smaller for thecombination of i and j than for i , since the momentum p h does not change. In typical calculations without fragmentation, all partons are clustered into jets and allobservables depend only on the kinematics of the partons indirectly through the kinematicsof jets. For a collinear limit, this means the relative magnitude of the momenta of thecollinear partons is irrelevant, as only their sum enters the observable. Because of this, itis sufficient for the kinematics of the subtraction term to correspond to the exact collinearconfiguration, replacing the collinear partons by a single parton carrying the appropriatecombination of their conserved quantities, such as flavour and momentum, as illustratedin fig. 1. If one of the collinear partons fragments, then the magnitude of its momentumdoes enter observables, as it is directly related to the momentum of the final hadron viathe rescaling by the momentum fraction. For the example shown in fig. 1, this implies therequirement x i p i = p h = x ij ( p i + p j ) , (3.3)where x i and x ij are the momentum fractions for the parton i and the combination of i and j , respectively. Similarly, the flavour of the fragmenting parton determines the sizeof the contribution through the fragmentation function. When moving to the subtractionkinematics, it is thus necessary to retain the information on the contribution from thefragmenting parton to the total momentum of the collinear partons, which for the exampleabove is e.g. the ratio between p i and p i + p j , and the flavour of the fragmenting particle.There is an important point to stress here concerning soft limits, as the situation isslightly different. A singular soft limit occurs when the total energy of a flavourless set ofpartons – containing gluons and equal numbers of quarks and anti-quarks of each flavour –becomes small. The standard observation is that a configuration containing a zero-energyflavourless set of partons cannot be distinguished from one where this set is removed fromthe final state. The kinematics of the subtraction term thus corresponds to the exact soft– 6 –onfiguration, removing the zero-energy flavourless set of partons from the final state. Thestatement that zero-energy, flavourless sets of partons can be removed from a final statewithout changing any observable is no longer true if one of those partons fragments, asthe hadron is assumed to always be observable on its own. One could in principle proceedas for the collinear case and construct the subtraction kinematics as usual, keeping theinformation about the flavour and momentum, the latter being zero by definition, of thefragmenting parton in the soft limit. However, this yields contributions where the hadronalways carries zero momentum. Not only is this an unphysical configuration, as the hadronhas a non-zero mass, the factorisation of the cross section into the hard process and afragmentation function only applies if the hard scale of the hadron, e.g. its transversemomentum, is much larger than its mass [23]. Additionally, fragmentation functions aredivergent as the momentum fraction goes to zero, so even if these soft limits of the partoniccross section were regulated, the hadronic cross section would still be divergent. Because ofthis, there are no (integrated) subtraction terms regulating the soft limit of a fragmentingparton.By considering exact singular limits, it has been explained that the kinematics ofsubtraction terms must be modified. The exact dependence of the kinematics on the fullphase space parametrisation is arbitrary, only in singular limits must the kinematics of thecross section and its corresponding subtraction terms match, i.e.: O m ( { y i } n , x ) → O ( { y i } n , x ) , (3.4)where the limit is any limit which is supposed to be regulated by dσ mn and eq. (3.4) shouldhold for all infrared-safe observables, where the momentum of the hadron is consideredan infrared-safe quantity within this framework. Without fragmentation, the analyticintegration performed to obtain the integrated subtraction term relies on the fact that O m does not depend on the parameters integrated over. As explained above, this maynot be the case when one of the partons fragments. An example would be a subtractionterm which regulates both a collinear and a soft singularity, which is a part of e.g. theCatani-Seymour dipole subtraction scheme [22], with one of the partons fragmenting. Inthis case, the energy of the hadron would depend on the energy of the soft parton, soit would depend on the parameter parameterising the soft limit, spoiling the ability toperform this integration fully analytically. An implementation of a subtraction schemecontaining such subtraction terms would therefore require laborious modifications beforegeneral fragmentation computations can be performed.Here a critical simplification exists for the sector-improved residue subtraction schemewith respect to many other subtraction schemes. If there are no subtraction terms whichregulate more than one type of singularity, i.e. every subtraction term is designed to countera singularity occurring as all elements of a single set of phase space parameters simulta-neously approach a singular point, then the kinematics of each subtraction term can bechosen to always match those of the cross section in the subtraction term’s characteristicsingular limit. Because this is a constant with respect to the variables which parameterisethe subtraction term and are integrated over analytically in the integrated subtractionterm, the integrated subtraction is unchanged by the introduction of fragmentation (aside– 7 –rom an overall factor given by the fragmentation function), avoiding the need to redo anyanalytic integration previously performed for the original subtraction scheme. For this rea-son, the sector-improved residue subtraction scheme is particularly suited for the extensionto fragmentation, since it does not contain any subtraction terms which regulate multiplesingularities [20].Aside from changes to the kinematics of the final state, the inclusion of fragmentationalso modifies the size of the contributions of different phase space points via multiplicationby the fragmentation function, as shown in eq. (3.2). There is a certain amount of freedomwhen it comes to the point at which the fragmentation function is evaluated for a subtrac-tion term, written in eq. (3.2) as the functions ˜ x m ( { y i } n , x ). The only strictly necessarycondition is that a subtraction term matches the singular behaviour of the cross section incertain singular limits. This requires that in any singular limit, the fragmentation functionis evaluated at the same point for both the cross section and its corresponding subtractionterms, i.e.: ˜ x m ( { y i } n , x ) → x , (3.5)where the limit is again any limit which is supposed to be regulated by dσ mn . The mostsimple choice ˜ x m ( { y i } n , x ) ≡ x (3.6)is made here. Note that in order to reuse the integrated subtraction terms from the casewithout fragmentation as explained above, ˜ x m ( { y i } n , x ) must fulfill an additional condition:it should not depend on the parameters parameterizing the singular limit regulated by dσ mn .This is trivially satisfied by the choice shown in eq. (3.6).The modifications discussed up until now are sufficient to perform calculations withfragmentation, but often lead to suboptimal numerical convergence. The reason for this isthat while the kinematics of the cross section and one of its subtraction terms match in thesingular limit, they do not in the remainder of the phase space. It is thus possible for bothcontributions to be large with opposite signs, but instead of mostly cancelling each other,the contributions are added to different bins of a calculated histogram. This missed-binningincreases the fluctuations within individual bins, increasing their Monte Carlo uncertaintyfor a given number of events and thus reducing the rate of numerical convergence. Tomitigate this, one can rescale the momentum fraction x for each contribution on an event-by-event basis, such that the value of an observable of choice is always identical for allcontributions for any given event. If this reference observable is now binned in a histogram,then missed-binning cannot occur by definition, potentially vastly improving the numericalconvergence.The final difference with respect to calculations without fragmentation is the intro-duction of collinear renormalisation counterterms for fragmentation functions. These areconceptually identical to those for PDFs and it is well-known how to obtain them in termsof splitting functions. The only difference with respect to the renormalisation of PDFs isthe need to use time-like splitting functions, which differ from the space-like ones startingat NLO [24, 25]. For completeness, in appendix B we present the explicit expressions for– 8 –he collinear counterterms while in appendix A we give in some detail the structure of thecross section for the process pp → t ¯ t + X → B + X . The fragmentation functions used in this paper are based on the perturbative fragmen-tation function approach [4], in which all fragmentation functions for the production ofheavy-flavoured hadrons can be related to a single non-perturbative fragmentation func-tion (NPFF) via convolutions with perturbatively calculable coefficients, called perturba-tive fragmentation functions (PFFs): D i → h ( x ) = (cid:0) D i → q ⊗ D NP h (cid:1) ( x ) , (4.1)where i can be any parton, h is the heavy-flavoured hadron and q is the heavy quark. Theheavy-quark PFFs were originally derived at NLO [4] and have since been computed atNNLO as well [26, 27]. The only ingredient required to compute FFs for the productionof heavy-flavoured hadrons is thus the NPFF. Typically, NPFFs are extracted from e + e − data, however, theoretically motivated ones also exist [28, 29].In the remainder of this work we will be interested in the case where the heavy quarkis the bottom, i.e. q = b , and the heavy-flavored hadron is a b -flavored one, i.e. h = B .At present, no such extraction at NNLO employing the PFF approach is available inthe literature. For this reason, three different sets of FFs were obtained from two differentextractions, each set corresponding to a different compromise. A third extraction, whichfollows an approximation of the PFF approach, was presented in ref. [30], but has not beenused here.The first two sets of FFs are based on the extraction of ref. [31]. The FF of that paperis not based on the PFF approach, instead relying on effective field theory calculations.Nonetheless, a NPFF was extracted at NLO and NNLO, including NNLL and N LL large- x resummation, respectively. Unfortunately, due to the different approach to the computationof FFs, there is no simple relation between the FF of that paper and one computed withinthe PFF approach. A reasonable conversion from one type of FF to the other has to bechosen. Another important point is that the extracted FF corresponds to the non-singlet(NS) combination, i.e. the difference between the bottom and anti-bottom FFs: D NS B ( µ F r , x ) = D b → B ( µ F r , x ) − D ¯ b → B ( µ F r , x ) . (4.2)The set of FFs used most centrally in this paper is labelled “FFKM”. Its initialconditions are obtained by taking the extracted non-singlet function of ref. [31] evaluatedat the initial scale µ F r = µ with µ = m b = 4 .
66 GeV, then calculating the FFs otherthan the bottom-quark FF from the PFFs and the extracted NPFF and, finally, addingthe anti-bottom FF to the non-singlet one to obtain the full bottom FF: D i → B ( µ , x ) = (cid:0) D i → b ⊗ D NP B (cid:1) ( µ , x ) , i (cid:54) = b , (4.3)– 9 –F set NPFF PFF Large- x DGLAPFFKM/FFKM(2) @ NLO NLO NLO NNLL NLLFFKM/FFKM(2) @ NNLO NNLO NNLO N LL NNLLCNO @ NLO NLO NLO NLL NLLCNO @ NNLO NLO NNLO NLL NNLL
Table 1 . The differences between the FFKM/FFKM(2) sets and the CNO set at NLO and NNLOin terms of perturbative and logarithmic orders. NPFF refers to the perturbative order at whichthe extraction of non-perturbative parameters was performed. PFF refers to the perturbativeorder of the PFFs used. The column “large- x ” shows the logarithmic order of the resummationof logarithms of 1 − x , while the column labelled “DGLAP” indicates the logarithmic order of theDGLAP resummation. D b → B ( µ , x ) = D NS B ( µ , x ) + D ¯ b → B ( µ , x ) . (4.4)The FFs at any other scale µ F r > µ are then obtained by evolving these initial conditionsusing the DGLAP evolution library APFEL [32].An alternative construction labelled “FFKM(2)” is to proceed as for the FFKM set,but as a final step the non-singlet contribution at each scale is replaced by the non-singletcontribution at that scale as provided by the authors of ref. [31]. This is not equivalent tothe FFKM set, since the FF of ref. [31] does not satisfy the non-singlet DGLAP evolutionequation.The third and final set of FFs, labelled “CNO”, is obtained by taking the extractionof ref. [33]. This extraction was performed using the PFF approach, but only at NLOincluding NLL large- x resummation. This time µ = m b = 4 .
75 GeV.NLO and NNLO versions of all three sets were constructed. The perturbative andlogarithmic orders of different components of the fragmentation functions are shown intable 1. All FFs are symmetrised with respect to particles and anti-particles. The scaleevolution is always performed using
APFEL , where the value and running of α s are alwayschosen to match those of the PDF set used at the same order. As an alternative toperforming the evolution with APFEL , the
MELA [34] library could have been used instead,as was e.g. done in ref. [35] to perform a detailed study of the evolution of heavy-quarkfragmentation functions. For simplicity, neither
MELA nor the results of ref. [35] have beenused here.In order to be able to estimate uncertainties due to the errors on the extracted FFs,multiple versions of all sets were constructed, corresponding to taking the extracted non-perturbative parameters and independently varying them by one standard deviation. Sincethere is only one parameter for the FFKM and FFKM(2) sets, this leads to three variationseach, while the CNO set involves two parameters, leading to 9 variations. For the CNOset, correlations between the parameters are ignored.All three FFs were found to be within reasonable agreement with each other, suggestingnone of the individual compromises are particularly significant.– 10 –
Applications b -fragmentation in top-quark decay As a first application we consider the process t → B + W + X with the subsequent decay W → (cid:96) + ν in NNLO in QCD. We work with top quark pole mass m t = 172 . µ R = µ F r = m t / / ≤ µ R /µ F r ≤
2. Perturbativecalculations for top decay at any accuracy (LO, NLO or NNLO) are always convolved withFF at NNLO. In all cases the value of the strong coupling α s is taken from the LHAPDF interface [36] as produced by the
NNPDF3.1
NNLO pdf set [37]. Further details about thisprocess and its setup can be found in appendix A as well as in ref. [38].In all observables discussed in this section we implement an energy cutoff of E ( B ) > x region of the FFs. Excluding this region is notconsequential for this work since in our implementation all power corrections ∼ ( m b ) n , n ≥
2, are neglected and our predictions are not valid in the very low x region anyway.As a check on our implementation we have verified that our calculation satisfies themomentum conservation sum rule, see appendix C.2 for details.We study the following observables: the invariant mass of the lepton and the hadron m ( B(cid:96) ) and the energy fraction of the B -hadron to its maximum energy E ( B ) /E ( B ) max ,where E ( B ) max = m t − m W m t . (5.1)The observables are shown in fig. 2. In both cases we show the absolute distributionsat different perturbative orders for the FFKM NNLO fragmentation function. The lowerpanel shows the ratio to the NLO result. The colored bands correspond to 7-point scalevariation. In fig. 3 we show a breakdown of the NNLO scale variation due to µ R and µ F r .Each one of these scales is varied (3-point variation) while the other scale is kept fixed atits central value. Similarly, fig. 4 shows the fragmentation function variation for the defaultFFKM fragmentation function at NNLO. Also shown are the central predictions at NNLObased on the other two FF sets: FFKM(2) and CNO.The invariant mass differential width m ( B(cid:96) ) is of particular interest since it is suitablefor extracting the top quark mass with high-precision [39]. It has previously been studiedwith NLO precision in ref. [40].The normalized energy spectrum is also interesting in top mass determinations since itdirectly exposes the fragmentation function. Therefore it allows one to directly assess thesensitivity of this observable to b -fragmentation and its potential for measuring NPFF’s.This observable has been studied in NLO+NLL QCD in [41–46]. The analytic expressionsof the coefficient functions for both m b = 0 and m b (cid:54) = 0 are known through NLO in QCD. b -fragmentation in top-quark pair-production and decay at the LHC In this section we present our predictions for the following B -hadron distributions in dilep-ton t ¯ t events at the LHC: the invariant mass of the B -hadron and charged lepton m ( B(cid:96) )– 11 – ]0 . . . . . . d Γ t / d m ( B ‘ ) [ pb / G e V ] t -decay m t = 172 . α s : NNPDF31 NNLOScale: µ R = µ Fr = m t / LO - FFKM NNLONLO - FFKM NNLONNLO - FFKM NNLO0 20 40 60 80 100 120 140 160 m ( B‘ ) [GeV] . . . . . . . . . r a t i o t o N L O []0 . . . . . . d Γ t / d ( E ( B ) / E ( B ) m a x ) [ pb ] t -decay m t = 172 . α s : NNPDF31 NNLOScale: µ R = µ Fr = m t / LO - FFKM NNLONLO - FFKM NNLONNLO - FFKM NNLO0 . . . . . . . . . . E ( B ) /E ( B ) max . . . . . . . . r a t i o t o N L O Figure 2 . Absolute differential top decay width as a function of the invariant mass m ( B(cid:96) ) (left)and the energy fraction E ( B ) /E ( B ) max (right). All curves are convoluted with the same FF: FFKMat NNLO. Shown is comparison for different perturbative orders: LO, NLO and NNLO. []0 . . . . . r a t i o t o NN L O t -decay m t = 172 . α s : NNPDF31 NNLOScale: µ R = µ Fr = m t / NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ R variation0 20 40 60 80 100 120 140 160 m ( B‘ ) [GeV] . . . . . r a t i o t o NN L O NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ Fr variation []0 . . . . . . . . . r a t i o t o NN L O t -decay m t = 172 . α s : NNPDF31 NNLOScale: µ R = µ Fr = m t / NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ R variation0 . . . . . . . . . . E ( B ) /E ( B ) max . . . . . r a t i o t o NN L O NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ Fr variation Figure 3 . As in fig. 2 but showing the scale variation of the NNLO prediction: µ R -only vs. totalscale variation (upper plot) and µ F r -only vs. total scale variation (lower plot). m ( B‘ ) [GeV] . . . . . . . . . r a t i o t o NN L O t -decay m t = 172 . α s : NNPDF31 NNLOScale: µ R = µ Fr = m t / NNLO-FFKMNNLONNLO-FFKM(2)NNLONNLO-CNONNLO 0 . . . . . . . . . . E ( B ) /E ( B ) max . . . . . . . . . r a t i o t o NN L O t -decay m t = 172 . α s : NNPDF31 NNLOScale: µ R = µ Fr = m t / NNLO-FFKMNNLONNLO-FFKM(2)NNLONNLO-CNONNLO
Figure 4 . As in fig. 2 but showing the fragmentation function variation of the default FF FFKMat NNLO. Shown also are the central predictions for the other two FF at NNLO: FFKM(2) andCNO. as well as B -hadron’s energy E ( B ). These two distributions are the t ¯ t equivalents of thedistributions discussed in sec. 5.1 in the context of top quark decay. The advantage ofworking with m ( B(cid:96) ) and E ( B ) is that they are defined in the detector frame and are,therefore, directly measurable without the need for reconstructing frames associated withthe top quark. Both m ( B(cid:96) ) and E ( B ) are of prime interest in the context of top quark– 12 –ass determination at the LHC and have been extensively studied in the past in NLOQCD [39, 40, 47, 48].The setup of the present calculation, which is closely related to the one in ref. [38], seealso appendix A, is as follows. We utilize the pdf set NNPDF3.1 . Its order is chosen in sucha way that it matches the order of the perturbative calculation. The value of the strongcoupling constant is obtained from the
LHAPDF library as provided by the
NNPDF3.1 pdfset. The order of the strong coupling constant evolution in the perturbative calculationis matched to the order of the pdf while the order of the coupling in the FF evolution ismatched to the order of the FF.The pdf variation utilizes the so-called reduced pdf set, see ref. [38] for details. Ourpredictions are based on fixed central scales µ R = µ F = µ F r = m t . (5.2)The reasons behind this scale choice are as follows. A fixed scale choice is well-justifiedin the kinematic ranges considered in this work. Furthermore, the use of fixed scales(instead of dynamic scales) can simplify the interpretation of the results especially whenthere are many scales and perturbative orders. The specific value of the central scale, m t / t ¯ t production. One may wonder if a central scale m t and not m t / b -fragmentation in t ¯ t production and decay inorder to make the interpretation of the scale variation of the prediction as transparent aspossible since this way all three scales appearing in this calculation have the same centralvalues.Scale variation is defined through a 15-scale variation, i.e. scaling up and down the common central scale by a factor of 2, subject to the constraints1 / ≤ µ R /µ F ≤ , / ≤ µ R /µ F r ≤ , (5.3)1 / ≤ µ F /µ F r ≤ . We use the G F scheme with the following parameters m W = 80 .
385 GeV , Γ W = 2 . ,m Z = 91 . , Γ Z = 2 . ,G F = 1 . · − GeV − ,α = √ G F π m W (cid:0) − ( m W /m Z ) (cid:1) . (5.4)Defining ξ = ( m W /m t ) , the leading order top-quark width is computed fromΓ (0) t = G F m t π √ − ξ ) (1 + 2 ξ ) = 1 . m t = 172 . . (5.5)– 13 –ur calculations are subject to typical phase space cuts: p T ( B ) ≥
10 GeV , | η ( B ) | ≤ . . (5.6)The predictions for the absolute distributions m ( B(cid:96) ) and E ( B ) through NNLO inQCD are shown in fig. 5. The bands around the three central predictions indicate their15-point scale variation. For both distributions we observe that the reduction of the scaleuncertainty when going both from LO to NLO and from NLO to NNLO is substantial. TheNNLO scale variation is about couple of percent in most bins. Notably, for the scale choice(5.2) the NNLO scale variation is asymmetric, unlike the LO and NLO ones. Because of thisasymmetry it is more useful to quantify the total width of the NNLO scale variation bandwhich never exceeds 10% and, in fact, in most bins is about half that value. This impliesthat the corrections due to missing higher order effects are probably at the one-percentlevel and thus rather small.We also observe that the size of the higher order corrections in both observables is mod-erate and in all cases the higher-order corrections are contained within the correspondinglower order scale variation band. The only exception is the lowest bin of the E ( B ) distri-bution however it is worth keeping in mind that this bin is strongly impacted by the cuts(5.6). The NNLO/NLO K-factor is rather small and tends to be within 5% for most binsin both distributions. It has a non-trivial shape relative to the NLO predictions once oneaccounts for the small size of the NNLO scale uncertainty band.The region of the m ( B(cid:96) ) distribution above about 150 GeV is impacted by correctionsbeyond the narrow width approximation which is utilized in this work (see ref. [38] fordetails). The monotonic increase in the shape of the NNLO/NLO K -factor of the E ( B )distribution suggests that at NNLO the maximum of that distribution is shifted towardshigher values of E ( B ) relative to NLO. Although in this paper we are not able to quantifythis shift with sufficient precision, we note that it may significantly affect any extractionof the top quark mass based on the proposal in refs. [47, 48]. A more precise estimate ofthis effect is possible but it will require a dedicated and more refined calculation which weleave for a future work.With the help of fig. 6 one can assess the origin of the scale variation in these twoobservables at NNLO. To that end we have shown a breakdown of the scale variation dueto one scale at a time (the other two being fixed at their central values) and compared tothe total scale variation eq. (5.3). It immediately becomes apparent that the bulk of thescale variation is due to the renormalization scale µ R . The second largest contribution isdue to the fragmentation scale µ F r while the contribution due to the factorization scalealone is tiny.In fig. 7 we compare at NNLO the three main sources of uncertainty for these twodistributions: scale, pdf and fragmentation uncertainties. The variations shown are for thedefault FFKM fragmentation function. As an alternative measure for the fragmentationuncertainty we show the central predictions based on the two alternative FFs: FFKM(2)and CNO. It is evident from this figure that scale variation is the dominant source ofuncertainty. This is true for all bins of both distributions. The second largest uncertaintyis the one due to fragmentation followed by the pdf uncertainty. The differences between the– 14 – ]0 . . . . . . . . . d σ / d m ( B ‘ ) [ pb / G e V ] LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / LO - FFKM NNLONLO - FFKM NNLONNLO - FFKM NNLO0 20 40 60 80 100 120 140 160 180 m ( B‘ ) [GeV] . . . . . . . r a t i o t o N L O []0 . . . . . . . . . d σ / d E ( B ) [ pb / G e V ] LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / LO - FFKM NNLONLO - FFKM NNLONNLO - FFKM NNLO0 25 50 75 100 125 150 175 200 E ( B ) [GeV] . . . . . . . r a t i o t o N L O Figure 5 . Absolute differential top-quark pair production and decay cross section as a function ofthe invariant mass m ( B(cid:96) ) (left) and the B -hadron energy E ( B ) (right). All curves are convolutedwith the same FF: FFKM at NNLO. Comparisons for LO, NLO and NNLO are shown. []0 . . . . . r a t i o t o NN L O LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ R variation []0 . . . . . r a t i o t o NN L O NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ F variation0 20 40 60 80 100 120 140 160 180 m ( B‘ ) [GeV] . . . . . r a t i o t o NN L O NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ Fr variation []0 . . . . . r a t i o t o NN L O LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ R variation []0 . . . . . r a t i o t o NN L O NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ F variation0 25 50 75 100 125 150 175 200 E ( B ) [GeV] . . . . . r a t i o t o NN L O NNLO - FFKM NNLO - Full scale variationNNLO - FFKM NNLO - µ Fr variation Figure 6 . As in fig. 5 but showing the scale variation of the NNLO prediction: µ R -only vs. total(upper plot), µ F -only vs. total (middle plot) and µ F r -only vs. total scale variation (lower plot). m ( B‘ ) [GeV] . . . . . . . . . r a t i o t o NN L O LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / NNLO-FFKMNNLONNLO-FFKM(2)NNLONNLO-CNONNLO 0 25 50 75 100 125 150 175 200 E ( B ) [GeV] . . . . . . . . . r a t i o t o NN L O LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / NNLO-FFKMNNLONNLO-FFKM(2)NNLONNLO-CNONNLO
Figure 7 . As in fig. 5 but showing the fragmentation and pdf variations of the default FFKM FF.Also shown are the central predictions for the other two FF at NNLO: FFKM(2) and CNO. – 15 –hree fragmentation functions tends to be consistent with the estimate of the fragmentationuncertainty although in some bins that difference is as large as twice the value of thefragmentation uncertainty estimate.In summary, the total uncertainty of the NNLO predictions for the m ( B(cid:96) ) and E ( B )distributions is within 5% for almost all bins and is dominated by the scale uncertainty.While in this first NNLO work on this subject we have considered the 15-point scale varia-tion eq. (5.3) around the central scale eq. (5.2) as the most straightforward generalizationof the usual restricted scale variation in processes involving a single factorization scale, itmay be beneficial to revisit this in the future and try to assess the impact and merits of amore restrictive scale variation and/or different dynamic or fixed scale choices. B -hadron FFs from t ¯ t events The focus of the previous discussions was on predictions for LHC observables given aset of fragmentation functions. Due to the limitations of the existing extractions from e + e − data one may naturally ask the question if LHC data can be used to improve theextraction of non-perturbative FFs. In this section we address this question in the contextof b -fragmentation in t ¯ t events. As it will become clear shortly, this study can easily beextended to other processes like direct b production.In principle, one can use any well-measured LHC B -hadron distribution to fit theNPFF. In order to increase the sensitivity to the NPFF it would be ideal if one usesdistributions that are as closely related to the FF’s as possible. An example for such adistribution is the B energy spectrum in top quark decay discussed in sec. 5.1. The onlydrawback of this distribution is that it requires the reconstruction of the decaying top quarkand, thus, cannot be measured directly. It is therefore preferable to have distributions withsimilar sensitivity to NPFF that are directly defined in the lab frame.In this work we propose one such distribution: the ratio p T ( B ) /p T ( j B ) of the transversemomentum of the identified B -hadron with respect to the transverse momentum of the jetthat contains it. We cluster jets with the anti- k T algorithm [49] with radius R = 0 .
8. Werequire that this jet fulfills p T ( B ) ≤ p T ( j B ) and | η ( j B ) | < .
4, consistent with the cuts ineq. (5.6). Note that both the B -hadron and its fragmentation remnants are included inthis jet-clustering, see the discussion around eq. (2.2).The differential p T ( B ) /p T ( j B ) distribution is shown in fig. 8. The shape and behaviorof this observable at different perturbative orders is fairly similar to the E ( B ) /E ( B ) max distribution in top decay shown in fig. 2. Higher order corrections are largely consistentwith the scale uncertainty bands of the lower perturbative order. The size of scale vari-ation at NNLO is below 5% except for large values of p T ( B ) /p T ( j B ) where it starts toincrease. We have checked that, just like in the case of m ( B(cid:96) ) and the B -hadron energy E ( B ) distributions shown in fig. 6, the scale variation in this observable is driven by therenormalization scale and in much smaller degree, by the fragmentation scale µ F r . Thevariation due to µ F alone is negligible.From fig. 8 one can also conclude that for intermediate and large values of p T ( B ) /p T ( j B )the uncertainty of this observable is driven by the uncertainty in the non-perturbative frag-mentation function. For values p T ( B ) /p T ( j B ) (cid:46) . ]0102030405060 d σ / d ( p T ( B ) / p T ( j B )) [ pb ] LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / LO - FFKM NNLONLO - FFKM NNLONNLO - FFKM NNLO0 . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . . . r a t i o t o N L O . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o NN L O LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / NNLO - FFKM NNLONNLO - FFKM (2) NNLONNLO - CNO NNLO
Figure 8 . Absolute differential cross section as a function of the transverse momentum ratio p T ( B ) /p T ( j B ) in top-quark pair production and decay. Comparison of the FFKM NNLO FFfor different perturbative orders showing scale variation (left) and comparison of FF, pdf andscale uncertainties (right). PDFs are matched to the corresponding perturbative order. The scalevariation bands are based on 15-point scale variation. the scale variation. The pdf uncertainty is negligible throughout the kinematic range.These observations imply that this observable has strong potential for constraining FF atNNLO in QCD at intermediate and large values of x .We next probe the sensitivity of the p T ( B ) /p T ( j B ) distribution to the following pa-rameters: the jet algorithm, the jet size and the B -hadron p T cut. Our aim is to determineoptimal values for these parameters which will facilitate the extraction of the fragmentationfunction.In fig. 9 we show the p T ( B ) /p T ( j B ) distribution for three different jet algorithms: anti- k T , k T [50, 51] and flavour- k T [52]. For ease of the comparison all jet algorithms have thesame jet size R = 0 .
4. For each jet algorithm we show the LO, NLO and NNLO corrections,including their scale variation. The pattern of higher-order corrections is almost identi-cal for the three jet algorithms. The three algorithms produce very similar distributions.This can be seen in the top left plot which shows a comparison of the three jet algorithmsat NNLO. There we see that the anti- k T and k T algorithms lead to almost identical be-havior. The flavour- k T algorithm also produces almost identical distribution for values of p T ( B ) /p T ( j B ) above about 0.6, but starts to deviate from the other two jet algorithmsfor lower values. Still the difference between the flavour- k T and the other two algorithmsis much smaller than the NNLO scale uncertainty. These comparisons indicate that fromthe viewpoint of this observable all three jet algorithms, anti- k T , k T and flavour- k T , aresuitable for the extraction of NPFF in t ¯ t events.Another comment about the use of the anti- k T and k T algorithms in this calculationis in order. It is well known [52] that starting from NNLO, flavorless jet algorithms are notautomatically infrared (IR) safe when applied to flavored problems. To achieve IR safetyof jets in the flavored context, dedicated jet algorithms are needed. One such proposal isthe flavour- k T algorithm of ref. [52]. Related ideas have been discussed in refs. [53, 54].The use of the anti- k T and k T algorithms is justified in the present work because of the– 17 – . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . . r a t i o t oa n t i - k T LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / NNLO-anti- k T -FFKMNNLONNLO- k T -FFKMNNLONNLO-flavour k T -FFKMNNLO 0 . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O anti- k T , R = 0.4 LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / LO-anti- k T -FFKMNNLONLO-anti- k T -FFKMNNLONNLO-anti- k T -FFKMNNLO0 . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O k T , R = 0.4 LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / LO- k T -FFKMNNLONLO- k T -FFKMNNLONNLO- k T -FFKMNNLO 0 . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O flavour k T , R = 0.4 LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / LO-flavour k T -FFKMNNLONLO-flavour k T -FFKMNNLONNLO-flavour k T -FFKMNNLO Figure 9 . As in fig. 8 but comparing different jet-algorithms: anti- k T , k T and flavour- k T . special nature of the observables computed here. Unlike a typical fixed order calculation,in this work we cluster not just partons but the B -hadron and its accompanying remnants.Since by construction all collinear singularities have been regulated at the level of thepartonic cross-section, a jet algorithm is no longer needed to ensure IR finiteness of thecalculation. In this sense our calculation is closer to an experimental setup than to a typicalfixed order partonic jet calculation. Since the fixed-order part of the B -hadron productioncross-section contains terms of the type log n ( m ) we expect that they will also be present inthe corresponding jet calculation. However due to the NNLL DGLAP resummation theyare likely to not play any role.We next consider the effect of the jet size R . In fig. 10 we compare predictions basedon the anti- k T algorithm with jet sized R = 0 . , . , . , .
8. We observe an expectedpattern of higher order corrections: as the jet size decreases, the observable becomes lessinclusive which results in decreased perturbative convergence. This is manifested throughthe increase of scale uncertainty at all orders considered in this calculation as well as larger K -factors. From this we concluded that from the viewpoint of theory, larger jet sizes arebetter for extracting fragmentation functions from the p T ( B ) /p T ( j B ) distribution.Finally, we consider the impact of the low p T ( B ) cut. To that end in fig. 11 we showthe p T ( B ) /p T ( j B ) distribution computed for three different values of this cut: 10 ,
20 and30 GeV. We show the LO, NLO and NNLO distribution for each p T ( B )-cut as well as acomparison of the three cuts at NNLO. In all cases we use same jet algorithm: anti- k T with R = 0 .
4. We observe that the intermediate-to-large p T ( B ) /p T ( j B ) region is not verymuch affected by the value of the low p T ( B ) cut which, in turn, means that the extractedfragmentation function at intermediate or large values of x is not very sensitive to this cut.From the top-left plot in fig. 11 we observe that in this region the NNLO scale variationfor all cut values is approximately the same.On the other hand, the value of the cut has a strong impact on the distribution atlow p T ( B ) /p T ( j B ). As the p T ( B ) cut is lowered, the distribution becomes divergent infixed order perturbation theory. This is consistent with the observed behavior of the– 18 – . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / anti- k T , R = 0.8 LO-FFKMNNLONLO-FFKMNNLONNLO-FFKMNNLO 0 . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / anti- k T , R = 0.6 LO-FFKMNNLONLO-FFKMNNLONNLO-FFKMNNLO0 . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / anti- k T , R = 0.4 LO-FFKMNNLONLO-FFKMNNLONNLO-FFKMNNLO 0 . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / anti- k T , R = 0.2 LO-FFKMNNLONLO-FFKMNNLONNLO-FFKMNNLO
Figure 10 . As in fig. 8 but comparing different jet sizes R = 0 . , . , . . . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . . . r a t i o t o p T ( B ) > G e V LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / NNLO- p T ( B ) > p T ( B ) > p T ( B ) > . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O anti- k T , R = 0.4, p T ( B ) >
10 GeV
LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / LO- p T ( B ) > p T ( B ) > p T ( B ) > . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O anti- k T , R = 0.4, p T ( B ) >
20 GeV
LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / LO- p T ( B ) > p T ( B ) > p T ( B ) > . . . . . . . . . . p T ( B ) /p T ( j B ) . . . . . . . r a t i o t o N L O anti- k T , R = 0.4, p T ( B ) >
30 GeV
LHC 13 TeV PDF: NNPDF31Scale: µ R = µ F = µ Fr = m t / LO- p T ( B ) > p T ( B ) > p T ( B ) > Figure 11 . As in fig. 8 but comparing different values of the p T ( B ) cut: p T ( B ) > , ,
30 GeV. distribution, which for smaller values of the p T ( B ) cut starts to show the typical signsof bad perturbative convergence: larger scale variation bands and increased K -factors.Finally, one should keep in mind that our calculation is performed with a massless b quarkand therefore misses corrections ∼ ( m b ) n for n ≥
2. For this reason it would be incompleteat low values of p T ( B ). For these reasons we conclude that if experimentally viable, alarger p T ( B ) cut would be preferable since it leads to more stable predictions and sinceany missing b mass corrections are automatically rendered negligible or at least significantlyreduced in importance. Heavy flavor production at hadron colliders has traditionally demanded improved theoret-ical precision which matches the large statistics accumulated at colliders like the Tevatronand the LHC. In processes like b and c production, identified b - or c -flavored hadrons are– 19 –opiously produced with transverse momenta much larger than their masses. For suchkinematics the heavy quark mass plays the role of an infrared regulator. In an appropri-ately defined formalism, like the perturbative fragmentation function one we utilize in thepresent work, such mass effects could be consistently neglected.In this work we extend for the first time the PFF formalism at hadron colliders toNNLO QCD. The novelty of the present work is that it develops a general, numeric,fully-flexible computational framework for perturbative cross sections for hadron colliderprocesses with identified hadrons in NNLO QCD. Our work also benefits from the fact thatall process-independent contributions needed for the description of heavy flavor fragmenta-tion in NNLO – like perturbative fragmentation functions, splitting functions and extractedfrom data non-perturbative fragmentation functions – are available in the literature. Ourframework is able to compute fully differential distributions with a single identified heavyhadron plus additional jets and non-strongly interacting particles. As a first applicationwe compute the NNLO QCD corrections to B -hadron production in t ¯ t production withdilepton decays. The predicted realistic differential distributions significantly benefit fromthe inclusion of the NNLO QCD corrections.There are a number of ways the current work can be extended and we plan to pursuethose in the near future. For example, one can compute open B production at high p T . Theframework developed here can be extended in a straightforward way to charm productionas well.One of the bottlenecks in this approach is the availability of high-quality non-perturbativefragmentation functions. These have previously been extracted from e + e − data but theprecision is not on par with current demand. In addition, the existing fragmentation func-tions are not fully compatible with our approach. To correct for this we intend to extract inthe future non-perturbative fragmentation functions from e + e − data within our framework.In this work we have also studied the prospect of using LHC data for extracting B -hadron fragmentation functions. To that end we have proposed, and studied in detail, adistribution which we find to be particularly well suited for this task: the ratio of the p T of the B hadron to the p T of the jet containing it. In the course of this study we have paidparticular attention to the thorny problem of flavored jets in NNLO QCD.Finally, an all-encompassing description of heavy flavor production in NNLO QCDwill require the merging of fixed order calculations at low p T with the high p T descriptionconsidered here. It is perhaps not too hard to envisage such a solution which, for example,builds on the FONLL approach at NLO. NNLO calculations with full mass dependenceare possible as was recently demonstrated in ref. [1]. While such a merging is beyond thescope of the present work it represents a natural future extension of the present work. Acknowledgments
The work of M.C. was supported by the Deutsche Forschungsgemeinschaft under grant396021762 - TRR 257. The work of T.G. was supported by the Deutsche Forschungsge-meinschaft (DFG) under grant 400140256 - GRK 2497: The physics of the heaviest particlesat the Large Hadron Collider. The research of A.M. and R.P. has received funding from– 20 –he European Research Council (ERC) under the European Union’s Horizon 2020 Researchand Innovation Programme (grant agreement no. 683211). A.M. was also supported bythe UK STFC grants ST/L002760/1 and ST/K004883/1.
A Structure of the cross section for t ¯ t production and top-quark decayincluding fragmentation In the narrow-width approximation for the top quark, the differential cross-section for t ¯ t production and decay factorizes into three sub-processes: the top-pair production differen-tial cross section and the differential widths for the top quark and antiquark dσ = dσ tt × d Γ t Γ t × d Γ t Γ t , (A.1)where × denotes the properly accounted for spin correlations between the various factorizedsub-processes. Through NNLO in QCD the three sub-processes can be expanded as follows dσ tt = dσ (0) tt + α s dσ (1) tt + α s dσ (2) tt , (A.2) d Γ t ( t ) = d Γ (0) t ( t ) + α s d Γ (1) t ( t ) + α s d Γ (2) t ( t ) . (A.3)Further details about the structure of the cross-section eq. (A.1) can be found in ref. [38].In the presence of fragmentation, i.e. for the process pp → t ¯ t + X → B + X , thecross section in eq. (A.1) is further split into contributions depending on the origin of thefragmenting parton: dσ = dσ tt × d Γ t → B Γ t × d Γ t Γ t + dσ tt × d Γ t Γ t × d Γ t → B Γ t + dσ ttB × d Γ t Γ t × d Γ t Γ t , (A.4)where the subscript B is introduced to explicitly label the sub-process which initiates thefragmentation into the hadron B .The fragmenting contributions have the following expansions through NNLO in QCD dσ ttB = α s dσ (1) ttB + α s dσ (2) ttB , (A.5) d Γ t ( t ) → B = d Γ (0) t ( t ) → B + α s d Γ (1) t ( t ) → B + α s d Γ (2) t ( t ) → B , (A.6)where: dσ ( n ) ttB = (cid:88) i dσ ( n ) tti ⊗ D i → B for n = 1 , , (A.7) d Γ ( n ) t ( t ) → B = (cid:88) i d Γ ( n ) t ( t ) → i ⊗ D i → B for n = 0 , , . (A.8)The type of parton i in the above equations that can fragment onto the observedhadron B depends on the perturbative order. At LO, for example, no additional partonsare present in dσ tt while the only parton present in the top quark (antiquark) decay is b (¯ b ). At higher orders also the gluon and other quark flavors start to contribute.– 21 – Collinear counterterms for processes involving fragmentation
Here we present the explicit expressions for the collinear counterterms required for thecalculation through NNLO QCD of any hadron collider process with fragmentation. Theresults below follow the conventions of ref. [20] and generalize the corresponding expressionsgiven in that reference to processes involving fragmentation.The NLO collinear renormalisation contribution readsˆ σ C ab → f ...f m [ ... ] ( p , p , k , ..., k m ) = α s π (cid:15) (cid:88) c (cid:90) dz (cid:34) (cid:18) µ R µ F (cid:19) (cid:15) P (0) ca ( z )ˆ σ B cb → f ...f m [ ... ] ( zp , p , k , ..., k m )+ (cid:18) µ R µ F (cid:19) (cid:15) P (0) cb ( z )ˆ σ B ac → f ...f m [ ... ] ( p , zp , k , ..., k m )+ 1 z (cid:88) i (cid:18) µ R µ F r (cid:19) (cid:15) P (0) f i c ( z )ˆ σ B ab → f ...c...f m [ ... ] ( p , p , k , ..., k i /z, ..., k m ) (cid:35) , (B.1)where [...] represents non-fragmenting final state particles, µ F is the PDF factorisationscale, µ F r is the fragmentation function factorisation scale and the relation P (0)T ab = P (0)S ba ≡ P (0) ba has been used. The superscripts S and T in the splitting functions stand for space-likeand time-like, respectively.The NNLO contributions readˆ σ C1 ab → f ...f m [ ... ] = α s π (cid:15) (cid:88) c (cid:90) dz (cid:34) (cid:18) µ R µ F (cid:19) (cid:15) P (0) ca ( z )ˆ σ R cb → f ...f m [ ... ] ( zp , ... )+ (cid:18) µ R µ F (cid:19) (cid:15) P (0) cb ( z )ˆ σ R ac → f ...f m [ ... ] ( p , zp , ... )+ 1 z (cid:88) i (cid:18) µ R µ F r (cid:19) (cid:15) P (0) f i c ( z )ˆ σ R ab → f ...c...f m [ ... ] ( ..., k i /z, ... ) (cid:35) , (B.2)ˆ σ C2 ab → f ...f m [ ... ] = α s π (cid:15) (cid:88) c (cid:90) dz (cid:34) (cid:18) µ R µ F (cid:19) (cid:15) P (0) ca ( z )ˆ σ V cb → f ...f m [ ... ] ( zp , ... )+ (cid:18) µ R µ F (cid:19) (cid:15) P (0) cb ( z )ˆ σ V ac → f ...f m [ ... ] ( p , zp , ... )+ 1 z (cid:88) i (cid:18) µ R µ F r (cid:19) (cid:15) P (0) f i c ( z )ˆ σ V ab → f ...c...f m [ ... ] ( ..., k i /z, ... ) (cid:35) + (cid:16) α s π (cid:17) (cid:15) (cid:88) c (cid:90) dz (cid:34) (cid:18) µ R µ F (cid:19) (cid:15) P (1)S ca ( z )ˆ σ B cb → f ...f m [ ... ] ( zp , ... )+ (cid:18) µ R µ F (cid:19) (cid:15) P (1)S cb ( z )ˆ σ B ac → f ...f m [ ... ] ( p , zp , ... )– 22 – 1 z (cid:88) i (cid:18) µ R µ F r (cid:19) (cid:15) P (1)T cf i ( z )ˆ σ B ab → f ...c...f m [ ... ] ( ..., k i /z, ... ) (cid:35) + (cid:16) α s π (cid:17) β (cid:15) (cid:88) c (cid:90) dz (cid:34) (cid:40)(cid:18) µ R µ F (cid:19) (cid:15) − (cid:18) µ R µ F (cid:19) (cid:15) (cid:41) P (0) ca ( z )ˆ σ B cb → f ...f m [ ... ] ( zp , ... )+ (cid:40)(cid:18) µ R µ F (cid:19) (cid:15) − (cid:18) µ R µ F (cid:19) (cid:15) (cid:41) P (0) cb ( z )ˆ σ B ac → f ...f m [ ... ] ( p , zp , ... )+ 1 z (cid:88) i (cid:40)(cid:18) µ R µ F r (cid:19) (cid:15) − (cid:18) µ R µ F r (cid:19) (cid:15) (cid:41) P (0) f i c ( z )ˆ σ B ab → f ...c...f m [ ... ] ( ..., k i /z, ... ) (cid:35) + (cid:16) α s π (cid:17) (cid:15) (cid:88) cd (cid:90) dz (cid:34) (cid:18) µ R µ F (cid:19) (cid:15) (cid:16) P (0) cd ⊗ P (0) da (cid:17) ( z )ˆ σ B cb → f ...f m [ ... ] ( zp , ... )+ (cid:18) µ R µ F (cid:19) (cid:15) (cid:16) P (0) cd ⊗ P (0) db (cid:17) ( z )ˆ σ B ac → f ...f m [ ... ] ( p , zp , ... )+ 1 z (cid:88) i (cid:18) µ R µ F r (cid:19) (cid:15) (cid:16) P (0) f i d ⊗ P (0) dc (cid:17) ( z )ˆ σ B ab → f ...c...f m [ ... ] ( ..., k i /z, ... ) (cid:35) + (cid:16) α s π (cid:17) (cid:15) (cid:88) cd (cid:90) (cid:90) dzdz (cid:34) (cid:18) µ R µ F (cid:19) (cid:15) P (0) ca ( z ) P (0) db ( z )ˆ σ B cd → f ...f m [ ... ] ( zp , zp , ... )+ 1 z (cid:88) i (cid:18) µ R µ F µ F r (cid:19) (cid:15) P (0) ca ( z ) P (0) f i d ( z )ˆ σ B cb → f ...d...f m [ ... ] ( zp , ..., k i /z, ... )+ 1 z (cid:88) i (cid:18) µ R µ F µ F r (cid:19) (cid:15) P (0) cb ( z ) P (0) f i d ( z )ˆ σ B ac → f ...d...f m [ ... ] ( p , zp , ..., k i /z, ... )+ 1 zz (cid:88) i In the following we detail two checks of our NNLO calculational setup defined in sec. 3. C.1 b -fragmentation in e + e − collisions As an important check of our numerical setup we calculate the coefficient functions in e + e − at NNLO QCD. In fig. 12 we show a typical e + e − observable, the normalized B -hadronenergy, computed at LO, NLO and NNLO QCD. It is compared at NNLO to a calculation ofthe same observable using the exact analytic form of the e + e − coefficient functions [55–58].We check separately the quark and gluon coefficient functions by comparing the B -hadronenergy including only the b and ¯ b contributions (left) or only the gluon one (right). The b + ¯ b and g contributions are separated according to eq. (2.1).– 23 – ]051015202530 d σ / d ( E ( B ) / E m a x ) [ A . U .] e + e − Q = m Z α s : NNPDF31 NNLOScale: µ R = µ Fr = m Z LO - FFKM NNLO b ¯ b onlyNLO - FFKM NNLO b ¯ b onlyNNLO - FFKM NNLO b ¯ b only []0 . . . . . . . . r a t i o t o N L O . . . . . . . . . E ( B ) /E max . . . . . . . S T R I PPE R / A n a l y t i c FFxMC unc.: NNLO - FFKM NNLO b ¯ b onlyNNLO Analytic - FFKM NNLO b ¯ b only [] − − − − − d σ / d ( E ( B ) / E m a x ) [ A . U .] e + e − Q = m Z α s : NNPDF31 NNLOScale: µ R = µ Fr = m Z LO - FFKM NNLO g onlyNLO - FFKM NNLO g onlyNNLO - FFKM NNLO g only[]0 . . . . . . . . . r a t i o t o N L O . . . . . . . . . E ( B ) /E max . . . . . . . . S T R I PPE R / A n a l y t i c FFxMC unc.: NNLO - FFKM NNLO g onlyNNLO Analytic - FFKM NNLO g only Figure 12 . Comparison of predictions for the B -hadron energy spectrum in e + e − collisions basedon our numerical calculation and on the exact analytic e + e − coefficient functions. A single partonicchannel is shown at a time: b + ¯ b (left) and gluon (right). The numerical setup is as follows. In all cases we use the FFKM fragmentation func-tions set at NNLO introduced in sec. 4. The calculations are performed for a fixed centralscale choice µ R = µ F r = m Z . The NNLO comparison between the two setups has beenperformed only for this central scale choice. The value of the strong coupling constant istaken from the LHAPDF interface as supplied with the NNPDF3.1 pdf set. The numericalvalues of all other parameters entering the calculation are given in eq. (5.4).As is evident from fig. 12 there is an excellent agreement between the two calculations,within the MC error of the numeric calculation, and in the full kinematic range considered.This agreement represent a very strong check on the correctness of our numerical setup forboth the quark and gluon coefficient functions. C.2 Sum rules in top decay It is well known [59] that heavy flavor production in e + e − collisions satisfies the followingsum rule: σ = 12 (cid:88) h (cid:90) dx x dσ h dx , (C.1)where h denotes a specific hadron that can be produced in the fragmentation of the heavyflavor and x = 2 E ( h ) /E had , with E had = Q being the energy available for hadronic radi-ation. The fragmentation functions which are implicit in the above equation satisfy thefollowing momentum conservation condition (cid:88) h (cid:90) dz z D i → h ( z ) = 1 . (C.2)– 24 –s an additional check of our computational setup we verify that a sum rule analogousto eq. (C.1) is satisfied in the case of b -production in top quark decay. To this end weconstruct a set of fake fragmentation functions that fulfill eq. (C.2): D b ( z ) = 60 z (1 − z ) ,D ¯ b ( z ) = 105 z (1 − z ) ,D g ( z ) = 30 z (1 − z ) ,D q ( z ) = 168 z (1 − z ) ,D ¯ q ( z ) = 504 z (1 − z ) . (C.3)For the purpose of checking the calculation of the coefficient functions it is sufficient toconsider the case of a single hadron species. The equivalent of eq. (C.1) for the case of topquark decay reads Γ = (cid:90) dx x d Γ dx , (C.4)with x = E ( B ) /E had , where E had = m t − E W is the energy available for hadronic radiationin top quark decay. The maximum value of E ( B ) is given in eq. (5.1).By comparing the RHS of eq. (C.4) with an independent direct calculation of the topquark width we have verified that the pure NNLO correction of order O ( α s ) satisfieseq. (C.4) with numerical precision of about 3%. References [1] S. Catani, S. Devoto, M. Grazzini, S. Kallweit and J. Mazzitelli, [arXiv:2010.11906 [hep-ph]].[2] M. Czakon, D. Heymes and A. Mitov, Phys. Rev. Lett. , no.8, 082003 (2016)[arXiv:1511.00549 [hep-ph]].[3] M. Czakon, D. Heymes and A. Mitov, JHEP , 071 (2017) [arXiv:1606.03350 [hep-ph]].[4] B. Mele and P. Nason, Nucl. Phys. B , 626-644 (1991) [erratum: Nucl. Phys. B ,841-842 (2017)][5] M. Cacciari, M. Greco and P. Nason, JHEP , 007 (1998) [arXiv:hep-ph/9803400 [hep-ph]].[6] M. Cacciari, S. Frixione, N. Houdeau, M. L. Mangano, P. Nason and G. Ridolfi, JHEP ,137 (2012) [arXiv:1205.6344 [hep-ph]].[7] M. Cacciari, M. L. Mangano and P. Nason, Eur. Phys. J. C , no.12, 610 (2015)[arXiv:1507.06197 [hep-ph]].[8] G. Aad et al. [ATLAS], Eur. Phys. J. C , no.7, 330 (2015) [arXiv:1503.05427 [hep-ex]].[9] V. Khachatryan et al. [CMS], JHEP , 123 (2016) [arXiv:1608.03560 [hep-ex]].[10] S. M. Berman, J. D. Bjorken and J. B. Kogut, Phys. Rev. D , 3388 (1971)[11] R. K. Ellis, W. J. Stirling and B. R. Webber, Camb. Monogr. Part. Phys. Nucl. Phys.Cosmol. , 1-435 (1996)[12] G. Altarelli and G. Parisi, Nucl. Phys. B , 298-318 (1977)[13] Y. L. Dokshitzer, Sov. Phys. JETP , 641-653 (1977) – 25 – 14] V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. , 438-450 (1972) IPTI-381-71.[15] A. Mitov, S. Moch and A. Vogt, Phys. Lett. B , 61-67 (2006) [arXiv:hep-ph/0604053[hep-ph]].[16] S. Moch and A. Vogt, Phys. Lett. B , 290-296 (2008) [arXiv:0709.3899 [hep-ph]].[17] A. A. Almasy, S. Moch and A. Vogt, Nucl. Phys. B , 133-152 (2012) [arXiv:1107.2263[hep-ph]].[18] M. Czakon, Phys. Lett. B , 259-268 (2010) [arXiv:1005.0274 [hep-ph]].[19] M. Czakon, Nucl. Phys. B , 250-295 (2011) [arXiv:1101.0642 [hep-ph]].[20] M. Czakon and D. Heymes, Nucl. Phys. B , 152-227 (2014) [arXiv:1408.2500 [hep-ph]].[21] M. Czakon, A. van Hameren, A. Mitov and R. Poncelet, JHEP , 262 (2019)[arXiv:1907.12911 [hep-ph]].[22] S. Catani and M. H. Seymour, Nucl. Phys. B , 291-419 (1997) [erratum: Nucl. Phys. B , 503-504 (1998)] [arXiv:hep-ph/9605323 [hep-ph]].[23] J. C. Collins, D. E. Soper and G. F. Sterman, Adv. Ser. Direct. High Energy Phys. , 1-91(1989) [arXiv:hep-ph/0409313 [hep-ph]].[24] G. Curci, W. Furmanski and R. Petronzio, Nucl. Phys. B , 27-92 (1980)[25] W. Furmanski and R. Petronzio, Phys. Lett. B , 437-442 (1980)[26] K. Melnikov and A. Mitov, Phys. Rev. D , 034027 (2004) [arXiv:hep-ph/0404143 [hep-ph]].[27] A. Mitov, Phys. Rev. D , 054021 (2005) [arXiv:hep-ph/0410205 [hep-ph]].[28] E. Braaten, K. m. Cheung, S. Fleming and T. C. Yuan, Phys. Rev. D , 4819-4829 (1995)[arXiv:hep-ph/9409316 [hep-ph]].[29] U. Aglietti, G. Corcella and G. Ferrera, Nucl. Phys. B , 162-201 (2007)[arXiv:hep-ph/0610035 [hep-ph]].[30] M. Salajegheh, S. M. Moosavi Nejad, H. Khanpour, B. A. Kniehl and M. Soleymaninia,Phys. Rev. D , no.11, 114001 (2019) [arXiv:1904.08718 [hep-ph]].[31] M. Fickinger, S. Fleming, C. Kim and E. Mereghetti, JHEP , 095 (2016)[arXiv:1606.07737 [hep-ph]].[32] V. Bertone, S. Carrazza and J. Rojo, Comput. Phys. Commun. , 1647-1668 (2014)[arXiv:1310.1394 [hep-ph]].[33] M. Cacciari, P. Nason and C. Oleari, JHEP , 006 (2006) [arXiv:hep-ph/0510032 [hep-ph]].[34] V. Bertone, S. Carrazza and E. R. Nocera, JHEP , 046 (2015) [arXiv:1501.00494 [hep-ph]].[35] G. Ridolfi, M. Ubiali and M. Zaro, JHEP , 196 (2020) [arXiv:1911.01975 [hep-ph]].[36] A. Buckley, J. Ferrando, S. Lloyd, K. Nordstr¨om, B. Page, M. R¨ufenacht, M. Sch¨onherr andG. Watt, Eur. Phys. J. C , 132 (2015) [arXiv:1412.7420 [hep-ph]].[37] R. D. Ball et al. [NNPDF], Eur. Phys. J. C , no.10, 663 (2017) [arXiv:1706.00428 [hep-ph]].[38] M. Czakon, A. Mitov and R. Poncelet, [arXiv:2008.11133 [hep-ph]].[39] A. Kharchilava, Phys. Lett. B , 73-78 (2000) [arXiv:hep-ph/9912320 [hep-ph]].[40] S. Biswas, K. Melnikov and M. Schulze, JHEP , 048 (2010) [arXiv:1006.0910 [hep-ph]]. – 26 – 41] G. Corcella and A. D. Mitov, Nucl. Phys. B , 247-270 (2002) [arXiv:hep-ph/0110319[hep-ph]].[42] M. Cacciari, G. Corcella and A. D. Mitov, JHEP , 015 (2002) [arXiv:hep-ph/0209204[hep-ph]].[43] G. Corcella and V. Drollinger, Nucl. Phys. B , 82-102 (2005) [arXiv:hep-ph/0508013[hep-ph]].[44] B. A. Kniehl, G. Kramer and S. M. Moosavi Nejad, Nucl. Phys. B , 720-736 (2012)[arXiv:1205.2528 [hep-ph]].[45] S. M. Moosavi Nejad, Phys. Rev. D , no.9, 094011 (2013) [arXiv:1310.5686 [hep-ph]].[46] S. M. Moosavi Nejad and M. Balali, Eur. Phys. J. C , no.3, 173 (2016) [arXiv:1602.05322[hep-ph]].[47] K. Agashe, R. Franceschini and D. Kim, Phys. Rev. D , no.5, 057701 (2013)[arXiv:1209.0772 [hep-ph]].[48] K. Agashe, R. Franceschini, D. Kim and M. Schulze, Eur. Phys. J. C , no.11, 636 (2016)[arXiv:1603.03445 [hep-ph]].[49] M. Cacciari, G. P. Salam and G. Soyez, JHEP , 063 (2008) [arXiv:0802.1189 [hep-ph]].[50] S. Catani, Y. L. Dokshitzer, M. H. Seymour and B. R. Webber, Nucl. Phys. B , 187-224(1993)[51] S. D. Ellis and D. E. Soper, Phys. Rev. D , 3160-3166 (1993) [arXiv:hep-ph/9305266[hep-ph]].[52] A. Banfi, G. P. Salam and G. Zanderighi, Eur. Phys. J. C , 113-124 (2006)[arXiv:hep-ph/0601139 [hep-ph]].[53] A. Buckley and C. Pollard, Eur. Phys. J. C , no.2, 71 (2016) [arXiv:1507.00508 [hep-ph]].[54] L. Dai, C. Kim and A. K. Leibovich, JHEP , 109 (2018) [arXiv:1805.06014 [hep-ph]].[55] P. J. Rijken and W. L. van Neerven, Phys. Lett. B , 422-428 (1996)[arXiv:hep-ph/9604436 [hep-ph]].[56] P. J. Rijken and W. L. van Neerven, Phys. Lett. B , 207-215 (1997)[arXiv:hep-ph/9609379 [hep-ph]].[57] P. J. Rijken and W. L. van Neerven, Nucl. Phys. B , 233-282 (1997)[arXiv:hep-ph/9609377 [hep-ph]].[58] A. Mitov and S. O. Moch, Nucl. Phys. B , 18-52 (2006) [arXiv:hep-ph/0604160 [hep-ph]].[59] P. Nason and B. R. Webber, Nucl. Phys. B , 473-517 (1994) [erratum: Nucl. Phys. B , 755 (1996)], 755 (1996)]