# Event engineering heavy flavor production and hadronization in high multiplicity hadron-hadron collisions

###### Abstract

Heavy flavor measurements in high multiplicity proton-proton and proton-nucleus collisions at collider energies enable unique insights into their production and hadronization mechanism because experimental and theoretical uncertainties cancel in ratios of their cross-sections relative to minimum bias events. We explore such event engineering using the Color Glass Condensate (CGC) effective field theory to compute short distance charmonium cross-sections. The CGC is combined with heavy-quark fragmentation functions to compute -meson cross-sections; for the , hadronization is described employing Nonrelativistic QCD (NRQCD) and an Improved Color Evaporation model. Excellent agreement is found between the CGC computations and the LHC heavy flavor data in high multiplicity events. Event engineering of the contributions of different quarkonium intermediate states in the CGC+NRQCD framework reveals that hadronization in rare events is dominated by fragmentation of the state.

###### pacs:

11.80.La, 12.38.Bx, 14.40.Lb, 14.40.PqIntroduction: The study of high multiplicity events in proton-proton () and proton-nucleus () collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) has focused attention on the spatial and momentum structure of rare parton configurations in the colliding projectiles obtained by variations in the multiplicity, energy and system size. Such “event engineering” first revealed the remarkable systematics of “ridge” like rapidity separated azimuthal angle hadron correlations, triggering debates regarding their initial state Dumitru:2010iy (); Dusling:2012iga (), or hydrodynamic origins Bozek:2011if (); Bozek:2012gr ().

Heavy flavor measurements add important elements to the discussion. The large quark masses provide a semi-hard scale to probe initial state dynamics. The possible diffusion of heavy quarks can help quantify final state interactions. A compelling example of event engineered heavy flavor measurements in and collisions at RHIC and the LHC are ratios of their yields in high multiplicity events relative to minimum bias events. When plotted versus event activity, the ratio of charged hadron multiplicity in rare relative to minimum bias events, many model dependencies cancel out. In particular, because nonperturbative features of hadronization are likely the same for both rare and minimum bias events, ratios of heavy flavor multiplicities are sensitive primarily to short distance interactions of intermediate states.

The exciting possibility that event engineering may help distinguish between intermediate states can be quantified in the Nonrelativistic QCD (NRQCD) Bodwin:1994jh () framework, wherein the inclusive differential cross-section of a heavy quarkonium state in and collisions is expressed as

(1) |

where are quantum numbers of the produced intermediate heavy quark pair, with , and denoting its spin, orbital, and total angular momenta, respectively. The symbol denotes a color singlet (CS, ) or color octet (CO, ) state. The are perturbative short distance coefficients for heavy quark pair production with quantum numbers and are universal nonperturbative long distance matrix elements (LDMEs) . The LDMEs can for instance be extracted from data on Quarkonium production at the Tevatron, and employed to make predictions for cross-sections at the RHIC and LHC. While NRQCD is successful, there are puzzling features in its comparisons with data. An important puzzle is that the magnitude of the linear combination of the and LDMEs extracted from hadroproduction data Ma:2010yw (); Butenschoen:2010rq () is larger than an upper bound set by BELLE data Zhang:2009ym (). This breaking of universality of the LDMEs brings into question NRQCD factorization.

In this work, we will show that the systematics of heavy flavor production in rare events in and collisions are sensitive to strongly correlated gluons in the colliding protons and nuclei. The strong correlations arise because rare events select hadron fluctuations with high parton occupancies. The dynamics of such configurations is controlled by an emergent semi-hard saturation scale , where is the longitudinal momentum fraction carried by a parton in the hadron Gribov:1984tu (); Mueller:1985wy (). Gluon modes with have the maximal occupancy possible in QCD, while modes are the hard partons of perturbative QCD (pQCD). Since grows with decreasing , and increasing nuclear size, the interplay of the dynamics of hard and soft modes evolves with the changing energy and centrality of the collision.

