D and B-meson production using kt-factorization calculations in a variable-flavor-number scheme

Within the framework of $k_t$-factorization, we compute the differential cross section for the production of $B$ and $D$ mesons, using a general-mass variable-flavor-number scheme. Our calculations include all relevant $2\to 2$ processes. We explain how to include the $2\to 1$ process in our calculations, but argue this is not (numerically) relevant at moderate transverse momentum due to its cancellation with the subtraction term. We apply this formalism to $pp$ collisions and compare our results with ALICE and LHCb data at central and forward rapidity.


Introduction
Heavy flavors, which play a particular role in perturbative Quantum Chro-moDynamics (pQCD), have been extensively studied. In the beginning, the observable was the total cross section, plotted as a function of the centerof-mass energy √ s. Later, with higher energies and statistics available, the differential cross section has been measured, for instance, at the Tevatron and Large Hadron Collider (LHC). We note y as the rapidity of the detected heavy flavor and p t as its transverse momentum. The theoretical description of this observable, requiring a more careful treatment compared to the total cross section, has known some troubles [1]. In the framework of collinear factorization, it is known that general-mass variable-flavor-number schemes (GM-VFNS) are more efficient than fixed-flavor-number schemes (FFNS) to describe the differential cross section at large p t . By more efficient, we mean that, at a given order, the former gives better results than the latter, which is particularly true at leading order (LO). This better efficiency is partially explained by the fact that the GM-VFNS resums to all orders some large logarithms ln p 2 t /m 2 Q , thanks to the heavy-quark distribution function.
Heavy-quark production has also been addressed by the first works on k t factorization [2][3][4][5]. During the last two decades, heavy-flavor data have been compared to k t factorization predictions by several groups, see, for instance, [6][7][8][9][10][11]. However, as pointed out in [12,13], some of the available calculations are not performed consistently with respect to the choice of the scheme. One reason for this is simply the lack of appropriate unintegrated parton distributions (uPDFs). Another reason is the complications coming from the use of off-shell matrix elements and the difficulty to compute the cross section at higher orders. The main goal of the present work is to provide for the first time k t -factorization calculations for D and B-meson production, using 2 → 2 processes evaluated in a GM-VFNS. The 2 → 1 process with an initial off-shell charm has already been addressed in [14].
In Sec. 2, we have a general discussion on the schemes used in calculations. Then, we present our GM-VFNS uPDFs in Sec. 3 and the event generator KaTie used for the evaluation of the off-shell cross section in Sec. 4. In Sec. 5, we compare our result to ALICE and LHCb data. Our calculations include all relevant 2 → 2 processes, which in k t -factorization are next-toleading order (NLO) contributions. In Sec. 6 we perform a comparison of our results for D and B mesons with other k t -factorization calculations. In Sec. 7, we discuss the implementation of the LO, 2 → 1, process and argue that it can be ignored at moderate transverse momentum, defined by α s ln p t /m Q ∼ O(1), due to its cancellation with the subtraction term. Finally, we give our conclusion in Sec. 8.

