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

Pietro Faccioli  Mariana Araújo  Valentin Knünz  Ilse Krätschmer  Carlos
Lourenço  João Seixas  Physics Department, Instituto Superior Técnico (IST),
Lisbon, Portugal  Laboratório de Instrumentação e Física Experimental de
Partículas (LIP), Lisbon, Portugal  European Organization for Nuclear Research
(CERN), Geneva, Switzerland  Institute of High Energy Physics (HEPHY), Vienna,
Austria

###### 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.

###### keywords:

Quarkonium , Polarization , NRQCD , QCD , Hadron formation

††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.

Figure 1: Mid-rapidity prompt quarkonium cross sections measured in pp
collisions at $\sqrt{s}=7$ TeV by ATLAS (red markers) [3, 4, 5] and CMS (blue
markers) [6, 7], as a function of $p_{\rm T}/M$. The curves represent a single
empirical function, with shape parameters determined by a simultaneous fit to
all data (of $p_{\rm T}/M>3$) and normalizations specific to each state (left
panel) or adjusted to the ${\rm J}/\psi$ points (right panel) to directly
illustrate the universality of the kinematic dependences.

##  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: Mid-rapidity $\chi_{c2}/\chi_{c1}$ and $\chi_{b2}/\chi_{b1}$ prompt
yield ratios, as well as the fraction of the prompt ${\rm J}/\psi$ yield
coming from $\chi_{c}$ decays, measured in pp collisions by ATLAS [5] and CMS
[9, 10], as a function of $p_{\rm T}/M$. The arrow indicates $5/3$, the value
predicted by heavy-quark spin-symmetry for the $\chi_{c2}/\chi_{c1}$ and
$\chi_{b2}/\chi_{b1}$ ratios in the hypothesis of pure colour-octet
production.

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.

Figure 3: Polar anisotropy parameter $\lambda_{\vartheta}$ measured by CMS in
pp collisions at $\sqrt{s}=7$ TeV, for prompt ${\rm J}/\psi$, $\psi{\rm(2S)}$
and $\Upsilon{\rm(1S)}$ dilepton decays [16, 17]. For improved clarity, values
corresponding to two or three rapidity bins were averaged, assuming
uncorrelated systematic uncertainties, and the very uncertain
$\Upsilon{\rm(2S)}$ and $\Upsilon{\rm(3S)}$ data are not shown.

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:

$\displaystyle\begin{split}\sigma(A+B\to
H+X)=\sum\limits_{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{split}$ (1)

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,

$\mathcal{L}({{}^{3}{\rm
S}_{1}^{[8]}}\to\chi_{c2,b2})\,/\,\mathcal{L}({{}^{3}{\rm
S}_{1}^{[8]}}\to\chi_{c1,b1})=5/3\;.$ (2)

Figure 4: SDCs (top) and corresponding helicity-frame $\lambda_{\vartheta}$
polarization parameters (bottom), for the main octet and singlet components of
directly-produced $\psi$ mesons (left), and for ${\rm J}/\psi$ mesons from
$\chi_{c1}$ (middle) and $\chi_{c2}$ (right) decays, calculated at NLO [20,
21, 22] for pp collisions at $\sqrt{s}=7$ TeV and a charmonium mass of 3 GeV.
The ${{}^{3}{\rm P}_{J}^{[8]}}$ SDC is defined as $\mathcal{S}({{}^{3}{\rm
P}_{0}^{[8]}})+3\mathcal{S}({{}^{3}{\rm
P}_{1}^{[8]}})+5\mathcal{S}({{}^{3}{\rm P}_{2}^{[8]}})$. The ${{}^{3}{\rm
P}_{J}^{[8]}}$ and ${{}^{3}{\rm P}_{1,2}^{[1]}}$ SDC functions are multiplied
by $m_{c}^{2}$, the squared charm-quark mass; they become negative at high
$p_{\rm T}/M$ and are plotted with flipped signs.

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]:

