Multiplicity dependence of $\chi_{c}$ and $\chi_{b}$ meson production

In this paper we analyze in detail the production of the $\chi_{c}$ and $\chi_{b}$ mesons in $pp$ collisions. Using the color dipole framework, we estimated the cross-sections in the kinematics of ongoing and forthcoming experiments, and found that our estimates are in reasonable agreement with currently available experimental data. We also analyzed the dependence on multiplicity of co-produced hadrons and found that it is significantly milder than that of $S$-wave quarkonia. We expect that the experimental confirmation of this result could constitute an important test of our understanding of multiplicity enhancement mechanisms in the production of different quarkonia states.


I. INTRODUCTION
The standard approach to the description of mesons containing heavy quarks is based on pomeron-pomeron fusion. The heavy quarks formed at the partonic level might hadronize into open heavy flavor mesons (e.g. D-or B-mesons) or form quarkonia states at later stages of the collision [1][2][3][4]. This framework provides reasonable estimates for the total and differential cross-sections, although it includes some uncertainties due to a rather limited knowledge of the fragmentation functions of open-flavor mesons or the Long Distance Matrix Elements (LDMEs) of quarkonia states [5][6][7][8][9][10][11]. Moreover, this approach recently encountered difficulties with the description of recent experimental data on multiplicity dependence of co-produced charged hadrons [12][13][14][15][16][17]. In fact, it was discovered by both the STAR and ALICE collaborations that the relative yields of the 1S quarkonia grow vigorously as function of multiplicity. This enhancement is seen in AA, pA [18,19] and even in pp collisions [20,21], which clearly signals that it is not related to collective effects. Similar enhancement was observed for D-and B-meson production [12]. As was mentioned in [22], these new findings cannot be easily accommodated in the framework of models based on the two-pomeron fusion picture and thus potentially could require introduction of new mechanisms both for AA and pp collisions.
Recently, in [23][24][25], it has been suggested that the experimentally found multiplicity dependence could indicate a sizable contribution from multipomeron mechanisms to charmonia production, and has been shown that the inclusion of the three-pomeron mechanism helps to describe the available data. Also, it was shown in [26] that in the D-meson production case the three-pomeron mechanism gives a sizable contribution, which might be responsible for up to 40 percent of all the produced D-mesons. The inclusion of this mechanism improves the overall description of the data, as well as allows to describe its multiplicity dependence. In order to understand better the role of the three-pomeron fusion mechanism in heavy charm production, it has been suggested in [27] to study the multiplicity dependence of diffractive production. Although the predicted cross-section is smaller than for the inclusive case, we expect that its multiplicity dependence could be studied during the High Luminosity Run 3 at the LHC (HL-LHC mode) [28][29][30].
In order to understand better the microscopic mechanisms of multiplicity enhancement in heavy quarkonia production, in this paper we suggest to study the multiplicity dependence of co-produced hadrons in the production of P -wave quarkonia, e.g. the lightest χ c and χ b -mesons. The production of χ c mesons has been extensively studied in the k T factorization approach in [31][32][33][34][35][36][37][38], where it was found that a two-pomeron (two-gluon) fusion mechanism provides a good description of the available data on its rapidity and transverse momentum dependence. Moreover, the color octet mechanism for P -wave quarkonia gives a small contribution due to smallness of the Long Distance matrix Elements (LDMEs) [31][32][33], so this fact minimizes the inherent uncertainty related to this mechanism. On the other hand, due to spin-orbital interactions, the P -wave quarkonia show up as a triplet of states with different angular momenta J = 0, 1, or 2. Independent studies of experimental cross-sections of each of these states provides a sensitive test of the underlying production mechanism. For this reason, we believe that P -wave quarkonia are ideally suited for the study of the multiplicity enhancement mechanisms. As we will show below, in the high energy limit the three-pomeron mechanism does not contribute to P -wave quarkonia production, so we expect that the multiplicity dependence for χ c and χ b mesons should be significantly milder than that of S-wave quarkonia. Since we are interested in the multiplicity dependence, instead of k T -factorization we will use a color dipole framework (also known as CGC/Saturation or CGC/Sat) [39][40][41][42][43][44][45][46][47]. The generalization of this framework to high-multiplicity events is well-known from the literature [48][49][50][51][52][53][54][55].
The paper is structured as follows. In the next section II we describe a framework for χ cJ and χ bJ quarkonia production. In Section III we make numerical estimates of the cross-sections, compare them with available experimental The leading order contribution to the cross-section of P -wave meson production via the two-pomeron fusion mechanism. The diagram includes two cut pomerons (upper and lower gluon ladders). Right plot: Possible contributions of the 3-reggeon mechanism (additional reggeon shown with gray color). As explained in the text, if the two lower gluons reggeize indepedendently (form two pomerons), then such contribution for P -wave quarkonia would vanish. In both plots the vertical dashed line stands for unitarity cuts. The produced meson M is shown with a double line and arrow. A summation over all possible permutations of gluon vertices in the heavy quark line/loop is implied. data and make predictions for future experiments. In Section IV we evaluate the dependence on multiplicity. Finally, in Section V we draw conclusions.

