Event engineering studies for heavy flavor production and hadronization in high multiplicity hadron-hadron and hadron-nucleus collisions

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 $D$-meson cross-sections; for the $J/\psi$, 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 in this CGC+NRQCD framework reveals a very rapid growth in the fragmentation of the $^3S_1^{[8]}$ state in rare events relative to minimum bias events.


I. INTRODUCTION
The study of high multiplicity events in proton-proton (p + p) and proton-nucleus (p + A) 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 [1,2], or hydrodynamic origins [3,4].
Heavy flavor measurements add important elements to the discussion because the large quark masses provide a semi-hard scale to probe initial state dynamics. A compelling example of event engineered heavy flavor measurements in p + p and p + A 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) [5] framework, wherein the inclusive differential crosssection of a heavy quarkonium state Q in p + p and p + A collisions is expressed as where κ = 2S+1 L [c] J are quantum numbers of the produced intermediate heavy quark pair, with S, L and J denoting its spin, orbital, and total angular momenta, respectively. The symbol c denotes a color singlet (CS, c = 1) or color octet (CO, c = 8) state. The dσ κ are perturbative short distance coefficients for heavy quark pair production with quantum numbers κ and O Q κ 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, an important puzzle is that the magnitude of the linear combination of the 1 S [8] 0 and 3 P [8] 0 LDMEs extracted from hadroproduction data [6,7] is larger than an upper bound set by BELLE e + e − data [8]. While this apparent breaking of universality may bring into question NRQCD factorization, we will show that event engineering offers a possible resolution to this puzzle.
In this work, we will show that the systematics of heavy flavor production in rare events in p+p and p+A collisions are sensitive to strongly correlated gluons in the colliding protons and nuclei. The dynamics of such configurations is controlled by an emergent semi-hard saturation scale Q s (x) in each of the colliding hadrons, where x is the longitudinal momentum fraction carried by a parton in the hadron [9,10]. Since Q s (x) grows with decreasing x, and increasing nuclear size, the interplay of the dynamics of hard and soft modes evolves with the changing energy arXiv:1803.11093v2 [hep-ph] 14 Oct 2018 and centrality of the collision.
A systematic framework to study gluon saturation is the Color Glass Condensate (CGC) effective field theory (EFT) [11][12][13]. The cross-sections for the production of heavy quarkonia in the CGC EFT for hadronhadron collisions were computed over a decade ago [14][15][16][17][18]. A more recent development is the CGC+NRQCD framework 1 [20], the novel element being that dσ κ 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 p + p and p + A collisions [21][22][23][24][25][26][27][28][29][30]. High multiplicity configurations are approximated by increasing the value of Q s (x) at the input large x scale in both protons and nuclei in multiples of Q 2 0 = 0.168 GeV 2 , the initial saturation scale at x = 0.01, determined from fits to the minimum bias e + p DIS data [31]. As also implemented in studies of ridge yields [32][33][34], increasing the saturation scale in this manner captures the fluctuations of protons and nuclei into larger numbers of color charges in rare events. More systematic treatments of high multiplicity "biased" color charge configurations are under development [35][36][37].
We will focus here 2 on measurements of D and J/ψ mesons in high multiplicity p + p and p + A collisions [42][43][44][45][46][47][48][49]. The striking feature of the data is that the production yields of D and J/ψ in high multiplicity events are significantly enhanced relative to minimum bias events. Interestingly, in p + p collisions, such growth is observed to be independent of collision energy. The models proposed to explain their systematics include percolation models [50,51], dipole models [52] and multiparton interaction models [53]. All these models approximate effects contained in the CGC EFT. Gluon saturation is included in the EPOS3 model [54], 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 and help resolve extant heavy flavor puzzles in collider experiments.

