Gluon fragmentation into $B_c^{(*)}$ in NRQCD factorization

The universal fragmentation functions of gluon into the flavored quarkonia $B_c$ and (polarized) $B_c^*$ are computed within NRQCD factorization framework, at the lowest order in velocity expansion and strong coupling constant. It is mandatory to invoke the DGLAP renormalization program to render the NRQCD short-distance coefficients UV finite in a point-wise manner. The calculation is facilitated with the sector decomposition method, with the final results presented with high numerical accuracy. This knowledge is useful to enrich our understanding toward the large-$p_T$ behavior of $B_c^{(*)}$ production at LHC experiment.

Like parton distribution functions (PDFs), fragmentation functions (FFs) are processindependent functions that encapsulate nonperturbative hadronization effect, and play central role in QCD phenomenology of collider physics. In accordance with the celebrated QCD factorization theorem [1], the inclusive production rate of an identified hadron H at large transverse momentum in high energy collision experiment is dictated by the fragmentation mechanism: (1) where A, B represent two colliding particles, dσ denotes the inclusive rate for producing the parton i, D i→H (z) is the fragmentation function for the parton i into H, which characterizes the probability for i to transition into any final state containing the hadron H specifying the fractional light-cone momentum z with respect to the parent parton i. The sum in (1) is extended over all species of partons, i = q,q, g. Similar to the PDFs, the scale dependence of fragmentation functions is also governed by the celebrated Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation. Taking the gluon fragmentation into the hadron H as example, the DGLAP equation reads where µ is interpreted as the renormalization scale, and P ig (ξ) are the corresponding splitting kernels. Once this FF is deduced at some initial scale µ 0 by any means, one then deduce its form at any other scale µ by solving the evolution equation (2). The fragmentation functions for light hadrons such as π, K, p, . . . , are hopelessly nonperturbative objects, which can only be extracted from experimental data. On the contrary, it was realized in the mid-90s that FFs for heavy quarkonia need not be genuinely nonpertubative entities, which nevertheless can be largely understood in perturbative QCD by the virtue of asymptotic freedom [2,3]. This philosophy is systematically embodied in the modern nonrelativistic QCD (NRQCD) factorization framework [4], that is, the quarkonium FFs can be expressed as the sum of products of short-distance coefficients (SDCs) and longdistance yet universal NRQCD matrix elements, with the series organized by the velocity expansion [2,3]. Since the nonperturbative NRQCD matrix elements are merely numbers rather than functions, the profiles of the quarkonia FFs are largely determined by perturbation theory, therefore NRQCD factorization approach is endowed with strong predictive power.
During the past three decades, a number of fragmentation functions for S-wave/P -wave charmonia/bottomonia have been investigated in the context of NRQCD factorization approach, typically at lowest order in α s (for a very incomplete list, see [2,3,. With the advance of higher-order calculational technique, the SDCs associated with some S-wave quarkonium FFs have recently been calculated through the next-to-leading order (NLO) in α s [14,19,[26][27][28][29].
Unlike J/ψ and Υ, the B c meson is the unique heavy quarkonium which are composed of two different heavy flavors: the b andc quarks. It is interesting to understand the production mechanism of this special heavy meson in hadron collision environment. The LO fragmentation functions for b/c → B in small x region, it is also of great phenomenological incentive to study the g → B ( * ) c fragmentation functions to predict their production rates at LHC. To produce a B c meson in this case, one has to first create a pair of cc and another pair of bb, therefore the perturbative order of this LO fragmentation process is already comparable with that of the NLO QCD correction to b/c → B ( * ) c . Thus, the computation of the g → B ( * ) c fragmentation functions at LO is already rather challenging technically, which have never been considered before. The aim of this work is to fill this gap, by computing the fragmentation functions for g → B ( * ) c at lowest order in velocity expansion and α s , by invoking the sector decomposition technique widely used in the area of multi-loop computation.
In literature there are several different strategies to extract the quarkonium FFs. Among them, the most systematic approach is to start from the gauge-invariant operator definition for the fragmentation functions pioneered by Collins and Soper long ago [30] (Note that this definition was first used by Ma to compute the quarkonium FFs in NRQCD [6]). One great virtue of this operator-based approach is to render renormalization program transparent. According to the operator definition given in [30], the g-to-B ( * ) c fragmentation function is expressed as (see also [15,17]): D = 4 − 2ε signifies the space-time dimensions, and µ is the renormalization scale associated with this nonlocal composite operator. G µν is the matrix-valued gluon field-strength tensor in the adjoint representation of SU(N c ). Here it is convenient to adopt the light-cone coordinates. Any four-vector A µ = (A 0 , A 1 , A 2 , A 3 ) can be recast in the light-cone format . We assume to work in a frame where the B c meson is moving along the z direction with the +-momentum P + , Moreover, k + = P + /z denotes the +-component momentum of injected by the gluon field strength operator. The symbol X indicates collectively those unobserved light hadrons accompanying the B c .
The eikonal factor Φ(0, x − , 0 ⊥ ) in (3) is a path-ordered exponential of the gluon field, inserted to ensure the gauge invariance of the FF: where A µ designates the matrix-valued gluon field in the adjoint representation, g s is the QCD coupling constant, and P denotes the path-ordering, The key observation underlying the NRQCD factorization that, prior to forming a B c state via soft interaction, the b andc quarks have to be first created with slow relative momentum, also within small distance ∼ 1/m c,b , which necessarily involves hard momentum transfer, thus can be studied in perturbation theory owing o asymptotic freedom of QCD. At the lowest order in velocity expansion, NRQCD factorization [4] allows one to refactorize the g → B ( * ) c FFs into the product of the short-distance coefficients and the long-distance NRQCD matrix elements: is the desired SDCs for H = B c , B * c , and the corresponding NRQCD production operators are defined by where ψ b and χ † c are the NRQCD field operators that annihilate a b quark and ac quark, respectively.
Using the vacuum saturation approximation, these vacuum NRQCD matrix elements can be appxomated by the radial wave function at the origin for B ( * ) c in phenomenological potential model: We can proceed to calculate the SDCs d H 1 (z) by the standard perturbative matching technique, i.e., by replacing the physical H state in (5) 1 ). Computing the left side of (5) using perturbative QCD, combined with the following vacuum matrix elements in perturbative NRQCD: we can directly solve the SDC d H 1 (z), order by order in α s . The fragmentation function defined in (3) is manifestly gauge-invariant. In practical calculation, we specialize to Feynman gauge for simplicity. We use a private Mathematica code to automatically generate the Feynman diagrams and the associated cut amplitudes that correspond to the perturbative fragmentation function defined in (3). The Feynman rules for the eikonal propagator and vertex [30], as well as those for conventional QCD propagators and vertices are also implemented with the aid of the package Qgraf [31]. There are seven diagrams on each side of the cut, so in total 49 cut diagrams at lowest order in α s . A simplifying feature is that the gluon propagator cannot be attached to the eikonal line at this perturbative order. For concreteness, in Fig. 1 we exhibit a typical LO cut diagram associated with the perturbative FF D g→bc .
The cut-amplitude structure of fragmentation function stemps from the insertion of the asymptotic out states in (3). Consequently, the corresponding cut-line phase space integration measure reads [15,17]  where k i (i = 1, 2) stands for the momentum of the i-th on-shell quark line (b and c) that pass through the cut. The integration over k + i can be transformed into a parametric integration in a finite interval, but the integration over the transverse momentum k i,⊥ are completely unbounded, i.e., from −∞ to +∞. This feature may persuade us that integration over k i,⊥ could be regarded as loop integration in D − 2-dimensional spacetime, with D = 4 − 2ε. Throughout this work, we adopt the dimensional regularization to regularize occurring UV divergences.
To project the bc pair onto the intended color/spin/orbital/color states, it is convenient to employ the covariant projector technique to expedite the calculation. At the lowest order in velocity expansion, it is legimate to partition the quark momenta inside the fictitious B ( * ) c state commensurate to their mass ratio: pc = rP , p b =rP , with r ≡ m c /(m b + m c ) and r ≡ m b /(m b + m c ). One then make the following substitutions in the quark amplitude to project out the desired contributions from the bc( 1 S (1) 0 ) and bc( 3 S With the aid of the above projectors (10), we employ the packages FeynCalc/FormLink [33,34] to conduct the Dirac/color trace operation. We also use the package Apart [35] to simplify the amplitude by the method of partial fraction, to make the integrand simpler.
It is appealing to use some modern multi-loop technique such as reverse unitarity method and integration by parts (IBP) to deal with phase space integration in (9). The gluon fragmentation into cc( 3 S (1) 1 ) and cc( 1 P (1) 1 ) have been analytically calculated in this manner [23,24,28]. Nevertheless, it turns out to be a quite subtle issue to use these techniques to deal with our case. In this work, similar to [27,29], we choose to apply the sector decomposition method [36,37] to evaluate the phase space integration as given in (9). Consequently, we will present our final results in an entirely numerical manner. In our opinion, the approach adopted in this work appears to be more amenable to automated calculation, and yield more accurate numerical predictions than the complicated subtraction approach first developted in [19] (see also [26,38,39]).
We first combine all the propagators in a cut diagram using Feynman parametrization, then accomplish two-loop integration over k 1,2 ⊥ in D−2-dimensional spacetime. We are then left with multi-fold integrals over Feynman parameters, which can be numerically calculated by the package FIESTA [40] which is based on the sector decomposition algorithm [36,37]. This method is typically useful with many finite multi-variable parametric integrals as output, with various UV poles explicitly isolated.
Upon summing all 49 cut diagrams, we find the total cut amplitude still contains a single uncanceled pole, whose coefficient varies with the momentum fraction z. This pole is clearly of UV origin, from the large k 1,2 ⊥ integration region in (9). This is a clear sign that the bare fragmentation function requires an additional operator renormalization, following the DGLAP paradigm [20,30]: where P qg (y) represents the splitting kernel for g → q: with T F = 1/2. The heavy-quark-to-B ( * ) c fragmentation functions at LO in α s are known long ago [5]. Recasting the results in the NRQCD factorization language, we have where the second equation implies the polarization-summed B * c , and the last one implies the transversely-polarized B * c (denoted by B * T c henceforth). For thec fragmentation into B ( * ) c , one simply makes the exchange r ↔r in above expressions.
The corresponding short-distance coefficients in (13) at LO in α s read [25]: For our purpose, we need compute these LO quark fragmentation functions to O(ǫ) in (11). Note in (11) the UV pole is subtracted in accordance with the MS procedure.
In Fig. 3, we also plot the various fragmentation functions D g→H (z) for H = B c , B * c and B * T c as function of z. We observe that these functions are regular near z = 0. In summary, in this work we have computed the g-to-B ( * ) c fragmentation functions within NRQCD factorization framework, at the lowest order in velocity expansion and α s . We start from the Collins and Soper's gauge-invariant operator definition, which makes the renormalization program transparent. We have employed an automated approach based on sector decomposition strategy to conduct the phase space integral in dimensional regularization. It turns out that, within tolerable amount of time, this method yields better numerical accuracy than the conventional NLO subtraction method. After implementing DGLAP renormalization procedure, we obtain the UV finite NRQCD short-distance coefficient functions in a pointwise manner. We hope that our results are helpful to strengthen our understanding of B ( * ) c production at large-p T at LHC experiment.
Note added. After this work was finished and while we were preparing the manuscript, a preprint has appeared very recently [42], which also computes the gluon-to-B ( * ) c fragmentation functions in NRQCD factorization approach, albeit using the conventional subtraction method.