Unintegrated PDFs: discussion on the scheme
In collinear factorization, the cross section for heavy-quark production in hadron-hadron collisions is given schematically by where f k/h is the collinear parton distribution for the parton k in the hadron h, andσ is the partonic cross section. Collinear PDFs are extracted by comparison of Eq. (1) with data. It is then clear that the output f k/h depends on the inputσ. A larger partonic cross section requires smaller PDFs, and the choice made forσ defines the scheme. Once the PDFs have been extracted, they can be used to make predictions. However, the new partonic cross section should be computed with the scheme used in the extraction of the PDFs.
The situation is similar in k t factorization. The cross section reads , depend on x, the fraction of the hadron longitudinal momentum carried by the parton, k t , the initial parton transverse momentum, and µ, the factorization scale.σ is the off-shell cross section. The uPDFs are generally not extracted from data 1 but built from the collinear PDFs by inverting the relation Note that different versions of this relation can be found in the literature. It is clear that the scheme of the uPDFs built from Eq. (3) should be identified to the scheme of the collinear PDFs appearing in this equation. Everything we said on the scheme is also true for the order of calculation. As a consequence, the different sets of uPDFs available on the TMDlib [16,17] cannot be compared by simply using them with the same cross section, as they have been obtained in different schemes and at different orders. For our present study, we use uPDFs and off-shell cross sections obtained in a GM-VFNS at order O(α 2 s ).
The numerical consequence of using a cross section at an order/scheme different from that of uPDFs depends on the case. Mixing the VFNS and FFNS could lead to the wrong estimation of the cross section for charm production by a factor of 4 [12,13]. Indeed in the VFNS, the cross section obtained at order O(α 2 s ) reads where FEP stands for flavor excitation processes, e.g., cg → cg, and FCP for flavor creation processes, e.g., gg → cc. 2 On the opposite, using the FFNS, we have We changed the subscript to k and l to point out that the sums in Eqs. (4) and (5) are not identical. We know thatσ (2) (FEP + FCP) ∼ 4σ (2) (FCP), due to the large contribution of the flavor excitation process cg → cg. Then, we deduce that Consequently, using the VFNS PDFs withσ(FCP) implies an underestimation of the cross section by a factor of 4. The (numerical) situation could improve at higher orders ifσ(FCP) andσ(FEP + FCP) are numerically closer, but mixing different schemes is always inconsistent and dangerous.
Strictly speaking, calculations performed in [12] are not completely consistent since the NLO PB uPDFs [15], being, in fact, next-to-next-to-leading order (NNLO) in k t -factorization, have been combined with a NLO cross section (O(α 2 s )). However, the numerical deviation from the calculations performed in the present paper is found to be small.
The main goal of the present work is to provide consistent GM-VFNS k t -factorization calculations. Consequently, we cannot use the PB uPDFs, and we build our own NLO (in k t -factorization) VFNS uPDFs in the next section.

The Watt-Martin-Ryskin unintegrated PDFs
Accordingly to Eq. (3), VFNS uPDFs of order O(α 2 s ) can be built from collinear PDF extracted in this scheme and at this order. In the Watt-Martin-Ryskin (WMR) approach [18], the uPDF for a parton a reads where P aa (z) are the usual unregulated splitting functions, except for P gg given by and T a (µ, k t ) is the Sudakov factor defined by A possible choice for ∆ is the strong ordering cutoff with µ ∼ p 2 t + m 2 the factorization scale. With this cutoff, However, the strong ordering cutoff does not stem from fundamental considerations, and the obtained uPDFs do not fulfill relation (3) with good precision. This is because ∆ is too large and cuts half of the phase space when k t = 0.5µ. Consequently, we keep the condition (11) and use ∆/d(x), with d(x) > 1 chosen such that the relation (3) is satisfied to a good approximation.
We built our uPDFs using the WMR formalism and the leading order CT14 PDFs [19]. These collinear PDFs have been extracted using the ACOT χ scheme [20][21][22]. In Fig. 1, we compare our unintegrated charm and gluon distributions with the PB and angular-ordering WMR results, obtained from the TMDlib. The main differences are encountered at small k t and at k 2 t > µ 2 . In particular, we observe the typical slow falloff of the AO WMR distribution, leading to an overestimation of the heavy-quark cross section [12]. In Figs