II. PRODUCTION MECHANISMS OF P -WAVE QUARKONIA
In high energy kinematics, the inclusive production gets its dominant contribution from the fusion of two pomerons, which for heavy quarkonia production is given by the diagram shown in the left panel of Figure 1. In the rest frame of one of the protons this process might be viewed as a fluctuation of the incoming virtual gluon into a heavyQQ pair, with subsequent scattering of theQQ dipole with the target proton. In the kinematics of LHC experiments the average light-cone momentum fractions x 1,2 carried by gluons are very small ( 1), and the gluon densities are enhanced. This enhancement implies that there could be sizable corrections from multiple pomeron exchanges between the heavy dipole and the target, which are formally suppressed in the heavy quark mass limit. For this reason, instead of a hard process on individual partons it is more appropriate to use the color dipole framework (also known as CGC/Sat) [39][40][41][42][43][44][45][46][47]. The color dipoles are eigenstates of interaction at high energies, and for this reason can be used as universal elementary building blocks automatically accumulating both the hard and soft fluctuations [56]. In fact, the light-cone color dipole framework has been successfully applied to phenomenological descriptions of both hadron-hadron and lepton-hadron collisions [57][58][59][60][61][62][63][64]. Another advantage of the CGC/Sat (color dipole) framework is that it allows a relatively straightforward extension for the description of high-multiplicity events, as discussed in [48][49][50][51][52][53][54][55]. The cross-section of quarkonia production process, shown in Figure 1, in the dipole approach is given by (see details in Appendix A)  [70,71]. In both plots a vertical dashed grey line stands for the unitarity cut, and 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 transverse size r.
where y and p T are the rapidity and transverse momenta of the produced quarkonia in the center-of-mass frame of the colliding protons; (z i , r i ) are the light-cone fractions of the quark and the transverse separation between quarks inside the dipole (with subindices i = 1, 2, standing for amplitude and its complex conjugate respectively); b 21 is the difference of impact parameters of the dipoles in the amplitude and its conjugate. We also use the notation Ψ M (r, z) for the light-cone wave function of quarkonium M (M = χ c , χ b ), and ΨQ Q for the quark-antiquark component of the gluon light-cone wave function (for the sake of completeness both are discussed in detail in Appendix B). The amplitude N M depends on a linear combination of forward dipole scattering amplitudes N (y, r) ≡ d 2 b N (y, r, b).
The notation x g g (x g , k T ) is used for the unintegrated gluon PDF. The expression for the p T -integrated cross-section has a similar structure and differs only by the replacement of the gluon uPDF x 1 g (x 1 , p T − k T ) by the integrated PDF x 1 g (x 1 , µ F ), taken at the scale µ F ≈ 2 m Q . The integrated gluon PDF x 1 g (x 1 , µ F ) in the CGC/Sat approach is closely related to the dipole scattering amplitude N (y, r) introduced earlier as [48, 65] where y = ln(1/x). Eq. (4) might be inverted and gives the gluon PDF in terms of the dipole amplitude, This result allows us to rewrite (1) entirely in terms of the dipole amplitude N . Now we would like to discuss the contributions of the three-pomeron mechanisms shown in the right panel of Figure 1. Usually it is expected that such contributions are suppressed in the heavy quark mass limit by α s (m Q ), and in certain cases additionally by Λ 2 QCD /m 2 Q . However, for charmonia the latter parameter might not be very small, and for this reason some corrections could be sizable. As was demonstrated in [72,73], this indeed happens in the case of 1S charmonia production, and these contributions are especially important in the events with large multiplicity of co-produced hadrons. For this reason we also need to estimate them in the case of χ c and χ b production. In the color dipole framework it is usually assumed that a universal dipole cross-section takes into account all such contributions. However, in phenomenological parametrizations usually such contributions are either taken into account using additional simplifying assumptions or disregarded altogether. For example, a widely used phenomenological parametrization "CGC", suggested in [66][67][68][69], was inspired by a solution of the Balitsky-Kovchegov (BK) equation and effectively resumms only fan diagrams shown in the left panel of Figure 2. This parametrization does not take into account the three-pomeron contributions at all. The alternative IP-Sat parametrization [70,71], inspired by a Glauber-like approach, resumms the set of BFKL ladder diagrams shown in the right panel of Figure 2. A central assumption of this approach is that the interaction of each 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 the general case it requires more careful treatment. For this reason, in general we cannot extract the contribution of the three-pomeron mechanism by just expanding the dipole amplitude, and instead we should evaluate it explicitly.
As was demonstrated in [23], the contribution due the to the three-pomeron mechanism has a structure similar to (1), although we need to replace the amplitude N M defined in (2) with a new dipole amplitudeÑ M , which explicitly takes into account two-Pomeron interaction of the dipole with the target proton. As was shown in [74], at high energies the dominant configuration (largest intercept) comes from the configuration when each pomeron reggeizes independently, which implies that the amplitudeÑ M has the form where κ is a numerical prefactor, which depends on the transverse profile (impact parameter dependence) of the dipole amplitude N , From the structure of (1) and the symmetry properties of the wave function of the P -wave quarkonia, we can see that such contribution vanishes. The configurations which might give nonzero contributions require antisymmetrization over color indices of both pomerons, thus forming a multireggeon state. It has been shown in the literature that such configurations have smaller intercepts and are suppressed at high energies [74]. For this reason we may consider that in LHC kinematics there is no three-pomeron contributions like those shown in the right panel of the Figure 1.
Finally, we will discuss briefly possible contributions of the color octet mechanism [5][6][7][8][9][10][11], which might be relevant in large-p T kinematics. For the P -wave quarkonia the relevant contribution is controlled by the Long Distance Matrix . The analyses available from the literature [31][32][33] conclude that the value of this LDME is very small, although the estimates of its exact value vary significantly, between 4.78 × 10 −5 and 2.01 × 10 −3 GeV 3 . In view of these findings, in what follows we will simply omit the contribution of the color octet mechanism.

