Production mechanisms of open-heavy flavor mesons

In this paper we discuss different mechanisms of open-heavy flavor meson production. Using the color dipole framework, we analyze in detail the contributions of the conventional two-pomeron fusion and the three-pomeron fusion correction. In a parameter-free way we found that the three-pomeron mechanism is significant for $D$-meson production in the small-$p_{T}$ kinematics, although it is less important at large $p_{T}$, as well as for $B$-mesons. The inclusion of the three-pomeron mechanism significantly improves the agreement of theoretical predictions with experimental data in the small-$p_{T}$ kinematics. We also consider the non-prompt charmonia production, and demonstrate that the theoretical results are in reasonable agreement with experimental data. Finally, we compare the theoretical predictions for the dependence on multiplicity of co-produced hadrons to experimental data recently measured by the ALICE collaboration. We found that, contrary to naive expectations, the contribution of the three-pomeron mechanism has only a mild effect on the self-normalized observables in the range of multiplicities studied at ALICE, and for this reason the two-pomeron fusion mechanism can describe reasonably well the experimentally observed multiplicity dependence.


I. INTRODUCTION
Hadrons containing heavy quarks present a widely used tool to test the predictions of Quantum Chromodynamics (QCD). In the heavy quark mass limit the dynamics of a heavy quark can be described perturbatively [1], which allows to test the perturbative QCD (pQCD) predictions. For this reason the production of heavy mesons has been extensively studied in the literature (see e.g. [2][3][4][5][6][7][8][9][10] for overview), and a reasonable description of the experimental data on inclusive production has been achieved. However, the existing theoretical models are constantly being challenged by the improving precision of the available data and the technical advances which make it possible to measure more complicated observables. In fact, the start of the High Luminosity Run 3 at LHC (HL-LHC mode) [11][12][13] will significantly enhance the available data and will give the possibility to analyze the mechanisms of different processes.
One of the observables which can be measured, thanks to the large luminosity, is the dependence of the cross-section on yields (multiplicity) of the charged particles co-produced together with a given heavy meson [14][15][16][17][18][19]. Since the charged particles are produced nonperturbatively, this observable allows to test an interplay of the soft and hard physics, while, as was explained in [20,21], the high multiplicity events allow to test the physics in a deeply saturated regime, which otherwise would require significantly larger energies. Recent experimental data [15][16][17][18][19] show that the yields of S-wave quarkonia (J/ψ, ψ(2S), Υ(1S)) grow rapidly as a function of the multiplicity of charged particles. Such behavior is at tension with conventional two-pomeron fusion mechanisms, and potentially might signal that there are sizable contributions from three-pomeron fusion mechanisms [20,21], which were previously disregarded as higher twist effects. For this reason it is important to revisit the analysis of open heavy-flavor (D-and B-) meson production and check if the conventional mechanisms can describe the multiplicity dependence. In the case of Dmesons such dependence was recently measured by the ALICE collaboration [14], while in B-meson production to the best of our knowledge there is no direct experimental data on the multiplicity dependence 1 , yet there are data on the multiplicity dependence of non-prompt J/ψ production, which proceeds via B → J/ψ decays [14]. These data allow to test the predicted multiplicity dependence for the case of heavier b-quarks. The range of multiplicites in the currently available experimental data is quite limited, because the statistics falls rapidly as a function of event multiplicity, but we expect that it will be significantly extended with data from Run 3 at LHC (HL-LHC mode) [11][12][13]. Another goal of this paper is to estimate the contribution of the three-pomeron mechanisms, which are usually disregarded in the heavy quark mass limit. While in the dipole picture it is believed that a universal dipole cross-section should take into account all such contributions, in phenomenological parametrizations of the dipole cross-section usually such contributions are taken into account only partially or not taken into account at all. For this reason in our explicit evaluation we estimate explicitly the role of such contributions. In particular, since the contribution of the threepomeron mechanism is expected to grow faster than that of the two-pomeron fusion, we pay special attention to the three-pomeron mechanism in large multiplicity events.
The paper is structured as follows. In Section II we discuss the framework used for the evaluation of the openheavy meson production. In Subsection II A we evaluate the contribution of the two-pomeron fusion mechanism and p Figure 1. A typical two-pomeron fusion diagram, taken into account in evaluation of the heavy quark production. In the dipole framework [22][23][24][25] the dipole cross-section is found as a solution of the Balitsky-Kovchegov (BK) equation, so effectively the BK pomeron includes additional fan-like contributions (shown by grey lines; resummation of all possible fan-like topologies is implied). The vertical dashed grey line stands for the unitarity cut, the blob in the lower part is the hadronic target (proton); the fermionic loop in the upper part of the figure includes a summation over all possible gluons. compare its predictions with experimental data. In Subsection II B we evaluate the contributions of the three-pomeron mechanisms and estimate numerically their relative contributions. Our major finding is that they are significantÄ for D-mesons for small p T 10 GeV, yet become negligible for large p T and for B-mesons. In Section III we develop the framework for the multiplicity dependence description in the dipole formalism and compare its predictions for the multiplicity dependence with available experimental data. Finally, in Section IV we draw conclusions.

