Quarkonium production at the LHC: a phenomenological analysis of surprisingly simple data patterns

TL;DR

LHC quarkonium production data exhibit a universal momentum scaling pattern across different states, suggesting a common underlying production mechanism independent of quantum numbers and mass.

Contribution

This study reveals a universal production pattern in quarkonium data, challenging existing models by showing no dependence on quantum numbers or mass.

Findings

Compatibility of various quarkonium cross sections with a universal momentum scaling

Absence of strong polarization signals in produced S-wave mesons

Data-driven analysis supports a universal production scenario

Abstract

The LHC quarkonium production measurements reveal a startling observation: the J/$\psi$, $\psi$(2S), $\chi_{c1,2}$ and $\Upsilon$(nS) $p_{\rm T}$-differential cross sections are compatible with one universal momentum scaling pattern. Considering also the absence of strong polarizations of directly and indirectly produced S-wave mesons, we are led to the conclusion that there is currently no evidence of a dependence of the partonic production mechanisms on the quantum numbers and mass of the final state. The experimental observations supporting this universal production scenario are remarkably significant, as shown by a new analysis approach, unbiased by specific theoretical calculations of partonic cross sections, which are only considered a posteriori, in comparisons with the data-driven results.

Equations

$$\begin{array}{rl}\sigma (A+B\to H+X)=\sum _{S,L,C}& \mathcal{S}(A+B\to Q\overline{Q}{[}^{2S+1}{L}_J^C]+X)\\ & \times \mathcal{L}(Q\overline{Q}{[}^{2S+1}{L}_J^C]\to H).\end{array}$$

$$\mathcal{L}({}^3{\mathrm{S}}_1^{[8]}\to {\chi }_{c2,b2})/\mathcal{L}({}^3{\mathrm{S}}_1^{[8]}\to {\chi }_{c1,b1})=5/3.$$

$$\begin{array}{rl}\psi (2S)\to {\chi }_{c2}\gamma /\psi (2S)\to {\chi }_{c1}\gamma & :0.95\pm 0.05,\\ {\rm Y}(2S)\to {\chi }_{b2}(1P)\gamma /{\rm Y}(2S)\to {\chi }_{b1}(1P)\gamma & :1.04\pm 0.08,\\ {\rm Y}(3S)\to {\chi }_{b2}(2P)\gamma /{\rm Y}(3S)\to {\chi }_{b1}(2P)\gamma & :1.04\pm 0.16.\end{array}$$

$${\sigma }_{\mathrm{dir}}(p_{\mathrm{T}}/M)={\sigma }_{\mathrm{dir}}^{\ast }[(1-f_{\mathrm{p}}^{\ast })g_{\mathrm{u}}(p_{\mathrm{T}}/M)+f_{\mathrm{p}}^{\ast }{g}_{\mathrm{p}}(p_{\mathrm{T}}/M)],$$

$$h(p_{\rm T}/M)=\frac{p_{\rm T}}{M}\cdot\bigg{(}1+\frac{1}{\beta-2}\cdot\frac{(p_{\rm T}/M)^{2}}{\gamma}\bigg{)}^{-\beta}\,.$$

$$\begin{array}{rl}{\lambda }_{\vartheta }^{{\chi }_0}& =0,\\ {\lambda }_{\vartheta }^{{\chi }_1}& ={\lambda }_{\vartheta }^{\mathrm{J}/\psi }/(4+{\lambda }_{\vartheta }^{\mathrm{J}/\psi }),\\ {\lambda }_{\vartheta }^{{\chi }_2}& =21{\lambda }_{\vartheta }^{\mathrm{J}/\psi }/(60+13{\lambda }_{\vartheta }^{\mathrm{J}/\psi }).\end{array}$$

$$\frac{⟨p^2⟩/m^2}{P^2/M^2}=1+\frac{(M^2-m^2{)}^2}{4m^2{P}^2}+\frac{(M^2-m^2{)}^2}{4m^2{M}^2}\times (1+⟨{\cos }^2\Theta ⟩),$$