III. NUMERICAL ESTIMATES
In the CGC/saturation approach, the dipole amplitude N y, r, b is expected to satisfy the non-linear Balitsky-Kovchegov equation [75] for the dipoles of small size r. In the saturation region this solution should exhibit a geometric scaling, being a function of one variable τ = r 2 Q 2 s , where Q s is the saturation scale [76][77][78][79]. Such behavior is implemented in different phenomenological parametrizations available from the literature. One of such parametrizations which we will use for our numerical estimates is the CGC parametrization [67], We would like to start our discussion of results from a comparison of the predicted p T -dependence of the crosssections with experimental data. The cross-section of χ c production is smaller than the cross-section of J/ψ, for this reason there is much less experimental data available from the literature. Since χ cJ is usually detected via the χ cJ → γ + J/ψ radiative decay channel, the experimental data are traditionally presented for the product of the cross-section onto the branching fraction B (χ cJ ) ≡ Br (χ cJ → γ + J/ψ) Br (J/ψ → µ + µ − ). For the χ bJ mesons we use a similar product of branching fractions with Υ(1S) instead of J/ψ. The values of the branching fractions are Table I: Values of the product of branching fractions B (χcJ ) ≡ Br (χcJ → γ + J/ψ) Br J/ψ → µ + µ − and B (χ bJ ) ≡ Br (χcJ → γ + Υ(1S)) Br Υ(1S) → µ + µ − , as given in [80].   Figure 3: Left: Comparison of the predicted pT -dependence for the χc1 and χc2 cross-sections. Experimental data are from ATLAS [81]. Right: Comparison of the model predictions for the ratio of χc1 and χc2 cross-sections at central rapidities, with experimental data from ATLAS [81], CMS [82] and LHCb [83] experiments. We added for comparison the LHCb data measured at off-forward rapidities, because we expect that the suppression of the cross-sections of χc1 and χc2 production at off-forward rapidities will be the same and thus will cancel in the ratio. For better visibility we use a logarithmic scale in the vertical axis.
known from [80] and for the sake of completeness are shown in Table I. As we can see, the values of B (χ c0 ) and B (χ b0 ) are extremely small compared to the other channels. For this reason observation of these states via radiative decays into 1S quarkonia is very difficult, and all the available data are given for χ c1 and χ c2 mesons. In Figures 3, 4 we compare the model predictions for the χ cJ production with available data from ATLAS [81], CMS [82], LHCb [83] and CDF [84]. We can see that the CGC/Sat model provides a reasonable description of the available data in a wide kinematic range and thus might be used for further analysis.   [86,87] and Power-like [88] parametrizations of the potential (see also Appendix B for more details).
In Figure 5 we show the p T -dependence of the cross-sections at different values of the collision energies √ s, which might be relevant for the future experimental data. We expect that the cross-section will proportionally increase as a function of energy, without changing its shape. In order to illustrate the dependence on the choice of wave function (∼potential model used for its evaluation), we have also shown in the same Figure the ratio of the crosssections evaluated with Cornell [86,87] and Power-like [88] parametrizations of the rest frame potential. While the cross-sections changes several orders of magnitude in the considered range of p T , the uncertainty due to choice of the potential does not exceed 15 per cent. We also got similar estimates for the parametrization [89]. The cross-sections of χ c0 and χ c2 get large contributions from the configuration with aligned spins of the quarks, whereas in the case of χ c1 there is also a sizable contribution from configurations where spins are anti-aligned, which explains the difference.
In Figure 6 we show our predictions for the p T -dependence of χ b1 and χ b2 mesons. The estimated cross-sections are smaller than for the χ c1 and χ c2 mesons, although within the reach of LHC experiments. So far there is no published data from LHC for the cross-sections of these mesons, but we hope that in the near future such measurements could be carried out.