The KaTie event generator
The calculations in the following have been performed with the help of the parton-level event generator KaTie [23]. It can generate events for which the partonic initial-state momenta are spacelike and have nonvanishing transverse components. The necessary uPDFs can be provided by TMDlib or as grid files in text format. In the latter case, KaTie takes care of the interpolation. KaTie operates at tree level and can deal with any processes within the standard model. It generates events employing importance sampling and numerically evaluates the matrix element using helicity amplitudes. The events are stored in event files which can be chosen to be in the LHEF format [24].  It is, however, also possible to make histograms directly using the provided tools.
The initial-state and final-state momenta are generated satisfying exact kinematics. Within k t factorization, this means that the initial-state momenta k µ i = x i P µ i + k µ i t are spacelike, with a longitudinal component along the lightlike momentum P µ i of either colliding hadron, plus transverse components k µ i t . The matrix elements are constructed as described in [25,26], as summed squares of helicity amplitudes. The essence of the method is that a spacelike external parton is represented as a pair of auxiliary lightlike partons satisfying eikonal Feynman rules, leading to manifestly gauge invariant amplitudes with exact kinematics.
In this paper, we consider parton-level processes that, within collinear factorization, would involve massive initial-state quarks with timelike momenta. The same construction of the amplitudes with spacelike initial-state momenta leads also for these processes to manifestly gauge invariant matrix elements. The only restriction is that the transverse momentum of the initial-state parton must not be much smaller than the mass. While gauge invariance is guaranteed, the correct on-shell limit cannot be reached for timelike momenta by naively taking very small transverse momentum. Note that the kinematics of the final-state massive quark is always exact, with its momentum timelike p 2 = m 2 Q . As mentioned earlier, KaTie operates at tree level. The tree-level matrix elements contain singularities as functions of the external momenta when these become soft or collinear with each other. In particular, the matrix elements behave singularly if a final-state gluon momentum becomes collinear with the momentum of an initial-state hadron. The latter also happens in k t factorization. We will consider partonic processes with a final-state gluon, which suffer from such a collinear singularity. In usual tree-level calculations, these singularities are avoided by phase space cuts that define the jet observables, for example, the demand of a minimal transverse momentum. We also use a small minimum on the transverse momentum to avoid the singularity. We see that the cross section depends only very mildly on this phase space cut if we vary the minimum between 0.5 and 2.0 GeV, the latter being of the order of the mass of the quarks we are considering.

D-meson production
Working with a GM-VFN scheme, the following processes are included where Q represents the heavy quark and q a light quark. The first and second lines correspond to flavor-creation and flavor-excitation processes, respectively. For consistency, we choose the charm mass equal to the one used for the CT14 PDFs, e.g., m c = 1.3 GeV. For the fragmentation of a charm quark into a D meson, we use the Peterson model of fragmentation function [27] with c = 0.05, and the fragmentation fractions given in table 1. The three first values are those used by FONLL [28].
The result for the pt distribution of D 0 mesons is compared to ALICE data [29] in Fig. 4. The error band corresponds to the factorization scale uncertainty, evaluated as usual by the variation of a factor of √ 2 above and below the central value, chosen to be where m t,i = p 2 t,i + m 2 c . The subscripts 1 and 2 label outgoing partons. One of these partons is a charm, since processes such as gg → gg with a final gluon fragmenting into a D meson have not been considered 3 . Note that we use m t,i even if p t,i corresponds to a gluon or light quark transverse momentum. In collinear factorization, the usual choice for the factorization scale is with m t the charm transverse mass. In the limit where the transverse momentum of the initial partons goes to zero, k 1t , k 2t → 0, Eqs. (14) and (15) coincide. In Fig. 5, we compare our calculations for the p t distribution of Figure 5: Same as Fig. 4 for D + mesons.
D + mesons with ALICE data. We observe a good agreement on the full p t range, and the central value alone provides a good description, except for the first bin of Fig. 5. A Comparison with figure 5 of Ref. [29] shows that the underestimation of D-meson data at small transverse momentum by theoretical calculations is usual. We will see in the next section that this is not the case for B mesons. Based on [29], where theoretical calculations are compared with the measured cross section, we observe that our work represents a significant improvement of k t -factorization calculations. We believe this is directly related to the consistent use of a GM-VFNS. To reach a similar result in an FFNS, it is probably necessary to include higher orders. We compare our results to other k t -factorization calculations in more detail in section 6. Note also that our central values are slightly better in comparison with FONLL calculations. Similar results obtained at 7 TeV for D * + and D + s are presented at the end of this paper in Figs. 17 and 18. We turn now our attention to ALICE measurement at 5 TeV [32]. The D 0 transverse momentum distribution is presented in Fig. 6. Within uncertainties, the agreement with our calcula- Figure 6: D 0 production at 5 TeV compared with ALICE data [32].
tions is excellent. We can also explore the production of D mesons at larger rapidities and energies. In Fig. 7, we compare our calculations with LHCb data [33] at 13 TeV in the rapidity range 2 < y < 2.5. Using the same set of parameters, we obtained a description of experimental data of similar quality compared to the central rapidity case. What has changed is the relative contributions of flavor excitation and creation processes. Indeed, in Fig. 8,  Figure 7: D 0 production at 13 TeV compared with LHCb data [33] in the rapidity range 2 < y < 2.5.
we observe that at forward rapidity the two contributions are closer (green and red curves). The interplay between flavor-excitation and flavor-creation processes being absent in an FFNS at order O(α 2 s ) is probably another good reason for using a GM-VFNS.
Finally, we quickly discuss other theoretical uncertainties related to our calculations. In Fig. 9, we show the result of varying the charm mass from 1.3 to 1.5 GeV. We observe only a small effect at small transverse momentum. Similarly, varying the p t cuts described in section 4 by a factor of 2 has only a small impact. The uncertainty related to the choice for the fragmentation function has been studied in [9]. Note that we use the same fragmentation function with the same value of the parameter c . We did not evaluate the uncertainty related to the choice of the hard scale discussed in [34], as it is not conventional.