$$\begin{array}{rl}{\sigma }^{{\chi }_{c(1,2)}}(p_{\mathrm{T}}/M)& ={\sigma }_{\mathrm{dir}}^{{\chi }_{c(1,2}}(p_{\mathrm{T}}/M)\\ & +B(\psi (2S)\to {\chi }_{c(1,2)}\gamma ){\sigma }^{\psi (2S)}(p_{\mathrm{T}}/M),\\ {\sigma }^{\mathrm{J}/\psi }(p_{\mathrm{T}}/M)& ={\sigma }_{\mathrm{dir}}^{\mathrm{J}/\psi }(p_{\mathrm{T}}/M)\\ & +B({\chi }_{c1}\to \mathrm{J}/\psi \gamma ){\sigma }^{{\chi }_{c1}}(p_{\mathrm{T}}/M)\\ & +B({\chi }_{c2}\to \mathrm{J}/\psi \gamma ){\sigma }^{{\chi }_{c2}}(p_{\mathrm{T}}/M)\\ & +B(\psi (2S)\to \mathrm{J}/\psi X){\sigma }^{\psi (2S)}(p_{\mathrm{T}}/M),\end{array}$$

$${\lambda }_{\vartheta }^{(1)}\oplus {\lambda }_{\vartheta }^{(2)}=\frac{f^{(1)}{\lambda }_{\vartheta }^{(1)}/(3+{\lambda }_{\vartheta }^{(1)})+f^{(2)}{\lambda }_{\vartheta }^{(2)}/(3+{\lambda }_{\vartheta }^{(2)})}{f^{(1)}/(3+{\lambda }_{\vartheta }^{(1)})+f^{(2)}/(3+{\lambda }_{\vartheta }^{(2)})}.$$

$$\begin{array}{rl}{\sigma }_{\mathrm{u}}& \to \sigma ({}^1{\mathrm{S}}_0^{[8]})+\frac{4-4{\lambda }_{\vartheta }({}^3{\mathrm{S}}_1^{[8]})}{4-{\lambda }_{\vartheta }({}^3{\mathrm{S}}_1^{[8]})}\sigma ({}^3{\mathrm{S}}_1^{[8]}),\\ {\sigma }_{\mathrm{p}}& \to \frac{3{\lambda }_{\vartheta }({}^3{\mathrm{S}}_1^{[8]})}{4-{\lambda }_{\vartheta }({}^3{\mathrm{S}}_1^{[8]})}\sigma ({}^3{\mathrm{S}}_1^{[8]}),\end{array}$$

Full text

††journal: Physics Letters B

1 Introduction, motivation and conceptual remarks

Most of the LHC analyses are devoted to “searches”; hundreds of papers report studies aimed at verifying if the LHC data agree with hypotheses put forward by theory models. Such analyses give the primary role to the theory, which makes predictions based on calculations with different levels of accuracy and/or precision: some calculations are performed under assumptions that often incorporate approximations and only up to a given fixed order in a perturbative series. In the case of heavy quarkonium production, non-relativistic quantum chromodynamics (NRQCD) [1] is the commonly-considered standard theory, to be probed by experimental analyses. This model is, indeed, the most sophisticated, complex, and conceptually profound presently available in this chapter of physics. It is however also true that the level of NRQCD calculations is not yet sufficiently satisfactory, and new improvements are continuously taking place. The absence of a reliable “standard model” is not, in itself, a problem, as long as we are aware of the transient nature of the calculations. Indeed, when experimental measurements are fitted (a suited word) within a given theoretical framework, we might be placing the data in a tight corset that moulds the patterns into desired shapes, preventing us from exploring a richer spectrum of options and, worse, potentially blinding us from the simplest and most natural interpretations.

In a previous publication [2], we showed how easy it is to create puzzles when conducting global fits of quarkonium production data in the framework of NRQCD. All we need is to believe that a certain superposition of the presently available next-to-leading order (NLO) perturbative QCD calculations of short-distance coefficients (SDCs) must be able to describe the measured (unpolarized) differential cross sections down to very low quarkonium transverse momentum ($p_{\rm T}$). The highest (statistical) precision of the lowest-$p_{\rm T}$ data drives the result of the fit and leads to an “inescapable prediction”: quarkonium production must be transversely polarized, already at not-so-high $p_{\rm T}$ values. The puzzling nature of the measured (absence of) quarkonium polarizations, which clearly contradict the predictions, can be trivially understood by simply realising that the existing NLO SDCs are not sufficiently accurate at low $p_{\rm T}$. Indeed, Ref. [2] shows that the existing SDCs are perfectly able to simultaneously describe, with a very good fit quality, the cross section and polarization data, provided we avoid the lowest $p_{\rm T}$ region. The puzzle was not in the data but rather in the belief that the current level of calculations already reached a sufficiently-high accuracy, even at low $p_{\rm T}$. If we assume that the theoretical calculations are accurate, we must conclude that the theory has been falsified by the LHC measurements (the corresponding global fits have astonishingly poor fit $\chi^{2}$ probabilities). Instead, we can reproduce all the existing data if we simply assume that the present calculations are inaccurate, a conclusion that keeps the door open regarding the validity of the theory at the fundamental level.

The study reported in the present Letter goes one step further in our data-driven approach: while respecting the NRQCD conceptual ideas, we reconsider its (hierarchical) boundaries and allow ourselves to explore a broader landscape, guided by the map that Nature placed at the disposal of those who give her the central stage and who humbly try to understand her language. We are fortunate to have access to the treasure map uncovered by the LHC experiments, which provides clear and easy-to-follow indications to those who embrace an unbiased vision. Indeed, the (high-$p_{\rm T}$) quarkonium measurements provided by the LHC experiments reveal strong model-independent indications regarding the mechanisms of prompt quarkonium production, somehow missed until now, presumably because the crystal-clear experimental patterns are obscured when seen through theoretically-driven perspectives. These considerations have physical and methodological consequences. On one hand, such indications must provide inspiration for developments in the theoretical description of quarkonium production. On the other hand, the quality reached by the experimental information suggests a new strategy for theory-data comparisons, where fits to measurements are performed with minimal theoretical ingredients and the results are compared only a posteriori with theory predictions based on fixed-order perturbative calculations.

In Section 2 we present our central observation, brought forth by a comparative analysis of prompt quarkonium cross section and polarization measurements in proton-proton collisions at the LHC, with no recourse to model considerations: the data are compatible with a surprisingly simple scenario, in which production and decay kinematics follow an almost universal trend, with no significant distinctions between states of different quantum numbers. This observation raises the central physics question of this paper: how different are the production mechanisms behind ${{}^{3}{\rm S}_{1}}$ and ${{}^{3}{\rm P}_{J}}$ quarkonia? In Section 3 we discuss how the seemingly universal data patterns compare to the NRQCD theory framework, with its relatively complex structure of hierarchies and constraints, and startling differences in the calculated kinematic behaviours of the participating processes.

We then discuss, in Section 4, a data-driven scenario. While the production of ${{}^{3}{\rm S}_{1}}$ quarkonia can be described by a completely general parametrization, exploiting the detailed experimental information available on cross sections and polarizations, parametrizing the production of ${{}^{3}{\rm P}_{J}}$ states still requires model assumptions, given the lack of corresponding polarization data. Our scenario is a direct formalization of the measured trends, assuming universal production and decay properties for the ${{}^{3}{\rm S}_{1}}$ and ${{}^{3}{\rm P}_{J}}$ quarkonia. This represents an important change with respect to the hierarchy of elementary mechanisms foreseen by NRQCD, where $\chi$ production has, a priori, a very different process composition with respect to the $\psi$ and $\Upsilon$ S-wave states. We test this hypothesis in a global fit of charmonium data, properly accounting for the relevant feed-down decays. We report the general kinematic relations used for the modelling of the momentum distributions and polarizations of the indirectly produced states. We also describe the details of our original data-fitting approach, where the individual physical contributions to quarkonium production are exclusively discriminated by their characteristic polarizations, while the $p_{\rm T}$ distributions and relative normalizations are parametrized by a flexible, empirical function. Perturbative calculations of the production kinematics are not used as ingredients anywhere in our analysis. This study shows that it is possible to entrust the measurements, given their current level of precision, with the full responsibility of determining the physical outcome of the fit, obtaining reliable and significant results. The theory calculations are only used a posteriori, in comparisons, for each individual process, with the distributions determined by the fit.

2 The surprising simplicity of the measured patterns

Figure 1 shows the first of the two interesting observations we want to discuss: a seemingly universal pattern in the shapes of the $p_{\rm T}$ distributions of all prompt quarkonia. Indeed, when presented as $p_{\rm T}/M$ distributions, where $M$ is the mass of the quarkonium state, the prompt ${\rm J}/\psi$, $\chi_{c1}$, $\chi_{c2}$, $\psi{\rm(2S)}$ and $\Upsilon{\rm(nS)}$ production cross sections are all compatible with a single kinematic dependence, at least for mid-rapidity ($|y|\lesssim 1$) and for not-too-small $p_{\rm T}/M$. This observation confirms, with higher-$p_{\rm T}$ data, a trend first noticed in Ref. [2]. The consistency of the $p_{\rm T}/M$ distributions is quantified in Section 4, while Ref. [8] offers a detailed evaluation of how precisely the current cross section measurements constrain the $\chi_{c}$ $p_{\rm T}/M$ trend to be similar to those of the ${\rm J}/\psi$ and $\psi{\rm(2S)}$.

Figure 2 presents a clearer view of the $\chi_{c}$ measurements, also adding $\chi_{b}$ data: the $\chi_{c2}/\chi_{c1}$ and $\chi_{b2}/\chi_{b1}$ yield ratios, as well as the ratio of ${\rm J}/\psi$ from $\chi_{c}$ decays to prompt ${\rm J}/\psi$ yield, show a flat dependence vs. $p_{\rm T}/M$. Since the prompt $\psi{\rm(2S)}$ mesons are fully directly produced while the ${\rm J}/\psi$ and $\Upsilon{\rm(nS)}$ states are significantly affected by feed-down contributions from $\chi_{c}$ and $\chi_{b}$ decays ($\simeq 25\%$ [5] and $\simeq 40\%$ [11], respectively), the perfect compatibility of their $p_{\rm T}/M$ shapes is another observation supporting the similarity of P- and S-wave quarkonium production, even in $p_{\rm T}/M$ ranges uncovered by the existing $\chi$ data.

The simplest explanation of such universality of the production kinematics is that one and the same parton-level process (or mixture of processes) describes the production of all states, irrespectively of their masses and quantum numbers. In fact, the kinematic dependence of the partonic cross section of a given process is invariant by simultaneous rescaling of all energy-related variables. This translates to an invariance by $p_{\rm T}/M$ rescaling, if the longitudinal momentum components can be neglected: such scaling is not expected to be valid in the low $p_{\rm T}/M$ and high-rapidity region of the LHCb data (see Refs. [12, 13, 14] and references therein). On the other hand, as shown in Ref. [2], current fixed-order perturbative calculations are seemingly unsuitable to describe cross section measurements for $p_{\rm T}/M\lesssim 3$, justifying that relatively high-$p_{\rm T}$ data remain preferable in data-theory comparisons.

The second experimental result guiding our analysis is the absence of significant polarizations in the measured quarkonium decay distributions, often considered a puzzling observation given that most theory predictions point to significant polarizations, increasing with $p_{\rm T}$ [15]. The most precise measurements of the $\lambda_{\vartheta}$ polar anisotropy parameter of the ${\rm J}/\psi$, $\psi{\rm(2S)}$ and $\Upsilon{\rm(nS)}$ dilepton decay distributions, shown in Fig. 3 (for the helicity frame) as a function of $p_{\rm T}/M$ and averaged over rapidity, are compatible with a small (even negligible) and constant polarization. While no direct $\chi$ polarization measurements exist, we remind that the three states shown in Fig. 3 have very different $\chi$ feed-down components, suggesting that weak polarizations are also to be expected for the $\chi_{c}$ and $\chi_{b}$ states.

We conclude that, today, there is no experimental evidence of differences in production and decay kinematics between quarkonium states of different masses and angular momentum properties, at least in the mid-rapidity region. Such a “universal production” scenario is a surprising experimental observation: in principle, conservation rules should make partonic production cross sections different for states of different quantum numbers. In the next section we discuss how the current models of quarkonium production relate to this observation.

3 The not-so-simple theory patterns

Current studies of quarkonium production phenomenology are based on NRQCD [1], a non-relativistic effective field theory whose pillar is the hypothesis of factorization of the long- and short-distance parts of the production process. Under this hypothesis, the inclusive prompt production cross section of a quarkonium state $H$, after a collision of initial systems $A$ and $B$, can be written as a linear combination of SDCs ($\mathcal{S}$) to produce heavy quark-antiquark pairs ($Q\overline{Q}$) of different colours (singlet or octet) and spin-angular momentum configurations:

[TABLE]

The coefficients $\mathcal{L}$ are the so-called long-distance matrix elements (LDMEs), representing the probabilities that the different $Q\overline{Q}$ states, with spin, angular-momentum and colour configuration described by the indices $S$, $L$, $J$ and $C$, evolve into the observable quarkonium. The LDMEs are assumed to be universal constants, for a given quarkonium, independent of the production mechanism and kinematics (collision system and energy, $p_{\rm T}$ and rapidity), while the SDCs are functions of the collision type and energy, as well as of the $Q\overline{Q}$ laboratory momentum, and are calculated perturbatively. The LDME coefficients can be determined by comparison with the quarkonium cross sections measured as a function of $p_{\rm T}$. It is important to remember, however, that the LDME and the SDC of a given term in the expansion are, individually, unobservable, as their separation depends on the NRQCD factorization scale. Also the intermediate $Q\overline{Q}$ pre-resonance state should not be considered a physical state. The fact that the LDMEs can be determined using experimental data does not mean that they are physical observables: their extraction presumes the knowledge of (equally unphysical) SDC functions. This also means that the LDME values, as determined from data, depend on the perturbative order of the SDC calculations and, effectively, lose their universal character, given that the evolution of perturbative calculations through successive orders usually follows different patterns for different collision systems and/or $p_{\rm T}$ and rapidity regions.

Assuming the limit of heavy component quarks strongly reduces and constrains, for each quarkonium, the wide and indiscriminate range of processes entering Eq. 1. In fact, for small values of the relative velocity, $v$, of the two heavy quarks (non-relativistic limit), a strong hierarchy arises in the LDME magnitudes and only a small number of dominating terms, depending on the final state, are kept in the linear combination. Moreover, some of these terms are ruled out by the perturbative calculations of the corresponding SDCs, because they lead to negligible yields in comparison to the observed prompt cross sections, a well known case being the ${{}^{3}{\rm S}_{1}^{[1]}}$  singlet contribution to ${\rm J}/\psi$ and $\psi{\rm(2S)}$ production [18]. As a result, the prompt production cross sections of the ${{}^{3}{\rm S}_{1}}$ quarkonia, ${\rm J}/\psi$, $\psi{\rm(2S)}$ and $\Upsilon{\rm(nS)}$, are expected to be dominated by the ${{}^{1}{\rm S}_{0}^{[8]}}$, ${{}^{3}{\rm S}_{1}^{[8]}}$ and ${{}^{3}{\rm P}_{J}^{[8]}}$ octet contributions, with comparable magnitudes (the ${{}^{3}{\rm S}_{1}^{[8]}}$ being suppressed at the short-distance level), the $\chi_{c}$ and $\chi_{b}$ production should mainly proceed through the ${{}^{3}{\rm S}_{1}^{[8]}}$ and ${{}^{3}{\rm P}_{1,2}^{[1]}}$ channels, and the $\eta_{c}$ production should reflect the ${{}^{1}{\rm S}_{0}^{[1]}}$ singlet and the ${{}^{3}{\rm S}_{1}^{[8]}}$ octet channels (the ${{}^{1}{\rm S}_{0}^{[8]}}$ and ${{}^{1}{\rm P}_{1}^{[8]}}$ terms being suppressed by the small SDC times LDME products [19]).

Besides $v$-scaling, heavy-quark spin-symmetry (HQSS) also plays a role in constraining the relative magnitudes of the NRQCD LDMEs. For example, the ratio between the octet LDMEs for $\chi_{c2,b2}$ and $\chi_{c1,b1}$ production is expected to be equal to the corresponding ratio of spin states,

[TABLE]

In the present study, we are mainly interested in the characteristic kinematic dependencies of the different contributing processes and on how they compare to the simplicity of the experimental scenario. Figure 4 shows the individual contributions of the dominating $Q\overline{Q}$ configurations, calculated at NLO [20, 21, 22]. The top panels show the $p_{\rm T}/M$ dependence of each SDC; the P-wave SDCs are shown (here and throughout the paper) multiplied by $m_{c}^{2}$. The bottom panels show the polarization parameter $\lambda_{\vartheta}=(1-3\,\beta_{0})/(1+\beta_{0})$ for S-wave dilepton decays, where $\beta_{0}$ is the longitudinal cross section fraction (angular momentum projection $J_{z}=0$) in the helicity frame, for directly-produced ${\rm J}/\psi$ or $\psi{\rm(2S)}$ mesons (left) and for mesons resulting from $\chi_{c1}$ (middle) or $\chi_{c2}$ (right) feed-down decays.

The P-wave contributions have rather peculiar kinematic behaviours, with cross sections becoming negative above certain threshold $p_{\rm T}/M$ values and unphysical polarization parameters ($|\lambda_{\vartheta}|>1$). In fact, the calculations of the singlet and octet P-wave SDCs present singularities that match (and can be cancelled by) the one appearing in the ${{}^{3}{\rm S}_{1}^{[8]}}$ LDME [1, 23]. Therefore, only the ${{}^{3}{\rm P}_{J}^{[1,8]}}+{{}^{3}{\rm S}_{1}^{[8]}}$ sum is required to be well-defined and physical, while the two “fragments” ${{}^{3}{\rm P}_{J}^{[1,8]}}$ and ${{}^{3}{\rm S}_{1}^{[8]}}$ are, effectively, an internal mathematical detail of the calculation, their relative contribution depending, for example, on the choice of the unobservable NRQCD factorization scale. In other words, the individual terms of the calculation are, on their own, unobservable. In particular, the ${{}^{3}{\rm P}_{J}^{[1,8]}}$ components at NLO contain large and negative transverse and/or longitudinal contributions, leading to the kinematic peculiarities mentioned above. While such characteristics can be acknowledged as technical details of a theoretical calculation involving an expansion over purely mathematical objects, they are, nevertheless, in striking contradiction with the universality of the $p_{\rm T}/M$ dependencies of the cross sections and polarizations measured for the different quarkonium states. The non-trivial $p_{\rm T}/M$ dependencies and sign-changing normalizations of the P-wave contributions require seemingly miraculous cancellations to yield results that, besides being physical, must reproduce the comparatively trivial experimental trends.

The different final states come from different pre-resonance mixtures, of characteristic kinematic trends: the ${{}^{3}{\rm S}_{1}^{[8]}}$ octet term, containing the contributions of processes where the $Q\overline{Q}$ comes from a single gluon, is (at high-enough $p_{\rm T}$) fully transversely polarized, while the ${{}^{1}{\rm S}_{0}^{[8]}}$ octet is unpolarized (by definition) and has a much steeper $p_{\rm T}/M$ shape. At high $p_{\rm T}/M$, all three P-wave components seem to asymptotically tend to the cross section and polarization shapes of their complementary-term ${{}^{3}{\rm S}_{1}^{[8]}}$. At lower $p_{\rm T}/M$, instead, in particular in the region covered by current $\chi$ measurements, the trends strongly deviate and, for example, the ${{}^{3}{\rm P}_{1}^{[1]}}$ and ${{}^{3}{\rm P}_{2}^{[1]}}$ terms sharply differentiate the $J=1$ and 2 states, possibly leading to strong polarizations. Comparing the SDCs of Fig. 4 with the experimental trends of Fig. 2, we see that, in particular, the different turnover points of the ${}^{3}P_{1}^{[1]}$, ${}^{3}P_{2}^{[1]}$ and ${{}^{3}{\rm P}_{J}^{[8]}}$ SDCs may prevent a smooth description of the flat $\chi_{2}/\chi_{1}$ and $\chi_{c}/{\rm J}/\psi$ ratios. It is true that the relative importance of the P-wave contributions with respect to the other participating processes is a priori unknown and has to be tuned using the data. In NRQCD, however, $v$-scaling rules predict an approximate equality in magnitude between the ${{}^{1}{\rm S}_{0}^{[8]}}$, ${{}^{3}{\rm S}_{1}^{[8]}}$ and ${{}^{3}{\rm P}_{J}^{[8]}}$ LDMEs for the production of ${{}^{3}{\rm S}_{1}}$ quarkonia, and between the ${{}^{3}{\rm P}_{1,2}^{[1]}}$ and ${{}^{3}{\rm S}_{1}^{[8]}}$ ones for the production of $\chi$ mesons. Moreover, the HQSS relation in Eq. 2 would lead to a $\chi_{2}/\chi_{1}$ ratio close to $5/3$ if the ${{}^{3}{\rm S}_{1}^{[8]}}$ octet term dominated $\chi$ production. The measured ratio (Fig. 2) violates the octet-only expectation by a factor of 2, indicating that the singlet components must be very important. On the other hand, large singlet terms would lead to a significant difference between the $J=1$ and 2 $p_{\rm T}/M$ dependencies, probably excluded by the measured ratio.

This picture gives the impression that the variety of pre-resonance contributions and their associated production and decay kinematics, as implied by $v$-scaling hierarchies, HQSS and current perturbative calculations, is redundant with respect to the observed universal $p_{\rm T}/M$ scaling and lack of polarization. The simple data patterns cannot be reproduced without invoking conspiring cancellations. In other words, the current state-of-the-art theory does not seem to naturally accommodate the experimental scenario, whose simplicity would be interpreted, in this framework, as a coincidence. On the other hand, the responsibility of achieving the necessary cancellations is left to the LDMEs (Eq. 1), which are not calculated and must be obtained from global fits to measurements. This procedure may lead to unstable results, given that the fits are driven by the necessity of cancelling kinematic SDC dependencies not seen in the data. Seemingly successful cancellations may be the result of a temporary coincidence, allowed by the experimental uncertainties on the observables most sensitive to the critical P-wave contributions: the $\chi$ cross sections are still measured with a rather low precision with respect to the ${\rm J}/\psi$ and $\psi{\rm(2S)}$ cross sections, and no $\chi$ polarization measurements exist so far. The resulting LDMEs may, therefore, be artificial and unphysical, a situation that cannot be immediately identified as a problem, given the absence of expected/calculated counterparts.

Furthermore, it is difficult to estimate the precision of the current SDC calculations, at least in the case of the P-wave terms. Such components show dramatic modifications from leading-order (LO) to NLO calculations, changing, in particular, from positive SDCs and reasonable polarizations at LO to mostly negative SDCs and unphysical polarizations at NLO [20, 21, 22, 24, 25, 26]. Moreover, leading-power fragmentation corrections [27], accounting for a subset of next-to-next-to-leading-order (NNLO) processes, still bring startling shape and normalization modifications to the ${{}^{3}{\rm P}_{J}^{[8]}}$ SDC (no corresponding calculations exist for the ${{}^{3}{\rm P}_{1,2}^{[1]}}$ terms), indicating a slow convergence of the perturbative expansion.

It is important to emphasise the potential fragility of the fit procedures currently adopted in NRQCD analyses of quarkonium production. These analyses rely on free LDMEs, representing unphysical quantities not comparable to direct observations or calculations, to obtain precise cancellations between the individually-unobservable SDCs of a mathematically-expanded physical cross section. Such cancellations must be accomplished over the entire kinematic phase space, while the LDMEs are independent of kinematics. The reliability of such procedures would require, in general, a highly over-constrained experimental scenario, not yet reached at present, especially in the $\chi_{c}$ and $\chi_{b}$ measurements and, to some extent, in all polarization measurements, which still lack precision and completeness. However, it is especially the precision of the SDC calculations that should be placed under control. Currently, no meaningful theoretical uncertainties are associated to the kinematic dependencies of the calculated SDCs, some of which show drastic shape modifications across perturbative orders. This is particularly worrisome given that significant deviations from ideal calculations may affect the results of the necessary term-to-term cancellations and even invalidate, in practice if not conceptually, the pivotal postulate that the coefficients of the factorization expansion are universal constants.

Motivated by these considerations, we adopt a completely data-driven analysis method, to ensure that the fit outcome is stable and susceptible to definite physical interpretations, even if at the cost of larger uncertainties in the results and predictions, truthfully reflecting the current experimental inputs. In this new fitting framework, only observable parameters and cross section contributions are considered, the perturbative kinematic calculations being replaced by empirical functions, to be determined by the data.

4 The universal almost unpolarized production scenario

Inspired by the scenario of maximum simplicity depicted by the LHC data (Section 2), in this section we discuss and test the correspondingly simple physical hypothesis that one common unpolarized mechanism dominates the production of both ${{}^{3}{\rm S}_{1}}$ and ${{}^{3}{\rm P}_{J}}$ quarkonium states, while allowing the existence of a “corrective”, transversely polarized term. This “universal almost unpolarized” (UAU) production mechanism does not obey the heavy-quark-limit hierarchies and constraints of NRQCD, in particular for what concerns the $\chi$ states. In the framework of NRQCD, universal production naturally translates into the existence of a common colour-octet channel leading, with suitable changes of quantum numbers, to any final quarkonium state. With $\chi_{1}$ and $\chi_{2}$ produced from the same octet pre-resonance state, HQSS would predict a $\chi_{2}/\chi_{1}$ ratio of $5/3$, very different from the measured values, which cluster around 0.8 (Fig. 2). On the other hand, it can be argued that this constraint on the ratio of octet components is only an approximation, since it neglects effects related to the $\chi_{2}-\chi_{1}$ mass difference and spin-orbital interactions, possibly leading to unequal wave functions of $\chi_{1}$ and $\chi_{2}$ [28, 29]. In fact, Eq. 2 has observable counterparts, for example, in the following ratios of branching fractions [30]:

[TABLE]

All these ratios are significantly different from $5/3$ and, interestingly, much closer to the measured $\chi_{2}/\chi_{1}$ ratios shown in Fig. 2. Incidentally, the ratio $B(\psi{\rm(2S)}\to\chi_{c2}\gamma)\,/\,B(\psi{\rm(2S)}\to\chi_{c0}\gamma)=1.10\pm 0.05$ [30] violates even more dramatically the corresponding spin-counting expectation of $5$. Furthermore, another HQSS relation, $\mathcal{L}({{}^{3}{\rm S}_{1}^{[8]}}\to\eta_{c})=\mathcal{L}({{}^{1}{\rm S}_{0}^{[8]}}\to{\rm J}/\psi)$, translates into the prediction of a large ${{}^{3}{\rm S}_{1}^{[8]}}$ contribution to $\eta_{c}$ production in pp collisions. But the pure ${{}^{1}{\rm S}_{0}^{[1]}}$ singlet cross section, a parameter-free NRQCD prediction, already adequately describes [19], alone, the measured $\eta_{c}$ cross section [31]. The addition of the ${{}^{3}{\rm S}_{1}^{[8]}}$ component significantly overshoots the data, a “puzzle” created by the HQSS constraint. We also note that spin-counting rules have been seen to fail even in the description of the production of simpler, heavy-light quark systems, predicting for example wrong yield ratios between charged and neutral $D$ mesons [32, 33].

The violation of HQSS-limit constraints is not the only source of conflict between the UAU hypothesis and NRQCD. The natural candidate for a single mechanism capable of producing any final state leading to unpolarized production is the ${{}^{1}{\rm S}_{0}^{[8]}}$ octet channel. In NRQCD, the role of this process in P-wave quarkonium production is disfavoured by the $v$-scaling rules, which favour the dominance of the combinations ${{}^{3}{\rm S}_{1}^{[8]}}+{{}^{3}{\rm P}_{1,2}^{[1]}}$, characterized by non-zero polarizations and much flatter high-$p_{\rm T}/M$ distributions (Fig. 4). In fact, a priori, the current SDC calculations predict a significant difference between the $\chi$ and $\psi$ production kinematics. Incidentally, the dominance of ${{}^{3}{\rm S}_{1}^{[8]}}$ as common production channel would fit much better the $v$-scaling hierarchies, such process being already a leading contribution for both $\chi$ and $\psi$ production, but this hypothesis leads to strong transverse polarizations of the ${{}^{3}{\rm S}_{1}}$ quarkonia and is excluded by the data. In summary, when seen as a “${{}^{1}{\rm S}_{0}^{[8]}}$ dominance model”, the data-inspired UAU scenario represents a stretching of the heavy-quark-limit hierarchies of NRQCD towards an unforeseen, stronger hierarchy.

While the UAU hypothesis may seem a rather specific physical assumption, its implementation in our analysis is general, as much as allowed by the existing data. Inspired by the pattern of Fig. 3, we consider that the production cross section can be decomposed into an “unpolarized” term and a “transversely polarized” term. These denominations are the only qualifications we use in our description of the participating processes. The two terms can be ideally associated, respectively, to the NRQCD ${{}^{1}{\rm S}_{0}^{[8]}}$ octet channel and to the net contribution of all other mechanisms, but any other scenario would fit in this framework, including one involving a superposition of (unobserved) longitudinally and transversely polarized sub-processes, provided that the observable net result is not longitudinal. This procedure could be generalized to scenarios including observed longitudinal polarizations by replacing the unpolarized contribution with a longitudinally polarized one and letting the data determine, in a model-independent way, the longitudinal and transverse cross-sections contributions.

In a production scenario characterized by non-longitudinal observable polarization, the cross section for the direct yield of a given ${{}^{3}{\rm S}_{1}}$ state can be parametrized as

[TABLE]

where $\sigma_{\rm dir}^{*}$ and $f_{\rm p}^{*}$ are, respectively, the total direct-production cross section and the fractional contribution of the polarized process, both calculated at a fixed reference value, which we choose to be $(p_{\rm T}/M)^{*}=2$. The $p_{\rm T}/M$ dependencies of the unpolarized and polarized yields are described by the shape functions $g_{\rm u}(p_{\rm T}/M)$ and $g_{\rm p}(p_{\rm T}/M)$, respectively, both normalized to unity at $(p_{\rm T}/M)^{*}$, $g(p_{\rm T}/M)=h(p_{\rm T}/M)/h((p_{\rm T}/M)^{*})$, with

[TABLE]

The parameter $\gamma$ defines the function in the low-$p_{\rm T}$ turn-on region and is only mildly sensitive to the data we are considering here; hence, in the fit we consider $\gamma$ as a common free parameter. The $\beta$ power-law exponent, instead, characterizes the high-$p_{\rm T}$ shape ($h\propto(p_{\rm T}/M)^{1-2\beta}$ for $p_{\rm T}/M\gg\sqrt{\gamma(\beta-2)}$); therefore, we distinguish the unpolarized and polarized cross sections with two different powers, $\beta_{\rm u}$ and $\beta_{\rm p}$, respectively.

The discrimination between the two physical contributions relies on the experimental polarization constraint, the unpolarized and polarized processes for ${{}^{3}{\rm S}_{1}}$ states being characterized by $\lambda_{\vartheta}=0$ and $\lambda_{\vartheta}=1$, respectively. The polarized yield fraction can be expressed, as a function of $p_{\rm T}/M$, as $f_{\rm p}=3\lambda_{\vartheta}/(4-\lambda_{\vartheta})$. This is a crucial difference of approach with respect to fits using SDC shapes fixed by calculations. In those analyses, the result of the fit is mainly or exclusively determined by the comparison of the calculated $p_{\rm T}$ dependencies with the cross section measurements, while the less precise polarization data are not included in the fits or have a relatively weak constraining effect. We adopt the opposite point of view, leaving to the polarization measurements, as functions of $p_{\rm T}/M$, the exclusive role of constraining both the relative normalizations and the relative differences between the momentum dependencies of the different process contributions. Given that the analysis becomes, therefore, fully driven by the experimental measurements, the precision of the results is no longer limited by unquantifiable theoretical uncertainties and will improve as the experimental panorama evolves.

Our previous analysis [2], using calculated SDCs, addressed the production of the $\psi{\rm(2S)}$ and $\Upsilon{\rm(3S)}$ quarkonia, the ones less affected by feed-down decays from P-wave quarkonia. We now concentrate on the charmonium system, complementing the $\psi{\rm(2S)}$ measurements with ${\rm J}/\psi$, $\chi_{c1}$ and $\chi_{c2}$ data. The bottomonium system, with its more articulated feed-down structure, not yet sufficiently constrained by data (no mid-rapidity $\chi_{b}$ measurements exist), is left for future consideration. We also profit from much improved $\psi{\rm(2S)}$ and ${\rm J}/\psi$ cross section measurements recently made available by ATLAS and CMS, extending to significantly higher $p_{\rm T}$. The similarity of the $p_{\rm T}/M$ dependencies shown by these precise measurements, discussed in Section 2, guides us to assuming that the direct production mechanism is the same for the ${\rm J}/\psi$ and $\psi{\rm(2S)}$ mesons. This translates, in our fit, in the use of a single “polarized fraction”, represented by the parameter $f_{\rm p}^{*}$, common to the two states.

No significant experimental indications exist, so far, regarding $\chi_{c}$ polarization. This observable may be indirectly constrained by the difference between the measured ${\rm J}/\psi$ and $\psi{\rm(2S)}$ polarizations, the former including about 25% of $\chi_{c}$ contribution while the latter is fully directly produced. Future measurements of the $\chi_{c}$ polarization and/or more precise ones of the ${\rm J}/\psi$ and $\psi{\rm(2S)}$ polarizations will allow us to perform a fully data-driven fit, without conjectures on differences between the $\chi_{c}$ production and the ${\rm J}/\psi$ or $\psi{\rm(2S)}$ production, that is, without constraints on the “universality” of S- and P-wave quarkonium production. At the present moment, however, the fit results are insensitive to different $\chi_{c}$ polarization hypotheses. To ensure that the fit is sufficiently constrained, given the precision of the existing data, we reduce the number of free parameters by adopting plausible relations between the properties of S- and P-wave states. It is only at this point, and only for the production of P-wave quarkonia, that our fit framework acquires model-dependent ingredients. While other hypotheses are considered in Ref. [8], here we compensate the missing $\chi_{c}$ polarization information by adopting the data-driven UAU conjecture: direct $\chi_{c}$ production is described by the same mixture of processes as the direct ${\rm J}/\psi$ and $\psi{\rm(2S)}$ production. In other words, the $p_{\rm T}/M$ shape functions of the polarized and unpolarized contributions, as well as the polarized fraction parameter $f_{\rm p}^{*}$, are assumed to be common to the P- and S-wave states.

As previously discussed, a possible realization of such a mechanism, not distinguishing between final states of different quantum numbers, is the production via colour-octet pre-resonances. This model inspires our definition of the possible $\chi_{c}$ polarizations. In general, while an intrinsically unpolarized pre-resonance will lead to unpolarized S- and P-wave final states, a transversely polarized one will yield different polarizations when it transforms to ${\rm J}/\psi$ ($L=0$, $J=1$) or $\chi_{c0,c1,c2}$ ($L=1$, $J=0,1,2$) mesons, because of different changes in orbital and total angular momenta. The following expressions, only reflecting angular momentum conservation, relate the decay anisotropy parameter $\lambda_{\vartheta}$ of $\chi_{c}$ mesons (in the ${\rm J}/\psi$ plus photon channel) to the corresponding one of a ${\rm J}/\psi$ meson (in the dilepton decay channel), when all transform from the same pre-resonance state via gluon emissions:

[TABLE]

This means, for example, that the “polarized” $\chi_{c}$ contribution, defined by $\lambda_{\vartheta}=+1$ for the ${\rm J}/\psi$ dilepton decay, is characterized by polarization parameters 0, $1/5$ and $21/73$, respectively, in the $\chi_{c0}$, $\chi_{c1}$ and $\chi_{c2}$ decays to ${\rm J}/\psi$ plus photon.

In our analysis we model the charmonium feed-down decay chains ($\psi{\rm(2S)}\to\chi_{c1,c2}\gamma$, $\psi{\rm(2S)}\to{\rm J}/\psi+X$ and $\chi_{c1,c2}\to{\rm J}/\psi\gamma$) taking into account momentum and polarization transformations. Calculating the decay kinematics gives the relation connecting the mother’s and daughter’s laboratory momenta,

[TABLE]

where $M$ ($m$) and $P$ ($p$) are, respectively, the mass and laboratory momentum of the mother (daughter) particle, and $\Theta$ is the emission angle of the daughter in the mother’s rest frame. The average symbols denote an integration over $\cos\Theta$, making linear terms in this variable disappear. The term $\langle\cos^{2}\Theta\rangle$ can vary between $1/5$ and $2/5$ for physical polarizations of the mother particle. In the calculation we only made the approximation that $m_{X}\ll\sqrt{M^{2}+m^{2}}$, where $m_{X}$ is the total mass of the accompanying decay particles (either a photon or pions, depending on the decay channel). The deviation from unity in the right term of the equation is never larger than a few percent for all relevant feed-down cases and for $P$ not much smaller than $M$: the ratio of laboratory momentum over mass remains practically unchanged, on average, in the transition from mother to daughter. Since the two particles have, for $p\gg(M-m)$, almost indistinguishable directions in the laboratory, the relation is also valid vectorially and we can assume $p_{\rm T}/m=P_{\rm T}/M$. In these conditions, we can model the kinematics of an indirectly produced particle using the $p_{\rm T}/M$ distribution of its mother particle. Therefore, we account for the momentum conversion along any feed-down chain by simply considering all observables as functions of $p_{\rm T}/M$.

The total prompt cross sections of the states affected by feed-down are obtained as

[TABLE]

while $\sigma=\sigma_{\rm dir}$ in the $\psi{\rm(2S)}$ case.

We now address the relevant polarization transfer rules. Measurements of the $\psi{\rm(2S)}\to{\rm J}/\psi\ \pi\pi$ decay angular distribution [34, 35, 36] reveal that, in very good approximation, the ${\rm J}/\psi$ retains the angular momentum alignment of the mother $\psi{\rm(2S)}$, being the $\pi\pi$ system dominated by its $J=0$ component. It can be assumed, therefore, that in the $\psi{\rm(2S)}$ to ${\rm J}/\psi$ feed-down the polarization is transmitted without modification, and that the shape parameters of the dilepton decay distribution (in particular, $\lambda_{\vartheta}$) do not change. Concerning the ${\rm J}/\psi$ mesons produced by $\chi_{c}$ feed-down decays, it has been shown [37] that, while the angular momentum projection changes from mother to daughter, the $\chi_{c}\to{\rm J}/\psi+\gamma$ decay and the ${\rm J}/\psi$ to leptons decay have surprisingly identical angular distribution shapes, the latter one having the advantage of being unaffected by the non-zero orbital angular momentum components of the photon. Therefore, also in this case (but for different reasons) the angular decay parameter $\lambda_{\vartheta}$ remains the same (even if obviously referring to two different decay channels) in the $\chi_{c}$ to ${\rm J}/\psi$ transition. Finally, the $\chi_{c}$ mesons produced by radiative decays of the $\psi{\rm(2S)}$ (as well as the ${\rm J}/\psi$ mesons from the subsequent $\chi_{c}$ radiative decays) have the same polarization parameters as in Eq. 6, after replacing $\lambda_{\vartheta}^{{\rm J}/\psi}$ with $\lambda_{\vartheta}^{\psi{\rm(2S)}}$.

To calculate the observable polarizations of the ${\rm J}/\psi$ and $\chi_{c}$ states, combining their direct and indirect contributions (Eq. 8), we use the general addition rule from Ref. [38], which, in the example of two processes contributing with fractions $f^{(1)}$ and $f^{(2)}$ of the cross section, reads as

[TABLE]

The fit has four state-independent (universal) free parameters: the shape parameters $\gamma$, $\beta_{\rm u}$ and $\beta_{\rm p}$ of the unpolarized and polarized cross sections and the fractional polarized contribution $f_{\rm p}^{*}$ (at the reference point $(p_{\rm T}/M)^{*}=2$). The other parameters are the normalizations of the direct production cross sections of the four charmonia.

Additional nuisance parameters are used to model the correlations induced by the luminosity uncertainties. The branching ratios used to extract the cross sections from the yields measured in specific decay channels, as well as those weighting the corresponding feed-down terms of the prompt cross sections, are also treated as constrained parameters, with central values and uncertainties taken from Ref. [30].

The cross section measurements depend on the polarization hypothesis assumed by the experiments for the determination of the detection acceptance. For each set of parameter values scanned during the fit and each $p_{\rm T}$-rapidity bin of the cross section measurement, we recalculate the average $\lambda_{\vartheta}$ parameter for the considered state and apply the corresponding polarization-dependent acceptance correction. The fit includes the charmonium measurements shown in Figs. 1, 2 and 3. The ATLAS $\chi_{c}$-to-${\rm J}/\psi$ feed-down fraction and $\chi_{c2}$-to-$\chi_{c1}$ yield ratio are not considered, to avoid correlations with the $\chi_{c}$ cross section data. The cross sections are only considered for $p_{\rm T}/M>2$. The global fit has 82 degrees of freedom and leads to a total $\chi^{2}$ of 29, an exceptionally good result. The agreement between the model and the data can be appreciated in Fig. 5.

The polarized fraction is quite small at low $p_{\rm T}/M$, reflecting the polarization data: $f_{\rm p}^{*}=(4.8\pm 2.1)\%$ (at $(p_{\rm T}/M)^{*}=2$). It increases almost linearly with increasing $p_{\rm T}/M$, becoming $(20\pm 5)\%$ at $p_{\rm T}/M=15$. The best values of the shape parameters are $\gamma=0.80\pm 0.17$, $\beta_{\rm u}=3.47\pm 0.05$ and $\beta_{\rm p}=2.97\pm 0.15$.

Figure 6 shows the universal unpolarized and polarized $p_{\rm T}/M$ distributions obtained from the fit, normalized to unity at $p_{\rm T}/M=2$ to illustrate the allowed spectrum of kinematic shapes, irrespectively of the normalization shift. The unpolarized curve (with a barely visible uncertainty band) is compared with the ${{}^{1}{\rm S}_{0}^{[8]}}$ octet SDC shape calculated at NLO with NNLO fragmentation contributions, represented by the dashed curve: the agreement is astonishing. Similarly, the ${{}^{3}{\rm S}_{1}^{[8]}}$ octet SDC calculation is in perfect agreement with the fitted polarized component, affected by larger uncertainties. This exercise supports the hypothesis that the universal unpolarized and polarized cross sections have a physical affinity with the ${{}^{1}{\rm S}_{0}^{[8]}}$ and ${{}^{3}{\rm S}_{1}^{[8]}}$ terms of the NRQCD factorization expansion. It is worth noting that the two SDCs have significantly different shapes, their ratio changing by a factor 10 in the range $3<p_{\rm T}/M<13$.

We conclude this section by discussing a new interesting hypothesis on long-distance polarization effects in quarkonium production and how it relates to the scenarios we are considering. A crucial ingredient of NRQCD is the “standard” conjecture that the ${{}^{3}{\rm S}_{1}^{[8]}}$ pre-resonance transfers its polarization (fully transverse at not-too-low $p_{\rm T}$) to the observable ${\rm J}/\psi$, $\psi{\rm(2S)}$ or $\Upsilon{\rm(nS)}$ without attenuation, if no intermediate feed-down transitions are involved. The process is modelled as a transition of the $Q\overline{Q}$ pair with the simultaneous emission of two gluons having total angular momentum $J=0$, in analogy with the observable process $\psi{\rm(2S)}\to{\rm J}/\psi\ \pi\pi$, where, as previously mentioned, the ${\rm J}/\psi$ inherits the $\psi{\rm(2S)}$ polarization. It has been recently observed [39] that another kind of process may be at play, consisting of two successive electric dipole transitions: ${{}^{3}{\rm S}_{1}^{[8]}}\to{{}^{3}{\rm P}_{J}^{[8]}}+g$, followed, for example, by ${{}^{3}{\rm P}_{J}^{[8]}}\ \to{\rm J}/\psi\ +g$. In this case, the resulting ${\rm J}/\psi$ and $\Upsilon$ polarization from ${{}^{3}{\rm S}_{1}^{[8]}}$  would be strongly attenuated with respect to the standard conjecture. In fact, from the polarization perspective the two steps are analogous to the radiative transitions $\psi{\rm(2S)}\to\chi_{c}+g$ and $\chi_{c}\to{\rm J}/\psi\ +g$: in the first one the polarization is attenuated, according to the relations in Eq. 6, while the second step preserves the value of $\lambda_{\vartheta}$. According to the calculations of Ref. [39], the ${{}^{3}{\rm S}_{1}^{[8]}}$ polarization may even be reduced to zero, if interference effects are important. This “double-transition” process should coexist with the one of the standard conjecture, in presently unknown proportions, helping NRQCD predictions to converge to an unpolarized scenario of ${\rm J}/\psi$, $\psi{\rm(2S)}$ and $\Upsilon{\rm(nS)}$ production. The new hypothesis does not affect $\chi_{c}$ polarizations, since the transition of ${{}^{3}{\rm S}_{1}^{[8]}}$ to $\chi_{c}$ is already modelled according to the minimal assumption of a single electric dipole transition. A further possible justification of an unpolarized ${{}^{3}{\rm S}_{1}^{[8]}}$ contribution has been provided in the $k_{\rm T}$-factorization approach [40], which takes into account the off-shellness of the initial gluons.

Our general decomposition of ${\rm J}/\psi$ and $\psi{\rm(2S)}$ cross sections into “unpolarized” and (fully transversely) “polarized” contributions can fit both scenarios, of transverse or unpolarized ${{}^{3}{\rm S}_{1}^{[8]}}$, since neither involves possible longitudinal contributions. However, whenever a parallelism to the NRQCD ${{}^{3}{\rm S}_{1}^{[8]}}$ term is made, as in Fig. 6, we choose to refer to the standard conjecture of fully transverse ${{}^{3}{\rm S}_{1}^{[8]}}$ polarization. In the new conjecture, instead, the ${{}^{3}{\rm S}_{1}^{[8]}}$ term would contribute to both the unpolarized ($\sigma_{\rm u}$) and polarized ($\sigma_{\rm p}$) cross sections, depending on the effective polarization ($0<\lambda_{\vartheta}({{}^{3}{\rm S}_{1}^{[8]}})<1$) transmitted by the ${{}^{3}{\rm S}_{1}^{[8]}}$ octet to the ${\rm J}/\psi$ or $\Upsilon{\rm(nS)}$ meson. The comparison in Fig. 6, for example, should be made with a different association of the unpolarized and polarized contributions to the two octet channels,

[TABLE]

where $\sigma({{}^{\rm 2S+1}{\rm S}_{J}^{[8]}})=\mathcal{S}({{}^{\rm 2S+1}{\rm S}_{J}^{[8]}})\times\mathcal{L}({{}^{\rm 2S+1}{\rm S}_{J}^{[8]}})$. The second equation shows that the polarized term should now be compared to a simply rescaled ${{}^{3}{\rm S}_{1}^{[8]}}$ function, leaving the pure shape comparison shown in the right panel unaltered. The comparison term for the unpolarized function, however, would become a combination of the ${{}^{1}{\rm S}_{0}^{[8]}}$ shape with a ${{}^{3}{\rm S}_{1}^{[8]}}$ contribution depending on the assumed $\lambda_{\vartheta}({{}^{3}{\rm S}_{1}^{[8]}})$ and ${{}^{3}{\rm S}_{1}^{[8]}}$/${{}^{1}{\rm S}_{0}^{[8]}}$ yield ratios. Independently of these assumptions, such a combination would be a less steep function of $p_{\rm T}/M$. It is interesting to notice that, given the very good agreement between the unpolarized shape and the ${{}^{1}{\rm S}_{0}^{[8]}}$ term alone (shown in Fig. 6), the present measurements already disfavour a non-negligible contribution from an additional unpolarized contribution having the shape of the ${{}^{3}{\rm S}_{1}^{[8]}}$ SDC. In other words, the data do not indicate that the ${{}^{3}{\rm S}_{1}^{[8]}}$ polarization, always considered to be $\lambda_{\vartheta}({{}^{3}{\rm S}_{1}^{[8]}})\approx+1$, should instead be closer to 0.

5 Discussion and conclusions

In this Letter we described a fully data-driven approach to the phenomenology of quarkonium production. We replaced perturbative calculations of the process kinematics, whose uncertainties are not reliably evaluated and in some cases are known to be large, with flexible empirical functions to be determined by the fit to the data. This choice assigns to the polarizations, mainly determined by basic physics considerations and, therefore, much more cleanly and stably predicted, the responsibility of discriminating between the different elementary processes of quarkonium production.

The proposed method implies an inversion of perspective in comparison to other analyses, where the imposition of the calculated individual cross-section contributions (the SDCs) results in a process discrimination based exclusively on the predicted $p_{\rm T}$ dependencies. In fact, polarization measurements are rarely included in global fits and, when included, they do not compete in constraining power with the cross section data, given the inferior statistical precision obtainable in multidimensional decay-angle analyses with respect to the much simpler yield counting. Our new fit procedure ignores the unobservable details of the theory calculations, in particular the existence of unphysical partial contributions (P-wave terms) and the factorization into the short- and long-distance matrix elements. The fit results exclusively reflect the data trends and uncertainties, remaining unbiased by theoretical ingredients. Comparisons with SDC calculations and LDME evaluations are made a posteriori. This is a safer procedure for data-theory comparisons, separating the data-driven constraints, with their rigorous statistical interpretation, from the model calculations, usually affected by unreliable uncertainties.

The experimental picture of quarkonium production in pp collisions is still incomplete. Once reasonably precise $\chi_{c}$ polarization measurements will become available, our procedure will no longer need any model input. At this moment, given the absence of $\chi_{c}$ polarization data, we must resort to some model assumptions regarding $\chi_{c}$ production. In this Letter, we assumed the UAU scenario, which simply reflects a full universality of charmonium production across all final states.

The disentangling of theory calculations from the fit to data allows us to face with awareness and devoted attention problems specific to one of the two heterogeneous worlds. For example, tensions originating from mutual experimental inconsistencies (such as those seen in Ref. [41]) become evident and can be addressed before and independently of the comparison to theory, without affecting the judgment of how well one or the other model describes the data. Most importantly, while the fit results only represent the data with their well defined and statistically-interpretable uncertainties, they can be confronted simultaneously to different models or multiple variants of the theory, the best way to appreciate the elusive concept of theory uncertainty. Instead, when model calculations are embedded in the fit, the quantification of theory uncertainties is either neglected or, in the best case, chosen rather arbitrarily (as in Ref. [2]), possibly leading to subjective fit results, not amenable to rigorous interpretations. Therefore, it is difficult to decide how well a model describes the data, since the goodness of the global fit (expressed by the $\chi^{2}$) has no statistical meaning when part of the uncertainties are lacking or are arbitrary.

In such cases, interesting indications can only come from macroscopic divergences, seemingly beyond any reasonable uncertainty, as in the case of the famous polarization puzzle. Even so, this puzzle persisted for so long exactly because of the entangled treatment of experimental and theoretical ingredients, which smeared the crucial (albeit a posteriori simple and predictable) evidence that the current perturbative SDC calculations, at NLO, are definitely steeper, at low $p_{\rm T}$, than the corresponding experimental distributions. This feature, rigorously studied in Ref. [2] by performing fits as a function of a selection cut progressively removing low-$p_{\rm T}$ data, is immediately visible when the results of our data-driven fits are compared to calculations, as done in Fig. 6. As also illustrated by the LDME determination presented in Ref. [8], this method gives a better insight into the limitations of the current calculations, allowing us to quantify the effects of the implied approximations.

It is especially important to realize that the present experimental and theoretical knowledge on quarkonium production is significantly more satisfactory than in the past. In the absence of spectacular inconsistencies to be deciphered and with the prospect of forthcoming LHC measurements filling the few remaining gaps and further improving the constraining power of the experimental picture, quarkonium studies finally deserve the prestigious ranking of precision QCD studies. Heading towards a detailed understanding of these phenomena, also through an accurate testing of NRQCD, a rigorous theory-data comparison method is essential.

Another central element in our analysis approach is the notion that particles characterized by the same partonic production mechanisms should have essentially the same $p/M$ distributions and that (for relatively small mass differences, as in the quarkonium spectrum) $p/M$ is also preserved in feed-down decays. This idea allows us to look for and identify physics patterns in a completely data-driven way. For this reason, we have limited our present analysis to mid-rapidity LHC data, for which the scaling variable $p/M$ is sufficiently well approximated by $p_{\rm T}/M$. The transformation of existing forward-rapidity $p_{\rm T}$ distributions into $p$ distributions requires a detailed analysis of rapidity dependencies and will be the subject of a future study.

The analysis reported in this Letter shows a surprising fact: even the latest CMS and ATLAS cross-section measurements, extending up to $p_{\rm T}$ of order 100 GeV, do not provide any evidence of differences in short-distance production mechanisms throughout the quarkonium spectrum. The $p_{\rm T}/M$ cross-section scaling is perfectly replicated among the seven measured states, from the ${\rm J}/\psi$ to the $\Upsilon{\rm(3S)}$, and all five S-wave states, for which the polarizations have been measured, show similar unpolarized behaviours. While it is true that the $\chi$ states still require more precise and complete direct measurements (in particular, first polarization measurements), significant indications can be derived from the comparison of the kinematic distributions of the several S-wave particles, which receive very different contributions from $\chi$ decays. In fact, the measured $p_{\rm T}/M$-differential ${\rm J}/\psi$ cross section, 25% of which arises from $\chi_{c}$ feed-down decays, is so similar to the $\psi{\rm(2S)}$ one, fully directly produced, that the $\chi_{c}$ $p_{\rm T}/M$ distribution is constrained to also have a very similar shape. The significance of this observation is discussed in Ref. [8], through a test where the $\chi_{c1}$ and $\chi_{c2}$ spectra are replaced with single (average) points, thereby removing the experimental constraints on the $\chi_{c}$ cross section shapes, which are then freely and independently determined by the global charmonium fit.

The idea that all measured quarkonia may have an undistinguishable partonic origin would imply that their final-state identities are acquired only in a subsequent, long-distance phase. This concept, a cornerstone of the NRQCD factorization framework, can be made extreme by assuming that all states, including the $\chi$ as much as the $\psi$ and $\Upsilon$ mesons, are predominantly produced from the same initial pre-resonance state, necessarily a colour-octet one. This is the hypothesis (UAU scenario) that we tested in this Letter, finding an excellent agreement with charmonium data. The existence of such a simple and satisfactory description of the data is quite remarkable and should stimulate theoretical considerations towards a natural explanation of the observed simplicity.

Interestingly, the UAU conjecture, as an exclusive alternative to the NRQCD heavy-quark-limit symmetry rules, can be tested in lighter quarkonium systems. For example, the observation of almost unpolarized $\phi$ production, with no significant dependence on the particle’s laboratory momentum, would support the existence of fundamental principles dictating that quark-antiquark bound-state formation mainly proceeds through unpolarized configurations. The general principles implicit in the UAU hypothesis would also naturally justify the seemingly surprising observation that the $\eta_{c}$ cross section only contains the ${{}^{1}{\rm S}_{0}^{[1]}}$ singlet term. When interpreted as the suppression of the “higher-order” members of an unexpectedly strong angular-momentum hierarchy, this observation becomes equivalent to the dominance of the ${{}^{1}{\rm S}_{0}^{[8]}}$ term in ${\rm J}/\psi$ production, another observation surpassing NRQCD expectations.

Finally, by condensing all the experimental information on quarkonium production (at least at not-too-low $p_{\rm T}$) in a simple analytical formula with few parameters, the UAU scenario provides a useful empirical tool for Monte Carlo simulations, as well as for the modelling of background or reference processes, for example in studies of nuclear modifications of quarkonium production in proton-nucleus and nucleus-nucleus collisions.

Acknowledgements

We are indebted to Hua-Sheng Shao, who kindly gave us the NLO calculations of the short distance coefficients. The work of P.F. is supported by FCT, Portugal, through the grant SFRH/BPD/98595/2013, while the work of I.K. is supported by FWF, Austria, through the grant P 28411-N36.