IV. MULTIPLICITY DEPENDENCE
As we found in the previous section, the CGC/Sat model (1) provides a reasonable description of the χ c1 and χ c2 production data, at Tevatron and LHC kinematics. 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 the number of produced charged particles N ch [90], therefore for the study of the multiplicity dependence it is more common to use a self-normalized ratio [21] dN M /dy where n = N ch / N ch is the relative enhancement of the charged particles in the bin, w (N M ) / w (N M ) and w (N ch ) / w (N ch ) are the self-normalized yields of quarkonium M and charged particles (minimal bias) events in a given multiplicity class; dσ M (y, √ s, n) is the production cross-sections for M , with rapidity y and N ch = ∆η dN ch /dη charged particles in the pseudorapidity window (η − ∆η/2, η + ∆η/2). If the inclusive cross-section of the process pp → M + X is proportional to probability to produce a meson M in a single pp collision, then the ratio (13) gives a conditional probability to produce a meson M in a pp collision in which N ch charged particles are produced. Due to the Local Parton-Hadron Duality (LPHD) hypothesis [91][92][93], the number of produced charged particles is directly proportional to the number of partons which stem from the individual pomerons and thus might 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 large 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 therefore we expect that the dipole amplitude (8) should maintain its form, although the value of the saturation scale Q s might be modified. As was demonstrated in [48][49][50], 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 (14) algebraically, 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 (14) allows to express Q 2 s as a function of n. Usually in the literature the logarithmic dependence on n, which stems from the running coupling in the denominator of (14) is disregarded, so (14) reduces to a simpler linearly growing dependence on n [48][49][50][51][52][53][54], The precision of this assumption was tested in [55], and it was found that a numerical solution of the running coupling Balitsky-Kovchegov (rcBK) equation differs from the approximate (15) by less than 10% in the region of interest (n 10). This correction is within the precision of current evaluations, and for this reason in what follows we will use (15) Figure 7: Left: multiplicity dependence for different χcJ states. While the cross-sections differ quite significantly due to spin structure, the self-normalized ratios are very close to each other. For the sake of reference we also added a dot-dashed grey curve for J/ψ production from our previous [23,72]. Right: Dependence of the multiplicity shapes on collision energies √ s. The plot is done for χc1 meson production, but the results for χc0, χc2 are almost identical. All evaluations are done assuming that charged particles and quarkonia are collected at central rapidities (|η, y| < 1), similar to what is available for J/ψ production from [16].
fall into the range 0.5-1 GeV, from (15) 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. Since increase of multiplicity and increase of energy (decrease of x) affect Q 2 s in a similar way, the study of the high-multiplicity events allows to study a deeply saturated regime which determines the dyanamics of all processes at significantly higher energies.
For phenomenological estimates of the multiplicity dependence it is very important that the rapidity bins used for collecting the co-produced charged hadrons overlaps with the rapidity bins used for collection of quarkonia. As was illustrated in our previous publications [23,26,72], the observed enhanced multiplicity should be shared equally between all cut pomerons which could contribute to co-production of hadrons in a given experimental setup. For this reason the strongest multiplicity dependence will show up in the case when both quarkonia and hadrons are collected at central rapidities. In what follows we will focus on this setup.
In Figures 7, 8 we show the multiplicity dependence of the χ cJ and χ bJ mesons for different energies. As we can see, the dependence is much milder than that of the 1S quarkonia (dot-dashed curve with label "J/ψ"), in agreement with our expectations based on dominance of two-pomeron mechanism. From (15) and the structure of the dipole cross-section (8) we can expect that each cut pomeron contributes to the multiplicity dependence factor ∼ n γ eff , where the parameter γ eff was defined in (10). Since χ b has smaller size than χ c , the typical values of γ eff are larger for the former than for the latter, and χ b has slightly faster dependence on multiplicity than χ c . Similarly, we can understand the change of multiplicity dependence with energy: due to prefactor 1/Y in (10), the average values of the parameter γ eff decrease as a function of √ s, and for this reason the dependence on multiplicity becomes milder for larger energies √ s pp .