A systematic framework to study of gluon saturation is the Color Glass Condensate (CGC) effective field theory (EFT) Gelis:2010nm (); Kovchegov:2012mbw (); Blaizot:2016qgz (). The cross-sections for the production of heavy quarkonia in the CGC EFT for hadron-hadron collisions were computed over a decade ago Gelis:2003vh (); Blaizot:2004wv (); Tuchin:2004rb (); Fujii:2005vj (); Fujii:2006ab (). A more recent development is the CGC+NRQCD framework^{1}^{1}1See Dominguez:2011cy () for a specialized discussion. Kang:2013hta (), the novel element being that in Eq. (1) is computed in the CGC EFT. There are several phenomenological studies of data from RHIC and LHC that employ these computations in and collisions Fujii:2013gxa (); Ma:2014mri (); Ma:2015sia (); Watanabe:2015yca (); Ducloue:2015gfa (); Fujii:2013yja (); Fujii:2015lld (); Ducloue:2016pqr (); Fujii:2017rqa (); Ma:2017rsu (). High multiplicity configurations are approximated by increasing the value of at the input large scale in both protons and nuclei in multiples of , the initial saturation scale determined from fits to the minimum bias DIS data Tribedy:2010ab (). Increasing the saturation scales at large is a simple way of representing the fluctuations of protons and nuclei into larger numbers of color charges in rare events. A similar procedure was followed in studies of ridge yields Dusling:2013qoz (); Dusling:2015rja (); Schenke:2016lrs (). The systematic treatment of high multiplicity “biased” color charge configurations has only just begun in the CGC EFT Dumitru:2017cwt (); Dumitru:2017ftq (); Dumitru:2018iko (): a desirable outcome of our study is to spur further effort.

We will focus here^{2}^{2}2The and open bottom computations require Sudakov resummation Mueller:2012uf (); Mueller:2013wwa (); Qiu:2013qka (); Qiu:2017xbx () and are beyond our scope here. on measurements of and mesons in high multiplicity and collisions Adam:2015ota (); Adam:2016mkz (); Abelev:2012rz (); Weber:2017hhm (); Khatun:2017yic (); Adamova:2017uhu (); Ma:2015xta (); Federicova:2017bmd (). These data demonstrate that the production yields of and in high multiplicity events are significantly enhanced relative to minimum bias events. The models proposed to explain their systematics include percolation models Ferreiro:2012fb (); Ferreiro:2015gea (), dipole models Kopeliovich:2013yfa () and multiparton interaction models Thakur:2017kpv (). All these models approximate effects contained in the CGC EFT. Gluon saturation is included in the EPOS3 model Werner:2013tya (), which also includes final state scattering effects. As we will show, the CGC+NRQCD EFT can address detailed differential questions regarding heavy flavor production mechanisms; further, it may allow us to use high multiplicity events to resolve extant heavy flavor puzzles in collider experiments.

Open flavor and quarkonium production: We first consider the spin and color averaged inclusive cross-section , which can be expressed in the CGC EFT as Blaizot:2004wv ()

(2) |

where , , ,with () and (), the transverse momentum and rapidity respectively of the produced charm (anti-charm) quarks. Further, and , where , denote the fractional longitudinal momentum fractions of the interacting gluons in the projectile and target respectively. The expression for the hard scattering matrix element is listed in Appendix A. The unintegrated gluon distribution function (UGDF) of the projectile proton is defined as Ma:2014mri ()

(3) |

Here () is the transverse area occupied by gluons in the proton (nucleus) and

(4) |

The fundamental Wilson line in the amplitude (complex conjugate amplitude) () represents multiple scattering of the quark with background fields at the position (). The expression corresponds to leading log x resummation in the CGC EFT and must not be confused with the LDMEs expectation value in Eq. (1).

The differential cross section for meson production is then given by

(5) |

where is the fragmentation function (FF) for , , mesons, with . It satisfies ; the branching ratio for the transition from to , in turn, satisfies with denoting all heavy flavor hadrons. We will employ the FFs introduced in Braaten:1994bz (); details of the implementation are discussed in Appendix B.

The color singlet () channel contribution of production cross-section in the CCG+NRQCD framework can be expressed as Ma:2014mri ()

(6) |

and the color octet (CO) intermediate states written as

(7) |

The hard matrix elements and are given in Appendix A. Note that in and where . Since Eq. (6) has a cubic dependence on , while Eq. (7) has only a quadratic dependence, it is evident that the short distance CS and CO cross-sections have different dependencies on the dynamics of saturated gluons in protons and nuclei.

We will compare the NRQCD results employing the above expressions with the cross-section computed in the Improved Color Evaporation Model (ICEM) Ma:2016exq (). The differential cross section for production in the CGC+ICEM framework is given by

(8) |