II. OPEN FLAVOR AND QUARKONIUM PRODUCTION
We first consider the spin and color averaged inclusive cross-section p + A(p) → c(p c ) +c(qc) + X, which can be expressed in the CGC EFT as [15] dσ cc 1 See [19] for a specialized discussion. 2 The Υ and open bottom computations require Sudakov resummation [38][39][40][41] and are beyond our scope here.
with p c⊥ (qc ⊥ ) and y c (yc), the transverse momentum and rapidity respectively of the produced charm (anti-charm) quarks. Further, y p = ln(1/x 1 ) and Y = ln(1/x 2 ), where x 1,2 = ( m 2 c + p 2 c⊥ e ±yc + m 2 c + q 2 c⊥ e ±yc )/ √ s, denote the 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 ϕ p,yp (k ⊥ ) is defined as [22] Here πR 2 p (πR 2 A ) is the transverse area occupied by gluons in the proton (nucleus) and where is the fundamental Wilson line in the amplitude (complex conjugate amplitude) representing multiple scattering of the quark with background fields at the position r ⊥ (0 ⊥ ). Note that · · · y here 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 D meson production is then given by where D c→D (z) is the fragmentation function (FF) for D 0 , D + , D * + mesons, with z = p D⊥ /p c⊥ . It satisfies dzD c→D (z) = Br(c → D); the branching ratio Br(c → D) for the transition from c to D, in turn, satisfies X Br(c → X) = 1 with X denoting all heavy flavor hadrons. We will employ here the BCFY [55] and KKKS [56] FFs; key details are discussed in Appendix B.
The color singlet (κ = 3 S 1 ) channel contribution of J/ψ production cross-section in the CGC+NRQCD framework can be expressed as [22] dσ κ cc,CS The hard matrix elements G κ 1 and Γ κ 8 are given in Appendix A. Note that x 1,2 = (2m c ) 2 + p 2 ⊥ e ±y / √ s in y p and Y where m c = m J/ψ /2. Since Eq. (6) has a cubic dependence on N Y , 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 J/ψ cross-section computed in the Improved Color Evaporation Model (ICEM) [57]. The differential cross section for J/ψ production in the CGC+ICEM framework is given by where J =q M 2 + p 2 ⊥ / [M ω c ωc| sinh(y c − yc)|] with ω c = m 2 c + p 2 c⊥ and ωc = m 2 c + q 2 c⊥ . Here M is the invariant mass of the cc.q and φ are respectively the relative momentum and angle between c andc in the cc pair rest frame [18]. F J/ψ represents the nonperturbative transition probability from the cc pair to the J/ψ meson. The principal difference between the ICEM and the conventional CEM [58][59][60] is that the J/ψ's transverse momentum differs from the pair's transverse momentum p ⊥ : p ⊥ = (m J/ψ /M )p ⊥ . In our computations, we will use m J/ψ = 3.1 GeV and 2 m D = 3.728 GeV.