B meson production
Changing Q by b in (12) and (13) gives the complete list of the processes considered. In particular, we did not take into account the cb → cb process.
In agreement with the CT14 PDFs, the bottom mass is set to m b = 4.75 GeV. We use the Peterson fragmentation function with b = 0.01 and the factorization scale Eq. (14) with m t,i = p 2 t,i + m 2 b . In agreement with the discussion in [35], we choose the fragmentation fraction The LHCb Collaboration measured the B ± double-differential cross sections at 7 and 13 TeV, in the rapidity range 2 < y < 4.5 [36]. In Fig. 10, we compare our results to LHCb data at 7 TeV and 2 < y < 2.5. We Figure 10: (B + + B − ) production at 7 TeV compared with LHCb data [36].
observe a global good agreement with our central predictions, with a slight overestimation at p t ∼ 5 GeV. All experimental data lay within theoretical uncertainties, estimated by varying the factorization scale. The uncertainty band is broader than for D mesons, which is due to the dependence of the bottom uPDF with µ, and to the fact that the main contribution is given by the flavor-excitation process gb → gb. Changing the Peterson for the Kartvelishvili et al. fragmentation function [37] with α = 7 gives a similar result. Our calculations at 13 TeV, presented in Fig. 11, show a similar agreement with LHCb data. Figure 11: (B + +B − ) production at 13 TeV compared with LHCb data [36].

