EUROPEAN LABORATORY FOR PARTICLE PHYSICS
CERN–PH
March 11, 2014
Quarkonium production in the LHC era:
[3mm] a polarized perspective
Pietro Faccioli, Valentin Knünz, Carlos Lourenço,
João Seixas and Hermine K. Wöhri
Abstract
Polarization measurements are usually considered as the most difficult challenge for the QCD description of quarkonium production. In fact, global data fits for the determination of the nonperturbative parameters of boundstate formation traditionally exclude polarization observables and use them as a posteriori verifications of the predictions, with perplexing results. With a change of perspective, we move polarization data to the centre of the study, advocating that they actually provide the strongest fundamental indications about the production mechanisms, even before we explicitly consider perturbative calculations.
Considering and measurements from LHC experiments and stateoftheart NLO shortdistance calculations in the framework of nonrelativistic QCD factorization (NRQCD), we perform a search for a kinematic domain where the polarizations can be correctly reproduced together with the cross sections, by systematically scanning the phase space and accurately treating the experimental uncertainties. This strategy provides a straightforward solution to the “quarkonium polarization puzzle” and reassuring signs that the theoretical framework is reliable. At the same time, the results expose unexpected hierarchies in the nonperturbative NRQCD parameters, that open new paths towards the understanding of boundstate formation in QCD.
Submitted to Phys. Lett. B
Laboratório de Instrumentação e Física Experimental de Partículas (LIP), Lisbon, Portugal
Physics Department, Instituto Superior Técnico (IST), Lisbon, Portugal
Institute of High Energy Physics (HEPHY), Vienna, Austria
European Organization for Nuclear Research (CERN), Geneva, Switzerland
1 Introduction, motivation and conceptual remarks
Up to the early 1990’s, quarkonium production was believed to be reasonably well described by the (leading order, LO) colour singlet model (CSM) [1], which basically assumes that the produced quarkantiquark () pair has, since its inception, the spin (), angular momentum () and colour quantum numbers of the observable quarkonium, and that no /changing transitions occur during the boundstate formation. Measurements of the J/ and production cross sections by the E789 fixedtarget experiment at Fermilab [2], which exceeded the predicted values by factors of 7 and 25, respectively, challenged this model. Nonetheless, since the data were collected at relatively low transverse momentum, , where nonperturbative effects not addressed by the CSM could dominate, this discrepancy was not immediately perceived as a potentially serious problem.
The results obtained by CDF at the Tevatron [3] covered a much higher range and showed an even larger difference between the calculated and measured prompt production cross sections (after subtracting the nonnegligible yield from B meson decays). The J/ discrepancy could be attributed to feeddown contributions from decays of the wave states and , poorlyknown experimentally and expected to dominate prompt J/ production because the directlyproduced J/ component was supposed to be phasespace suppressed by the need of an extra gluon. The clearcut perception that there was a problem in the understanding of directlyproduced charmonium production was unravelled by the measurement, insensitive to feeddown decays and a factor 50 higher than calculated in the CSM, a surprisingly large discrepancy nicknamed “the CDF anomaly”.
Meanwhile, significant progress was being made on the theory side, with the birth of the nonrelativistic quantum chromodynamics (NRQCD) factorization approach [4]. While in the CSM the observed boundstate meson can only result from quark pairs produced in a singlet state, NRQCD includes terms where the original quark pairs are in colouroctet states. In this effective field theory, the nonperturbative evolution that converts the (coloured) into a physicallyobservable bound meson, possibly changing and/or , is described by longdistance matrix elements (LDMEs), factorized from the partonlevel contributions. In NRQCD the LDMEs are supposed to be constant (i.e. independent of the momentum) and universal (i.e. processindependent). They cannot be calculated within the theoretical framework and need to be determined by comparisons to experimental data.
The NRQCD formalism was greeted with encouraging words given its seeming success in reproducing the CDF charmonia differential cross sections. One should bear in mind, however, that the colouroctet terms have free normalisations (essentially determined by the LDMEs), so that we should not be very surprised if the measured cross sections can be well described (the CSM, instead, had zero adjustable parameters).
A very reasonable way of evaluating if a given theory provides a suitable representation of reality is to fix its free parameters through fits to a given set of measurements and then check how well it predicts other physical observables, not previously considered. This procedure has been followed over the last years in quarkonium production studies: first the NRQCD LDMEs are determined by fitting the crosssection measurements to a superposition of singlet and octet terms; then the resulting model is used to predict the quarkonium polarizations. The outcome (see Ref. [5] and references therein) is that the predicted quarkonium polarizations are very different from the measured ones, a situation dubbed “the quarkonium polarization puzzle”.
It is worth emphasizing at this point that the detailed NRQCD modelling of quarkonium production crucially depends on the LDMEs, which are determined by the quality (and variety) of the available experimental measurements. This, by itself, is not an uncommon situation in highenergy physics. For instance, most QCD calculations require the use of datadriven parton distribution functions (PDFs) and fragmentation functions (FFs) to translate partonlevel calculations into predictions that can be compared to experimental (hadronlevel) data. Assuming that the PDFs and FFs are universal, one can measure them using suitable processes (like deepinelastic scattering and interactions) and then use them as inputs for calculating other measured processes.
But why is it that, in past quarkonium production studies, only the crosssection measurements were used to fit the NRQCD LDMEs? This is obviously not the only viable strategy. One could start by fitting the polarization measurements and then predict the differential cross sections, apart from their absolute normalizations; or, more democratically, one could make a global fit of both sets of measurements, cross sections and polarizations. The answer is that quarkonium polarization measurements are very complex and require exceptional care in the corresponding data analyses. Most of the quarkonium polarization results published before 2011 are incomplete and ambiguous [6]. Results obtained by the CDF and D0 Tevatron experiments, in particular, have been plagued by a series of suspicious observations, with at least two cases (CDF Run 1 versus CDF Run 2 for the J/ [7, 8] and CDF versus D0 for the [9, 10]) where two measurements mutually excluded each other. In these conditions, it is not surprising to see the fundamental role of polarization measurements purposely downgraded to an a posteriori crosscheck of the predictions.
The experimental situation has dramatically improved in the last years, first with the most recent CDF measurement of the three polarizations [11], and then with the LHC measurements, made by CMS for the five wave quarkonium states (two charmonia [12] and three bottomonia [13]) and by LHCb for the J/ [14]. All of these studies were made following the muchimproved methodologies proposed in a series of recent publications [6, 15, 16, 17]. The good mutual consistency of all these recent results reflects their vastly improved robustness with respect to the previous measurements.
These new results allow for a change of strategy. We can now proceed with “global fits” of quarkonium data, considering the polarization measurements at the same level as the cross sections. Actually, polarization is much more straightforwardly related to the variables of the theory than the momentum distributions, the different colour channels for production being characterized by simple and distinctive polarization patterns. This consideration guides our analysis, allowing us to make immediate simplifications and improve the robustness of the results.
To fully benefit from the improved quality and constraining power of the new measurements, efforts must be devoted to a careful treatment of the experimental and theoretical uncertainties. Correlations induced by systematic uncertainties due to luminosity measurements must be correctly taken into account. A more complex type of correlation is the one induced by the strong dependence of the acceptance determinations on the shape of the dilepton decay distributions. Previouslyreported NRQCD global fits use differential cross sections measured with acceptance corrections evaluated assuming unpolarized production, ignoring the large uncertainty that the experiments assign to the lack of prior knowledge about the quarkonium polarization. This imprecision must be corrected, at least to ensure logical consistency in cases leading to predictions of significantlypolarized quarkonia. As a further improvement, theoretical uncertainties must be modeled in the fitting algorithm, in order to allow for a realistic evaluation of the goodness of fit. This crucial aspect of the verification of the theory has not been properly addressed in previous analyses.
Furthermore, we should also revisit the very spirit and motivation of these fits. In past studies, measurements made at rather low , even lower than the mass of the quarkonium state, have been included in the NRQCD fits. At first sight, this might seem a good idea, because the lowest data points are usually the ones with the best statistical accuracy and, hence, are the ones that most strongly constrain the free parameters of the fits. However, the NRQCD factorization approach [4] requires that the shortdistance and longdistance phases of the quarkonium production process can be factorized, i.e. there must be a sharp separation between the typical distance scales of the hard scattering process, , and of the boundstate, , where is the heavyquark mass and its velocity in the rest frame. Effectively, this means that we cannot expect that the NRQCD calculations reproduce the measurements at low (especially given that the presently available calculations are limited to nexttoleading order in , NLO) and, hence, using those data to constrain the fitted parameters is a priori unjustified. Even worse, their high statistical accuracy might lead the fit into strongly biased results. When we compare data to theory, the pertinent question is not “Is NRQCD a valid theory for heavy quarkonium production?” but rather “Is there a kinematic domain in which NRQCD is a valid theory for heavy quarkonium production?”. To search for possible domains of validity of NRQCD factorization, at the present status of the perturbative calculations, in the description of charmonium and bottomonium production, we will apply progressively changing kinematic thresholds to the data and study the corresponding variations in the fit results.
The paper is organized as follows. After a preliminary description of the theory ingredients in Section 2, we present in Section 3 some basic considerations that drive our analysis, inspired by recurring patterns in the data and in particular in the polarization measurements. The study of the possible domain of validity of NLO NRQCD factorization based on LHC data is described in Section 4. The main results of the global analysis are presented in Section 5 and discussed in Section 6.
2 Theory ingredients
In the hypothesis of factorization of short and longdistance effects, the cross section for the inclusive production of the bound meson (plus unobserved particles, ) in a collision of initial systems and is expressed by the formula
(1) 
In each term of the sum, the (kinematicsdependent) shortdistance coefficient (SDC), , is proportional to the partonlevel cross section for the production of the preresonance in a given angular momentum (,) and colour () configuration, while the corresponding (constant) LDME, , is proportional to the probability of the boundstate formation. The theory ingredients of our analysis of and production in pp collisions at the LHC are the perturbative calculations of the SDCs for colouroctet and coloursinglet pairs, and the corresponding polarizations. The LDMEs are fit parameters determined from data. According to NRQCD scaling rules, the dominating contributions to the production of wave vector quarkonia are the coloursinglet () and three colouroctet (, and ) channels. We use the calculations made at NLO reported in Ref. [18], provided for a rest energy of the coloursinglet or colouroctet preresonance state GeV. Figure 1 shows the individual contributions (products of dependent SDCs times constant LDMEs) of the four colour configurations, and their sum, compared to the J/ cross section measured by CDF [19]. The LDMEs multiplying the octet SDCs have been obtained from a global fit of hadro and photoproduction data [20, 21].
Figure 2 illustrates how the individual SDCs change from LO to NLO. It is worth noting that the SDC, for above 7.5 GeV, changes from positive at LO to negative at NLO. Figure 3 shows the dependence of the polarization parameters , calculated at NLO for vector quarkonia produced in different colour configurations, where and () is the transverse (longitudinal) short distance cross section, in the helicity frame (HX). At LO, except for small deviations at low , vector quarkonia have either equal to (from or ) or to (from or ).
Among the colouroctet contributions, the shortdistance cross section raises attention for several peculiarities, calling for extra efforts in improved calculations. Firstly, both the yield and its polarization change drastically from LO to NLO. Secondly, at NLO they have unphysical behaviours, the yield being negative at low or high , depending on the sign of the corresponding LDME, and the polarization parameter reaching values higher than (and even diverging for a certain value of , when ). One might be tempted to argue that this is not a conceptual problem, since the cross section and polarization are not observable for each individual subprocess; only the sum over all subprocesses must lead to physicallymeaningful observables. However, at least in principle and at least in some phase space corner, the exact cancellation of the unphysical effects may be affected by approximations in the models (including the use of the model outside its domain of validity), by a not sufficiently accurate treatment of the experimental constraints on the theoretical parameters, or even simply by their statistical/systematic fluctuations. Fits relying on delicate compensations clearly demand special care.
3 Datadriven considerations
Our analysis is inspired and guided by two main datadriven considerations. The first is illustrated by Fig. 4, which shows the differential cross sections for the production of seven different quarkonium states, as measured by the ATLAS and CMS experiments [23, 24, 25, 26, 27]. We applied a mass rescaling to the variable in order to equalize the kinematic effects of different average parton momenta and phase spaces. When transformed to distributions, the shapes of the differential cross sections of these seven states are well described (at least for ) by a simple empirical function [22], with common values of its two shape parameters (the normalized of the global fit is 1.1 with 77 degrees of freedom).
The easiest conjecture explaining this common behaviour is that a very simple composition of processes, probably dominated by one single mechanism, is responsible for the production of all quarkonia. If several mechanisms were simultaneously at play, we would expect to see variations of their mixture because the differences in the masses of the component quarks and in the binding energy of the observed hadrons should induce changes in the nonperturbative effects. We must also keep in mind that the production kinematics addressed by these measurements differ from each other in that they contain almost pure wave ( and ) or wave ( and ) contributions or, because of feeddown effects, a mixture of the two (J/, , ). If confirmed with higher precision, the observed scaling would provide a strong physical indication without relying on explicit theoretical calculations. In fact, since the kinematics of coloursinglet processes is necessarily dependent on the angular momentum of the observed state, observing that states of different angular momentum quantum numbers are produced perturbatively with identical kinematics directly implies that coloursinglet processes play a negligible role.
The second, even stronger, hint comes from the quarkonium polarization measurements. As shown in Fig. 5, the polarizations of the wave quarkonia recently measured by CDF [11] and at the LHC [28, 12, 13, 14] cluster around the unpolarized limit, with no significant dependencies on or rapidity, no strong changes from directlyproduced states to those affected by wave feeddown decays, and no evident differences between charmonium and bottomonium. This observation strengthens the conjecture that, in “zeroorder” approximation, all quarkonia are dominantly produced by a single mechanism. Naturally, the polarization observable has an immediate interpretation in terms of angular momentum properties, especially strong given the peculiarity of the unpolarized result: the dominating channel must be the one leading to the “groundstate” preresonance object .
One may wonder whether this conclusion is in contradiction with NRQCD, given that, as seen in the previous section, current stateoftheart analyses, based on the fit of only distributions, point to a mixture of the , and channels definitely leading to transverse polarization [20, 21]. The answer is that the process mixture resulting from the fit depends in a dramatic way on the range where the fit is performed. We can have opposite physical indications when the data are fitted down to the lowest measured (the results being dominated by the more precise low data points) or when we assume a “validity domain” of the theory starting from a higher value. To illustrate this statement, Fig. 6 shows how the and SDCs for J/ production compare to the data when they are normalized to the lowest or highest point. While the low points are well described by the contribution and determine, therefore, a prediction of transverse polarization if included in the fit, at higher the data are closer in shape to the SDC: a fit starting at a value in the range 10–15 GeV would lead to a dominance of the unpolarized contribution. This observation illustrates the crucial importance of performing a scan of kinematic thresholds to search for a possible domain of validity of the theory. This procedure will be the subject of Section 4.
Despite oftenheard claims to the contrary, a careful look at Fig. 1 reveals that the fitted curve is a very unsatisfactory description of the measurements, given their rather small uncertainties. It is usually argued that theoretical uncertainties, not included in the fit, can cover the observed discrepancy, reconciling theory and data. However, as shown in Figs. 1 and 3, the octet is the only component that significantly changes from LO to NLO, in polarization and in the shape of the distribution (changes in normalization are absorbed in the LDMEs and do not affect the fit quality). Actually, judging from the difference between the LO and NLO calculations, the current theoretical uncertainty in the term is so large that by considering it in the fit we would introduce an excessive freedom, running the risk that this undetermined contribution would artificially absorb the datatheory discrepancy: the fit would improve its “mathematical quality” at the expense of losing all its physical impact.
We should also mention that, particularly in cases where a model does not describe faithfully the data, the fit can lead to meaningless and unstable results. It is helpful, at least as an initial step — in our case, the kinematic domain scan — to reduce the freedom of the fit to a minimum of essential parameters, with the aim of obtaining stable and univocal results in each tested condition. Besides its large uncertainty, the mathematical peculiarities of the SDC, mentioned in Section 2, represent a further danger to the robustness of the fit. Therefore, we will perform our domain scan considering only the , and components. More than a practical solution, this choice emerges from the previously discussed datadriven expectation that, within the domain of validity of the theory, the octet must be the dominating contribution; the term, with its unphysical polarization, can only represent a relatively small correction. In any case, the impact of this initial assumption will be tested a posteriori.
4 Kinematic domain scan
Our analysis considers a total of 121 data points, measured in pp collisions at 7 TeV by three LHC experiments: ATLAS ( cross sections [25]), CMS ( [24] and [26] cross sections; and [12] and [13] polarizations) and LHCb ( [29] and [30] cross sections). They correspond to data of GeV (43 points) and to data of GeV (78 points), including differential cross sections (99 points) and polarizations (, 22 points). We only consider and production because these states are not significantly affected by feeddown effects and can be treated as being directly produced. The description of the production of the remaining wave quarkonia contains a considerably higher number of free parameters, the LDMEs of the and states, and will be addressed more appropriately when additional constraints from detailed measurements of the wave cross sections and polarizations will become available.
The experiments have provided a thorough account of the dependence of each crosssection data point on the polarization hypothesis assumed in the acceptance determination. It is crucial to take this effect properly into account because it induces a strong correlation between the crosssection data points and the actual polarization prediction tested in the fit, thereby correlating the data points themselves. In our fit procedure, for each value of the explored parameter space (i.e., the and LDMEs of the and : four free parameters) we start by calculating the polarization prediction for each crosssection measurement; then we use this polarization value to recalculate the cross sections, using an acceptance correction taken from the tables provided by the experiments (with suitable interpolations). We also explicitly treat the pointtopoint correlations induced by the luminosity uncertainties in the crosssection measurements. To do this, for each crosssection data set we introduce a nuisance parameter representing a global rescaling of all points, constrained according to the relative luminosity uncertainty.
Concerning the theoretical ingredients, we use SDCs (and their longitudinal and transverse components) calculated for the production of a pair of GeV rest energy [18]. To obtain the shape of the SDC for the production of a object of rest energy equal to the mass of the considered quarkonium state, , we must rescale the variable by . It can be objected that this is not sufficient, because the rest energy of the pair, , is not necessarily equal to . However, we must also consider what happens, from the kinematic point of view, in the transition from the to the observable quarkonium, because of the emission or absorption of soft gluons. It can be shown that the average quarkonium threemomentum and the threemomentum (both in the laboratory) are related by the approximate expression . The approximation is excellent (corrections ) for of the order of the energy splitting between the radial and orbital angular momentum excitations of quarkonia. Towards midrapidity we can assume that, on average, . Even if, assuming that were known, we rescaled the by to obtain the observed quarkonium kinematics, we should then also scale the by . The net result of the two scalings is equivalent to one overall scaling by , which is, therefore, an asmuchaspossible accurate representation of the quarkonium production kinematics.
We complement the rescaling of the SDCs with a normalization rescaling exponentially depending on the quarkonium mass, approximately reflecting the normalization dependence of the measured distributions (Fig. 4). Since any normalization shift in the SDC is effectively reabsorbed in a rescaling of the corresponding LDME (except for the singlet term, which gives a negligible contribution), this choice has no influence on the fit quality nor on the results for the cross sections of the individual octet processes, . Only these cross sections, denoted by in what follows, can be directly compared among analyses; the LDMEs fitted from the data, for instance, differ by about a factor of two if we use the unscaled SDCs.
From a physical point of view, our redefinition of the SDCs equalizes the meaning of one given LDME among different states: two states of different mass but same value of are characterized, with this convention, by approximately the same probability of the transition; using the GeV convention for both and J/, for example, the two probabilities would differ by about a factor of two.
As explained in the previous sections, our central question is whether it is possible to define a kinematic domain, at sufficiently high , where current NRQCD calculations give a statistically satisfying description of the available data. The first step in our analysis is, therefore, a series of fits performed by selecting only the data points (for cross sections and polarizations, from all considered experiments) satisfying the selection , for a progressively changing choice of . Following the datadriven motivations given in Section 3, we consider the LDMEs and as the only two parameters of interest, assuming that is negligible. For stability reasons, we perform the scan without including a modelling of the theoretical uncertainties (which will be discussed in Section 5). In the absence of theoretical uncertainties, the fits to and data are essentially decorrelated and can be treated as two independent procedures.
Figure 8 shows how drastically the quality of the two fits change with varying , while Fig. 8 shows the corresponding behaviour of the fitted parameters.
Let us first consider the case. Above the normalized of the fit stops showing a decreasing trend and reaches a value of order 1 (the exact value being lower than one simply indicates the presence of correlations in the published pointtopoint systematic uncertainties and is not relevant for our studies). Correspondingly, the LDMEs cease to show any systematic trend above GeV and start a statisticallike oscillation around a common value, with obviously increasing uncertainty. This behaviour is typical of the stabilization of the fit results, when tensions between data and theory disappear and further rejection of data already removes points inside the domain of the theory. With this criterium, we can consider the results as well described by the theory for above GeV. Also in the scan we reach a small and rather stable normalized value, even if the trend of the LDMEs still leaves open the possibility that a complete stabilization may only happen at slightly higher values of .
While future data, extending with better precision towards higher , are needed for a conclusive statement, we do not expect a significant change of the physical conclusions, as can be appreciated from Fig. 9. The relative importance of the colouroctet cross section with respect to the total contribution of colouroctet processes, calculated at an arbitrary reference and midrapidity, saturates close to unity in the case, clearly indicating that the octet state dominates production, whereas the trend points to a more democratical share between the and contributions, at such high values.
Figure 9 also shows, for both quarkonia, that a fit performed by undiscriminatingly including all available data down to the lowest would lead, with high significance, to the opposite physical conclusion: that production is negligible and the wave cross sections are dominated by the (transversely polarized) contribution. As anticipated in Section 3, this conclusion, “traditionally” presented as a prediction of NRQCD, is in reality a result completely determined by the use of data not belonging to the domain of validity of the theory calculations.
We conclude this section by clarifying that our considerations would not be modified by the inclusion of photoproduction data, given that all such measurements are presently restricted to the low region, excluded by our study. Therefore, the hypothesis that the LDMEs are universal cannot be tested until precise photoproduction data will become available at high .
5 Results and predictions
Given the results shown in the previous section, we continue our analysis only using the 44 data points (30 cross sections and 14 polarizations) that belong to the kinematic domain GeV and GeV. These numerical values are clearly affected by some degree of arbitrariness and might have to be adjusted, at least in the case, when more precise high data will become available.
We start by addressing our datadriven assumption that we can neglect the contributions in the description of wave quarkonium production. Given the very good quality of the fits performed with , the current and are far from indicating the necessity of a nonvanishing wave octet component. Despite the caveats exposed in Sections 2 and 3 (in particular, we must be very critical regarding fits including this octet component, given its overwhelming theoretical uncertainty), we have repeated the fit with the additional free parameter . In the case, the central values of the fit results do not change, the fractional contributions to the octet cross sections at being (), () and (), fully consistent with our hypothesis and with the result shown in Fig. 9. The large and uncertainties are strongly anticorrelated. The fit becomes strongly underconstrained, still favouring, nevertheless, the octetcrosssection component (at ).
For a realistic evaluation of the LDMEs, we have included theoretical uncertainties in the fit procedure. For the and SDCs and polarizations, we assume as uncertainty (corresponding to a variation) the magnitude of the difference between NLO and LO calculations (shown in Section 2). In the case of the coloursinglet cross section and polarization, we keep the NLO calculation as central model but define the uncertainty as the difference with respect to the partial NNLO calculation of Ref. [31], except for the uncertainty in the cross section case, taken to be the LO model, to constrain the cross section to positive values. The nuisance parameters describing these allowed variations of the shortdistance ingredients are kept common to the and in one global fit, accounting for the correlation that the theoretical uncertainties induce between them.
Figure 10 shows the fitted data and the bestfit curves for the cross sections and polarizations, including the individual coloursinglet and colouroctet contributions. Uncertainty bands are also shown for the and cross section and polarization components, while the coloursinglet contribution is represented by the LO (dashed), NLO (dotdashed) and partial NNLO (dotted) calculations, together with the corresponding bestfit curve (solid), which is lower than the NLO calculation.
Figure 12 shows the and probability densities of the fitted and LDMEs, in the form of twodimensional contours, while Fig. 12 shows the corresponding distributions of the LDME ratio. Remarkably, the magnitudes of the two matrix elements are very different, in contradiction with usual expectations, as discussed in the next section.
Overall, in both the and cases, the octet is the dominating production channel. However, at very high the contribution seems to start prevailing in the case, as can be observed in the top left panel of Fig. 10, indicating that veryhigh mesons might be produced with a strong transverse polarization. This is clearly illustrated in Fig. 13, which shows the cross sections and polarizations corresponding to the fitted LDMEs, extrapolated to values well beyond the ranges probed by the existing measurements. These predictions were calculated using the twodimensional probabilitydensity function for the fit parameters, thereby taking into account parameter correlations and the modelled theoretical uncertainties, besides the experimental ones.
To put our results in the context of the existing literature, we stress that this is the first analysis specifically dedicated to the strategy for the comparison of existing theory calculations to measurements, with a datadriven attitude and a focus on the treatment of the experimental results. Previous “globalfit” analyses were, instead, byproducts of studies centred in the calculation of the shortdistance ingredients of NRQCD. In fact, remarkable efforts by a few groups triggered a great progress on this front over the last years, leading to full NLO crosssection and polarization calculations for different collision systems, energies and kinematic domains, for all relevant coloursinglet and colouroctet processes, including waves. The comparisons with data, however, have not followed detailed strategies and rigorous reproducibility criteria, so that it is often difficult to appreciate the consequences and prospects of the different fit approaches. It is sometimes impossible to understand the reasons of the differences in the fit results, which in some cases are very significant, giving the wrong overall impression that the NRQCD framework is either very unstable with respect to variations in the inputs or that it can at most give orderofmagnitude evaluations of cross sections and qualitative estimates of polarizations.
Some analyses include crosssection measurements down to GeV, others apply a fixed threshold GeV. With different data sets from analysis to analysis, it is difficult to quantify exactly how the choice affects the results and the quality of the theorydata agreement. In fact, the quantification of the agreement only addresses the cross sections or is not even reported. Moreover, since the polarization uncertainty correlations and luminosity uncertainties are never mentioned, one has to assume that they are neglected or assumed to be uncorrelated among different kinematic intervals, a choice that introduces an artificial freedom in the predicted shapes.
We take as examples the three recent analyses of prompt charmonium production reported in Refs. [20, 21] (A1), [32] (A2) and [33] (A3). A1 considers a large amount of J/ data from Tevatron, LHC, RHIC and photoproduction, using only distributions as constraints and neglecting the feeddown from decays. A2 only considers J/ data from CDF, for GeV, including the polarization as constraint and neglecting the feeddown. A3 uses CDF and LHCb data for J/, and , with GeV, excluding polarizations and including the modelling of the feeddown for the J/. The outcomes of A1 and A3 are, despite the very different strategies, substantially similar: a strong transverse polarization is predicted for directly produced wave charmonium in the range covered by the CMS measurements. A2 reproduces the unpolarized scenario by allowing a mutual cancellation of the transverse polarizations of the and crosssection terms, which are found to have opposite signs. However, the study suggests that no unique scenario can describe at the same time the CDF distributions ( GeV) and the LHC ones ( GeV). The latter are shown to lie entirely, with their error bars, between two curves, corresponding to the octet cross section term containing either 0% or 100% contribution of , which exclude the CDF bestfit result. It is concluded that the current LHC measurements lack constraining power on the parameter space, especially on the LDME. No quantification of the goodness of agreement is reported for any of the considered scenarios.
6 Discussion on the observed LDME hierarchies
The results of our study point to the existence of a much stronger hierarchy of LDMEs than the one predicted by the usual powercounting scheme of NRQCD, based on elegant and very general considerations on the formal structure of the NRQCD Lagrangian and operators. In NRQCD the three octet contributions , and are expected to scale in the same way for small values of the heavyquark velocity and, therefore, to have similar magnitudes. Two basic considerations compete in the determination of this result. On one hand, independently of the observed particle (wave or wave quarkonium), transitions from octet (or singlet) states with nonzero orbital angular momentum are suppressed, reflecting the fact that the perturbative state must be produced at shortdistance and small relative momentum. In particular, transitions from preresonance wave octet (and singlet) states are suppressed by a factor (coming from two additional spatial derivatives in the structure of the respective operators). On the other hand, the probability of softgluon emission depends on the process. The process is a chromoelectric transition (, ), while is a chromomagnetic transition (, ): their probabilities scale, respectively, like and . The process is predominantly a doublechromoelectric transition, its probability scaling like . These two considerations alone lead to the prediction that the three processes should have comparable probabilities (all scaling like ).
In order to understand why, instead, the data indicate the hierarchy , these rules must apparently be integrated with further conjectures on the mechanism of quarkonium formation. For example, the dependence of the interaction potential on the colour state of the quarkantiquark pair may play a role. With colour neutralization, the shortdistance potential changes from weakly repulsive, , to attractive, . Therefore, in the octettosinglet transition the pair undergoes a significant decrease in potential energy, , of the order of the kinetic energy of the bound state, i.e., of the energy splitting between radial and orbital angular momentum excitations of the quarkonium, –0.6 GeV (very similar for charmonium and bottomonium). Transitions in which the kinetic energy decreases () should therefore be disfavoured, because they require that the emitted soft gluons have comparatively high energy: . In particular, the transition, with GeV, should be suppressed, while the transition from , with GeV, would be the least subjected to the energy requirement on the gluon radiation. This sort of threshold effect may explain why and productions are dominated by the and octets.
Another fact worth of attention is that the measured suppression with respect to is not as strong for the as for the . This may possibly reflect the fact that the quark in a bottomonium state, having larger average momentum than the quark in a charmonium state, can emit higherenergy gluons. Clearly, this is only one possible conjecture. Alternative velocityscaling schemes [34, 35] also go in the direction of a better qualitative description of the measured LDME hierarchy, by reducing or eliminating the relative suppression of the chromomagnetic octettosinglet transition with respect to the chromoelectric one, therefore favouring the singleemission transition over the doubleemission transition . Incidentally, the different quality of the interaction potential for singlet and octet quarkantiquark pairs may also have a role in the observed dominance of octet processes over the singlet ones, given that the expansion of the initial “pointlike” towards boundstate sizes is energetically favoured when the shortdistance potential is repulsive rather than attractive.
These reasonings do not pretend to represent univocal explanations of the measured effects. They should be considered as illustrations of how the observation of definite scaling hierarchies for the LDMEs as a function of quarkonium mass, binding energy and quark flavour can have strong implications concerning the longdistance processes at play. Clarifying such hierarchies is one of the most stimulating reasons justifying accurate quarkonium production measurements at high, to be made at the LHC, so that we can pave the way towards a clearcut understanding of boundstate formation in QCD.
Concerning wave quarkonium production, the double chromoelectric transition is disfavoured with respect to the single one, besides being suppressed because of the higher angular momentum of the colouroctet state. Also for wave quarkonium production we expect, therefore, that the wave octet contribution is negligible. Among the remaining ones, the single chromoelectric transition should be favoured with respect to the double chromomagnetic+chromoelectric transition, but the relative importance of the two may be influenced by the enhancement of the latter and, a priori, both should be taken into consideration.
The possibility of a nonnegligible role of the contribution in production, in analogy with and production, is suggested by the two experimental facts discussed in Section 3: 1) the approximate universality of the scaling of wave quarkonium cross sections, indicating that states with a significant feeddown behave similarly to the others; 2) the absence of a clear polarization pattern differentiating directly produced states from those affected by a large feeddown. Future polarization measurements will be crucial to distinguish between the and contributions, respectively characterized by lack of polarization ( from ) or by moderate transverse polarizations (high and from have and , respectively, in the centreofmass helicity frame ^{1}^{1}1Values calculated following the method of Ref. [36], applied to the decay chain , , where has angular momentum projection or (transverse polarization), and assuming electricdipole gluon and photon radiations.).
7 Summary and conclusions
Nonrelativistic QCD, a rigorous and consistent effective theory based on QCD, should provide an accurate description of heavy quarkonium production. However, the efforts to validate NRQCD as a working framework have brought to light serious and persistent mismatches between data and calculations, especially concerning polarization. Recent CMS measurements of the polarizations of (directly produced) and have seemingly removed any residual ambiguity in this evidence.
We have addressed the “quarkonium production puzzle” through a deep reconsideration of the strategy for theorydata comparison. While the polarization data are traditionally excluded from global NRQCD analyses of quarkonium production (and used only as a posteriori verifications of the predictions), we argue that they are actually the most stringent and straightforward constraints in discriminating the underlying fundamental processes and we move them from the periphery to the centre of the study.
In fact, the measured unpolarized scenario points to a straightforward Occamrazor interpretation: the different colouroctet contributions to the wave quarkonium yield follow a magnitude hierarchy reflecting their degree of polarization. The unpolarized channel should dominate, while the one, with a polarization more transverse than what is physically allowed, should at most be a tiny correction. A small contribution, characterized by a fully transverse (but physical) polarization, would be sufficient to explain the possible tendency of the measured polarizations towards slightly transverse values at higher .
The data show another interesting pattern: the differential cross sections of seven quarkonium states are compatible with a common scaling, at least for . Given that these quarkonia include two essentially pure wave states ( and ), three wave states affected by a significant feeddown from wave states (J/, and ) and two wave states ( and ), their very similar behaviour suggests that quarkonium production is the result of a simple mixture of processes, stable with varying mass and quantum numbers. This observation clearly favours a scenario where one single process dominates and, together with the polarization argument, makes it even less reasonable to consider that the unpolarized measurements could be the result of a delicate cancellation of strongly polarized processes. Furthermore, given that both the and octets are transversely polarized, their mutual cancellation implies that one of them needs to contribute with a negative cross section.
These datadriven considerations guide us in our global fit of LHC measurements of and cross sections and polarizations. Having a prior expectation of what a reasonable result will be helps us avoiding the pitfalls of illposed, underconstrained or unstable fits. By excluding polarization data from the fits, previous analyses have effectively chosen to restrict the safe domain of the theory to the description of the unpolarized crosssection observables. We propose a different definition of field of validity, including polarization observables as crucial players while possibly excluding the lowest data, knowing that fixedperturbativeorder factorization calculations are supposed to work only at sufficiently high . The systematic search for the domain of validity of the theory through a scan of the kinematic phase space is a crucial step in our analysis.
For the first time in this kind of studies, we perform a rigorous treatment of correlated experimental uncertainties, including the dependence of experimental acceptances on the polarizations. Once a candidate domain of validity is defined, we also include in the fit a modelling of the theoretical uncertainties, so that they are reflected in the output parameters. This effectively introduces a partial correlation between the and systems: charmonium and bottomonium fits become one global quarkonium fit.
Bringing the polarization data to the centre of the stage and decreasing the (statistically strongest) weight of the low data is a “Copernican revolution” that seems to provide a straightforward solution to the puzzle: the cross sections and the polarizations are both perfectly fitted by the theory in a domain approximately defined by the selection cut . Confirming our initial expectation, no wave component is needed to describe the data. We also find that the data favour a coloursinglet component smaller than the NLO calculation and even ten times smaller than the partial NNLO calculation.
These facts, together with the further hierarchy , are physically intriguing and are to be interpreted as strong indications for the understanding of the mechanisms of boundstate formation (an example of such interpretations being presented in Section 6). Furthermore, finding that the octet term gives a negligible contribution to quarkonium production is also extremely interesting from another perspective. In fact, contrary to what happens in the wave octet case, the SDCs of the dominant wave octet components have a very stable shape from LO to NLO, indicating that, at the present status of the perturbative calculations, the theoretical uncertainties in the framework are relatively small. This points to a great potential of NRQCD as a precision instrument to address and isolate the intriguing aspect of the process, the formation of the bound state, as described by the nonperturbative LDMEs. If these observations are confirmed by future data, the LHC measurements will provide precise determinations of the LDMEs of all quarkonium states in a consistent framework. On the other hand, we must call attention to the fact that the existing photoproduction data belong to the kinematic domain that our study has excluded. Therefore, a test of the universality of the LDMEs must wait for precise high measurements in processes different from direct production in pp collisions.
We must also mention that, while the distributions for cannot be described at NLO using the same process mixture implied by higher data, they are, nevertheless, still compatible with the zeropolarization pattern, smoothly continuing the high trend (see Fig. 5). In other words, there is no indication from data alone of a change in production mechanism from high to low . The implied dominance of quarkonium production via an intermediate isotropic wave function (presumably ) finds its simplest explanation in one of the crucial aspects of the factorization concept: the quantum numbers of the produced change during the boundstate formation, making it possible that, for example, a quarkonium exhibits a distinctive polarization pattern. At the same time, the indication comes invariably from low and high data and is, therefore, more “universal” than the validity of the current factorized NLO calculation, as established by the results of our high fits. Furthermore, the polarization is zero at all perturbative orders and the factorization prediction for the production via is, obviously, unpolarized when resummed to all orders in any kind of perturbative expansion, in agreement with data down to low . This leaves open the possibility that factorized calculations may describe simultaneously high and low data, if higher perturbative orders improve the description (specifically, by reducing the steepness of the distribution at low ). Also in this case, polarization data show their power in driving us towards encouraging indications on the reliability of the NRQCD factorization framework.
Finally, we have also extrapolated the fitted cross sections and polarizations to very high , providing predictions for future LHC measurements; the term seems to become more important, at least for the , increasing the fraction of transversely polarized mesons.
We congratulate our colleagues from the LHC collaborations for their highquality quarkonium production measurements, without which this study could not have been made. We acknowledge very interesting discussions with Geoff Bodwin and Sergey Baranov. We are indebted to Mathias Butenschön and Bernd Kniehl, who kindly gave us their NLO calculations of the short distance coefficients. The work of P.F. is supported by FCT, Portugal, through the grant SFRH/BPD/42343/2007, while the work of V.K. is supported by FWF, Austria, through the grant P24167N16.
References