$\displaystyle\begin{split}\psi{\rm(2S)}\to\chi_{c2}\,\gamma\;/\;\psi{\rm(2S)}\to\chi_{c1}\,\gamma\quad~{}~{}~{}&:0.95\pm
0.05\;,\\\
\Upsilon{\rm(2S)}\to\chi_{b2}{\rm(1P)}\,\gamma\;/\;\Upsilon{\rm(2S)}\to\chi_{b1}{\rm(1P)}\,\gamma&:1.04\pm
0.08\;,\\\
\Upsilon{\rm(3S)}\to\chi_{b2}{\rm(2P)}\,\gamma\;/\;\Upsilon{\rm(3S)}\to\chi_{b1}{\rm(2P)}\,\gamma&:1.04\pm
0.16\;.\end{split}$ (3)

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

$\sigma_{\rm dir}(p_{\rm T}/M)=\sigma_{\rm dir}^{*}[(1-f_{\rm p}^{*})\,g_{\rm
u}(p_{\rm T}/M)+f_{\rm p}^{*}\,g_{\rm p}(p_{\rm T}/M)]\;,$ (4)

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

$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}\,.$ (5)

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:

$\displaystyle\begin{split}\lambda_{\vartheta}^{\chi_{0}}&=0\;,\\\
\lambda_{\vartheta}^{\chi_{1}}&=\lambda_{\vartheta}^{{\rm
J}/\psi}/(4+\lambda_{\vartheta}^{{\rm J}/\psi})\;,\\\
\lambda_{\vartheta}^{\chi_{2}}&=21\lambda_{\vartheta}^{{\rm
J}/\psi}/(60+13\lambda_{\vartheta}^{{\rm J}/\psi})\;.\end{split}$ (6)

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,

$\frac{\langle
p^{2}\rangle/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+\langle\cos^{2}\Theta\rangle)\;,$
(7)

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

$\displaystyle\begin{split}\sigma^{\chi_{c(1,2)}}(p_{\rm
T}/M)&=\sigma^{\chi_{c(1,2}}_{\rm dir}(p_{\rm T}/M)\\\
&+B(\psi{\rm(2S)}\to\chi_{c(1,2)}\gamma)\,\sigma^{\psi{\rm(2S)}}(p_{\rm
T}/M)\;,\\\ \sigma^{{\rm J}/\psi}(p_{\rm T}/M)&=\sigma^{{\rm J}/\psi}_{\rm
dir}(p_{\rm T}/M)\\\ &+B(\chi_{c1}\to{\rm
J}/\psi\gamma)\,\sigma^{\chi_{c1}}(p_{\rm T}/M)\\\ &+B(\chi_{c2}\to{\rm
J}/\psi\gamma)\,\sigma^{\chi_{c2}}(p_{\rm T}/M)\\\ &+B(\psi{\rm(2S)}\to{\rm
J}/\psi X)\,\sigma^{\psi{\rm(2S)}}(p_{\rm T}/M)\;,\end{split}$ (8)

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

$\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)})}\;.$
(9)

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.

Figure 5: Comparison between the data points and the best fit curves obtained
in the UAU scenario: ${\rm J}/\psi$ (top left) and $\psi{\rm(2S)}$ (top right)
cross sections; $\chi_{c}$ cross sections and $\chi_{c2}/\chi_{c1}$ ratio
(bottom left); and $\lambda_{\vartheta}$ for the ${\rm J}/\psi$ and
$\psi{\rm(2S)}$ (bottom right).

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: Unpolarized (left) and polarized (right) cross sections in the UAU
scenario, normalized to unity at $p_{\rm T}/M=2$, with uncertainty bands
reflecting correlated variations in the fit parameters. The ${{}^{1}{\rm
S}_{0}^{[8]}}$ and ${{}^{3}{\rm S}_{1}^{[8]}}$ SDCs are also shown,
arbitrarily normalized, calculated at NLO with NNLO fragmentation
contributions [27].

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,