B-meson production at very large rapidity
Keeping the same fragmentation function, we observe a deviation between our calculations and LHCb data in the rapidity range 4 < y < 4.5, see Fig.  12. Changing parameter b of the fragmentation function restores the agreement between theory and experiment. However, this change is unwanted as it implies a rapidity-dependent fragmentation function which is not part of the k t -factorization formalism, at least in its simplest formulation. Our understanding of Fig. 12 is that it shows the limit of uPDFs built from collinear PDFs. Indeed, observables at large rapidities trigger smaller and larger values of x compared to central rapidity. From Figs. 2 and 3, we observe that the agreement between integrated uPDFs and collinear PDFs Figure 12: (B + + B − ) production at 7 TeV compared with LHCb data [36] in the rapidity range 4 < y < 4.5. at x = 0.01 is not as good as at x = 0.001, and the situation is even worse at larger x. The agreement between integrated uPDFs (built from collinear PDFs) and collinear PDFs is, in general, not perfect because relation (3) holds only approximately. This relation presents other issues discussed in Ref. [13]. We expect that uPDFs extracted directly from data will give a better description, and it seems to be the case. In Fig 13, we show a comparison between LHCb data for B mesons and our calculations, with our uPDFs replaced by the PB uPDFs [15]. The latter are obtained in a GM-VFNS by fitting experimental data. Note that the use of the PB uPDFs, extracted at the NNLO, with our NLO off-shell cross section is not consistent. However, in a VFNS, the numerical impact of this inconsistency is small. With the PB uPDFs, the agreement between theory and experiment is good, including in the rapidity range 4 < y < 4.5. 4 Figure 13: Results obtained with the PB uPDFs [15]. The LHCb data [36] are described satisfactorily in both rapidity ranges 2 < y < 2.5 and 4 < y < 4.5.