III. RESULTS FOR D-MESON AND J/ψ PRODUCTION
With the expressions in Eqs. (5), (6), (7) and (8), we can simultaneously study D-meson and J/ψ production with increasing event activity, as represented by the inclusive charged hadron multiplicity. The latter is computed in a k ⊥ factorized approximation to the CGC EFT [31,61,62] as shown in Appendix C. The dynamical ingredients in all the computations are the UGDs in the projectile and the target. Therefore fixing these, and their energy evolution (see Appendix D) from single inclusive production provides significant predictive power.
In Appendix E, we present numerical results for the charged hadron multiplicity. As shown there, these initial scales Q 2 sp,0 (Q 2 sA,0 ) at x = 0.01 for protons (nuclei) that enter into the UGDs are well constrained by the data on p ⊥ versus dN ch /dη of charged hadrons. For the event engineering studies, the UGDs are obtained by varying Q 2 sp,0 (Q 2 sA,0 ) within a range of 1-3 (4-12) times their corresponding minimum bias values (Q 2 0 = 0.168 GeV 2 ). With Q sp,0 and Q sA,0 thereby constrained, the UGDs can be used to compute the isospin averaged D-meson cross-section. Figure 1 compares our model prediction to the mid-rapidity LHC high multiplicity data in both p + p and p + A collisions, normalized to the minimum bias value, versus dN ch /dη 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 α s . Likewise, the ratio on the x-axis minimizes nonperturbative uncertainties from geometry effects in both protons and nuclei. The agreement with p + p data at √ s = 7 TeV shown in Fig. 1 (a) is remarkably good for both p ⊥ windows. The experimental error bars are however large for the rarest events. Figure 1 (b) shows the model comparison to LHC p + A data at √ s = 5.02 TeV/nucleon.
While model agreement with data in the 1 < p ⊥ < 2 GeV window is quite good, it overshoots data for 2 < p ⊥ < 4 GeV though it has the same qualitative trend. Because one varies both Q sp,0 and Q sA,0 , there is room for finetuning. Appendix F shows that D-meson p ⊥ distributions for minimum bias events are well reproduced out to p ⊥ ∼ 5 GeV in both p + p and p + A collisions.
The very same UGDs are used to compute J/ψ production. Remarkably, the relative contribution of dσ κ for each κ changes with increasing event activity. Figure 2 shows that relative yield of 3 S [8] 1 is larger than the other channels for all dN ch /dη, and it increases significantly with increasingly rare events. This implies that a very rapid growth in J/ψ production in rare events in the 3 S [8] 1 channel relative to minimum bias. The growth in the contributions of the 1 S [8] 0 and 3 P [8] J channels are relatively much smaller. This enhanced contribution of the short distance contributions in the 3 S [8] 1 channel suggests the LDMEs of the 1 S [8] 0 and 3 P [8] J channels could potentially be smaller. This may provide a way forward in reconciling the LDMEs extracted from hadroproduction with the universality requirement extracted from BELLE e + e − data, hence providing a possible resolution of the NRQCD puzzle mentioned previously.
The relative large 3 S contribution suggests that the simpler ICEM model where, gluon fragmentation through this channel dominates, may be sufficient to describe J/ψ production and we will do so in the following. In future, we will study rare events directly in the CGC+NRQCD framework. Figure 3 (a) shows that the data on ratios of the J/ψ cross-section in p + p collisions is √ s-independent. In the CGC, as seen previously for ridge yields [33], the energy dependence of cross-sections is controlled by Q s (x), which also governs the charged hadron multiplicity; events at different energies with the same Q s are therefore identical. Figure 3 (a) predicts that RHIC p+p data at √ s = 0.5 TeV will conform to this expectation. In Fig. 3 (b), we compare the CGC+ICEM model to data in p + A collisions. Since many nonperturbative uncertainties cancel in these ratios, the agreement with both p + p and p + A data demonstrates that the CGC EFT captures key features of the short distance cross-sections.

IV. SUMMARY
We outlined the potential of event engineered heavy flavor measurements to uncover the dynamics of rare parton configurations at collider energies. Our CGC EFT studies 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 hand, the treatment of rare multiplicity biased configurations is significantly more complex than computations developed to study minimum bias configurations and demands further theoretical development.
Our work further illustrates the potential of event engineering to distinguish between intermediate states with differing quantum numbers that contribute to the hadronization of quarkonia. The finding that the hadronization contribution of the 3 S 1 state to J/ψ production grows rapidly suggests the growing importance of hard gluon fragmentation in J/ψ hadronization. As noted, this result may provide an important clue in resolving the universality requirements on LDMEs from BELLE e + e − data, thereby possibly resolving a puzzle between the magnitudes of the LDMEs extracted from hadron collision data relative to e + e − data.
A systematic theoretical uncertainty is that the dilutedense approximation to CGC EFT we employ is valid only when Q s,proj. /k ⊥,proj. < Q s,target /k ⊥,target . The full "dense-dense" EFT computation is beyond the scope of present computations; these are beginning to be quantified [63]. This systematic uncertainty is reduced at forward rapidities in p + p collisions and at both central and forward rapidities in p + A collisions. The ratios considered mitigate these uncertainties; further, the requirement that we reproduce charged particle multiplicities is a powerful constraint. Our results for the J/ψ ratios at forward rapidities are presented in Appendix F. Within the uncertainties noted, we find good agreement with data. The model, with the parameters thus fixed, can for example be compared to data on J/ψ-hadron correlations at the LHC [64]. Finally, a source of systematic uncertainty in our computation we have not discussed is the possible role of higher twist fragmentation contributions at low p ⊥ . The short distance hard matrix elements ensure any such contribution is suppressed by α s (m Q ). Such higher order contributions, as well as other α s suppressed contributions to the matrix elements are not included in our treatment. Our framework however can be systematically improved in future to include such effects. Ξ = Ξ qq,qq + Ξ qq,g + Ξ g,g , where In the above, The Lipatov vertex C µ that appears here, can be written in component form as