where with and . Here is the invariant mass of the . and are respectively the relative momentum and angle between and in the pair rest frame. represents the nonperturbative transition probability from the pair to the meson. The principal difference between the ICEM and the conventional CEM Fritzsch:1977ay (); Gluck:1977zm (); Barger:1979js () is that the ’s transverse momentum differs from the pair’s transverse momentum : . In our computations, we will use and .

Results for -meson and production: With the expressions in Eqs. (5), (6), (7) and (8), we can simultaneously study -meson and production with increasing event activity, as represented by the inclusive charged hadron multiplicity. The latter is computed in a factorized approximation to the CGC EFT Tribedy:2010ab (); Tribedy:2011aa (); Albacete:2012xq () as shown in Appendix C. The dynamical ingredients in this computation, as in the case of heavy-flavor production, are the UGDs in the projectile and the target. Therefore fixing these from single inclusive production provides significant predictive power. In Appendix D, we outline the ingredients in computing the energy evolution of the UGDs. In Appendix E, we present numerical results for the charged hadron multiplicity. As shown there, these initial scales are well constrained by data on the average transverse momentum of charged hadrons versus .

With the and thereby constrained, the corresponding UGDs can be used to compute the the isospin averaged -meson cross-section. Figure 1 compares our model prediction to the LHC data on rare events in both and collisions, normalized to the minimum bias value, versus likewise normalized to its minimum bias value. As is clear from Eqs. (2)–(5), the ratio plotted on the y-axis is fairly insensitive to uncertainties arising from choice of fragmentation functions, proton and nuclear size, and the coupling constant . Likewise, the ratio on the -axis minimizes nonperturbative uncertainties from geometry effects in both protons and nuclei. The agreement with data at TeV shown in Figure 1 (a) is remarkably good for both windows. The experimental error bars are however large for the rarest events. Figure 1 (b) shows the model comparison to LHC data at TeV/nucleon. While model results overshoot data for GeV, the same qualitative trend is seen as in data. Note that because one varies both and (for values given in the figure caption) there is room for finetuning; we have not attempted to do so. Appendix F also shows that -meson distributions for minimum bias events are well reproduced out to GeV in both and collisions.

The same UGDs can be used to compute production in rare events. Remarkably, the relative contribution of for each set of quantum numbers changes with increasing event activity. The relative contribution of different channels plotted in Fig. 2 shows that is larger than the other channels for all , and it increases significantly with increasingly rare events. Most importantly, the LHC data are in between curves of and other channels. This may imply that production with high event activity is dominated by the channel but has only small contributions from and channels. This is consistent with the universality requirement extracted from BELLE data. The above discussion also shows that production in rare events can provide useful constraints for LDMEs.

Because of dominance, we can use the simpler ICEM model to describe production. This is reasonable because in the ICEM production is dominated by gluon fragmentation, and thus mainly produces intermediate states. In Fig. 2 we indeed find that ICEM curve is also in between curves of and other channels. We will use the CGC+ICEM model in the following. Figure 3 (a) shows that the data on ratios of the cross-section in collisions is –independent. In the CGC, as seen previously for ridge yields Dusling:2015rja (), this is natural because the energy dependence of cross-sections is controlled by , which also governs the charged hadron multiplicity; events at different energies with the same are therefore identical. Figure 3 (a) also provides a prediction that RHIC data at TeV will conform to this expectation. In Fig. 3 (b), we show the comparison of the CGC+ICEM model to data on ratio in collisions. The agreement with data is quite good; since many of the nonperturbative uncertainties cancel in the ratios, this agreement with the and data demonstrates that the CGC EFT captures key features of the short distance cross-sections. It also demonstrates the reasonability of dominance, which agrees with the universality requirement.

A systematic theory uncertainty is that the dilute-dense approximation to CGC EFT we employ is valid only when . The full “dense-dense” EFT computation is cumbersome and beyond the scope of present computations. This systematic uncertainty is reduced at forward rapidities in collisions and at both central and forward rapidities in collisions, which however necessitates introducing nonperturbative scales for both projectile and target. The ratios considered will mitigate these uncertainties; further, the requirement that we reproduce charged particle multiplicities is a powerful constraint. Our results for the ratios at forward rapidities are presented in Appendix F. Within the uncertainties noted, we find good agreement with data. Conversely, the and -meson data can be employed to constrain variations in values with increasing event activity. The model, with the parameters thus fixed, can for example be compared to data on -hadron correlations at the LHC Acharya:2017tfn ().