V. CONCLUSIONS
In this paper we analyzed in detail the production of χ c and χ b mesons in the CGC/Sat approach. We found that the model predictions for the p T -dependent cross-section are in agreement with available experimental data for χ c1 and χ c2 mesons in LHC kinematics. We also made predictions for χ b mesons, which might be checked in the ongoing and future experiments, both at RHIC and at LHC. We also studied the dependence of the cross-sections on multiplicity of co-produced hadrons, and found that it is significantly milder than that of 1S quarkonia (J/ψ, Υ...). This happens because the dominant production mechanism of for the P -wave quarkonia is the two-pomeron fusion, whereas the three-pomeron contributions are strongly suppressed at high energies. Our evaluation is largely parameter-free and relies only on the choice of the parametrization for the dipole cross-section (8) and the wave function of the meson.
The explanation of multiplicity dependence in the CGC/Sat approach differs from other approaches suggested for the description of multiplicity dependence, like e.g. the percolation approach [94] Figure 8: Left: multiplicity dependence for different χ bJ states. While the cross-sections differ quite significantly due to spin structure, the self-normalized ratios are very close to each other. Right: Dependence of the multiplicity shapes on collision energies √ s. The plot is done for χ b1 meson production, but results for χ b0 , χ b2 are almost identical. All evaluations are done assuming that charged particles and quarkonia are collected at central rapidities (|η, y| < 1), similar to what is available from [16].
(a) (b) (c) Figure 9: 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 might be also viewed as gluon-gluon fusion gg → g, with subsequent gluon fragmentation g →QQ. In CGC parametrization of the dipole cross-section approach each "gluon" is replaced with reggeized gluon (BK pomeron), which satisfies the Balitsky-Kovchegov equation and corresponds to a fan-like shower of soft particles.
elastic amplitude [95]. While for 1S quarkonia all approaches give comparable descriptions, this is not so for P -wave quarkonia. For this reason we expect that the measurement of the multiplicity dependence of χ c and χ b would be an important litmus test for all the models which describe production of quarkonia, and we hope that it will be done both at LHC and RHIC.

VI. ACKNOWLEDGEMENTS
We thank our colleagues at UTFSM university for encouraging discussions. This research was partially supported by the project ANID PIA/APOYO AFB180002 (Chile) and Fondecyt (Chile) grant 1180232. Also, we thank Yuri Ivanov for technical support of the USM HPC cluster, where some evaluations were performed.

