Quarkonium production in the LHC era: a polarized perspective

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 non-perturbative parameters of bound-state 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 psi(2S) and Y(3S) measurements from LHC experiments and state-of-the-art NLO short-distance calculations in the framework of non-relativistic 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 non-perturbative NRQCD parameters, that open new paths towards the understanding of bound-state formation in QCD.

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 high-energy physics. For instance, most QCD calculations require the use of datadriven parton distribution functions (PDFs) and fragmentation functions (FFs) to translate parton-level 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 deep-inelastic scattering and e + e − interactions) and then use them as inputs for calculating other measured processes.
But why is it that, in past quarkonium production studies, only the cross-section 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 Υ(1S) [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 Υ(nS) polarizations [11], and then with the LHC measurements, made by CMS for the five S-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 much-improved 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 QQ 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. Previously-reported 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 significantly-polarized 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 p T , 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 p T 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 short-distance and long-distance phases of the quarkonium production process can be factorized, i.e. there must be a sharp separation between the typical QQ distance scales of the hard scattering process, ∼ max(1/p T , 1/m Q ), and of the bound-state, ∼ 1/(m Q v), where m Q is the heavy-quark mass and v its velocity in the QQ rest frame. Effectively, this means that we cannot expect that the NRQCD calculations reproduce the measurements at low p T (especially given that the presently available calculations are limited to next-to-leading order in α s , 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.

Theory ingredients
In the hypothesis of factorization of short-and long-distance effects, the cross section for the inclusive production of the bound meson H (plus unobserved particles, X) in a collision of initial systems A and B is expressed by the formula In each term of the sum, the (kinematics-dependent) short-distance coefficient (SDC), S, is proportional to the parton-level cross section for the production of the pre-resonance QQ in a given angular momentum (L,S) and colour (C) configuration, while the corresponding (constant) LDME, L, is proportional to the probability of the bound-state formation. The theory ingredients of our analysis of ψ(2S) and Υ(3S) production in pp collisions at the LHC are the perturbative calculations of the SDCs for colour-octet and colour-singlet QQ pairs, and the corresponding polarizations. The LDMEs are fit parameters determined from data. According to NRQCD v-scaling rules, the dominating contributions to the production of S-wave vector quarkonia are the colour-singlet ( 3 S J ) channels. We use the calculations made at NLO reported in Ref. [18], provided for a rest energy of the colour-singlet or colour-octet pre-resonance state E QQ = 3 GeV. Figure 1 shows the individual contributions (products of p T -dependent SDCs times constant LDMEs) of the four QQ 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 photo-production data [20,21].  [18] for J/ψ production, scaled by LDMEs fitted [20,21] to the CDF data [19].  J SDC, for p T above 7.5 GeV, changes from positive at LO to negative at NLO. Figure 3 shows the p T dependence of the polarization parameters λ ϑ , calculated at NLO for vector quarkonia produced in different QQ colour configurations, where λ ϑ = (S T − S L )/(S T + S L ) and S T (S L ) is the transverse (longitudinal) short distance cross section, in the helicity frame (HX). At LO, except for small deviations at low p T , vector quarkonia have λ ϑ either equal to +1 (from 3 S       Figure 3: Polarizations parameters, λ ϑ , calculated at NLO for the colour-singlet term and for the three colour-octet terms [18].
Among the colour-octet contributions, the 3 P

[8]
J short-distance cross section raises attention for several peculiarities, calling for extra efforts in improved calculations. Firstly, both the QQ 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 p T , depending on the sign of the corresponding LDME, and the polarization parameter λ ϑ reaching values higher than +1 (and even diverging for a certain value of p T , when S T = −S L ). 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 physically-meaningful 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.

Data-driven considerations
Our analysis is inspired and guided by two main data-driven 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 p T variable in order to equalize the kinematic effects of different average parton momenta and phase spaces. When transformed to p T /M distributions, the shapes of the differential cross sections of these seven states are well described (at least for p T /M > 3) by a simple empirical function [22], with common values of its two shape parameters (the normalized χ 2 of the global fit is 1.1 with 77 degrees of freedom).  Figure 4: Mid-rapidity quarkonium p T /M distributions measured at √ s = 7 TeV by the ATLAS and CMS experiments [23,24,25,26,27]. The solid curve is a fit to the J/ψ data of CMS (for p T /M > 3), using a power-law function [22], while the dashed curves are replicas with normalizations adjusted to the individual datasets.
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 non-perturbative 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 S-wave (ψ(2S) and Υ(3S)) or P -wave (χ c1 and χ c2 ) contributions or, because of feed-down effects, a mixture of the two (J/ψ, Υ(1S), Υ(2S)). If confirmed with higher precision, the observed p T /M scaling would provide a strong physical indication without relying on explicit theoretical calculations. In fact, since the kinematics of colour-singlet 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 colour-singlet processes play a negligible role.  [11], ALICE [28], CMS [12,13] and LHCb [14].
The second, even stronger, hint comes from the quarkonium polarization measurements. As shown in Fig. 5, the polarizations of the S-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 p T or rapidity, no strong changes from directly-produced states to those affected by P -wave feed-down decays, and no evident differences between charmonium and bottomonium. This observation strengthens the conjecture that, in "zero-order" 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 "ground-state" pre-resonance object 1 S [8] 0 . One may wonder whether this conclusion is in contradiction with NRQCD, given that, as seen in the previous section, current state-of-the-art analyses, based on the fit of only p T distributions, point to a mixture of the 1 S 1 and 1 S [8] 0 SDCs of J/ψ production [18] normalized to the first data point (dashed lines) or to the last data point (solid lines) of the CDF measurement [19].
transverse polarization [20,21]. The answer is that the process mixture resulting from the fit depends in a dramatic way on the p T range where the fit is performed. We can have opposite physical indications when the data are fitted down to the lowest measured p T (the results being dominated by the more precise low-p T data points) or when we assume a "validity domain" of the theory starting from a higher p T value. To illustrate this statement, Fig. 6 shows how the 3 S [8] 1 and 1 S [8] 0 SDCs for J/ψ production compare to the data when they are normalized to the lowest-or highest-p T point. While the low-p T points are well described by the 3 S [8] 1 contribution and determine, therefore, a prediction of transverse polarization if included in the fit, at higher p T the data are closer in shape to the 1 S [8] 0 SDC: a fit starting at a p T value in the range 10-15 GeV would lead to a dominance of the unpolarized 1 S [8] 0 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 often-heard 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 3 P [8] J octet is the only component that significantly changes from LO to NLO, in polarization and in the shape of the p T 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 3 P [8] J 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 data-theory 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 3 P [8] J 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 3 S 1 and 1 S [8] 0 components. More than a practical solution, this choice emerges from the previously discussed data-driven expectation that, within the domain of validity of the theory, the 1 S octet must be the dominating contribution; the 3 P [8] J 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.

Kinematic domain scan
Our analysis considers a total of 121 data points, measured in pp collisions at 7 TeV by three LHC experiments: ATLAS (Υ(3S) cross sections [25]), CMS (ψ(2S) [24] and Υ(3S) [26] cross sections; and ψ(2S) [12] and Υ(3S) [13] polarizations) and LHCb (ψ(2S) [29] and Υ(3S) [30] cross sections). They correspond to ψ(2S) data of p T > 4 GeV (43 points) and to Υ(3S) data of p T > 10 GeV (78 points), including p T -differential cross sections (99 points) and polarizations (λ ϑ , 22 points). We only consider ψ(2S) and Υ(3S) production because these states are not significantly affected by feed-down effects and can be treated as being directly produced. The description of the production of the remaining S-wave quarkonia contains a considerably higher number of free parameters, the LDMEs of the χ 1 and χ 2 states, and will be addressed more appropriately when additional constraints from detailed measurements of the P -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 cross-section 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 1 S [8] 0 and 3 S [8] 1 LDMEs of the ψ(2S) and Υ(3S): 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 point-to-point correlations induced by the luminosity uncertainties in the cross-section 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 QQ pair of E 0 = 3 GeV rest energy [18]. To obtain the shape of the SDC for the production of a QQ object of rest energy equal to the mass of the considered quarkonium state, M , we must rescale the p T variable by M/E 0 . It can be objected that this is not sufficient, because the rest energy of the QQ pair, E QQ , is not necessarily equal to M . However, we must also consider what happens, from the kinematic point of view, in the transition from the QQ to the observable quarkonium, because of the emission or absorption of soft gluons. It can be shown that the average quarkonium threemomentum p and the QQ three-momentum P (both in the laboratory) are related by the approximate expression p 2 /P 2 M 2 /E 2 QQ . The approximation is excellent (corrections ≤ 2%) for |E QQ − M | of the order of the energy splitting between the radial and orbital angular momentum excitations of quarkonia. Towards mid-rapidity we can assume that, on average, p T /p QQ T M/E QQ . Even if, assuming that E QQ were known, we rescaled the QQ p T by E QQ /E 0 to obtain the observed quarkonium kinematics, we should then also scale the p T by M/E QQ . The net result of the two scalings is equivalent to one overall scaling by M/E 0 , which is, therefore, an as-much-as-possible accurate representation of the quarkonium production kinematics.
We complement the p T rescaling of the SDCs with a normalization rescaling exponentially depending on the quarkonium mass, approximately reflecting the normalization dependence of the measured p T /M distributions (Fig. 4). Since any normalization shift in the SDC S is effectively reabsorbed in a rescaling of the corresponding LDME L (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, σ(A + B → QQ[ 2S+1 L [8] J ] → H + X) = S(A + B → QQ[ 2S+1 L [8] J ] + X) × L(QQ[ 2S+1 L [8] J ] → H). Only these cross sections, denoted by σ( 2S+1 L [8] J ) in what follows, can be directly compared among analyses; the LDMEs fitted from the ψ(2S) 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 L(QQ[ 2S+1 L [8] J ] → H) are characterized, with this convention, by approximately the same probability of the QQ[ 2S+1 L [8] J ] → H transition; using the E 0 = 3 GeV convention for both ψ(2S) 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 p T , 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 p T > p min T , for a progressively changing choice of p min T . Following the data-driven motivations given in Section 3, we consider the LDMEs L( 1 S [8] 0 → H) and L( 3 S [8] 1 → H) as the only two parameters of interest, assuming that L( 3 P [8] J → H) is negligible. For stability reasons, we perform the p min T 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 ψ(2S) and Υ(3S) data are essentially decorrelated and can be treated as two independent procedures.  Let us first consider the Υ(3S) case. Above p min T /M ∼ 3 the normalized χ 2 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 point-to-point systematic uncertainties and is not relevant for our studies). Correspondingly, the LDMEs cease to show any systematic trend above p min T ∼ 30 GeV and start a statistical-like 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 Υ(3S) results as well described by the theory for p T above ∼ 30 GeV. Also in the ψ(2S) scan we reach a small and rather stable normalized χ 2 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 p min T . While future data, extending with better precision towards higher p T , 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 σ( 1 S [8] 0 ) colour-octet cross section with respect to the total contribution of colour-octet processes, calculated at an arbitrary reference p T /M = 6 and mid-rapidity, saturates close to unity in the ψ(2S) case, clearly indicating that the 1 S [8] 0 octet state dominates ψ(2S) production, whereas the Υ(3S) trend points to a more democratical share between the 1 S [8] 0 and 3 S [8] 1 contributions, at such high p T values. Figure 9 also shows, for both quarkonia, that a fit performed by undiscriminatingly including all available data down to the lowest p T would lead, with high significance, to the opposite physical conclusion: that 1 S [8] 0 production is negligible and the S-wave cross sections are dominated by the (transversely polarized) 3 S [8] 1 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-p T 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 p T .

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 p Υ(3S) T > 30 GeV and p ψ(2S) T > 12 GeV. These numerical values are clearly affected by some degree of arbitrariness and might have to be adjusted, at least in the ψ(2S) case, when more precise high-p T data will become available.
We start by addressing our data-driven assumption that we can neglect the 3 P

[8]
J contributions in the description of S-wave quarkonium production. Given the very good quality of the fits performed with L( 3 P [8] J ) = 0, the current Υ(3S) and ψ(2S) are far from indicating the necessity of a non-vanishing P -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 L( 3 P  1 SDCs and polarizations, we assume as uncertainty (corresponding to a ±1σ variation) the magnitude of the difference between NLO and LO calculations (shown in Section 2). In the case of the colour-singlet 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 −1σ 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 Υ(3S) and ψ(2S) in one global fit, accounting for the correlation that the theoretical uncertainties induce between them. Figure 10 shows the fitted data and the best-fit curves for the cross sections and polarizations, including the individual colour-singlet and colour-octet contributions. Uncertainty bands are also shown for the 1 S 1 cross section and polarization components, while the colour-singlet contribution is represented by the LO (dashed), NLO (dot-dashed) and partial NNLO (dotted) calculations, together with the corresponding best-fit curve (solid), which is lower than the NLO calculation. Figure 11 shows the ψ(2S) and Υ(3S) probability densities of the fitted 1 S LDMEs, in the form of two-dimensional contours, while Fig. 12 shows the corresponding distributions of the 3 S 1 / 1 S [8] 0 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 ψ(2S) and Υ(3S) cases, the 1 S [8] 0 octet is the dominating production channel. However, at very high p T the 3 S [8] 1 contribution seems to start prevailing in the Υ(3S) case, as can be observed in the top left panel of Fig. 10   p T Υ(3S) 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 p T values well beyond the ranges probed by the existing measurements. These predictions were calculated using the two-dimensional 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 data-driven attitude and a focus on the treatment of the experimental results. Previous "global-fit" analyses were, instead, byproducts of studies centred in the calculation of the short-distance 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 cross-section and polarization calculations for different collision systems, energies and kinematic domains, for all relevant colour-singlet and colour-octet processes, including P -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 order-of-magnitude evaluations of cross sections and qualitative estimates of polarizations. Some analyses include cross-section measurements down to p T = 3 GeV, others apply a fixed threshold p T > 7 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 theory-data 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 re-ported in Refs. [20,21] (A1), [32] (A2) and [33] (A3). A1 considers a large amount of J/ψ data from Tevatron, LHC, RHIC and photo-production, using only p T distributions as constraints and neglecting the feed-down from χ c decays. A2 only considers J/ψ data from CDF, for p T > 7 GeV, including the polarization as constraint and neglecting the χ c feeddown. A3 uses CDF and LHCb data for J/ψ, ψ(2S) and χ c , with p T > 7 GeV, excluding polarizations and including the modelling of the feed-down 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 S-wave charmonium in the p T range covered by the CMS measurements. A2 reproduces the unpolarized scenario by allowing a mutual cancellation of the transverse polarizations of the 3 S J cross-section terms, which are found to have opposite signs. However, the study suggests that no unique scenario can describe at the same time the CDF p T distributions (7 < p T < 20 GeV) and the LHC ones (7 < p T < 70 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 1 S [8] 0 , which exclude the CDF best-fit result. It is concluded that the current LHC measurements lack constraining power on the parameter space, especially on the 1 S 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 power-counting 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 1 S [8] 0 , 3 S [8] 1 and 3 P [8] J are expected to scale in the same way for small values of the heavy-quark velocity v and, therefore, to have similar magnitudes. Two basic considerations compete in the determination of this result. On one hand, independently of the observed particle (S-wave or P -wave quarkonium), transitions from octet (or singlet) states with non-zero orbital angular momentum are suppressed, reflecting the fact that the perturbative QQ state must be produced at short-distance and small relative momentum. In particular, transitions from pre-resonance P -wave octet (and singlet) states are suppressed by a factor v 2 (coming from two additional spatial derivatives in the structure of the respective operators). On the other hand, the probability of soft-gluon emission depends on the process. The 3 P [8] J → ψ/Υ + g process is a chromoelectric transition (∆L = ±1, ∆S = 0), while 1 S [8] 0 → ψ/Υ + g is a chromomagnetic transition (∆L = 0, ∆S = ±1): their probabilities scale, respectively, like v 2 and v 4 . The 3 S [8] 1 → ψ/Υ + gg process is predominantly a double-chromoelectric transition, its probability scaling like v 4 . These two considerations alone lead to the prediction that the three processes should have comparable probabilities (all scaling like v 4 ).
In order to understand why, instead, the data indicate the hierarchy 3 P 0 , 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 quark-antiquark pair may play a role. With colour neutralization, the short-distance potential changes from weakly repulsive, V 8 ≥ 0, to attractive, V 1 < 0. Therefore, in the octet-to-singlet transition the QQ pair undergoes a significant decrease in potential energy, ∆V (8 → 1) V 1 −T , of the order of the kinetic energy T of the bound state, i.e., of the energy splitting between radial and orbital angular momentum excitations of the quarkonium, ∼ 0.4-0.6 GeV (very similar for charmonium and bottomonium). Transitions in which the QQ kinetic energy decreases (∆T < 0) should therefore be disfavoured, because they require that the emitted soft gluons have comparatively high energy: E g = |∆V | − ∆T . In particular, the 3 P 1 suppression with respect to 1 S [8] 0 is not as strong for the Υ(3S) as for the ψ(2S). This may possibly reflect the fact that the b quark in a bottomonium state, having larger average momentum than the c quark in a charmonium state, can emit higher-energy gluons. Clearly, this is only one possible conjecture. Alternative velocity-scaling 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 octet-to-singlet transition with respect to the chromoelectric one, therefore favouring the single-emission transition 1 S [8] 0 over the doubleemission transition 3 S [8] 1 . Incidentally, the different quality of the interaction potential for singlet and octet quark-antiquark pairs may also have a role in the observed dominance of octet processes over the singlet ones, given that the expansion of the initial "point-like" QQ towards bound-state sizes is energetically favoured when the short-distance 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 long-distance processes at play. Clarifying such hierarchies is one of the most stimulating reasons justifying accurate quarkonium production measurements at high-p T , to be made at the LHC, so that we can pave the way towards a clear-cut understanding of bound-state formation in QCD.
Concerning P -wave quarkonium production, the double chromoelectric transition 3 P

[8]
J → χ+gg is disfavoured with respect to the single 3 S [8] 1 → χ+g one, besides being suppressed because of the higher angular momentum of the colour-octet state. Also for P -wave quarkonium production we expect, therefore, that the P -wave octet contribution is negligible. Among the remaining ones, the single chromoelectric transition 3 S [8] 1 → χ + g should be favoured with respect to the double chromomagnetic+chromoelectric 1 S [8] 0 → χ + gg transition, but the relative importance of the two may be influenced by the ∆T > 0 enhancement of the latter and, a priori, both should be taken into consideration.
The possibility of a non-negligible role of the 1 S [8] 0 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 p T /M scaling of S-wave quarkonium cross sections, indicating that states with a significant χ feed-down behave similarly to the others; 2) the absence of a clear polarization pattern differentiating directly produced states from those affected by a large χ feed-down. Future χ polarization measurements will be crucial to distinguish between the 1 S [8] 0 and 3 S [8] 1 contributions, respectively characterized by lack of polarization (χ from 1 S [8] 0 ) or by moderate transverse polarizations (high-p T χ 1 and χ 2 from 3 S [8] 1 have λ ϑ = +1/5 and +21/73, respectively, in the centre-of-mass helicity frame 1 ).

Summary and conclusions
Non-relativistic 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) ψ(2S) and Υ(3S) have seemingly removed any residual ambiguity in this evidence.
We have addressed the "quarkonium production puzzle" through a deep reconsideration of the strategy for theory-data 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 Occam-razor interpretation: the different colour-octet contributions to the S-wave quarkonium yield follow a magnitude hierarchy reflecting their degree of polarization. The unpolarized 1 S [8] 0 channel should dominate, while the 3 P [8] J one, with a polarization more transverse than what is physically allowed, should at most be a tiny correction. A small 3 S [8] 1 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 p T .
The data show another interesting pattern: the differential cross sections of seven quarkonium states are compatible with a common p T /M scaling, at least for p T /M > 3. Given that these quarkonia include two essentially pure S-wave states (ψ(2S) and Υ(3S)), three S-wave states affected by a significant feed-down from P -wave states (J/ψ, Υ(1S) and Υ(2S)) and two P -wave states (χ c1 and χ c2 ), 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 3 P 1 octets are transversely polarized, 1 Values calculated following the method of Ref. [36], applied to the decay chain 3 S [8] 1 → χ 1,2 + g, χ 1,2 → ψ/Υ + γ, where 3 S [8] 1 has angular momentum projection J z = +1 or −1 (transverse polarization), and assuming electric-dipole gluon and photon radiations. their mutual cancellation implies that one of them needs to contribute with a negative cross section.
These data-driven considerations guide us in our global fit of LHC measurements of ψ(2S) and Υ(3S) cross sections and polarizations. Having a prior expectation of what a reasonable result will be helps us avoiding the pitfalls of ill-posed, under-constrained 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 cross-section observables. We propose a different definition of field of validity, including polarization observables as crucial players while possibly excluding the lowest-p T data, knowing that fixed-perturbative-order factorization calculations are supposed to work only at sufficiently high p T . 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 ψ(2S) and Υ(3S) 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-p T 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 p T /M > 3. Confirming our initial expectation, no P -wave component is needed to describe the data. We also find that the data favour a colour-singlet component smaller than the NLO calculation and even ten times smaller than the partial NNLO calculation.
These facts, together with the further hierarchy L( 3 S 0 ), are physically intriguing and are to be interpreted as strong indications for the understanding of the mechanisms of bound-state formation (an example of such interpretations being presented in Section 6). Furthermore, finding that the 3 P [8] J octet term gives a negligible contribution to quarkonium production is also extremely interesting from another perspective. In fact, contrary to what happens in the P -wave octet case, the SDCs of the dominant S-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 non-perturbative 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 photo-production 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-p T measurements in processes different from direct production in pp collisions.
We must also mention that, while the p T distributions for p T /M < 3 cannot be described at NLO using the same process mixture implied by higher-p T data, they are, nevertheless, still compatible with the zero-polarization pattern, smoothly continuing the high-p T trend (see Fig. 5). In other words, there is no indication from data alone of a change in production mechanism from high to low p T . The implied dominance of quarkonium production via an intermediate isotropic wave function (presumably 1 S [8] 0 ) finds its simplest explanation in one of the crucial aspects of the factorization concept: the quantum numbers of the produced QQ change during the bound-state formation, making it possible that, for example, a J = 1 quarkonium exhibits a distinctive J = 0 polarization pattern. At the same time, the indication comes invariably from low-and high-p T data and is, therefore, more "universal" than the validity of the current factorized NLO calculation, as established by the results of our high-p T fits. Furthermore, the 1 S [8] 0 polarization is zero at all perturbative orders and the factorization prediction for the production via 1 S [8] 0 is, obviously, unpolarized when resummed to all orders in any kind of perturbative expansion, in agreement with data down to low p T . This leaves open the possibility that factorized calculations may describe simultaneously high-and low-p T data, if higher perturbative orders improve the p T description (specifically, by reducing the steepness of the 1 S [8] 0 p T distribution at low p T ). 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 p T , providing predictions for future LHC measurements; the 3 S [8] 1 term seems to become more important, at least for the Υ(3S), increasing the fraction of transversely polarized mesons.