Summary: We outlined in this letter the potential of event engineered heavy flavor measurements to uncover the dynamics of rare parton configurations at collider energies. Our studies within the framework of the CGC+NRQCD EFT, suggest that the short distance dynamics in such events requires saturation scales that are an order of magnitude greater than those in minimum bias events. On the one hand, these harder scales suggest that the weak coupling CGC framework is more reliable for rare events. On the other, the treatment of multiplicity biased configurations is significantly more complex than computations developed to study minimum bias configurations and demands further theoretical development. For instance, analyses of the relative contribution of initial versus final state effects, for the most part, do not incorporate these ideas.

Within the systematic uncertainties outlined, the event engineering studies show that the CGC does an excellent job of describing the short distance physics in the production of both open flavor -mesons and the onium state. This is particularly compelling because the nonperturbative scales are fixed by comparisons to data on charged particle multiplicities in rare events and several nonperturbative parameters cancel out in the ratios. Our work also illustrates the potential of event engineering to distinguish between intermediate states with differing quantum numbers that contribute to the hadronization of mesons. We find that with increasing event activity the state dominates production, consistent with an interpretation of the dominance of hard gluon fragmentation in hadronization. Remarkably, this is also consistent with the universality requirement from BELLE data. An interesting possibility to explore in future is to use event engineering to extract the LDMEs using the CGC+NRQCD framework at low to moderate values of ; this is especially promising now that that next-to-leading order CGC computations are becoming available. For the reasons outlined, these developments have the potential to resolve puzzles in quarkonium production and will be addressed in follow-up work.

Acknowledgements:

KW is supported by Jefferson Science Associates, LLC under U.S. DOE Contract No. DE-AC05-06OR23177 and U.S. DOE Grant No. DE-FG02-97ER41028. P.T. and R.V. are supported by the U. S. Department of Energy Office of Science, Office of Nuclear Physics, under contracts No. DE-SC0012704. This work is part of and supported by the DFG Collaborative Research Centre “SFB 1225 (ISOQUANT)”. We would like to thank Anton Andronic, Juergen Berges, Peter Braun-Munzinger, Rongrong Ma, Silvia Masciocchi, Jianwei Qiu, Lijuan Ruan, and Johanna Stachel for useful comments. R.V. would in particular like to thank Adrian Dumitru and Vladimir Skokov for sharing their deep insights on a theoretical framework for multiplicity biased events.

## Appendix A Appendix A: Hard Matrix Elements

### a.1 Hard Matrix element in production

### a.2 Nrqcd

For the color singlet channel, reads Ma:2014mri ()

(12) |

where , , with and . Note here due to momentum conservation at LO.

For the color octet channels, reads Kang:2013hta ()

(13) | ||||

(14) | ||||

(15) |

## Appendix B Appendix B: -meson fragmentation functions

The Braaten, Cheung, Fleming, and Yuan (BCFY) Braaten:1994bz () heavy-quark fragmentation functions (FF) provide different -distributions for pseudoscalar mesons and vector mesons. Following Refs. Cacciari:2012ny (); Cacciari:2003zu (), we will set the different FF for , , and production to be

(16) | ||||

(17) | ||||

(18) |

where the original BCFY FFs are given by Braaten:1994bz ()

(19) | ||||

(20) |

is determined analytically from . Here describes production involving the effect of the decay into , and reads

(21) |

We shall fix GeV and GeV. is a single nonperturbative parameter and can be interpreted as the ratio of the constituent mass of the light quark to the mass of the heavy meson like . One can easily estimate .

## Appendix C Appendix C: Inclusive hadron production

We will review here charged hadron production in and collisions in the CGC framework Tribedy:2010ab (); Tribedy:2011aa (); Albacete:2012xq (). The differential cross section for inclusive gluon production in collisions () in the -factorization formula at LO Kovchegov:2001sc (); Blaizot:2004wu () is given by

(22) |

where . Now in and , one should read . The impact parameter dependence is encoded in the saturation scale of the proton and nucleus for simplicity. is a normalization factor which takes account of information about a transverse area for overlap region between the projectile proton and the target nucleus. However, throughout this paper, we leave it an arbitrary constant, since we shall consider the ratio of the hadron multiplicity in rare events to that in minimum bias events.

For inclusive hadron production at finite transverse momentum, a light hadron FF () is involved with the gluon production cross section, as usual. However, it is unclear whether fragmentation function is applicable to low hadron production. Nevertheless, we shall take into account gluon fragmentation function because such a fragmenting process can play a significant role to provide us with reliable predictive power to describe data of charged hadron production. We shall go though this further below.