NRQCD
For the color singlet 3 S 1 channel, G 1 reads [22] For the color octet channels, Γ κ 8 reads [20] Appendix B: D-meson fragmentation functions We will discuss here heavy-quark fragmentation functions (FF) that provide different z-distributions for pseudoscalar mesons and vector mesons. We consider specifically the Braaten-Cheung-Fleming-Yuan (BCFY) FF [55] and the Kneesch-Kniehl-Kramer-Schienbein (KKKS) FF [56]. Consider the BCFY FF first, following Refs. [67,68], we will set the different FF for D 0 , D + , and D * production to be where the original BCFY FFs are given by [55] D (P ) N is determined analytically from  Figure 4 also displays the DGLAP evolution of the BCFY FFs by setting (B1)-(B3) as initial conditions and evolving µ from 1.5 GeV to 10.5 GeV. Clearly, the DGLAP evolution significantly modifies the initial BCFY FFs.
Turning now to the KKKS FF, in the KKKS set 3 , 3 Numerical points of the KKKS FF as well as other FF set are available online thanks to [69]. the µ dependence of the FFs for D-mesons was again taken into account through DGLAP evolution. As to initial conditions, the functional form D c→D (z, µ 0 ) = N z −(1+γ 2 ) (1 − z) a e −γ 2 /z is set at µ 0 = 1.5 GeV. All the input parameters N , a, γ are determined by global fitting of all available e + e − data. In Fig. 4, the KKKS FFs at µ = 10.5 GeV are compared to the BCFY FFs together with CLEO e + e − data [65]. The data comparisons obviously prefer the KKKS FFs to describe e + e − data, although one must keep in mind that the data are normalized cross-sections for D-mesons production, not heavy quark FFs themselves. Indeed, both the BCFY FF with the DGLAP evolution and the KKKS FF overshoot the data points at lower z because we do not convolute hard scattering part with the FFs here for simplicity. If we take into account hard scattering part correctly, the KKKS FFs should agree with the data [56].

Appendix C: Inclusive hadron production
We will review here charged hadron production in p+p and p + A collisions in the CGC framework [31,61,62]. The differential cross-section for inclusive gluon production in p + A collisions (p + A → g(p g ) + X) in the k ⊥ -factorization formula at LO [70,71] is given by where k ⊥ ≤ p g⊥ . Now in y p and Y , one should read x 1,2 = p g⊥ e ±y / √ s. The impact parameter dependence is encoded in the saturation scale of the proton and nucleus for simplicity.K b 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 (D h ) is involved with the gluon production cross-section, as usual. However, it is unclear whether fragmentation function is applicable to low p ⊥ 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 D h (z) = 6.05z −0.714 (1 − z) 2.92 which corresponds to the NLO parametrization of the Kniehl-Kramer-Potter (KKP) FF for charged hadron production at µ = 2 GeV [72]. Now charged hadron multiplicity at pseudorapidity η can be written as where J y→η = p g⊥ cosh η/ p 2 g⊥ cosh 2 η + m 2 h is the Jacobian for transforming the expression in y-space to that in η-space. We have assumed that y = y h = y g and defined p ⊥ ≡ zp g⊥ for simplicity. σ inel is an inelastic cross-section in p + A collisions. We will put a cut off p max = 10 GeV and p min = 0.1 GeV in Eq. (C2) in our numerical calculations. z min is determined from the kinematical condition, x 1,2 ≤ 1. The rapidity in dσ g /d 2 p g⊥ dy is replaced with where we assumed that hadron's transverse momentum is strongly correlated with the gluon's transverse momentum p g⊥ so that we use p g⊥ in the Jacobian and Eq. (C3).
With regard to the mass scale of the charged hadron, we fix m h as 300 MeV. One must keep in mind that the rapidity of the produced gluon is shifted by ∆y = 0.465 as y → y−∆y in Eq. (C3) to perform numerical calculations in p + A collisions at the LHC.