Comparison with theory
In this section, we discuss other k t -factorization calculations and perform a comparison with our results. In Fig. 14 [29]) with our calculations and ALICE data [29].
and [30]. The calculations by Nefedov et al. include the R + R → g and R + R → c +c processes, where R is a reggeized gluon. The outgoing gluon is transformed into a D meson by a scale-dependent fragmentation function D g→D (z, µ 2 ). We observe a good agreement of Nefedov et al. calculations with data but on a limited p t range. It should be mentioned that the description of low-p t data is the hardest part. Szczurek et al. calculations include the g + g → c +c process and do a fair description of the data. However, the situation has been improved by our calculations, which show a better accuracy and stability.
The situation is quite different at forward rapidity, where we compare Nefedov et al. [30] and Szczurek et al. [9] results with LHCb data [38], see Fig. 15. Indeed, the results obtained by these two groups are less satisfying, while our calculations show the same accuracy. They use the WMR approach Figure 15: Comparison of Nefedov et al. [30] and Szczurek et al. results [9] with our calculations and LHCb data [38] at 7 TeV and 2 < y < 2.5.
with the angular ordering cutoff to build uPDFs from collinear PDFs. The main point is that they choose a set of PDFs determined in a VFN scheme. In section 2 we reminded that once the order and scheme of PDFs have been fixed, there is no choice left on the partonic processes to be included. Considering their choice of PDFs, we believe that Szczurek et al. would improve their computation of the partonic cross section by applying the VFN scheme, and therefore include the non-negligible c + g → c + g process. The absence of this process explains the observed underestimation. In fact, the underestimation should be even worse, but it is partially compensated by the use of the too large angular-ordered WMR uPDFs [12,13]. 5 The solution proposed in [40] is not in agreement with factorization, since the authors increase the D-meson cross section by adding some O(α 3 s ) and O(α 4 s ) corrections to the partonic cross section, but keep working with O(α 2 s ) PDFs (for default calculations). A similar discussion could apply to Nefedov et al.
calculations since they use the same uPDFs, but the situation is unclear as they include a 2 → 1 process. 6 Finally, the comparison of LHCb data for B + mesons with calculations by Nefedov et al. [31] is presented in Fig. 16. Note that they obtained  [31] with our calculations and LHCb data [36] at 7 TeV and 2 < y < 2.5. We performed a linear interpolation of theory data points in order to compute the ratio to experiment. The ratio between Nefedov et al. and LHCb data at p t = 20 GeV is not shown because the distance between the two last pink circles is too large and linear interpolation not accurate enough. better results at 3 < y < 3.5, with only a reasonable overestimation of the experimental cross section at p t > 15. In conclusion, k t -factorization results have been improved by our O(α 2 s ) GM-VFNS calculations. We believe this is directly related to the consistent use of factorization. Insisting on the importance of working fully either in a VFNS or in a FFNS was the main goal of this work. Our results can be improved further by including the 2 → 1 process, discussed quickly in section 7, and scale-dependent fragmentation functions.
7 Treatment of the 2 → 1 process It is sometimes believed that there is a double counting if we take into account both the gg → QQ and gQ → gQ processes. It is not the case since double counting happens between two different orders, for instance, when one adds the NLO contribution to a LO calculation. There is consequently no double counting between two LO or NLO processes. On the opposite, the leading-order contribution in k t -factorization is given by the 2 → 1 process, and we can expect double counting between this process and the 2 → 2 processes used in this work. The goal of this section is to discuss this in more detail.
It is interesting to note that the LO and NLO graphs for deep inelastic scattering are similar to the 2 → 1 and 2 → 2 processes in k t -factorization. At leading order, the virtual photon scatters from a (heavy) quark producing an on-shell quark 7 . At next-to-leading order, we find 2 → 2 processes such as γ * +g → Q+Q. As explained in [20], there is a double counting between the LO and NLO contributions, because the production of the second particle in the NLO process is also accounted for by the evolution equation of the LO process. A subtraction term is required, and the cross section reads In principle, finding the subtraction term is automatic. One should apply first collinear factorization at the partonic level, allowing one to determinê σ, where the hat means that this partonic cross section is free of infrared divergences. In DIS at order α s , the partonic cross section for heavy-quark production with the process initiated by the parton of flavor a is where (0) and (1)  is associated to the partonic process γ * +b → Q+X. Ref. [20] gives an explicit example for a = g. Using the fact thatσ where we replaced b by Q in the last equality, the LO process being γ * +Q → Q. This procedure automatically generates a term with a minus sign: the subtraction term.σ (1) a is free of infrared divergences due to the cancellation between σ a . In a second step, factorization is applied at the hadronic level, where f a/h gives the distribution of a parton of flavor a in the hadron h. We believe that the fully consistent treatment of heavy-quark production in k t -factorization is given by Eq. (17), and we plan to work on this soon. However, it is possible to understand why the 2 → 2 process alone provides a good description of the experimental data. On the one hand, in DIS, there is a cancellation between the subtraction term and the LO term in the region µ > m Q and α s ln(µ/m Q ) ∼ O(1). It is the region studied in this work, and we expect a similar situation in k t -factorization. On the other hand, we can argue that the 2 → 1 process at NLO, with a loop in the conjugate amplitude, is negligible with respect to the 2 → 2 process. Indeed, both contributions are of order O(α 2 s ), but the latter is enhanced by the divergences discussed in Sec. 4. Consequently, we can expect that the full contribution, Eq. (17), will be numerically close to the result given by the 2 → 2 process alone.

Conclusion
We have presented for the first time k t -factorization calculations for heavyquark production, using a GM-VFN scheme. We use uPDFs and an off-shell cross section defined in this scheme at the order O(α 2 s ), sometimes misleadingly called leading order. Indeed, the leading order corresponds to the 2 → 1 process, and Eq. (17) shows how the LO and NLO contributions should be put together. We expect a nearly complete cancellation between the LO contribution and the subtraction term at moderate p t , as well as a negligible role of the 2 → 1 process at one loop, justifying the use of the 2 → 2 contributions alone. However, the implementation of Eq. (17) could be important in the region µ m Q , with potential application to jet physics.
The agreement between our calculations and experimental data is excellent, for both D and B mesons. In particular, the agreement has been improved compared to older k t -factorization calculations, performed either in an FFN scheme or in a mix of VFN and FFN schemes. However, we have seen that at very large rapidity our calculations fail if we use uPDFs built from collinear PDFs. In conclusion, accurate and consistent k t -factorization calculations should be performed following the example of collinear factorization: the uPDFs should be exctrated from data, and the cross section computed in a scheme identical to the one used for the uPDFs. Figure 18: Results for D + s compared to ALICE data at 7 TeV [29].