In our numerical computations, we employ which corresponds to the LO parametrization of the KKP FF Kniehl:2000fe (). Now charged hadron multiplicity at pseudorapidity can be written as

(23) |

where is the Jacobian for transforming the expression in -space to that in -space. We have assumed that and defined for simplicity. is an inelastic cross section in collisions. We will put a cut off GeV and GeV in Eq. (23) in our numerical calculations. is determined from the kinematical condition, . The rapidity in is replaced with

(24) |

where we assumed that hadron’s transverse momentum is strongly correlated with the gluon’s transverse momentum so that we use in the Jacobian and Eq. (24). With regard to the mass scale of the charged hadron, we fix as 300 MeV. One must keep in mind that the rapidity of the produced gluon is shifted by as in Eq. (24) to perform numerical calculations in collisions at the LHC.

### c.1 Appendix D: Small- evolution

The rapidity or energy dependence of the dipole amplitude, to leading accuracy in , is given by the nonlinear Balitsky-Kovchegov (BK) equation Balitsky:1995ub (); Kovchegov:1996ty ():

(25) |

where the running coupling evolution kernel in Balitsky’s prescription Balitsky:2006wa () is given by

(26) |

with being the size of the parent dipole size prior to one step in evolution. The one loop coupling constant in coordinate space is employed to solve the rcBK equation. We can use the initial dipole amplitude at or to be of the form given by the McLerran-Venugopalan (MV) model McLerran:1993ni (); McLerran:1993ka ():

(27) |

where is an anomalous dimension, is the initial saturation scale in the proton at . The infrared cutoff is chosen by freezing . For the initial input parameters in the rcBK equation, we set , , , , and . These parameters in this initial condition are obtained from global data fitting at HERA-DIS and given in Ref. Albacete:2012xq (); rcbk (). For the target nucleus, where for minimum bias events in collisions is obtained from fitting the New Muon Collaboration data on the nuclear structure functions Dusling:2009ni (). For the purpose of our discussion, we shall fix simply for heavy nuclei such as Pb and Au in our numerical calculations. Indeed, several previous studies Ma:2015sia (); Fujii:2015lld (); Fujii:2017rqa (); Ma:2017rsu () adopting the smaller value of succeeded in describing nuclear modification factor of and meson at RHIC and the LHC.

At large values of , we need to extrapolate the parametrization of the dipole amplitude to these values. In Refs. Ma:2014mri (); Ma:2017rsu (), the adjoint dipole distribution in Eq. (3) at is determined to be where the coefficient can be determined by matching the UGDF to collinear gluon distribution function. However, it is unclear whether the above matching procedure is applicable to high multiplicity events. In lieu, at large , we adopt the simple extrapolation ansatz for (3) Gelis:2006tb ():

(28) |

We also apply the same procedure on the target side.

### c.2 Appendix E: Numerical results for inclusive hadron production

We first clarify our setup for numerical calculations in this paper. Assuming the CGC framework is yet applicable to collisions at collider energies, the only deference between collisions and collisions is the initial saturation scale for the target modulo the geometrical transverse size of the target. Regarding input parameters, we do not set , , and to specific values here and leave these factors arbitrary in our numerical computations, since those parameters are irrelevant to the relative yield of . With regard to strong coupling constant in Eqs. (2)(6)(7)(22), we fix it as a constant value like because all the differential cross sections in this paper have been derived at leading order in .

Figure 4 shows relative in collisions at the LHC at mid rapidity by varying the initial saturation scale . We take the saturation scales of the projectile proton and the target proton to be symmetrical; . The averaged is obtained by setting with . The solid line is the result obtained by using the KKP FF, while the dashed lines correspond to the result without using the KKP FF. It is clear that the relative grows almost linearly as increases when the KKP FF is used.

The computation of the multiplicity in collisions is generally more complicated because it depends on the combination of the saturation scale of the projectile proton and that of the target nucleus. In Fig. 4 (b), several combinations of and are depicted in different lines. We set the averaged in collisions as the result with and . In contrast to collisions, the relative in collisions does not show a rapid growth with increasing and , even if we employ the KKP FF.

The mean transverse momentum of hadrons produced in high multiplicity events in and collisions is an important observable to check whether the CGC framework describes bulk data. The definition of is given by

(29) |

Figure 5 shows dependence of