Appendix D: Small-x evolution
The rapidity or energy dependence of the dipole amplitude, to leading accuracy in N c , is given by the nonlinear Balitsky-Kovchegov (BK) equation [73,74]: where the running coupling evolution kernel in Balitsky's prescription [75] is given by with r ⊥ = r 1⊥ + r 2⊥ being the size of the parent dipole size prior to one step in Y evolution. The one loop coupling constant in coordinate space α s (r 2 ⊥ ) = 1/ 9 4π ln 4C 2 r 2 ⊥ Λ 2 +â is employed to solve the rcBK equation. We can use the initial dipole amplitude at x = x 0 = 0.01 or Y 0 = ln 1/x 0 to be of the form given by the McLerran-Venugopalan (MV) model [76,77]: where γ is an anomalous dimension, Q sp,0 is the initial saturation scale in the proton at x = x 0 . The infrared cutoffâ is chosen by freezing α s (r → ∞) ≡ α fr . For the initial input parameters in the rcBK equation, we set Q 2 sp,0 = 0.168 GeV 2 , γ = 1.119, C = 2.47, Λ = 0.241 GeV, and α fr = 1.0. These parameters in this initial condition are obtained from global data fitting at HERA-DIS and given in Ref. [62,78]. For the target nucleus, Q 2 sA,0 = cA 1/3 Q 2 sp,0 where c 0.5 for minimum bias events in p + A collisions is obtained from fitting the New Muon Collaboration data on the nuclear structure functions F 2,A (x, Q 2 ) [79]. For the purpose of our discussion, we shall fix simply Q 2 sA,0 = 2 Q 2 sp,0 for heavy nuclei such as Pb and Au in our numerical calculations. Indeed, several previous studies [23,27,29,30] adopting the smaller value of Q 2 sA,0 succeeded in describing nuclear modification factor of J/ψ and D meson at RHIC and the LHC.
At large values of x ≥ x 0 = 0.01, we need to extrapolate the parametrization of the dipole amplitude to these x values. In Refs. [22,30], the adjoint dipole distribution in Eq.
where the coefficient a(x) 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 x ≥ x 0 , we adopt the simple extrapolation ansatz for (3) [81]: We also apply the same procedure on the target side.

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 p + p collisions at collider energies, the only deference between p + p collisions and p + A collisions is the initial saturation scale for the target modulo the geometrical transverse size of the target. Regarding input parameters, we do not setK b ,K ch , and σ inel to specific values here and leave these factors arbitrary in our numerical computations, since those parameters are irrelevant to the relative yield of N ch . With regard to strong coupling constant α s in Eqs. (2)(6)(7) and Eq. (C1), we fix it as a constant value like α s ∼ 0.2 because all the differential cross-sections in this paper have been derived at leading order in α s . Figure 5 shows relative dN ch /dη in p + p collisions at the LHC at mid rapidity by varying the initial saturation scale Q 2 sp,0 . We take the saturation scales of the projectile proton and the target proton to be symmetrical; Q 2 sp1,0 = Q 2 sp2,0 . The averaged N ch is obtained by setting Q 2 sp1,0 = Q 2 sp2,0 = Q 2 0 with Q 2 0 = 0.168 GeV 2 . 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 N ch grows almost linearly as Q 2 sp,0 increases when the KKP FF is used. The computation of the multiplicity in p + A 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. 5 (b), several combinations of Q 2 sp,0 and Q 2 sA,0 are depicted in different lines. We set the averaged N ch in p + A collisions as the result with Q 2 sp,0 = Q 2 0 and Q 2 sA,0 = 2Q 2 0 . In contrast to p + p collisions, the relative N ch in p + A collisions does not show a rapid growth with increasing Q 2 sp,0 and Q 2 sA,0 , even if we employ the KKP FF. The mean transverse momentum p ⊥ of hadrons produced in high multiplicity events in p + p and p + A collisions is an important observable to check whether the CGC framework describes bulk data. The definition of p ⊥ is given by (E1) Figure 6 shows N ch dependence of p ⊥ for single hadron production in p + p and p + A collisions at the LHC. We fix normalization of dN ch /dη in p + p and p + A collisions to fit minimum bias data respectively. Using the KKP FF, one can obtain a reasonable description of the data in p + p collisions at the LHC. In p + A collisions, numerical results with larger saturation scales for the projectile proton and the target nucleus show a nice agreement with data at the highest multiplicity. These comparisons clearly substantiate the robustness of the CGC framework in describing bulk data.   m c = 1.5 ≈ m J/ψ /2 is used in Eqs. (6) and (7). As noted in [22], some of the dependence of on the quark masses in the short-distance cross-sections, is canceled out by the dependence of the LDMEs on quark mass. In Fig. 7 (a), differential cross-sections for D 0 , D + , D * production in minimum bias p + p collisions at the LHC are shown. As showed in Fig. 4, the KKKS FFs agree quite well with e + e − data relative to the BCFY FFs even after DGLAP evolution is taken into consideration. However, both these FF sets are in agreement with data on D meson production in p + p collisions for p ⊥ > 1 GeV. Specifically, for the region in p ⊥ of interest, from 1 GeV to 4 GeV, the BCFY curves and the KKKS curves are indistinguishable. Indeed, for the double ratio of minimum bias result to high multiplicity result, it makes little difference for our results. We, of course, anticipate that better data at high multiplicity can help us to confirm whether the tension with e + e − data for the BCFY FFs is also seen in hadron-hadron collisions.
In Fig. 7 (a), K-factor of 2.5 is required to describe data if we set the effective transverse area as R p = 0.6 fm. However, a smaller value of R p can be also taken and is compatible with matching of unintegrated gluon distributions to gluon collinear PDFs at x = 0.01 [22]. A smaller   Ref. [85]. Data at forward rapidity are from [47].
transverse area can therefore bring 50% uncertainties to K since higher order NLO effects cannot be distinguished from uncertainties in the transverse area. Indeed, this is a strong motivation for considering double ratios as we do, because the K-factor cancels out in the ratio.
In p + A collisions, we determine the effective transverse area of the target nucleus R A by imposing that nuclear modification factor R pA = dσ pA /(Adσ pp ) for cc production should approach unity at asymptotically high p ⊥ . This condition leads to R A = A/N γ R p with N = Q 2 sA,0 /Q 2 sp,0 . Now the initial condition for Q 2 sA,0 = 2Q 2 sp,0 with γ = 1.119 for the rcBK equation gives R A = 9.79R p . Using this value of R A with the same K-factor, Fig. 7 (b) shows a nice agreement with data in minimum bias p + A collisions. Figure 8 (a) shows that for minimum bias Q 2 sp,0 = Q 2 0 , the relative contributions of dσ κ /dp ⊥ for κ = 3 S  J are similar to that of the ICEM at low p ⊥ and differs from the 3 S [8] 1 . In contrast, the p ⊥ distribution of the latter is harder than the other channels at large p ⊥ , a trend similar to that of the ICEM. This is understand-able because high p ⊥ J/ψ are likely to be produced via gluon fragmentation with the quantum numbers of the 3 S [8] 1 channel. In contrast, Fig. 8 (b) shows that for rare Q 2 sp,0 = 5 Q 2 0 configurations, the normalized cc differential cross-section for the 3 S [8] 1 channel is close to that of the ICEM over the entire p ⊥ range. The other channels are relatively harder at low p ⊥ and softer at higher p ⊥ .
We show in Fig. 9 comparisons of the ICEM with data on N ch dependence of J/ψ production in p + p and p + A collisions at the LHC at forward rapidity. In contrast to mid rapidity, at forward rapidity, the symmetrical treatment; Q 2 sp1,0 = Q 2 sp2,0 overshoots the data slightly in p+p collisions. Data point at dN ch / dN ch ∼ 4 seems to favor the asymmetrical treatment; Q 2 sp1,0 < Q 2 sp2,0 . This is consistent with a naive expectation that a phase space for the gluon distribution of the projectile proton can shrink at forward rapidity (x 1 ∼ O(1)) where a dilute-dense approximation is robust. One can find the similar trend for forward J/ψ production in p + A collisions.
Predictions for mean transverse momentum of J/ψ production in p + A collisions at the LHC are given in Fig. 10. At mid rapidity, only the minimum bias data is available. The CGC prediction shows that p ⊥ of J/ψ depends on the change of the Q 2 sA,0 largely but does not change rapidly as N ch increases. On the other hand, at forward rapidity, our numerical results overestimate J/ψ's p ⊥ at high N ch . The comparable results for p ⊥ of average D (D 0 , D + , D ) production in p + A collisions at mid-rapidity using the BCFY FFs with r = 0.1 is shown in Fig. 11, showing a relatively flat dependence on event activity compared to the J/ψ.