II. EVALUATION OF THE INCLUSIVE CROSS-SECTION
In this section we will focus on the production of open heavy-flavor mesons (D-and B-mesons). The cross-section for heavy meson production can be related to the cross-section for heavy quark production as [4,5,9,10].
where y is the rapidity of the heavy meson (D-or B-meson), y * = y − ln z is the rapidity of the heavy quark, p T is the transverse momentum of the produced D-meson, D i (z) is the fragmentation function which describes fragmentation of the parton i into a heavy meson, and dσ pp→QiQi+X /dy * is the cross-section of heavy quark production with rapidity y * . The fragmentation functions for the Dand B-mesons, as well as non-prompt J/ψ production, are known from the literature (see the Appendix B for details). Since the dominant contribution to all mentioned states stems from the heavy cand b-quarks (prompt and non-prompt mechanisms respectively), their production can be evaluated in the heavy quark mass limit, and for this reason in what follows we will focus on the evaluation of the cross-section dσ pp→QiQi+X /dy * d 2 p * T , which appears in the integrand of (1).

A. Two-pomeron contribution
The conventional mechanism widely used for description of the heavy meson production is a pomeron-pomeron fusion (see Figure 1). The corresponding cross-section in the dipole approach is given by [9,10] dσ pp→QiQi+X (y, √ s) where y and p T are the rapidity and transverse momenta of the produced heavy meson in the center-of-mass frame of the colliding protons; k T is the transverse momentum of heavy quark, g (x 1 , p T ) in the first line of (2) is the unintegrated gluon PDF; Ψ g→QQ (r, z) is the light-cone wave function of theQQ pair with transverse separation between quarks r and the light-cone fraction carried by the quark z, and we use for it the standard perturbative expressions [26] The meson production amplitude N M depends on the mechanism of QQ pair formation. For the case of the twopomeron fusion, it is given in leading order by [10] (see also Appendix A) For the p T -integrated cross-section the gluon uPDF x 1 g (x 1 , p T − k T ) must be replaced with the integrated gluon PDF x g G (x g , µ F ) , which should be taken at the scale µ F ≈ 2 m Q . In the LHC kinematics at central rapidities (our principal interest) this scale significantly exceeds the saturation scale Q s (x), which justifies the use of two-pomeron approximation. However, in the kinematics of small-x g (large energies) there are sizable non-perturbative (nonlinear) corrections to the evolution in the dipole approach, and therefore in this kinematics the corresponding scale µ F should be taken at the saturation momentum Q s . In this approach the gluon PDF x 1 G (x 1 , µ F ) is closely related to the dipole scattering amplitude N (y, r) =´d 2 b N (y, r, b) as [27,28] where y = ln(1/x). Eq. (9) can be inverted and it gives the gluon uPDF in terms of the dipole amplitude, The corresponding unintegrated gluon PDF can be rewritten as [29] x g x, which allows to rewrite the result in a symmetric and self-consistent form, which in turn permits a straightforward generalization of the high-multiplicity events.
Here and in what follows, for our numerical evaluations we will we use the "CGC" parametrization of the dipole cross-section [23] N  dσ/dp T [μb/GeV] Figure 2. The pT -dependence of the cross-section dσ/dpT for D + -mesons, evaluated with the two-pomeron fusion mechanism and integrated over the rapidity bin. Left plot: comparison with data in the LHC kinematics, at central and forward rapidities.
The experimental data are from [32][33][34]. Right plot: Comparison with experimental data from the Tevatron at central rapidities.
The experimental points are from the CDF and D0 collaborations [30,31]. For other mesons the pT -dependence has a similar shape, although it differs by a numerical factor of two (a more detailed comparison with data can be found in [10,35]   The experimental data are from CMS [36](" √ s=7 TeV, |y| < 2.1" data points), ATLAS [37]( √ s=7 TeV, |y| < 0.5 data points) and CDF [38] ( √ s=1.96 TeV, |y| < 0.6 data points). Right plot: The pT -dependence of the cross-section dσ/dy dpT for nonprompt J/ψ. Comparison with experimental data from CMS [39] ( √ s = 5 TeV data) and CDF [40]( √ s = 1.96TeV data) at central rapidities. For ψ(2S) the pT -dependence has a similar shape and differs only by normalization. In both plots, for some experimentally measured bin-integrated cross-sections dσ/dpT , it was converted into dσ/dpT dy dividing by the width of the rapidity bin (this is justified since in the LHC kinematics at central rapidities y ≈ 0 the cross-section is flat).
In Figures 2, 3 we show the p T -dependence for both D-meson and B-meson production, as well as for the case of non-prompt J/ψ mesons. We can see that in the large p T region the two-pomeron mechanism describes very well all the available data. At small-p T 5 GeV there are no direct measurements for B-mesons, although there are data for non-prompt J/ψ (from decays of the B-mesons), and we can see that the model describes the data available from Tevatron [30,31]. However, for D-mesons the agreement is marginal in this kinematics, and the two-pomeron mechanism systematically overestimates the experimental data by more than 2 σ. Such behavior is not related to technical details of our evaluation (like the choice of the dipole cross-sections or fragmentation functions) and was also observed by other authors (see e.g. [9,10]). Since this small-p T region gives the dominant contribution to the p Tintegrated cross-section, the two-pomeron mechanism will also overestimate this observable. As we will demonstrate in the next section, the agreement with data in the small-p T kinematics improves after the inclusion of the multigluon contributions.   [22][23][24][25] of the color singlet dipole cross-section N (z, r) (resummation of all possible tree-like topologies is implied). Right plot: The BFKL ladder diagrams resummed in the IP-Sat (b-Sat) parametrization [26,41]. In both plots a vertical dashed grey line stands for the unitarity cut, the blob in the lower part is the hadronic target (proton); two fermionic lines in the upper part of the blob stand for the dipole of the transverse size r.

B. Three-pomeron contribution
As we discussed in the introduction, for the c-quarks potentially there could be a sizable contribution from the 3-gluon fusion mechanism. While usually it is believed that such contributions are suppressed by α s (m Q ), and in certain cases additionally by Λ 2 QCD /m 2 Q , we have seen from [20,21] that potentially such contributions might give a sizable correction for charmonia production, especially in large-multiplicity events. In the framework of the dipole model it is usually assumed that the universal dipole cross-section takes into account all such contributions. However, in phenomenological parametrizations usually such contributions either are taken into account with additional simplifying assumptions or some of the contributions are disregarded. For example, a widely used phenomenological parametrization "CGC" suggested in [22][23][24][25], was inspired by a solution of the Balitsky-Kovchegov (BK) equation and effectively resums only fan diagrams shown in the left panel of Figure 4. This parametrization does not take into account the three-pomeron contributions at all. The competing IP-Sat parametrization [26,41], which was inspired by Glauber-like approach, resums thefor set of BFKL ladder diagrams shown in the right panel of Figure 4. A central assumption which allows for numerical simplifications is that the interaction of the BFKL ladder (pomeron) with a dipole of size r is given by ∼ α s µ 2 r 2 x g(x), which might work for small color-singlet dipoles, but in general cases requires a more careful treatment. Although for sufficiently small dipoles the predictions of both approaches agree with each other, for the subleading terms the CGC and IP-Sat dipole cross-sections might differ significantly. For this reason in general we cannot extract the contribution of the three-pomeron mechanism from (2, 8) and need to evaluate it explicitly. However, we should take into account that in contrast to J/ψ production at the same order of perturbation theory, we may get also interference terms of the leading-order with subleading order contributions. Since we work in the eikonal approximation, these diagrams will differ only by a numerical (combinatorial) factor. Due to these interference contributions the correction is not positively defined.
As was demonstrated in the Appendix A, for the three-pomeron contribution we can show that the corresponding cross-section is given by where and σ eff ≈ 20 mb is an effective cross-section discussed in detail in A. Both functions N ± (z, r 1 , r 2 ) are invariant with respect to permutation r 1 ↔ r 2 . For the p T -integrated cross-sections it is possible to show that the integration reduces to r 1 = r 2 = r, so the cross-sections N ± simplify tõ For the interference term we get in a similar way In general the contribution of the interference term (20) is negative, and larger by magnitude than the direct contribution (15), and for this reason the total correction of the three-pomeron mechanism in general is negative. This contribution is strongly suppressed at large p T because in this kinematics a typical dipole size r ∼ 1/p T , and the contributions (15,20) have an additional suppression ∼ O r 2 ∼ O p −2 T compared to (2,8). In order to illustrate the relative size of the three-pomeron mechanism (15) and the interference term (20), in Figure 6 we plotted the ratio of the cross-sections evaluated with three-pomeron and two-pomeron mechanisms, We can see that for the c-quarks at small p T both contributions might be substantial and constitute up to a factor of two correction. For the b-quarks it does not exceed ten per cent even for p T ≈ 0, in agreement with heavy quark mass s =7 TeV, y=0 c→D + , All Figure 6. The relative contribution of the 3-gluon to the 2-gluon mechanism, as defined in (21). The curves with labels "c → D + " and "b → D + " correspond to prompt and non-prompt contributions to D + -meson production (for other D-mesons the results are similar). The additional label "|3 IP | 2 " in some curves implies that for the contribution of the 3-pomeron fusion cross-section only the contribution (15) was taken into account, whereas for the curves with label "All" we also took into account the contribution of the interference term (20). We can see that for c-quarks the contribution of the 3-gluon mechanism in the small-pT kinematics is significant and changes the result by a factor of two, whereas for b-quarks it is just a minor correction which does not exceed 10% even for pT ≈ 0. For large pT the relative contribution decreases for all quark flavors, and for pT 10 GeV it becomes negligible.
limit expectations. In the large-p T kinematics the relative weight of the three-pomeron contribution is suppressed and does not exceed a few per cent for p T 10 GeV. In Figure 7 we show the p T -dependence of the cross-section, taking into account both two-and three-pomeron mechanisms. We can see that in the region of small-p T the agreement with data is much better than with just the two-pomeron mechanism shown in Fig 3. For this reason in what follows we will take into account both the the two-and three-pomeron mechanisms.

A. Theoretical framework
As was illustrated in the previous section, the dipole approach (2,8,15,20) with a CGC dipole parametrization provides a very reasonable description of the inclusive Dand B-meson production. Nevertheless, the description of the multiplicity dependence presents more challenges at the conceptual level, because there are different mechanisms to produce enhanced number of charged particles N ch . The probability of multiplicity fluctuations decreases rapidly as a function of number of produced charged particles N ch [42], and for this reason in the study of the multiplicity dσ/dp T [μb/GeV] Figure 7. The pT -dependence of the cross-section dσ/dpT for D + -mesons, evaluated taking into account both the three-pomeron and interference contributions. We use the same notations as in the Figure (2); the data are integrated over the rapidity bin.
Left plot: comparison with data in the LHC kinematics at central and forward rapidities. The experimental data are from [32][33][34]. Right plot: Comparison with experimental data from the Tevatron at central rapidities. The experimental points are from the CDF and D0 collaborations [30,31]. For other mesons the pT -dependence has a similar shape, although it differs by a numerical factor of two (a more detailed comparison with data might be found in [10,35]).
dependence it is more common to use a normalized ratio [43] dN M /dy dN M /dy where n = N ch / N ch is the relative enhancement of the number of charged particles in the pseudorapidity window (η − ∆η/2, η + ∆η/2); N ch = ∆η dN ch /dη is the average number of charged particles in the pseudorapidity window (η − ∆η/2, η + ∆η/2);w (N M ) / w (N M ) and w (N ch ) / w (N ch ) are the self-normalized yields of heavy meson M (M = D, B) and charged particles (minimal bias events) in a given multiplicity class; dσ M (y, √ s, n) is the production cross-sections for heavy mesons M with rapidity y and N ch = n N ch charged particles in the pseudorapidity window (η − ∆η/2, η + ∆η/2), whereas dσ ch (y, √ s, n) is the production cross-sections for N ch = n N ch charged particles in the same pseudorapidity window. If the inclusive cross-section of the process pp → M + X is proportional to the probability to produce a meson M in a single pp collision, then the ratio (22) gives a conditional probability to produce a meson M in a pp collision in which N ch charged particles were produced. Due to Local Parton-Hadron Duality (LPHD) hypothesis [44][45][46] the number of produced charged particles is directly proportional to the number of partons which stem from the individual pomerons and thus can be studied using perturbative methods.
In the color dipole approach analyzed in this paper, we expect that the multiplicity dependence is enhanced due to a larger average number of particles produced from each pomeron. Nevertheless, we still expect that each such cascade ("pomeron") should satisfy the nonlinear Balitsky-Kovchegov equation, and for this reason we expect that the dipole amplitude (11) should maintain its form, although the value of the saturation scale Q s might be modified. As was demonstrated in [27,47,48], the observed number of charged multiplicity dN ch /dy of soft hadrons in pp collisions is given by the so-called KLN-style formula where c is a numerical coefficient, and N I P is the number of BK pomerons. Solving algebraic Eq.(23), we could extract Q 2 s as a function of dN ch /dy. Taking into account that the distribution dN ch /dy is almost flat, we may approximate n = N ch / N ch ≈ (dN ch /dy)/ dN ch /dy , so (23) allows to express Q 2 s as a function of n. Frequently in the literature the logarithmic dependence on n, which stems from the running coupling in denominator of (23), is disregarded, and therefore (23) reduces to a simpler linearly growing dependence on n [27,[47][48][49][50][51][52], The precision of this assumption was tested in [9], and it was found that a numerical solution of the running coupling Balitsky-Kovchegov (rcBK) equation differs from the approximate (24) by less than 10% in the region of interest (n 10). Since this correction is within the precision of current evaluations, in what follows we will use (24) for our estimates. While at LHC energies it is expected that the typical values of saturation scale Q s (x, b) fall into the range 0.5-1 GeV, from (24), we can see that in events with enhanced multiplicity this parameter might exceed the values of heavy quark mass m Q and lead to an interplay of large-Q s and large-m Q limits. Thus the study of the high-multiplicity events gives us access to a new regime which otherwise would require significantly higher energies.
As was shown in [27,47,[53][54][55], for dilute systems the saturation scale Q s is closely related to the gluon density of the target, This qualitative relation just reflects the fact that for dilute system the saturation scale Q 2 s is proportional to the density of partons (gluons), which is described by the gluon density G. It is tempting to extrapolate the relation (25) to study the multiplicity dependence of the gluon density. However, such interpretation might be useful only for the case when n is not very large, while for the events with very large multiplicity (n 1) the concept of the gluon density becomes quite obscure, since in this case the twist expansion does not work (it is heavily broken by higher order terms). The n-dependence of (24) hints to the fact that the contributions with large number of cut pomerons should be enhanced compared to the n = 1 case. Indeed, in the heavy quark mass limit and for not very large n, the typical dipole size in (2) is given by r ∼ m −1 Q , so from the structure of (11) we can see that this enhancement is given by a factor ∼ n γ eff . However, in the deeply saturated regime (n 1), when Q 2 s (x, b; n) m 2 Q , the typical dipole size is controlled by the saturation scale and thus the n-dependence should be the same for all multipomeron contributions. We would also like to mention that the uncut pomerons do not contribute to the observed enhancement of charged particles and thus should not be taken into account in the multiplicity evaluation.
To conclude, the suggested mechanism introduces a dependence on multiplicity of soft produced particles, and it is quite different from other approaches such as the percolation approach [56] or the modification of the slope of the elastic amplitude [57]. Moreover, it can be applied both to the production of soft and hard particles. In the following subsection we will use this approach for analyzing the multiplicity dependence of quarkonia production.

B. Phenomenological estimates
In a typical experimental setup the detector used to collect charged particles N ch usually covers some small rapidity bin ∆η, in which a relative enhancement of multiplicity n = N ch / N ch is observed. This enhancement is a result of superpositions of increased multiplicities from individual pomerons. Since each pomeron hadronizes independently, we may expect that we should apply the modification of the dipole cross-section discussed in the previous section to each of the pomerons, modifying the corresponding dipole amplitude. However, the number of pomerons which can participate in the observed multiplicity enhancement depends crucially on the details of how the experiment is done, as explained in Figure 8.
For the simplest case when the bins used for the collection of mesons and charged particles are well-separated in rapidity (left panel in Figure 8), it is clear that all charged particles can stem only from cut pomerons at a given rapidity (lower pomerons in Figure 8), so the cross-section to produce heavy meson M and N particles can be found as a mathematical expectation, convoluting the probability of a given partition (N 1 , ...N k ) with the value of the corresponding cross-section, viz.: In (27) which implies that a contribution of a pomeron to the observed number of charged particles equals the sum of all possible contributions of its parts. The exact evaluation of (27) requires knowledge of the function P (N, N ), which is an essentially nonperturbative object, studied in the literature in the context of different models (see e.g. [58] for details), and which has been suggested that it might be described by the Poisson distribution. Fortunately, for phenomenological estimates we may minimize the sensitivity to the choice of the model for P (N, N ), taking into account the following facts: • The dependence on n k , which stems from the cross-section in the second line of (27), is very mild and saturates (becomes constant) in the region of large n. In contrast, the functions P (N k , N k ) decrease exponentially for n k ≡ N k / N k 1, and for this reason in the evaluation of (27) we may replace each n k with some average value n k (which in general depends on n).
• Since all active pomerons have the same average number of particles N k , taking into account a very mild dependence of the cross-section on n k , we may apply iteratively (28) and rewrite (27) as n i ≡ n/k, i = 1, ..., k.
Thus effectively we come to the conclusion that the observed multiplicity is shared equally between all the pomerons at a given rapidity window. The convolution of the distributions P (n i ) just cancels in the ratio (22), so the latter can be rewritten as p p (y − ∆y 2 , y + ∆y 2 ) (η 1 − ∆η 2 , η 1 + ∆η 2 ) Figure 9. (color online) Case of "minimal bias" measurement done by the V0 detector at ALICE. The charged particles are collected in both backward and forward directions, and only their sum is used to measure N ch . For the sake of definiteness we consider only the three-pomeron fusion mechanism and for simplicity disregard the fragmentation of heavy quark (convolution with fragmentation function), which slightly smears the distribution over rapidity. The rapidity bin used for the collection of heavy mesons (blue box) does not overlap with the bins used for the collection of charged particles (red boxes); however, the elevated multiplicity cannot be unambiguously attributed to be (shared between) the upper or lower pomerons.
where we use notation dσ instead of dσ for the cross-section, to emphasize that we took out the normalization to probability distribution of charged particles (factor P (N, N )) and share multiplicity equally between all pomerons in a given rapidity window. The situation becomes more complicated for the setup when the heavy meson and charged particles bins partially overlap, as shown in the right panel of Figure 8. In this case in the intersection region the enhanced multiplicity is due to either the upper or to the lower cut pomerons, depending on the position of the heavy quark inside the bin. For the sake of simplicity we will consider the case of complete overlap of both bins, which is realized in the case of ALICE measurements at central rapidities. In this case we have to average over the rapidity of heavy quark in the relation (27), taking into account that for the upper pomerons N k ∼ (dN/dy) (y − y min ), whereas for the lower pomerons N k ∼ (dN/dy) (y max − y). Following the assumptions formulated after Eq. (28), we can see that instead of (30) we should use dN M /dy dN M /dy =´y max y min dy dσ pp→QQ+X (y, η, Q 2 , n)/dy ∆y dσ pp→QQ+X (Y η,, Q 2 , n =1)/dy (31) with values n k = n (y − y min ) / (k up ∆y) and n k = n (y max − y) / (k down ∆y) where k up , k down is the number of pomerons in the the upper and lower parts of the diagram.
Finally, we would like to stop briefly on the case of the so-called minimum bias configuration studied by the V 0 detector at ALICE. In this case the total charge is accumulated in two rapidity bins, in very forward and very backward directions, as seen from Figure 9. In this configuration we do not have overlap of yand η-bins, yet still we cannot assign the enhanced multiplicity either to upper or to lower pomerons. Instead of this, we should use a sum of all possible partitions of the total number of charged particles. The corresponding cross-section is given by (30), although we should use n k = n/k up and n k = n/k down , where k up , k down is the number of pomerons in the upper and lower parts of the diagram.
We would like to start the discussion of the numerical results from estimates of the role of the three-pomeron mechanism. In Figure 10 we plotted the ratio of the three-pomeron and two-pomeron contributions (21), which was discussed in Section II B. We can see that in the large-multiplicity events the relative contribution of the threepomeron mechanism is sizable even for n = 1. As a function of multiplicity this ratio grows, and starting from n ≈ 5 for D-mesons it becomes a dominant contribution, as can be seen from the right panel of Figure 10. Such increase s =7 TeV, y=0 n=1 c→D + , All c→D + , |3 IP 2 b→D + , All b→D + , |3 IP 2  (21). The curves with labels "c → D + " and "b → D + " correspond to prompt and non-prompt contributions to D + -meson production (for other D-mesons the result is similar). The additional label "|3 IP | 2 " in some curves implies that for the contribution of the 3-pomeron fusion cross-section only the contribution (15) was taken into account, whereas for the curves with label "All" we also took into account the contribution of the interference term (20). We can see that for c-quarks the contribution of the 3-pomeron mechanism is substantial and constitutes up to 50% of the total result at small pT , whereas for b-quarks it does not exceed 10% even for pT ≈ 0. For large pT the relative contribution decreases both for the c-and b-quarks, and for pT 10 GeV becomes negligible. Right plot: the same evaluation for high-multiplicity events with n = 5. The total result (sum of direct and interference terms) decreases as a function of multiplicity and nearly vanishes for the high-multiplicity events near n ≈ 5, as explained in the text.
can be understood from a comparison of the structures of (8,15): in the heavy quark mass limit the size of the dipole is given by r Q ∼ 1/max (m Q , p T ), and as we explained in Section III A, in this regime each cut pomeron yields an additional factor ∼ n γ eff . However, such grow cannot continue up to infinity. When the saturation scale Q s becomes significantly larger than all the other scales, the relevant dipole size is given by r Q ∼ 1/Q s and the system saturates, i.e, having a very mild dependence on n. The interference term (20) has the same number of cut pomerons as the two-pomeron mechanism, and for this reason we expect that the ratio (21) should be largely independent on n. Since at n = 1 the contribution of the interference term is larger than (15) and has opposite sign, their sum decreases by absolute value in the range 1 n 5, and changes sign near n ≈ 5. For this reason a combined contribution of the three-pomeron mechanisms (15,20) does not affect the multiplicity dependence significantly in the multiplicity range studied in [14]. We also may notice that due to the increase of the saturation scale Q 2 s , the transition to the large-p T regime happens at significantly larger p T .
Currently the data on the multiplicity dependence of open charm mesons are available from the ALICE experiment [14]. As we can see from Figure 11, our results can perfectly describe the available data on D-meson production. For the sake of definiteness we make the comparison with experimentally available averaged contribution of the D + , D 0 and D * 0 mesons. In the same figure we have shown the contribution of the non-prompt mechanism (dashed lines). As expected, this contribution has the same dependence on n, though numerically its relative size in the inclusive cross-section varies depending on the kinematics (values of p T ).
Finally, in Figure 12 we have shown the multiplicity dependence of the non-prompt J/ψ, which are formed from the decays of the b-mesons. The experimental data clearly show that the multiplicity dependence grows faster than linear. The dipole approach provides a very reasonable description of the available multiplicity dependence.

IV. CONCLUSIONS
In this paper we studied the mechanisms of the open-heavy flavor meson production. We took into account both the standard two-pomeron mechanism, as well as estimated the contribution of three-pomeron fusion. We found that the latter correction is important for D-mesons for small-p T data, where it changes the result by a factor of two and allows to improve considerably the agreement of theoretical predictions with data. The correction is less relevant for B-mesons, where it does not exceed ten per cent. As a function of transverse momentum p T , the relative dN ch /dη 〈dN ch /dη〉 Figure 11. Theoretical multiplicity dependence of the prompt and non-prompt mechanisms of the D-meson production in different bins in pT at central rapidities. The experimental data are from ALICE [14] and correspond to prompt production mechanism for averaged contribution of D 0 , D + and D * + mesons. For the sake of reference we have also shown a dotted line which corresponds to a linear dependence. The charged particles and D-mesons are collected in rapidity window |y| < 0.5, |η| < 1.
weight of the three-pomeron contribution decreases, and for p T 10 GeV the correction does not exceed one per cent. Such behavior agrees with general expectations based on large-p T and heavy quark mass limit evaluations. Since the p T -integrated cross-section is dominated by the small-p T -region, we conclude that it is sensitive to the contributions of the three-pomeron mechanism. Our evaluation is largely parameter-free and relies only on the choice of the parametrization for the dipole cross-section (11).
The suggested approach is able to describe the multiplicity dependence measured by ALICE [14]. Contrary to naive expectations, the relative contribution of the three-pomeron mechanism has a rather complicated dependence on multiplicity. Due to interplay of direct and interference contributions, shown in Figure 5, the relative contribution of the three-pomeron correction decreases as a function of multiplicity for small n, changes sign near n ≈ 5 and starts growing at larger values of n. For this reason this contribution does not lead to a pronounced multiplicity dependence pp, s =7 GeV  [14]. For the sake of reference we have also shown a short-dashed line which corresponds to a linear dependence. It is expected that the charged particles are collected at rapidity window |y| < 0.9, |η| < 1 .
(a) (b) (c) Figure 13. The diagrams which contribute to the heavy meson production cross-section in the leading order perturbative QCD. The contribution of the last diagram (c) to the meson formation can be also viewed as gluon-gluon fusion gg → g with subsequent gluon fragmentation g →QQ. In the CGC parametrization of the dipole cross-section approach, each "gluon" is replaced with a reggeized gluon (BK pomeron), which satisfies the Balitsky-Kovchegov equation and corresponds to a fan-like shower of soft particles.
for the range of multiplicities available from LHC [14]. This result differs dramatically from quarkonia production, where the measurement of multiplicity dependence was suggested as a means to estimate the role of the three-pomeron contributions [20,21]. This difference occurs due to lack of the interference contributions shown in the right panel of Figure 5, in the quarkonia production case. reason it is very instructive to discuss different contributions in parallel with the perturbative k T -factorization-style approach, tacitly assuming that each gluon should be understood as a parton shower ("pomeron"). The rules which allow to express the cross-sections of hard processes in terms of the color singlet dipole cross-section can be found in [59,60]. In the high-energy eikonal picture, the interaction of the quarks and antiquark with a target are given by ±ig t a γ (x ⊥ ), where x ⊥ is the transverse coordinate of the quark, and the function γ (x ⊥ ) is related to a gluonic field of the target. This function is related to a dipole cross-section σ(x, r) as where r is the transverse size of the dipole, and z is the light-cone fraction of the dipole momentum carried by the quarks. The equation (A1) can be rewritten in the form For very small dipoles, the dipole cross-section is related to the gluon uPDF 2 σ (x, r) = 4πα s 3ˆd so the functions γ (x, r) can be also related to the unintegrated gluon densities. For many high energy processes dominated by the pomeron-pomeron fusion mechanism it is possible to express the exclusive amplitude or inclusive cross-section as a sum of the contributions which have the same structure as the left-hand side of (A2). For some processes the last term in (A2) eventually cancels after summation over all possible diagrams, so the color dipole density matrix becomes expressed in terms of the linear combination of the color singlet dipole cross-sections σ(x, r) with different arguments. While in the deeply saturated regime we can no longer speak about individual gluons (or pomerons), we expect that the relations between the dipole amplitudes and color singlet cross-sections should be valid even in this case. This is a crucial assumption which essentially constitutes one of the elements of multiplicity dependence discussed in Section II B, and which gives a good description of the multiplicity dependence [20,21].
For the case of D-meson production, the leading-order contribution is given by the diagrams shown in Figure (13) and yields for the cross-section the result given in (2,8) (see [59,60] for details). In the evaluation of the p T -dependence, we should project the amplitude in coordinate space (state with definite quark coordinate r Q ) onto the state with constant momentum p T , by taking an additional Fourier transform´d 2 p T exp (ip T · r). After squaring the amplitude in momentum space, this implies the inclusion into (A1) of the additional factor ∼´d 2 r 1 d 2 r q e ip T ·(r1−r2) , where r 1,2 are the coordinates of the quark in the amplitude and its conjugate. In the frame where the momentum of the primordial gluon is not zero, we get an additional convolution with the p T -distribution of the incident ("primordial") gluons, as shown in (2), and was demonstrated in [10] 3 .
For the three-pomeron contribution the above-given approach can be extended. However, for the description of the interaction with the target we need to model the multigluon interactions, which are not taken into account by (A1). In the perturbative limit, the corresponding interactions are described in terms of the so-called Double Parton Distribution Functions (DPDFs) (see [61,62] for a review and discussion). In general these objects have a complicated structure, and are not related to the gluon uPDFs. However, at high energies the correlations between the partons are negligible [63][64][65], so the DPDFs can be expressed as products of independent uPDFs. Thanks to this property, the cross-sections of the so-called Double Parton Scattering processes can be represented as products of Single Parton Scattering (SPS) processes. In the Color Glass Condensate model [66][67][68] this assumption is fulfilled automatically. In the color dipole approach, we expect that the cross-section will be given by a linear combination of structures 4 k=1 γ( r k ), which in view of (A1) could be separated into a sum of products of color singlet dipole amplitudes.
Taking into account all the diagrams shown in Figure 14, we obtain for the amplitude of the three-pomeron process a is the color index of the incident (projectile) gluon, b and c are the color indices of the gluons attached to the target (verticalt-channel gluons in Figure 14), r Q , rQ are the coordinates of the quarks. For the evaluation of the cross-section we should square the amplitude, and potentially could get different structures with bc = b c . Indeed, as was demonstrated in [61,69]. For the Double Parton Distribution Functions (DPDFs) the corresponding cross-section is described by 6 different color structures, which take into account possible color state of the gluon pairs in the amplitude and its conjugate, viz.: δ bb δ cc , f bb k f cc k , d bb k d cc k , t bb , cc However, as was illustrated in [70], the largest intercept has a configuration when the two gluons are in a relative color singlet state (two cut pomerons), which corresponds to the first term in (A4). This configuration dominates at high energies, and for this reason in what follows we will take into account only this contribution.
For the evaluation of the p T -dependent cross-section we need to project the coordinate space quark distribution onto the state with definite transverse momentum p T , and so we have for the square of the amplitude For the p T -integrated cross-section these formulas simplify since we will have to put x i ≡ y i , i = Q,Q. As discussed earlier, at high energies the correlations between the partons are negligible, and the two gluons reggeize independently and are in color singlet state with respect to each other [64], so the four-pomeron configuration, up to numerical factor ∼ σ −1 eff ≈ (20 mb) −1 , can be found from (A5) applying iteratively the relation (A1). It is possible to demonstrate that after such procedure we can express the three-pomeron dipole amplitude in terms of the color singlet dipole cross-sections as given in (15). The evaluation of the contribution (20) follows a similar algorithm, although we have to take into account that one of the pomerons is uncut, and this reason it does not contribute to the growth of the multiplicity.  Table I. The values of parameters used for evaluation of the D-meson fragmentation function with parametrization (B4) (see [72] for details). evaluated in detail in [5]. Due to space limitations we do not write out the full expressions for this function and instead in Figure 15 we compare the fragmentation functions D b→B and D b→J/ψ . These two functions differ by more than two orders of magnitude, and for this reason in order to facilitate comparison, we plotted the distributions normalized to unity,D(z) = D(z)/´1 0 dz D(z). As we can see, the distribution D b→J/ψ is significantly wider and has a peak near smaller values of z ≈ 0.5.
For the case of D-mesons we should take into account that there are two complementary mechanisms, a direct (prompt) production, and indirect (non-prompt) mechanism from decays of B-quarks. In both cases we use a fragmentation function taken from [72], with parameters given in Table I. Despite of the significant difference between the values of constants between D + and D 0 mesons, the two parametrizations have very similar shapes and differ only by a factor of two in normalization.