[1]
R. Baier and R. Rückl,
Z. Phys. C19 (1983) 251.
E.W.N. Glover, A.D. Martin, W.J. Stirling, Z. Phys. C38 (1988) 473; Erratumibid. C49 (1991) 526. 
[2]
M.H. Schub et al. (E789 Coll.),
Phys. Rev. D52 (1995) 1307;
Erratumibid. D53 (1996) 570.
P.L. McGaughey, Nucl. Phys. A610 (1996) 394c. 
[3]
A. Sansoni et al. (CDF Coll.),
Nucl. Phys. A610 (1996) 373c.
F. Abe et al. (CDF Coll.), Phys. Rev. Lett. 79 (1997) 572.  [4] G.T. Bodwin, E. Braaten, G.P. Lepage, Phys. Rev. D51 (1995) 1125; Erratumibid. D55 (1997) 5853.
 [5] N. Brambilla et al. (QWG Coll.), Eur. Phys. J. C71 (2011) 1534.
 [6] P. Faccioli, C. Lourenço, J. Seixas and H.K. Wöhri, Eur. Phys. J. C69 (2010) 657.
 [7] T. Affolder et al. (CDF Coll.), Phys. Rev. Lett. 85 (2000) 2886.
 [8] A. Abulencia et al. (CDF Coll.), Phys. Rev. Lett. 99 (2007) 132001.
 [9] D. Acosta et al. (CDF Coll.), Phys. Rev. Lett. 88 (2002) 161802.
 [10] V.M. Abazov et al. (D0 Coll.), Phys. Rev. Lett. 101 (2008) 182004.
 [11] T. Aaltonen et al. (CDF Coll.), Phys. Rev. Lett. 108 (2012) 151802.
 [12] S. Chatrchyan et al. (CMS Coll.), Phys. Lett. B727 (2013) 381.
 [13] S. Chatrchyan et al. (CMS Coll.), Phys. Rev. Lett. 110 (2013) 081802.
 [14] R. Aaij et al. (LHCb Coll.), Eur. Phys. J. C73 (2013) 2631.
 [15] P. Faccioli, C. Lourenço, J. Seixas and H.K. Wöhri, Phys. Rev. Lett. 102 (2009) 151802.
 [16] P. Faccioli, C. Lourenço and J. Seixas, Phys. Rev. Lett. 105 (2010) 061601.
 [17] P. Faccioli, C. Lourenço and J. Seixas, Phys. Rev. D81 (2010) 111502(R).
 [18] M. Butenschön, B.A. Kniehl, Phys. Rev. Lett. 108 (2012) 172002 and private communication.
 [19] D. Acosta et al. (CDF Coll.), Phys. Rev. D71 (2005) 032001.
 [20] M. Butenschön, B.A. Kniehl, Nucl. Phys. B (Proc. Suppl.) 222–224 (2012) 151.
 [21] M. Butenschön, B.A. Kniehl, Mod. Phys. Lett. A 28 (2013) 1350027.
 [22] I. Abt et al. (HERAB Coll.), Eur. Phys. J. C60 (2009) 525.
 [23] G. Aad et al. (ATLAS Coll.), Nucl. Phys. B850 (2011) 387.
 [24] S. Chatrchyan et al. (CMS Coll.), JHEP 02 (2012) 011.
 [25] G. Aad et al. (ATLAS Coll.), Phys. Rev. D87 (2013) 052004.
 [26] CMS Coll., CMSPASBPH12006.
 [27] ATLAS Coll., ATLASCONF2013095.
 [28] B. Abelev et al. (ALICE Coll.), Phys. Rev. Lett. 108 (2012) 082001.
 [29] R. Aaij et al. (LHCb Coll.), Eur. Phys. J. C72 (2012) 2100.
 [30] R. Aaij et al. (CMS Coll.), Eur. Phys. J. C72 (2012) 2025.
 [31] J.P. Lansberg, Phys. Lett. B679 (2009) 340.
 [32] K.T. Chao, Y.Q. Ma, H.S. Shao, K. Wang and Y.J. Zhang, Phys. Rev. Lett. 108 (2012) 242004.
 [33] B. Gong, L.P. Wan, J.X. Wang and H.F. Zhang, Phys. Rev. Lett. 110 (2013) 042002.
 [34] G.A. Schuler, Int. J. Mod. Phys. A12 (1997) 3951.
 [35] S. Fleming, I.Z. Rothstein and A.K. Leibovich, Phys. Rev. D64 (2001) 036002.
 [36] P. Faccioli, C. Lourenço, J. Seixas and H.K. Wöhri, Phys. Rev. D83 (2011) 096001.