Appendix A: Evaluation of the dipole amplitudes
In this Appendix, for the sake of completeness, we explain the main technical steps and assumptions used for the derivation of the cross-section (1). The general rules, which allow to express the cross-sections of hard processes in terms of the color singlet dipole cross-section, might be found in [39][40][41][42][43][44][45][46][47]. In the heavy quark mass limit the strong coupling α s (m Q ) is small, so the interaction of a heavyQQ dipole with gluons might be considered perturbatively. At the same time, we tacitly assume that each such gluon should be understood as a parton shower ("pomeron").
In the high-energy eikonal picture, the interaction of the quarks and antiquark with a t-channel gluon are described by a factor ±ig t a γ (x ⊥ ), where x ⊥ is the transverse coordinate of the quark, and the function γ (x ⊥ ) is related to a distribution of gluons in the target. This function is related to a dipole scattering amplitude N (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) might be rewritten in the form The value of the constant is related to the infrared behavior of the theory, and for the observables which we consider in this paper, it cancels exactly. For very small dipoles, the dipole scattering amplitude is related to the gluon uPDF as [108] N (x, r) = 4πα s 3 so the functions γ (x, r) might be also related to the unintegrated gluon densities in coordinate space. With the help of (A2), it is possible to express the exclusive amplitudes and inclusive cross-sections as linear combinations 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.
For the case of P -wave production, we should take into account that the quark and antiquark transverse coordinates are given by b i − z i r i and b i +z i r i respectively, where subindex i takes values i = 1, 2 for the amplitude and its conjugate. 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 , so we have to evaluate the additional convolution is the difference of impact parameters of the dipole center of mass in the amplitude and its conjugate, and b = (b 1 + b 2 )/2. The integral over d 2 b is evaluated using (A2), and after simple algebraic simplifications we can get (1). For the three-pomeron we consider that in a target rest frame we have a scattering of a dipole in the field of two pomerons. Since at high energies each pomeron reggeizes independently, we expect that the amplitude will be a mere product of the scattering amplitudes on each pomeron, thus yielding (6).
All evaluations of this appendix were done in the frame where the momentum of the primordial gluon is zero. In any other frame we should take into account an additional convolution with the momentum distribution of the incident ("primordial") gluons, as appears in the first line of (1).

Appendix B: Wave functions
For evaluations of 1, we will need explicit parametrizations for theQQ component of the light cone gluon wave function ΨQ Q and the P -wave quarkonia wave function Ψ M . We may expect that in the heavy quark mass limit for ΨQ Q we may use the well-known perturbative expressions available from the literature [96,97], where the superscript index (±1) of ΨQ Q refers to helicity of the projectile gluon, and h,h in the r.h.s. are the helicities of the quark and antiquark respectively. Most previous evaluations of quarkonia production [31][32][33][34][35][36][37] were done in the heavy quark mass limit, assuming that the wave function Ψ M (z, r) might be approximated with its small-r Taylor expansion, which starts with a linear term for P -wave. In this case the dependence on the wave function reduces to dependence on the value of the slope |R (0)|. This scheme is justified in the heavy quark mass limit. However, it is clear that for charmonia the heavy quark mass limit might not work very well, and for this reason we will maintain the full r-dependence of the wave function.  For our evaluations we construct a light-cone wave function from the rest frame wave functions evaluated in the potential models using the Brodsky-Huang-Lepage (BHL) prescription [98] (see also a similar scheme [99]). It is known that for heavy quarkonia the results of the BHL prescription are close to the wave functions evaluated in Covariant Spectator Theory [100] as well as lattice evaluations [101][102][103]. According to BHL prescription, the light-cone wave function related to the rest frame wave function as [98,104,105] Φ LC (z, r) = d 2 k ⊥ e ik ⊥ ·r ⊥ k 2 ⊥ + m 2 where ψ RF (k) is the Fourier image of the rest frame wave function. Due to spin-orbital interaction, we expect that the spinorial structure of the wave functions of χ cJ , χ bJ will be depend crucially on the angular momentum J, and for this reason the production cross-sections of χ cJ , χ bJ mesons will acquire dependence on J [109]. The rest frame quarkonium wave function might be written using the Clebsch-Gordan coefficients as where {χ s } are spinors corresponding to definite projection of spin s of the quark and antiquark, M is the projection of orbital angular momentum, R n1 (r) is the radial wave function of the P -wave quarkonium, and Y 1M is an ordinary spherical harmonic. A set of useful relations between Clebsch-Gordan coefficients, which facilitate the summations, might be found in [32,38]. For the evaluations of radial part R n1 we used expressions found in potential models. For the sake of comparison in evaluation we considered three different choices of the potential: • The Cornell potential which was introduced in [86,87], • The potential with logarithmic large-r behaviour suggested in [89], V Log (r) = α + β ln r.
• The power-like potential introduced in [88], We checked that the wave functions obtained with all three potentials have similar shapes. As we mentioned earlier, in the heavy quark mass limit we expect that the wave function might be approximated by its behavior near the r ≈ 0. The wave functions of the P -wave quarkonia have a node at r ≈ 0, so the relevant parameter, which determines the cross-sections in this limit, is the value of the slope |R (0)| 2 . To facilitate comparison with other approaches, in Table II we provide the values of the latter parameter. For the evaluation of the overlap with gluon wave function (B1), it is convenient to rewrite the quarkonium wave functions in a helicity basis. Conventionally, this is done applying the Melosh-Wigner spin rotation operators defined in [106,107].