$\displaystyle\begin{split}\sigma_{\rm u}&\rightarrow\sigma({{}^{1}{\rm
S}_{0}^{[8]}})+\frac{4-4\lambda_{\vartheta}({{}^{3}{\rm
S}_{1}^{[8]}})}{4-\lambda_{\vartheta}({{}^{3}{\rm
S}_{1}^{[8]}})}\,\sigma({{}^{3}{\rm S}_{1}^{[8]}})\;,\\\ \sigma_{\rm
p}&\rightarrow\frac{3\lambda_{\vartheta}({{}^{3}{\rm
S}_{1}^{[8]}})}{4-\lambda_{\vartheta}({{}^{3}{\rm
S}_{1}^{[8]}})}\,\sigma({{}^{3}{\rm S}_{1}^{[8]}})\;,\end{split}$ (10)

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.

## References

## References

  * [1] G. T. Bodwin, E. Braaten, and P. Lepage  Phys. Rev. D51 (1995) 1125, doi:10.1103/PhysRevD.55.5853,10.1103/PhysRevD.51.1125, arXiv:hep-ph/9407339. [Erratum: Phys. Rev.D55,5853(1997)].
  * [2] P. Faccioli et al.  Phys. Lett. B736 (2014) 98, doi:10.1016/j.physletb.2014.07.006, arXiv:1403.3970.
  * [3] ATLAS Coll.  JHEP 09 (2014) 079, doi:10.1007/JHEP09(2014)079, arXiv:1407.5532.
  * [4] ATLAS Coll.  Phys. Rev. D87 (2013) 052004, doi:10.1103/PhysRevD.87.052004, arXiv:1211.7255.
  * [5] ATLAS Coll.  JHEP 07 (2014) 154, doi:10.1007/JHEP07(2014)154, arXiv:1404.7035.
  * [6] CMS Coll.  Phys. Rev. Lett. 114 (2015) 191802, doi:10.1103/PhysRevLett.114.191802, arXiv:1502.04155.
  * [7] CMS Coll.  Phys. Lett. B749 (2015) 14, doi:10.1016/j.physletb.2015.07.037, arXiv:1501.07750.
  * [8] P. Faccioli et al.  Submitted to Phys. Lett. B arXiv:XXX.YYY.
  * [9] CMS Coll.  Eur. Phys. J. C72 (2012) 2251, doi:10.1140/epjc/s10052-012-2251-3, arXiv:1210.0875.
  * [10] CMS Coll.  Phys. Lett. B743 (2015) 383, doi:10.1016/j.physletb.2015.02.048, arXiv:1409.5761.
  * [11] LHCb Coll.  Eur. Phys. J. C74 (2014) 3092, doi:10.1140/epjc/s10052-014-3092-z, arXiv:1407.7734.
  * [12] LHCb Coll.  Phys. Lett. B714 (2012) 215, doi:10.1016/j.physletb.2012.06.077, arXiv:1202.1080.
  * [13] LHCb Coll.  Phys. Lett. B718 (2012) 431, doi:10.1016/j.physletb.2012.10.068, arXiv:1204.1462.
  * [14] LHCb Coll.  Eur. Phys. J. C74 (2014) 2872, doi:10.1140/epjc/s10052-014-2872-9, arXiv:1403.1339.
  * [15] P. Faccioli, C. Lourenço, J. Seixas, and H. Wöhri  Eur. Phys. J. C69 (2010) 657, doi:10.1140/epjc/s10052-010-1420-5, arXiv:1006.2738.
  * [16] CMS Coll.  Phys. Lett. B727 (2013) 381, doi:10.1016/j.physletb.2013.10.055, arXiv:1307.6070.
  * [17] CMS Coll.  Phys. Rev. Lett. 110 (2013) 081802, doi:10.1103/PhysRevLett.110.081802, arXiv:1209.2922.
  * [18] CDF Coll.  Phys. Rev. Lett. 79 (1997) 572, doi:10.1103/PhysRevLett.79.572.
  * [19] M. Butenschön, Z.-G. He, and B. Kniehl  Phys. Rev. Lett. 114 (2015) 092004, doi:10.1103/PhysRevLett.114.092004, arXiv:1411.5287.
  * [20] K.-T. Chao et al.  Phys. Rev. Lett. 108 (2012) 242004, doi:10.1103/PhysRevLett.108.242004, arXiv:1201.2675.
  * [21] H.-S. Shao, Y.-Q. Ma, K. Wang, and K.-T. Chao  Phys. Rev. Lett. 112 (2014) 182003, doi:10.1103/PhysRevLett.112.182003, arXiv:1402.2913.
  * [22] H.-S. Shao  Comput. Phys. Commun. 198 (2016) 238, doi:10.1016/j.cpc.2015.09.011, arXiv:1507.03435.
  * [23] A. Petrelli et al.  Nucl. Phys. B514 (1998) 245, doi:10.1016/S0550-3213(97)00801-8, arXiv:hep-ph/9707223.
  * [24] M. Butenschön and B. Kniehl  Nucl. Phys. Proc. Suppl. 222 (2012) 151, doi:10.1016/j.nuclphysbps.2012.03.016, arXiv:1201.3862.
  * [25] M. Butenschön and B. Kniehl  Mod. Phys. Lett. A28 (2013) 1350027, doi:10.1142/S0217732313500272, arXiv:1212.2037.
  * [26] M. Butenschön and B. Kniehl  Phys. Rev. Lett. 108 (2012) 172002, doi:10.1103/PhysRevLett.108.172002, arXiv:1201.1872.
  * [27] G. T. Bodwin, H. S. Chung, U.-R. Kim, and J. Lee  Phys. Rev. Lett. 113 (2014) 022001, doi:10.1103/PhysRevLett.113.022001, arXiv:1403.3612.
  * [28] S. P. Baranov  Phys. Rev. D83 (2011) 034035, doi:10.1103/PhysRevD.83.034035.
  * [29] S. P. Baranov, A. V. Lipatov, and N. P. Zotov  Phys. Rev. D93 (2016) 094012, doi:10.1103/PhysRevD.93.094012, arXiv:1510.02411.
  * [30] K. A. Olive  Chin. Phys. C40 (2016) 100001, doi:10.1088/1674-1137/40/10/100001.
  * [31] LHCb Coll.  Eur. Phys. J. C75 (2015) 311, doi:10.1140/epjc/s10052-015-3502-x, arXiv:1409.3612.
  * [32] C. Lourenço and H. Wöhri  Phys. Rept. 433 (2006) 127, doi:10.1016/j.physrep.2006.05.005, arXiv:hep-ph/0609101.
  * [33] A. David  Phys. Lett. B644 (2007) 224, doi:10.1016/j.physletb.2006.10.071.
  * [34] BES Coll.  Phys. Rev. D62 (2000) 032002, doi:10.1103/PhysRevD.62.032002, arXiv:hep-ex/9909038.
  * [35] CLEO Coll.  Phys. Rev. D30 (1984) 1433, doi:10.1103/PhysRevD.30.1433.
  * [36] CLEO Coll.  Phys. Rev. D58 (1998) 052004, doi:10.1103/PhysRevD.58.052004, arXiv:hep-ex/9802024.
  * [37] P. Faccioli, C. Lourenço, J. Seixas, and H. Wöhri  Phys. Rev. D83 (2011) 096001, doi:10.1103/PhysRevD.83.096001, arXiv:1103.4882.
  * [38] P. Faccioli, C. Lourenço, and J. Seixas  Phys. Rev. D81 (2010) 111502, doi:10.1103/PhysRevD.81.111502, arXiv:1005.2855.
  * [39] S. Baranov  Phys. Rev. D93 (2016) 054037, doi:10.1103/PhysRevD.93.054037.
  * [40] S. Baranov, A. Lipatov, and N. Zotov  Eur. Phys. J. C75 (2015) 455, doi:10.1140/epjc/s10052-015-3689-x, arXiv:1508.05480.
  * [41] ATLAS Coll.  Eur. Phys. J. C76 (2016) 283, doi:10.1140/epjc/s10052-016-4050-8, arXiv:1512.03657.