Hadroproduction at Large Hadron Collider Challenges NRQCD Factorization

We report on our analysis [1] of prompt ηc meson production, measured by the LHCb Collaboration at the Large Hadron Collider, within the framework of nonrelativistic QCD (NRQCD) factorization up to the sub-leading order in both the QCD coupling constant αs and the relative velocity v of the bound heavy quarks. We thereby convert various sets of J/ψ and χc,J long-distance matrix elements (LDMEs), determined by different groups in J/ψ and χc,J yield and polarization fits, to ηc and hc production LDMEs making use of the NRQCD heavy quark spin symmetry. The resulting predictions for ηc hadroproduction in all cases greatly overshoot the LHCb data, while the color-singlet model contributions alone would indeed be sufficient. We investigate the consequences for the universality of the LDMEs, and show how the observed tensions remain in follow-up works by other groups.


Introduction
Despite great efforts from both experimental and theoretical sides ever since the J/ψ meson was discovered in 1974, the genuine mechanism describing the production of heavy quarkonium states, i.e. heavy quark-antiquark pairs Q Q = cc, b b, has remained mysterious.The non-relativistic QCD (NRQCD) factorization formalism [2], which is built on the non-relativistic field theory [3], is now arguably the most probable candidate theory in town.In this approach the production of heavy quarkonium factorizes into short-distance coefficients (SDCs), which are process dependent and can be calculated perturbatively as expansions of the QCD coupling constant α s , and long-distance matrix elements (LDMEs), which are supposed to be universal, and can be obtained through fitting to experimental data.The LDMEs are organized as powers of the small velocity v of the heavy (anti-)quark in the meson rest frame.The main difference between NRQCD factorization formalism and the colorsinglet (CS) model is that in CS model (CSM) the Q Q pair can only be in a CS state and have the same quantum numbers (total spin S , orbital angular momentum L, and total angular momentum J) as the quarkonium considered, whereas the NRQCD factorization formalism accommodates for contributions of intermediate Q Q [n] pairs in all possible Fock states n = 2S +1 L [a]  J , where a = 1, 8 stands for CS and color octet (CO), respectively.On the theory side, the NRQCD factorization formalism resolves a conceptual problem of the CSM, namely the existence of uncanceled infrared (IR) divergences when L > 0 [4].
Despite its theoretical advancement, NRQCD factorization looks less compelling when confronted with J/ψ data measured in various environments, like in e + e − annihilation, γγ collisions, photoproduction, or hadroproduction.The SDCs for all of these processes have been calculated at NLO in the α s expansion and some of them even at the v 2 subleading order (see Ref. [5] for a recent review and references therein for details).By a global fit to the worldwide J/ψ yield data using SDCs at QCD NLO [6], three CO LDMEs O J/ψ (1 S 8 0 ) , O J/ψ ( 3 S 8 1 ) , and O J/ψ ( 3 P 8 0 ) , have been successfully pinned down in compliance with the velocity scaling rules.However, the resulting predictions of the J/ψ polarization in hadroproduction strikingly disagree with the measurements at the Fermilab Tevatron and the CERN LHC [7].On the other hand, reasonably well fits are obtained when only J/ψ yield and polarization data is considered in hadroproduction with a cut on the J/ψ transverse momentum p T ≈ 7 GeV [8][9][10], but they all largely overshoot the world's data from other colliders than hadron colliders [11].
Recently, the LHCb Collaboration measured the production of the η c meson (J PC = 0 −+ ), the spin-singlet partner of the J/ψ meson, via its decay to p p [12].The measurement was performed at center-of-mass energies √ S = 7 and 8 TeV in the rapidity range 2.0 < y η c < 4.5 in bins of p η c T .This measurement could provide another independent test of the NRQCD factorization conjucture and shed light on the J/ψ polarization puzzle, since the LDMEs of these two states are related by the heavy quark spin symmetry (HQSS), another pillar of NRQCD factorization.According to the NRQCD velocity scaling rules, the leading CS and CO Fock states of direct η c production are 1 S [1]   0 at O(v 3 ), and 1 S [8]  0 , 3 S [8]  1 , and 1 P [8]  1 at order v 7 .For prompt η c production there may also be a large feeddown contribution from h c 1 , similar to the J/ψ case, where the χ cJ feeddown contributes up to 20 − 30% [9,13,14].The leading CS and CO Fock states of h c production are 1 P [1]  1 , 1 S [8]  0 at order v 5 .So far, calculations of the SDCs of direct η c hadroproduction have been done only at leading order (LO) in α s and v 2 and were missing the 1 S [8]  0 channel [15].From J/ψ production, we learned that the NLO QCD corrections can be orders of magnitude larger than LO contributions [6,[8][9][10] and that the relativistic corrections can not be ignored either [16][17][18].Therefore, it is of great interest to do a fully-fledged NRQCD analysis of prompt η c production at NLO in α s and subleading order in v 2 .

NRQCD Factorization Formalism
We work in the QCD collinear parton model implemented in the fixed-flavor scheme with n f = 3 active quark flavors in the colliding protons, which are represented by the parton density functions (PDFs) at the factorization scale µ f .Together with the NRQCD factorization, at α s NLO and v 2 subleading order, the cross sections can be expressed as [1]  0 , 1 S [8]  0 , 3 S [8]  1 , 1 P [8]   1 where O h (n) ( P h (n) ) are the v 2 LO (subleading) LDMEs with h = η c , h c , and dσ cc [n] (dσ cc [n]  v 2 ) the corresponding SDCs including their NLO QCD corrections.To account for the mass difference between η c and h c mesons, we use p h c T = p η c T m h c /m η c for the transverse h c momentum p h c T .This is a fair approximation.The definitions of the O and P operators for S -wave production and the O operators for the P-wave production are given in [2,17].We define the P production operators for the P-wave Fock states in a similar way as The HQSS implies that the η C (h c ) and J/ψ (χ c0 ) LDMEs are up to corrections of order v 4 related via [2] Q η c ( 1 S [1]  0 / 1 S [8]  0 ) = with Q = O, P.
We calculate the QCD corrections and the v 2 relativistic corrections using the techniques developed in Refs.[17,19,20].Because of the limitation of length, we will not go into any details but instead just highlight some some key points.When we calculate the QCD corrections to P-wave production there are uncancelled infrared divergences in the real radiative corrections.These match the product of LO S -wave SDCs and the QCD corrections to their LDMEs after appropriate MS operator renormaliztion.The latter are given by where O h (n) 0 are the tree-level LDMEs and µ r and µ λ are the QCD and NRQCD renormalization scales, respectively.The µ λ dependence of the renormaliztion scales can then be determined by solving the operator evolution equation The computation of QCD corrections to 1 P [1]  1 state hadroproduction was carried out in [22].Within the uncertainties due to phase-space-slicing method, our results are consistent with theirs.However, when we vary µ λ around its default value, we find significant differences to [22] and we can trace the difference back to the fact that they perform the variation of µ λ only in the SDCs but not in the scale dependent LDMEs.The QCD corrections to the SDCs for the 1 S [1]  0 and 1 P [8]  1 channels as well as the SDCs for the v 2 subleading order operators of 1 S [1]  0 , 1 P [1]  1 and 1 P [8]  1 are calculated by us for the first time.
to 2 times their default values in order to estimate the scale uncertainties.The values m η c = 2983.6MeV, m h c = 3525.28MeV and Br(h c → η c + γ) = 51% are taken from Ref. [24].
To investigate how large the NLO QCD corrections and v 2 corrections to the contributing SDCs are for unit LDMEs, we define the K-factor and the R-factor as K(n) and R(n) as functions of p η c T for the 7 TeV LHCb kinematic setup are shown in Fig. 1.We observe that the QCD corrections turn the SDC of 1 P [1]  1 negative with a huge K-factor, while the R-factor is tiny.This feature is similar, for example, to the 3 P [8]  J SDC of direct J/ψ hadroproduction.The 1 P [8]   1 SDC, however, remains positive at QCD NLO.As for the R-factors, they are almost independent of p T when p η c T > 7GeV, and they all are of order unity except for n = 1 P [1]  1 , which indicates that the relativistic corrections are indeed of relative order v 2 .
We adopt two approaches to compare our NRQCD predictions with the LHCb data [12].In the first one, we implement the HQSS via Eq.( 3) to obtain η c and h c LDMEs from the J/ψ and χ c0 LDME sets determined by four independent fits at QCD NLO, which are listed in Table 1.In the cases where no χ cJ and CS J/ψ LDMEs are available, we simply omit those contributions.In the 1 η c : -1 S [8]   0 η c : -1 P [8]   1 Total h c : -1 P [1] 1 η c : -1 S [8]   0 η c : -1 P [8]   1 Total h c : -1 P [1]   1 h c : 1 S [8]    following we will see that η c production almost exclusively happens through direct production, and in particular mainly through the 3 S [8]  1 channel.This observation will provide a retroactive justification to omit the other LDMEs where not available.In Fig. 2, we compare the NRQCD and CSM default predictions evaluated with the four LDME sets in Table 1, including the NLO QCD corrections but not yet the relativistic corrections, with the LHCb data.We consider three sources of errors, which are added quadratically to give our uncertainty bands.The first comprises of the scale uncertainties of µ λ , µ r , and µ f .The second one is made up by the fit errors of the LMDEs given in Table 1.The third one is due to the fact that HQSS relations (Eq.( 3)) only hold at v 2 LO and that the values of the LDMEs P h (n) are unknown.These effects are estimated by setting P h (n) = ξm 2 c O h (n) and varying ξ between −0.5 and 0.5 so as to be consistent with the result v 2 ≈ 0.23 from potential model calculations [25].The effect of the relativistic corrections thus enters through this third error.
We also show the contribution of each Fock state separately choosing the default values of the parameters in Fig. 2. We see now that the h c feeddown contributions are indeed negligible because of the small SDCs of 1 P [1]  1 and 1 S [8]  0 .The smallness of the feeddown contributions could not have been anticipated without this explicit calculation, having in mind the large χ cJ feeddown to prompt J/ψ.The most interesting feature shown in Fig. 2 is however that the CSM prediction, made up basically by the 1 S [1]  0 channel alone, can describe the LHCb data almost perfectly, leaving little room for CO contributions after all.While the contributions from 1 S [8]  0 and 1 P [8]  1 channels are indeed tiny for all four J/ψ LDME sets considered, the 3 S [8]  1 contribution, however, turns out to be dominating the NRQCD prediction, and lets it overshoot the LHCb data by up to one order of magnitude.Even the LDME set from Ref. [6], which has the least discrapency to the LHCb data, still yields an unacceptable value of lower bounds of the respective error bands in Fig. 2, the χ 2 /d.o.f descends to 36.7/7, which is still too large to be acceptable.
In our second approach, we determine the LDMEs of η c and h c by directly fitting the LDMEs to the LHCb data without making use of the HQSS.As shown in the previous paragraph, we can neglect the h c feeddown contributions, and also the contributions of the 1 S [8]  0 and 1 P [8]  1 states, from the beginning.Their SDCs are simply not large enough to compensate the O(v 4 ) suppression of the respective CO LDMEs relative to the CS LDME O η c ( 1 S [1]  0 ): The SDCs of 1 S [8]  0 are only of the same order as those of 1 S [1]  0 , while the 1 P [8]  1 SDCs are even smaller.Therefore we are only left with the contributions from the 1 S [1]  0 and 3 S [8]  1 states.Our fitting procedure is that we first determine the CS LDMEs from the η c inclusive decay width [26] and the partial decay width of η c → γγ, and then with those CS LDMEs as inputs fit the CO LDMEs to the LHCb data.We are allowed to do the CS LDME extraction from the decay rates since in the CS case the differences between the decay and production LDMEs are only of order v 4 .In the determination of the CS LDMEs, we set m c = 1.5 GeV, α = 1/137, α s (2m c ) = 0.26 and take the values Γ(η c ) = (32.3± 1.0) MeV and Br(η c ) → γγ = (1.57± 0.12) × 10 −4 from [24].We perform the fitting once excluding and once including the O(v 2 ) corrections.When excluding the O(v 2 ) corrections, we obtain 2 O η c ( 1 S [1]  0 ) = (0.24 ± 0.02) GeV 3 and where the first error is the fit error, and the second one is due to varying the factorization and renormalization scales, yielding a very nice description of the LHCb data with χ 2 /d.o.f = 1.4/6.The results are shown in the upper two panels in Fig. 3.When including the O(v 2 ) corrections, we obtain O η c ( 1 S [1]  0 ) = (0.50 ± 0.009) GeV 3 , P η c ( 1 S [1]  0 ) /m 2 c = (0.14 ± 0.005) GeV 3 , and O η c ( 3 S [8]  1 ) = (−6.6 ± 2.3 ± 4.4) × 10 −3 GeV 3 , where O η c ( 3 S [8]  1 ) = O η c ( 3 S [8]  1 ) − 1.1 P η c ( 3 S [8]  1 ) .Note that since R( 3 S [8]  1 ) does almost not change at all with p T (see Fig. 1), we cannot fit O η c ( 3 S [8]  1 ) and P η c ( 3 S [8]  1 ) separately, but instead fit their linear combination O η c ( 3 S [8]  1 ) .The result, again looking good, is shown in the lower two panels of Fig. 3. Now setting the CS contribution to zero and assuming that the contributions of the other Fock states are positive, we obtain an upper bound of O η c ( 3 S [8]  1 ) < (1.5 ± 0.2 +0.3 −0.4 ) × 10 −2 GeV 3 .If we further set P η c ( 3 S [8]  1 ) /(m 2 c O η c ( 3 S [8]  1 ) ) = v 2 = 0.23, we obtain O J/ψ ( 1 S [8]  0 ) = O η c ( 3 S [8]  1 ) < (1.9 ± 0.3 +0.4 −0.5 ) × 10 −2 GeV 3 .
We would like to point out that if we similarly fit the J/ψ CS LDMEs to its decay rate into e + e − and three gluons, we obtain that O J/ψ ( 3 S [1]  1 ) /3 = 0.63 ± 0.003 GeV 3 , P J/ψ ( 3 S [1]  1 ) /(3m 2 c ) = (7.7 ± 0.01) × 10 −2 GeV 3 , where the decay widths Γ(J/ψ) → e + e − and Γ(J/ψ) → ggg are again taken from [24].By comparing these values to the η c CS LDME values obtained through the fit to the η c decay rates in the previous paragraph, we find that the HQSS holds nicely among these two fits of the J/ψ and η c CS LDMEs.The absence of large α s and v 2 corrections in direct η c production (see Fig. 1) and the smallness of h c feeddown contributions (see Fig. 2) are further arguments in favour of this η c LDME determination, since it shows that the theoretical predictions of prompt η c production are well under control, in contrast to the J/ψ case.We are thus confident about the above constraint (6).

Tensions in Different Fits and Consequences for LDME Universality
Our work shows that it seems not possible to accomodate for all charmonium data with a single set of LDMEs, even if one restricts the analysis only to middle-and high-p T hadroproduction.The basic conflict is that the LHCb η c data demands a low O J/ψ ( 1 S [8]  0 ) value respecting (6), while the observed unpolarized J/ψ hadroproduction demands a dominant 1 S [8]  0 channel at high p T and thus a large value of O J/ψ ( 1 S [8]  0 ) , as can explicitly be seen from the fit values of [8,10].This tension is also evident from the two works [27,28], which appeared shortly after our analysis and similarly deal with the LHCb η c data within the NRQCD factorization framework.In [27], LDME values are used which nicely describe both the J/ψ and η c hadroproduction yield.Here the tension lies again with the J/ψ polarization: At high-p T the helicity frame polarization value λ θ does not drop below 0.4, in clear conflict to LHC and Tevatron data, see e.g. the first panel of Fig. 8 in [29] and the second panel of Fig. 2 in [27].Another route is chosen in [28].Here, the CS LDME O η c ( 3 S [8]  1 ) is fitted along with the CO LDMEs to J/ψ and η c yield and polarization data to obtain a reasonably well description.Although there are still discrapencies with the Tevatron polarization data (see low panel of Fig. 3 in [28]), the main tension here is between the CS decay and production LDMEs.Their value of O J/ψ ( 3 S [1]  1 ) /3 = O η c ( 1 S [1]  0 ) = (0.16 ± 0.08) GeV 3 is incompatible with the values extracted from decay rates and potential model calculations of both J/ψ and η c , see e.g.Table 1, or the corresponding values of the previous section 3 .It is needless to say that the LDME values of both works [27] and [28] fail to describe J/ψ data from e + e − collisions, photoproduction, as well as hadroproduction data where p T <7 GeV.

Summary
We calculated the α s corrections to the SDCs of the 1 S [1]  0 and 1 P [8]  1 cc production channels as well as the SDCs of the v 2 subleading order operators of the 1 S [1]  0 , 1 P [1]  1 , and 1 P [8]  1 channels.We then compared the NRQCD predictions to recent LHCb prompt η c production data using η c and h c LDMEs 3 For completeness we add that in [28] two LDME fits were done with two different choices for the renormalization and factoriztion scales, the value quoted here is the one using µ r = µ f = m T , which is the reasonable scale choice used in basically all other quarkonium production calculations and fits.For µ r = µ f = E ηc the value O ηc ( 1 S [1]  0 ) = (0.23 ± 0.12) GeV derived from up-to-date J/ψ production LDMEs with the help of HQSS.We demonstrated that the CS contribution alone can already nicely describe the LHCb measurement, while the NRQCD predictions overshoot the data by far, resulting in an unacceptably large χ 2 /d.o.f value of 5.24 or above.We found that the 3 S [8]  1 contribution is the dominating CO channel and that the h c feeddown contributions are negligable.Therefore we were allowed to determine the CS 1 S [1]  0 LDMEs from η c decays, then fit O η c ( 3 S [8]  1 ) directly to the LHCb data and thus to infer a conservative upper bound on O J/ψ ( 1 S [8]  0 ) = O η c ( 3 S [8]  1 ) .The value obtained is much smaller than is to be expected from the NRQCD velocity scaling rules and much smaller than the values obtained in the previous J/ψ production fits.We then carried on to describe how the tensions between the different data remain in two later η c data analyses performed by other groups.
HQSS is a direct consequence of NRQCD.Furthermore, the LHCb data was taken at moderate p T values, where neither large logarithms ln(p 2 T /m 2 c ) nor a breakdown of NRQCD factorization is anticipated by anyone.Therefore, assuming that the LHCb η c measurement will be confirmed by other experiments, we have to conclude that either the universality of the NRQCD LDMEs is in question or that another important ingredient to current NLO NRQCD analyzes has been overlooked so far.

1 TotalFigure 2 .
Figure 2. The LHCb[12] measurements of dσ/d p T for prompt η c hadroproduction at √ s = 7 TeV (upper panel) and 8 TeV (lower panel) are compared with the default predictions of NRQCD (solid lines) and the CSM (dot-dashed lines) at NLO, but without relativistic corrections, evaluated with the four LDME sets in Table1.The theoretical errors as explained in the text are indicated by the yellow and blue bands, respectively.For comparison, also the default contributions due to the individual Fock states are shown.Red color (minus sign in the legend) indicates negative values.

Figure 3 .
Figure 3.Comparison between NRQCD predictions calculated in the second approach as described in the text and LHCb measurements [12] of dσ/d p T for prompt η c production at √ s = 7TeV (left panel) and 8TeV (right panel).The results shown in the lower(upper) panel are calculated with(out) the O(v 2 ) corrections.