11institutetext: LIP, Lisbon, Portugal 22institutetext: CERN, Geneva,
Switzerland 33institutetext: DESY, Hamburg, Germany

# From prompt to direct ${\rm J}/\psi$ production:
new insights on the $\chi_{{c}1}$ and $\chi_{{c}2}$ polarizations and feed-
down contributions from a global-fit analysis of mid-rapidity LHC data

Pietro Faccioli 11 Carlos Lourenço 22 Thomas Madlener 33

(Received: May 12, 2020 / Revised version: date)

###### Abstract

While the prompt ${\rm J}/\psi$ cross section and polarization have been
measured with good precision as a function of transverse momentum, $p_{\rm
T}$, those of the directly produced ${\rm J}/\psi$ are practically unknown,
given that the cross sections and polarizations of the $\chi_{{c}1}$ and
$\chi_{{c}2}$ mesons, large indirect contributors to ${\rm J}/\psi$
production, are only known with rather poor accuracy. The lack of precise
measurements of the $\chi_{{c}J}$ polarizations induces large uncertainties in
the level of their feed-down contributions to the prompt ${\rm J}/\psi$ yield,
because of the polarization-dependent acceptance corrections. The experimental
panorama of charmonium production can be significantly improved through a
consistent and model-independent global analysis of existing measurements of
${\rm J}/\psi$, $\psi{\rm(2S)}$ and $\chi_{c}$ cross sections and
polarizations, faithfully respecting all the correlations and uncertainties.
In particular, it is seen that the $\chi_{{c}J}$ polarizations and feed-down
fractions to ${\rm J}/\psi$ production have a negligible dependence on the
${\rm J}/\psi$ $p_{\rm T}$, with average values
$\lambda_{\vartheta}^{\chi_{{c}1}}=0.55\pm 0.23$,
$\lambda_{\vartheta}^{\chi_{{c}2}}=-0.39\pm 0.22$, $R^{\chi_{{c}1}}=(18.8\pm
1.4)\%$ and $R^{\chi_{{c}2}}=(6.5\pm 0.5)\%$. The analysis also shows that
$(67.2\pm 1.9)$% of the prompt ${\rm J}/\psi$ yield is due to directly-
produced mesons, of polarization constrained to remarkably small values,
$\lambda_{\vartheta}^{{\rm J}/\psi}=0.04\pm 0.06$.

###### pacs:

12.38.AwGeneral properties of QCD and 12.38.QkExperimental tests of QCD and
13.20.GdDecays of ${\rm J}/\psi$, $\Upsilon$, and other quarkonia

##  1 Introduction

The study of the production of heavy quarkonia, bound states of charm or
beauty quarks and antiquarks, offers the best experimental laboratory to
understand how quarks combine into hadrons, the least understood sector of
quantum chromodynamics (QCD), the theory of strong interactions bib:NRQCD ;
bib:QWGYR ; bib:QWG2011 . Indeed, the heaviness of these quarks implies that
they have relatively small velocities in the $Q\overline{Q}$ system, so that
their production can be theoretically described in two steps: a short-distance
regime where the $Q\overline{Q}$ is produced (typically through gluon fusion),
calculable using perturbative QCD, followed by an intrinsically non-
perturbative (long-distance) transition, particularly challenging to
understand at the present moment.

While it is known since long that experimental measurements of quarkonium
cross sections and polarizations should, in principle, lead to significant
progress in our QCD-based understanding of hadron formation, that prospect has
faced serious hurdles for a long time, first because of several challenges in
the execution of the measurements and the poor reliability of the resulting
data bib:EPJC69 , and second because most theoretical analyses of the data
have not properly taken into consideration the correlations and uncertainties
affecting the measurements, as explained in Refs. bib:FaccioliPLB736 ;
bib:Faccioli:PLB773 .

The high-quality measurements made over the last decade at the LHC, with a
remarkable level of detail and precision, provided a much improved
experimental situation. In particular, double-differential cross sections, in
transverse momentum, $p_{\rm T}$, and rapidity, $y$, have been measured in pp
collisions at $\sqrt{s}=7$, 8 and 13 TeV for the ${\rm J}/\psi$,
$\psi{\rm(2S)}$ and $\Upsilon$(nS) vector states, by ATLAS bib:ATLASdimuon ;
bib:ATLASpsi2S ; bib:ATLASYnS , CMS bib:BPH-14-001 ; bib:BPH-12-006 ;
bib:BPH-15-005  and LHCb bib:LHCb_psi_cs ; bib:LHCb_psip_cs ;
bib:LHCb_Upsilon_7_8 ; bib:LHCb_psi_Upsilon_8 ; bib:LHCb_jpsi_13 ;
bib:LHCb_psip_7_13 ; bib:LHCb_Upsilon_13 . Also the polarizations have been
measured for these states, at mid-rapidity by CMS bib:BPH-13-003 ;
bib:BPH-11-023  and at forward rapidity by LHCb bib:LHCb_psi_pol ;
bib:LHCb_psip_pol ; bib:LHCb_Upsilon_pol . In comparison with the very
significant experimental progress made at the LHC regarding the differential
cross sections and polarizations of vector quarkonia, the corresponding
knowledge of the $\chi_{c}$ and $\chi_{b}$ states has remained rather poor,
limited until recently to cross sections or cross section ratios affected by
relatively large uncertainties bib:ATLASchic ; bib:BPH-11-010 ; bib:BPH-13-005
; bib:LHCb_chic12_7 ; bib:LHCb_chic_7 . First experimental measurements on the
polarizations of the $\chi_{{c}1}$ and $\chi_{{c}2}$ states have recently been
reported by CMS bib:BPH-13-001 , effectively constraining the _difference_
between the polarizations of the two states but leaving their individual
values mostly unknown.

The polarization of the produced particle is not only, by itself, a crucial
element for the understanding of the underlying formation mechanisms, but is
also an essential element for the determination of the acceptance correction
of the corresponding cross section measurement. In fact, all the LHC
collaborations published their cross section measurements bib:ATLASdimuon ;
bib:ATLASpsi2S ; bib:ATLASYnS ; bib:BPH-14-001 ; bib:BPH-12-006 ;
bib:BPH-15-005 ; bib:LHCb_psi_cs ; bib:LHCb_psip_cs ; bib:LHCb_Upsilon_7_8 ;
bib:LHCb_psi_Upsilon_8 ; bib:LHCb_jpsi_13 ; bib:LHCb_psip_7_13 ;
bib:LHCb_Upsilon_13 ; bib:ATLASchic ; bib:BPH-11-010 ; bib:BPH-13-005 ;
bib:LHCb_chic12_7 ; bib:LHCb_chic_7  together with tables providing correction
factors that reflect how the central values of each measurement change when
the detection acceptance (computed, by default, for unpolarized production) is
recomputed for other (extreme) polarization scenarios. One cannot
underestimate the crucial importance of this knowledge, as the corrected
production yields can vary by more than 50% depending on the polarizations
assumed in the evaluation of the acceptances, an effect that vastly dominates
over the statistical and systematic measurement uncertainties. A consequence
of the incomplete experimental information on the $\chi_{{c}1}$ and
$\chi_{{c}2}$ polarizations is, therefore, that also their cross sections (and
cross section ratios) continue to carry a large associated uncertainty.

This lack of knowledge also has a strong effect on the understanding of ${\rm
J}/\psi$ production. In fact, while the _prompt_ ${\rm J}/\psi$ differential
cross section and polarization are the most precisely measured observables
among all measurements in the field of quarkonium production, their
interpretation in terms of properties of the physically-relevant _directly_
produced ${\rm J}/\psi$ mesons remains obscured by a large uncertainty, given
the very significant and not well known fraction of _indirect_ production from
$\chi_{c}$ feed-down decays.

This uncertainty also blurs the pattern of how the production of quarkonia of
different masses, binding energies and quantum numbers is modified by the QCD
medium produced in high-energy heavy-ion collisions bib:Faccioli:Fate .
Indeed, the feed-down fractions from heavier states are a crucial ingredient
in the observation of signatures of the sequential suppression mechanism,
according to which the production rate of quarkonium states should be
progressively suppressed, as the temperature of the medium increases,
following a hierarchy in the binding energy of the state bib:DPS ; bib:KKS .

In this paper we report the results of a global fit of mid-rapidity charmonium
measurements made by ATLAS and CMS at 7 TeV, including the recent $\chi_{c}$
polarization measurement, to derive the best possible determinations of the
$\chi_{c}$-to-${\rm J}/\psi$ feed-down fractions and $\chi_{c}$ polarizations,
and also of the properties of direct ${\rm J}/\psi$ production. The analysis
is completely independent of any quarkonium production theoretical model. It
only relies on the published measurements, which include the indirect
constraints that the differences between the ${\rm J}/\psi$ and
$\psi{\rm(2S)}$ data impose on the $\chi_{c}$ cross sections and
polarizations, through the feed-down contributions present in the ${\rm
J}/\psi$ case and absent in the $\psi{\rm(2S)}$ case. The transition
probabilities from heavier to lighter states needed in this work are all well
known and listed in the PDG tables bib:PDG . Our results are, therefore, fully
data driven.

Besides reporting the $\chi_{c}$ polarizations and feed-down contributions to
${\rm J}/\psi$ production determined by the global fit, we also discuss how
improved constraints on the $\chi_{c}$ polarizations could be obtained with
new charmonium measurements (other than direct measurements of the $\chi_{c}$
polarizations themselves).

##  2 Experimental data and fit parametrization

The data considered in our analysis are the ${\rm J}/\psi$, $\psi{\rm(2S)}$,
$\chi_{{c}1}$ and $\chi_{{c}2}$ differential cross sections measured by ATLAS
and CMS at 7 TeV bib:ATLASpsi2S ; bib:BPH-14-001 ; bib:ATLASchic , as well as
the $\chi_{{c}2}$ over $\chi_{{c}1}$ cross section ratio bib:BPH-11-010  and
the ${\rm J}/\psi$ and $\psi{\rm(2S)}$ polarizations bib:BPH-13-003  measured
by CMS at the same energy. All of these measurements have been reported as
functions of $p_{\rm T}$. To constrain the $\chi_{{c}J}$ polarizations we also
include the recent CMS measurement of the $\chi_{{c}2}$ over $\chi_{{c}1}$
yield ratio versus $|\cos\vartheta|$ ($\vartheta$ being the lepton emission
angle in the rest frame of the daughter ${\rm J}/\psi$) in three ${\rm
J}/\psi$ $p_{\rm T}$ bins. The total number of independent data points is 108
and they are all shown in Figs. 1, 3 and 3.

Figure 1: Mid-rapidity prompt quarkonium cross sections measured in pp
collisions at $\sqrt{s}=7$ TeV by ATLAS (red markers) bib:ATLASpsi2S ;
bib:ATLASchic  and CMS (blue markers) bib:BPH-14-001 . The inset shows the
$\chi_{{c}2}$ to $\chi_{{c}1}$ cross section ratio bib:BPH-11-010 . The curves
represent the result of the fit described in the text.

Figure 2: Polar anisotropy parameter $\lambda_{\vartheta}$, in the helicity
frame, measured by CMS in pp collisions at $\sqrt{s}=7$ TeV, for prompt ${\rm
J}/\psi$ and $\psi{\rm(2S)}$ dimuon decays bib:BPH-13-003 . Values
corresponding to two (${\rm J}/\psi$) or three ($\psi{\rm(2S)}$) rapidity bins
were averaged. The curves represent the result of the fit described in the
text.

Figure 3: The $\chi_{{c}2}/\chi_{{c}1}$ yield ratio vs. $|\cos\vartheta|$ (in
the helicity frame), for three ${\rm J}/\psi$ $p_{\rm T}$ bins (8–12, 12–18
and 18–30 GeV), as measured by CMS in pp collisions at $\sqrt{s}=8$ TeV
bib:BPH-13-001 . The curves represent the result of the fit described in the
text.

The ${\rm J}/\psi$ and $\psi{\rm(2S)}$ $p_{\rm T}$-differential cross sections
measured by ATLAS in the dimuon decay channel bib:ATLASdimuon  have not been
included in our global-fit analysis because the data, reported in eight
equidistant $|y|$ bins in the range $|y|<2$, show shapes as a function of
$p_{\rm T}$ that vary quite strongly among the $|y|$ bins. These variations,
in particular between the $|y|<0.25$ and $0.25<|y|<0.5$ bins, are clearly not
statistical fluctuations and significantly exceed what one could expect from
the reported (systematic) uncertainties, pointing to an internal inconsistency
affecting these two data sets.

To derive global observables from the analysis, a parametrization of yields
and polarizations is necessary, because the kinematic binning of the reported
distributions is not identical among all data sets and quarkonium states.
Moreover, for the comparison/combination of results concerning objects of
different mass scales the absolute transverse momentum, in which the data are
binned, is not the best variable: it is, in fact, preferable to use a
relative, dimensionless variable, in our case chosen as $p_{\rm T}/M$, the
ratio between the $p_{\rm T}$ and the mass of the quarkonium state. The
convenience of this variable will become apparent in the next paragraphs,
where we sometimes use the definition $\xi\equiv p_{\rm T}/M$ to simplify the
notation.

The $p_{\rm T}/M$ dependences of the considered particle yields are
parametrized using a shape function $g(\xi)$, normalized to unity at the
arbitrary reference point $\xi^{*}=5$, a value close to the centre of gravity
of the data: $g(\xi)=h(\xi)/h(\xi^{*})$, with

$h(\xi)=\xi\cdot\bigg{(}1+\frac{1}{\beta-2}\cdot\frac{\xi^{2}}{\gamma}\bigg{)}^{-\beta}\,.$
(1)

This functional shape describes very well the quarkonium transverse momentum
distributions in different kinematic domains bib:HERAb ; bib:Faccioli:PLB773 ;
bib:Faccioli:EPJC78p118 . The parameter $\gamma$ (having the meaning of the
average $p_{\rm T}/M$ squared) defines the function in the low-$p_{\rm T}$
turn-on region and is only mildly sensitive to the data we are considering;
hence, only one $\gamma$ parameter will be considered, common for all states
and polarization configurations and treated as global free parameter of the
fit. The power-law exponent $\beta$ describes the asymptotic high-$p_{\rm T}$
behaviour: $g\propto\xi^{-(2\beta-1)}$ for $\xi\gg\sqrt{\gamma\,(\beta-2)}$.
It is, therefore, the shape parameter actually characterising each considered
(sub-)process.

To ensure correct feed-down relations between the charmonium family members,
we include a detailed account of how the mother’s momentum and polarization
are transferred to the daughter in the relevant decays:
$\psi{\rm(2S)}\to\chi_{c1,2}\;\gamma$; $\psi{\rm(2S)}\to{\rm J}/\psi\;X$;
$\chi_{c1,2}\to{\rm J}/\psi\;\gamma$. The rule for the momentum propagation
from mother to daughter is, on average, $p_{\rm T}/m=P_{\rm T}/M$, where $M$
($m$) and $P_{\rm T}$ ($p_{\rm T}$) are, respectively, the mass and laboratory
transverse momentum of the mother (daughter) particle bib:Faccioli:PLB773 .
The polarization transfer rules were calculated in the electric dipole
approximation and precisely account for the observable dilepton distribution
with no need of higher-order terms bib:chiPol . In particular, the
$\lambda_{\vartheta}^{\chi_{{c}1}}$ and $\lambda_{\vartheta}^{\chi_{{c}2}}$
(polar anisotropy) parameters refer to the shapes of the corresponding
daughter-${\rm J}/\psi$’s dilepton decay distributions, which are the ones
directly measured and fully reflect the $\chi_{c}$ polarization state, while
being insensitive to the uncertain contributions of higher-order photon
multipoles bib:chiPol . Consequently, the terms _longitudinal_ and
_transverse_ will here always refer to the yields of events where the daughter
${\rm J}/\psi$ has, respectively, angular momentum projection $J_{z}=0$
($\lambda_{\vartheta}=-1$) and $J_{z}=\pm 1$ ($\lambda_{\vartheta}=+1$), even
if the mother $\chi_{{c}J}$ has a very different correspondence between
angular momentum configuration and polar anisotropy parameter: for example, a
$\chi_{{c}1}$ state with $J_{z}=0$ or $J_{z}=\pm 1$ leads to
$\lambda_{\vartheta}=+1$ and $-1/3$, respectively. All polarizations are
considered and defined in the centre-of-mass helicity frame.

The polarizations are parametrized as functions of $p_{\rm T}/M$ by
considering for each directly produced state its _total_ and _longitudinal_
cross sections, both parametrized as described by Eq. 1. The ratio of
longitudinal to total cross section, the longitudinal fraction, is calculated
from the polar anisotropy parameter for the two-body decay distribution of the
directly produced state, $\lambda_{\vartheta}^{\mathrm{dir}}(\xi)$, as
$f_{\mathrm{long}}(\xi)=\left[1-\lambda_{\vartheta}^{\mathrm{dir}}(\xi)\right]/\left[3+\lambda_{\vartheta}^{\mathrm{dir}}(\xi)\right]$.
With no further input or prior information than the charmonium data
themselves, a complete parametrization of cross sections and polarizations of
the four considered charmonium states (${\rm J}/\psi$, $\chi_{{c}1}$,
$\chi_{{c}2}$ and $\psi{\rm(2S)}$) would require eight $\beta$ parameters,
describing the $p_{\rm T}/M$ dependences of the four total direct-production
cross sections and those of the corresponding longitudinal cross sections. The
remaining eight “parameters of interest” are the normalizations of the four
direct-production cross sections and the four corresponding polar anisotropy
parameters, all conventionally considered at the reference point $(p_{\rm
T}/M)^{*}=5$.

However, not all of the shape parameters are varied independently in the fit.
Some of them are treated as common to several states and/or polarized
subprocesses, on the basis of data-driven considerations or basic physics
considerations. The picture of mid-rapidity differential cross sections and
polarizations shows, in fact, a characteristic simplicity. As discussed in
Refs. bib:Faccioli:PLB773 ; bib:Faccioli:EPJC78p118 , the production cross
sections of the ${{}^{3}{\rm S}_{1}}$ and ${{}^{3}{\rm P}_{J}}$ quarkonium
states measured by ATLAS and CMS, at both 7 and 13 TeV, follow remarkably
uniform patterns as a function of $p_{\rm T}/M$. Such scaling is expected,
from dimensional analysis considerations bib:Faccioli:EPJC78p118 , if the
fundamental production processes are identical, in quality and relative
contributions, for all states, which is a conceivable scenario when we only
consider states of identical quantum numbers. However, an interesting, albeit
unexpected, aspect of such “universal” picture of mid-rapidity production is
that there are currently no indications of a difference between the $p_{\rm
T}/M$ distribution shapes of the P-wave states and those of the S-wave states
bib:Faccioli:EPJC78p268 . Moreover, the measured charmonium and bottomonium
decay distributions indicate similar polarizations ($\lambda_{\vartheta}$ in
the helicity frame) for all vector states, independently of their different
feed-down contributions from $\chi_{c}$ and $\chi_{b}$ states. Finally, all
polarizations are perfectly compatible with being independent of $p_{\rm
T}/M$.

These indications of uniform kinematic behaviours are particularly significant
when we consider the large mass variation between the charmonium and
bottomonium families. Moreover, the uniformity of such scaling patterns has
been verified at both 7 and 13 TeV bib:Faccioli:EPJC78p118 . It is, therefore,
reasonable to adopt such quarkonium-wide observations as simplifying
assumptions for the parametrization of our charmonia-only fit.

We start by imposing that the _directly-produced_ vector mesons ${\rm J}/\psi$
and $\psi{\rm(2S)}$ have identical $p_{\rm T}/M$-dependent kinematic patterns,
while keeping, obviously, two independent yield normalizations at $(p_{\rm
T}/M)^{*}$. This assumption, suggested and supported by the observed
universality of the $p_{\rm T}/M$-differential cross sections and
polarizations, is also adopted by construction in theoretical models based on
the factorization hypothesis (such as NRQCD), where the kinematics-dependent
short distance cross section terms depend only on the heavy quark mass but not
on the final bound state bib:NRQCD . Therefore, in the fit, a single
parameter, $\beta_{\mathrm{total}}^{\psi}$, describes the $p_{\rm T}/M$
dependences of the total ${\rm J}/\psi$ and $\psi{\rm(2S)}$ direct-production
cross sections, while another one, $\beta_{\mathrm{long}}^{\psi}$, represents
the shapes of the two corresponding longitudinal cross sections. The
assumption that both states have the same _direct-production_ polarization is
reflected in the choice of one common free parameter for the polar anisotropy
parameter at the reference point,
$\lambda_{\vartheta}^{\psi,\mathrm{dir}}(\xi^{*})$. These constraints ensure
that any difference observed between the ${\rm J}/\psi$ and $\psi{\rm(2S)}$
$p_{\rm T}/M$ distributions and polarizations is attributed to the
$\chi_{{c}1}$ and $\chi_{{c}2}$ feed-down contributions: ${\rm J}/\psi$ and
$\psi{\rm(2S)}$ measurements become, therefore, indirect constraints on the
$\chi_{{c}1}$ and $\chi_{{c}2}$ cross sections and polarizations.

Given the relatively large uncertainties affecting some of the data sets used
in the analysis, most notably the $\chi_{{c}J}$ cross sections and the
$\psi{\rm(2S)}$ polarization, it seems judicious to refrain from having many
free parameters in the fit model, at least in the “default” analysis.
Therefore, we further require that the production cross sections of all _four_
charmonium states follow explicitly the $p_{\rm T}/M$ universality suggested
by the ensemble of mid-rapidity quarkonium data: one common parameter
describes the asymptotic power-law behaviour of their total cross sections,
$\beta_{\mathrm{total}}^{\chi_{{c}1}}=\beta_{\mathrm{total}}^{\chi_{{c}2}}=\beta_{\mathrm{total}}^{\psi}$.
In order not to limit the range of possible physical outcomes, no further
constraint is imposed on the polarizations: three independent parameters
represent the polar anisotropies of the directly produced states at the
reference $(p_{\rm T}/M)^{*}$:
$\lambda_{\vartheta}^{\psi,\mathrm{dir}}(\xi^{*})$,
$\lambda_{\vartheta}^{\chi_{{c}1},\mathrm{dir}}(\xi^{*})$, and
$\lambda_{\vartheta}^{\chi_{{c}2},\mathrm{dir}}(\xi^{*})$. Furthermore, three
independent $\beta$ exponents characterize the shapes of the longitudinal
cross sections of the directly produced states: $\beta_{\mathrm{long}}^{\psi}$
(common to the ${\rm J}/\psi$ and $\psi{\rm(2S)}$ states),
$\beta_{\mathrm{long}}^{\chi_{{c}1}}$ and
$\beta_{\mathrm{long}}^{\chi_{{c}2}}$. The set of “parameters of interest” is
completed by the normalizations of the four (total) direct production cross
sections, defined in the fit model as the ${\rm d}\sigma/{\rm d}p_{\rm T}$
values (in nb/GeV) at $(p_{\rm T}/M)^{*}$. It is worth noting that the
measured cross sections and cross section ratios have been published in the
form of products of the production cross sections and branching fractions into
the detected decay channel, while our analysis always considers the pure
production cross sections, obtained dividing the measured values by the
relevant branching fractions, taken from Ref. bib:PDG .

The baseline fit model we have just described has four free $\beta$ exponents
and, hence, will be referred to in the remaining of this article by the
“$4\beta$” label. We have also repeated the analysis in two alternative fit
configurations, analogously labelled as the “$6\beta$” and “$1\beta$” models.
In the more unconstrained $6\beta$ scenario, the $\chi_{{c}1}$ and
$\chi_{{c}2}$ total cross sections are free to have power-law $p_{\rm T}$
trends different from each other and from that of the ${\rm J}/\psi$ and
$\psi{\rm(2S)}$ mesons, so that the results will provide a test of the
importance of the $p_{\rm T}/M$-universality we have assumed in the baseline
model. The more constrained $1\beta$ variant imposes a common power-law
exponent on all total _and longitudinal_ cross sections, so that the four
charmonia are assumed to have identical (fully universal) $p_{\rm T}/M$
shapes, not only in the differential cross sections but also in the
polarizations. In other words, in this scenario only the _magnitudes_ of the
cross sections and polarizations can be different among the four (directly-
produced) states, being therefore an effective way to directly obtain $p_{\rm
T}/M$-averaged values of the polarizations, of the yield ratios, and of the
feed-down fractions. The number of free power-law exponents is the only
difference between the $4\beta$, $6\beta$ and $1\beta$ variants, the three
longitudinal fractions and four direct cross section normalizations at
$(p_{\rm T}/M)^{*}$ remaining free parameters in all options.

While the ${\rm J}/\psi$ and $\psi{\rm(2S)}$ polarization measurements impose,
as discussed, indirect constraints on the $\chi_{{c}1}$ and $\chi_{{c}2}$
polarizations (mainly on their sum), direct constraints are provided (mainly
on their difference) by the three $\chi_{{c}2}$ over $\chi_{{c}1}$ yield
ratios versus $|\cos\vartheta|$, measured by CMS in three ranges of ${\rm
J}/\psi$ transverse momentum. Those data points are parametrized with the
expression

$R_{i}\;\frac{1+\lambda_{\vartheta}^{\chi_{{c}2}}(\xi_{i})\cos^{2}\vartheta}{1+\lambda_{\vartheta}^{\chi_{{c}1}}(\xi_{i})\cos^{2}\vartheta}\,,$
(2)

where $\xi_{i}$ ($i=1,2,3$) are the three average-$p_{\rm T}/M$ values of the
measurement, $R_{i}$ are the three ratio normalizations, treated as
independent fit parameters, and $\lambda_{\vartheta}^{\chi_{{c}J}}(\xi_{i})$
($J=1,2$) are the polar anisotropy parameters of the ${\rm J}/\psi$ from
$\chi_{{c}1}$ or $\chi_{{c}2}$, calculated at the three $p_{\rm T}/M$ values
using the parametrized longitudinal and total $\chi_{{c}1}$ or $\chi_{{c}2}$
cross sections, defined above, also including the small $\psi{\rm(2S)}$ feed-
down contribution.

Correlations between the data points are taken into account by defining a
number of nuisance parameters. First, independently for ATLAS and CMS, all the
cross sections are scaled by a global factor that, while being a free
parameter in the fit, is constrained by a Gaussian function of mean unity and
width equal to the relative uncertainty of the integrated luminosity, reported
in the experimental publications. In other words, the fit quality incurs a
penalty reflecting the difference between the best-fit scale factor and unity,
normalized by the uncertainty. By equally scaling all the cross sections of a
given experiment, these two nuisance parameters induce a correlation between
the $\psi{\rm(2S)}$, $\chi_{{c}1}$ and $\chi_{{c}2}$ cross sections measured
by ATLAS and between the ${\rm J}/\psi$ and $\psi{\rm(2S)}$ cross sections
measured by CMS. Second, also the branching factions needed to convert the
measured values to production cross sections are analogously scaled by
Gaussian-constrained nuisance parameters, the Gaussian widths being the
relative uncertainties reported in Ref. bib:PDG . This second set of nuisance
parameters induces correlations between the ATLAS and CMS data. It turns out
that all the post-fit nuisance parameters are identical to unity, except for
the one related to the $\psi{\rm(2S)}\to\mu\mu$ branching fraction, which
deviates from unity by 1%, a negligible departure given that the uncertainty
is six times larger.

Another, very important, source of correlations between all data points is the
dependence of the detection acceptances on the polarization. For each set of
parameter values considered in the fit scan, the expected values of the
polarizations and cross sections are calculated, for all states, as functions
of $p_{\rm T}$, using the shape-parametrization functions described above. The
expected $\lambda_{\vartheta}$ values can be immediately compared to the
measured ones, for the determination of the corresponding $\chi^{2}$ terms,
while for the calculation of the cross section $\chi^{2}$ terms we first scale
the measured cross sections by acceptance-correction factors calculated for
the $\lambda_{\vartheta}$ value under consideration. These correction factors
are computed, for each data point, using the tables published by the
experiments (for exactly this purpose) for the cross sections of particles
produced with fully transverse or fully longitudinal polarizations.

##  3 Results of the global fit analysis

As can be appreciated from the information presented in Table 1, the fit
quality is excellent in all three fit variants. No tension or difference in
trends is visible between the 108 data points and the best fit curves, as
shown in Figs. 1, 3 and 3 for the baseline $4\beta$ case. The smallness of the
$\chi^{2}$ per degree of freedom, $\chi^{2}/{\rm ndf}=40/93$, corresponding to
an exceptionally (and suspiciously) good fit $\chi^{2}$ probability, points to
the existence of unaccounted correlations between systematic uncertainties in
the data points. Indeed, it is very likely that a fraction of the systematic
uncertainties assigned in the experimental publications to each of the $p_{\rm
T}$ bins actually reflects an effect that commonly affects a broad region of
the distribution, leaving its shape essentially unchanged, so that the true
point-to-point uncorrelated uncertainties are somewhat smaller than those we
have used. In any case, it is certainly informative to compare the
$\chi^{2}/{\rm ndf}$ value of the baseline analysis with those of the two
variants mentioned before: $\chi^{2}/{\rm ndf}=39/91$ ($6\beta$) and 43/96
($1\beta$). We see that, given the precision of the presently-available
experimental inputs, there is no advantage in using the fit model with two
more free parameters. Indeed, according to the Akaike information criterion
(AIC) bib:AIC , the likelihood of the $6\beta$ model is much smaller than that
of the $4\beta$ model.

Table 1: Fit quality information. In all three fit variants there are 108 data points and 10 nuisance parameters. | $6\beta$ | $4\beta$ | $1\beta$
---|---|---|---
No. of free parameters | 27 | 25 | 22
No. of degrees of freedom | 91 | 93 | 96
Fit $\chi^{2}$ | 39 | 40 | 43

Interestingly, the more constrained $1\beta$ fit model, which imposes a common
value to all the six $\beta$ exponents, provides a description of the data
that is essentially as good as that of the baseline option, despite having
three less free parameters. The slightly worse fit $\chi^{2}$ is compensated
by the extra simplicity of the model, leading to a large increase in the AIC
relative likelihood. We will refrain from using this observation to highlight
the implication that all charmonium states are seemingly produced with
identical kinematical patterns, both in terms of cross sections and in terms
of polarizations, so that a rather straightforward model is able to faithfully
reproduce all the data points considered in our study. Instead, we simply
argue that this remarkable observation should trigger further experimental
measurements of quarkonium cross sections and polarizations, with
significantly improved precision, so that the validity of the $1\beta$ model
can be scrutinised much more accurately. Only then we will be able to conclude
if this strongly constrained fit is merely a very effective and economic
description of the presently existing data, providing a reliable computation
of $p_{\rm T}/M$-averaged results, or if we are seeing a smoking-gun signature
of a fully-universal scenario, reflecting a deeper symmetry at the core of
hadron formation than assumed in today’s theories of quarkonium production.

Table 2: Values of the fitted parameters of interest, in the three considered scenarios. | $6\beta$ | $4\beta$ | $1\beta$
---|---|---|---
Direct-production ${\rm d}\sigma/{\rm d}p_{\rm T}$ (nb/GeV), at $p_{\rm
T}/M=5$
${\rm J}/\psi$ | $1.274\pm 0.059$ | $1.283\pm 0.059$ | $1.281\pm 0.058$
$\psi{\rm(2S)}$ | $0.230\pm 0.010$ | $0.232\pm 0.010$ | $0.233\pm 0.010$
$\chi_{{c}1}$ | $1.008\pm 0.097$ | $0.966\pm 0.092$ | $1.008\pm 0.085$
$\chi_{{c}2}$ | $0.604\pm 0.071$ | $0.577\pm 0.068$ | $0.617\pm 0.063$
Direct $\lambda_{\vartheta}$ at $p_{\rm T}/M=5$
$\lambda_{\vartheta}^{\psi,\rm dir}$ | $-0.005\pm 0.072$ | $\hskip 5.00002pt0.022\pm 0.062$ | $\hskip 5.00002pt0.040\pm 0.060$
$\lambda_{\vartheta}^{\chi_{c1},\rm dir}$ | $\hskip 5.00002pt0.504\pm 0.303$ | $\hskip 5.00002pt0.371\pm 0.268$ | $\hskip 5.00002pt0.521\pm 0.247$
$\lambda_{\vartheta}^{\chi_{c2},\rm dir}$ | $-0.402\pm 0.276$ | $-0.533\pm 0.250$ | $-0.392\pm 0.233$
Kinematical-dependence parameters
$\gamma$ | $0.642\pm 0.176$ | $0.643\pm 0.157$ | $0.601\pm 0.148$
$\beta_{\mathrm{total}}^{\psi}$ | $3.358\pm 0.032$ | $\left.\rule[6.45831pt]{0.0pt}{15.0694pt}\right\\}3.379\pm 0.022$ | $\left.\rule[6.45831pt]{0.0pt}{37.88878pt}\right\\}3.385\pm 0.021$
$\beta_{\mathrm{total}}^{\chi_{c1}}$ | $3.444\pm 0.122$
$\beta_{\mathrm{total}}^{\chi_{c2}}$ | $3.562\pm 0.171$
$\beta_{\mathrm{long}}^{\psi}$ | $3.494\pm 0.142$ | $3.431\pm 0.097$
$\beta_{\mathrm{long}}^{\chi_{c1}}$ | $3.178\pm 0.715$ | $3.530\pm 0.599$
$\beta_{\mathrm{long}}^{\chi_{c2}}$ | $3.429\pm 0.322$ | $3.674\pm 0.273$

The fitted values of all the parameters of interest are collected in Table 2,
for the three variants: $6\beta$, $4\beta$ and $1\beta$. The obtained
$\chi_{c}$-to-${\rm J}/\psi$ and $\psi{\rm(2S)}$-to-${\rm J}/\psi$ feed-down
fractions are shown in Fig. 4. The results of the baseline ($4\beta$) fit are
represented by the $p_{\rm T}/M$-dependent central values (solid lines),
enveloped by filled bands of widths equal to the 68.3% confidence level
uncertainties, obtained by integrating the multivariate normal distribution
representing the joint probability distribution of all parameters over the
physical domains of all the variables not shown in the figure. The results of
the $1\beta$ fit option, independent of $p_{\rm T}/M$ by construction, are
also shown in the figure, as dashed lines (central values) surrounded by empty
rectangles (uncertainties); the corresponding numerical values are collected
in Table 3.

Figure 4: The fractions of the total prompt ${\rm J}/\psi$ production rate due
to feed-down decays from the $\chi_{{c}1}$, $\chi_{{c}2}$ and $\psi{\rm(2S)}$
mesons, as a function of $p_{\rm T}/M$. Table 3: $p_{\rm T}/M$-averaged values
of the feed-down fractions, as determined in the $1\beta$ global-fit analysis.
Virtually identical values are obtained in the $4\beta$ fit option, for
$p_{\rm T}/M=5$. The derived direct ${\rm J}/\psi$ fraction is $(67.2\pm
1.9)$%. Feed-down fractions (%)
---
$\chi_{{c}1}\to{\rm J}/\psi$ | $18.8\pm 1.4$
$\chi_{{c}2}\to{\rm J}/\psi$ | $\hskip 5.00002pt6.5\pm 0.5$
$\chi_{{c}1}+\chi_{{c}2}\to{\rm J}/\psi$ | $25.3\pm 1.8$
$\psi{\rm(2S)}\to{\rm J}/\psi$ | $\hskip 5.00002pt7.5\pm 0.3$
$\psi{\rm(2S)}\to\chi_{{c}1}$ | $\hskip 5.00002pt2.2\pm 0.2$
$\psi{\rm(2S)}\to\chi_{{c}2}$ | $\hskip 5.00002pt3.4\pm 0.3$
Figure 5: The polarization parameter $\lambda_{\vartheta}$ of the
$\chi_{{c}1}$ and $\chi_{{c}2}$ mesons (left), as well as of the ${\rm
J}/\psi$ mesons directly produced (same as of the $\psi{\rm(2S)}$) and
produced in $\chi_{{c}1}$ plus $\chi_{{c}2}$ decays (right), as a function of
$p_{\rm T}/M$.

The corresponding results for the polarizations ($\lambda_{\vartheta}$ in the
helicity frame) are presented in Fig. 5, in the left panel for the
$\chi_{{c}1}$ and $\chi_{{c}2}$ mesons and in the right panel for the
directly-produced ${\rm J}/\psi$ mesons. The right panel also shows the
polarization of the ${\rm J}/\psi$ mesons produced in decays of both
$\chi_{c}$ mesons, an observable determined with a better precision than each
of the individual (anti-)correlated polarizations shown on the left panel. As
in Fig. 4, the dashed lines and empty rectangles represent the $p_{\rm
T}/M$-independent results obtained in the $1\beta$ fit, effectively
representing averages over $p_{\rm T}/M$ of the $4\beta$ fit results, shown as
solid lines and filled bands. The respective numerical values are collected in
Table 4, which also shows the derived polarization of promptly produced ${\rm
J}/\psi$ mesons, naturally intermediate between the values of the directly
produced mesons and of those emitted in the $\chi_{{c}J}$ decays.

Table 4: Polarizations determined in the $4\beta$ fit variant (for $p_{\rm
T}/M=5$) and in the $1\beta$ option (averaged over $p_{\rm T}/M$).
$\lambda_{\vartheta}$
---
| $4\beta$ | $1\beta$
$\chi_{{c}1}$ | $\hskip 5.00002pt0.46\pm 0.23$ | $\hskip 5.00002pt0.55\pm 0.23$
$\chi_{{c}2}$ | $-0.52\pm 0.24$ | $-0.39\pm 0.22$
${\rm J}/\psi$ from $\chi_{{c}1}+\chi_{{c}2}$ | $\hskip 5.00002pt0.19\pm 0.18$ | $\hskip 5.00002pt0.27\pm 0.19$
Direct ${\rm J}/\psi$ = $\psi{\rm(2S)}$ | $0.022\pm 0.062$ | $\hskip 5.00002pt0.040\pm 0.060$
Prompt ${\rm J}/\psi$ | $\hskip 5.00002pt0.048\pm 0.037$ | $\hskip 5.00002pt0.087\pm 0.024$

Our global-fit analysis provides significant improvements in the determination
of several interesting observables. To start with, individual (purely data-
driven) values of the $\chi_{{c}1}$ and $\chi_{{c}2}$ polarizations are
extracted, as reported in Table 4. Equally important, the feed-down fractions
(shown in Table 3) are determined with a rather good precision, of around 10%,
even for the small fractions of $\chi_{{c}1}$ and $\chi_{{c}2}$ production
yields due to radiative decays of $\psi{\rm(2S)}$ mesons. Finally, the ratio
between the $\chi_{{c}2}$ and the $\chi_{{c}1}$ cross sections (times the
corresponding $\chi_{{c}J}\to{\rm J}/\psi\,\gamma$ branching fractions),
becomes much more precisely determined:
$B\sigma(\chi_{{c}2})/B\sigma(\chi_{{c}1})=0.343\pm 0.024$. Figure 7 provides
a graphical illustration of the improvement reached in the precision of the
$\chi_{{c}2}$ to $\chi_{{c}1}$ cross section ratio. The left panel shows the
very strong dependence of the original ATLAS and CMS measurements (and of the
level of their mutual compatibility) on the unknown $\chi_{{c}J}$
polarizations. The $\chi_{{c}J}$ polarization constraints contributing
(directly and indirectly) to our global fit strongly reduce the uncertainty
associated with the polarization dependence of the acceptance correction,
leading to the rather well aligned points shown on the right panel, where the
acceptance corrections reflect the best-fit polarization scenario (curiously,
very close to the $J_{z}=0$ extreme). The filled band represents the final
result and the difference between the widths of the filled and dashed bands
reflects the residual polarization uncertainty, a rather small effect in
comparison with the impact seen in the left panel.

Figure 6: Ratio between the $\chi_{{c}2}$ and $\chi_{{c}1}$ cross sections,
times the corresponding $\chi_{{c}J}\to{\rm J}/\psi\,\gamma$ branching
fractions, as a function of $p_{\rm T}/M$, for two extreme polarization
scenarios (left) and for the polarizations determined in our global-fit
analysis (right).

Figure 7: Two-dimensional distributions showing the correlations between the
$\lambda_{\vartheta}^{\chi_{{c}2}}$ and $\lambda_{\vartheta}^{\chi_{{c}1}}$
polarizations (left) and between the $\chi_{{c}2}$ and $\chi_{{c}1}$ feed-down
contributions to the prompt ${\rm J}/\psi$ production rate (right). The
results of the three fit variants are represented by the solid ($4\beta$),
dashed ($1\beta$) and dotted-dashed ($6\beta$) lines.

Among the remaining physical results, particularly interesting is the
polarization of the $\psi{\rm(2S)}$ and of the directly produced ${\rm
J}/\psi$, shown by the pink band in Fig. 5, constrained to be “zero” with a
previously unseen precision and no signs of momentum dependence. This is a
unique result for a vector state; both Drell–Yan dileptons Lam:1978pu ;
bib:Faccioli:PRLgeneralizeLT ; Guanziroli:1987rp ; Conway:1989fs ; Zhu:2006gx
; Zhu:2008sj  and vector boson Mirkes:1994eb ; Aaltonen:2011nr ;
Khachatryan:2015paa ; Abbott:2000aj ; Acosta:2003jm ; Acosta:2005dn ;
Chatrchyan:2011ig ; ATLAS:2012au  polarizations are known to be significantly
non-zero and momentum dependent, as are those of low-$p_{\rm T}$ quarkonia
Abt:2009nu ; Brown:2000bz . In fact, there are only two ways to obtain a
vector particle in an angular momentum state having zero observable
polarization. One is to prepare a mixture of two (or more) very different,
strongly polarized states: as demonstrated in Ref.
bib:Faccioli:PRLgeneralizeLT , for a _single_ , individually produced angular
momentum state there is always a polarization axis with respect to which
$\lambda_{\vartheta}=+1$. The _exact_ compensation of two strongly polarized
production processes, leaving no margin for a residual momentum-dependent
deviation of $\lambda_{\vartheta}$ from zero, would be an astonishing
coincidence; in fact, it could only be attributed to the existence of unknown
symmetries governing charmonium production, at least in the mid-rapidity limit
bib:Faccioli:NRQCD2019 . The other possibility is that the ${\rm J}/\psi$
originates from the decay of a $J=0$ state, as expected to happen in the
production from the ${}^{1}{\rm S}_{0}$ colour octet term in NRQCD (and in the
feed-down from $\chi_{{c}0}$). In this possible subprocess, while the
polarization continues to be naturally fully transverse along the direction of
the recoil gluon (or photon), it undergoes a complete rotational smearing when
seen in the experimental polarization frame, whose $z$ axis is fully
decorrelated from such natural direction when the mother-daughter mass
difference is small. While the production via ${}^{1}{\rm S}_{0}$ octet is
foreseen, it is not naturally predicted to be the only, dominating mechanism.
In either case, a precise confirmation of a $p_{\rm T}$-independent
unpolarized scenario has strong and rather remarkable physical implications.

Figure 8: Two-dimensional distributions showing the correlations between the
$\lambda_{\vartheta}$ polarization parameters of the $\chi_{{c}1}$ or
$\chi_{{c}2}$ mesons and the corresponding values for the ${\rm J}/\psi$
(left), the $\psi{\rm(2S)}$ (middle) and the difference between the two
(right). The results of the three fit variants are represented by the solid
($4\beta$), dashed ($1\beta$) and dotted-dashed ($6\beta$) lines.

The variations represented by the bands in Figs. 4 and 5 are generally
correlated. Figure 7 shows the correlations between the $\chi_{{c}1}$ and
$\chi_{{c}2}$ polarizations (left) and between the $\chi_{{c}1}$ and
$\chi_{{c}2}$ feed-down contributions to ${\rm J}/\psi$ production (right),
for the three fit variants. Particularly interesting are the correlations
shown in Fig. 8, where we can see that a significantly improved knowledge of
the $\chi_{{c}J}$ polarizations will derive from new, precise measurements of
the $\psi{\rm(2S)}$ polarization and, above all, of the difference between the
${\rm J}/\psi$ and $\psi{\rm(2S)}$ polarizations. This latter measurement can
be performed, for example, by determining the ratio between the ${\rm J}/\psi$
and $\psi{\rm(2S)}$ angular distributions in a given $p_{\rm T}/M$ interval,
with the cancellation of a large part of the important systematic
uncertainties related to acceptance and efficiency descriptions. It would,
therefore, represent a particularly clean constraint on the sum of the
$\chi_{{c}1}$ and $\chi_{{c}2}$ polarizations. Finally, Fig 9 shows that the
$\chi_{{c}1}$ and $\chi_{{c}2}$ polarizations are also correlated with the
corresponding feed-down fractions to ${\rm J}/\psi$ production, so that
precise measurements of those feed-down fractions will also reduce the
$\lambda_{\vartheta}^{\chi_{{c}1}}$ and $\lambda_{\vartheta}^{\chi_{{c}2}}$
uncertainties.

Figure 9: Two-dimensional distributions showing the correlation between the
$\lambda_{\vartheta}^{\chi_{{c}1}}$ ($\lambda_{\vartheta}^{\chi_{{c}2}}$)
observable and the $\chi_{{c}1}$ ($\chi_{{c}2}$) feed-down contribution to the
prompt ${\rm J}/\psi$ production rate. The results of the three fit variants
are represented by the solid ($4\beta$), dashed ($1\beta$) and dotted-dashed
($6\beta$) lines.

##  4 Summary

We have performed a global study of charmonium production at mid-rapidity and
LHC energies, with the aim to improve the current knowledge of $\chi_{{c}J}$
polarizations and feed-down fractions to ${\rm J}/\psi$ production, and to
extract the kinematic properties of the directly produced ${\rm J}/\psi$. The
analysis is fully data driven, not relying on any theoretical inputs. It uses
LHC data on ${\rm J}/\psi$, $\psi{\rm(2S)}$, $\chi_{{c}1}$ and $\chi_{{c}2}$
differential cross sections and decay angular distributions, measured by ATLAS
and CMS at mid-rapidity, in pp collisions at $\sqrt{s}=7$ TeV.

A first result of the analysis is that all polarizations and cross section
ratios are found to be perfectly compatible with being $p_{\rm
T}/M$-independent. When the $p_{\rm T}/M$ scaling of the cross sections —
whose evidence is further strengthened by the particularly significant
comparison between charmonium and bottomonium data at 7 and 13 TeV — is
imposed as direct constraint in the fit, the feed-down fractions,
polarizations and cross section ratios are determined with good precision. We
note that the ${\rm J}/\psi$ feed-down fractions are perfectly compatible with
values obtained in a global analysis of low-$p_{\rm T}$ data from fixed-target
experiments bib:Faccioli:FD2009 , an observation that confirms the
independence of such ratios on $p_{\rm T}/M$ and collision energy.

While the $\chi_{{c}1}$ and $\chi_{{c}2}$ individual polarizations remain the
least well known observables, the global-fit of all available data provides a
first determination of their individual values. The significant improvement in
our knowledge can be seen in Fig. 10, where the result of the global fit (in
the $1\beta$ fit variant) reported in this paper (pink contour) is compared to
the almost orthogonal results obtained with two complementary subsets of
constraints: the direct $\chi_{{c}J}$ polarization data shown in Fig. 3 (blue
line, from Ref. bib:BPH-13-001 ) and all the indirect experimental information
(red line).

Figure 10: Two-dimensional distributions showing the correlation between the
$\lambda_{\vartheta}^{\chi_{{c}2}}$ and $\lambda_{\vartheta}^{\chi_{{c}1}}$
polarizations, in the $1\beta$ fit variant, when using three sets of
measurements: only the direct $\chi_{{c}J}$ polarization data bib:BPH-13-001
(blue), which essentially constrains the
$\lambda_{\vartheta}^{\chi_{{c}2}}-\lambda_{\vartheta}^{\chi_{{c}1}}$
difference, all the other data (red), mostly constraining the
$\lambda_{\vartheta}^{\chi_{{c}1}}+\lambda_{\vartheta}^{\chi_{{c}2}}$ sum, and
the result of the global fit presented in this paper (pink).

Further improvements do not need to come from future $\chi_{{c}J}$
polarization measurements, which are notoriously challenging; precise data on
the ${\rm J}/\psi$ and $\psi{\rm(2S)}$ polarizations, as well as on the
$\chi_{{c}J}$ feed-down fractions, can also lead to better determinations of
the $\chi_{{c}J}$ polarizations. In particular, a significant improvement can
be obtained through the measurement of the _difference_ between the prompt
${\rm J}/\psi$ and $\psi{\rm(2S)}$ polarizations, potentially very precise
given the cancellation of most systematic uncertainties. New measurements of
the $\psi{\rm(2S)}$ polarization, especially towards higher $p_{\rm T}$, are
also a top priority. In fact, the polarization of _directly produced_ ${\rm
J}/\psi$ and $\psi{\rm(2S)}$ states, accessible for the first time as a result
of our global fit analysis, is found to be very small
($\lambda_{\vartheta}^{{\rm J}/\psi}=0.04\pm 0.06$) and $p_{\rm
T}$-independent. It is important to clarify if this fine-tuned balance between
transverse and longitudinal yields is only attained within the relatively
narrow $p_{\rm T}$ window covered by the presently available data, in which
case it can be seen as a mere coincidence, or remains unbroken up to higher
$p_{\rm T}$ values, in which case it can be seen as a clear sign of a highly
peculiar underlying production mechanism, probably involving a not yet
understood symmetry.

Acknowledgement: The work reported in this paper has been performed using an
extended version of the analysis framework previously used to obtain the
results reported in Ref. bib:Faccioli:EPJC78p268 , a study made in
collaboration with M. Araújo, J. Seixas, I. Krätschmer and V. Knünz.

## References

  * (1) G. T. Bodwin, E. Braaten, and P. Lepage  Phys. Rev. D 51 (1995) 1125, doi:10.1103/PhysRevD.55.5853,10.1103/PhysRevD.51.1125, arXiv:hep-ph/9407339. [Erratum: Phys. Rev.D55,5853(1997)].
  * (2) Quarkonium Working Group Coll. arXiv:hep-ph/0412158.
  * (3) N. Brambilla et al.  Eur. Phys. J. C 71 (2011) 1534, doi:10.1140/epjc/s10052-010-1534-9, arXiv:1010.5827.
  * (4) P. Faccioli, C. Lourenço, J. Seixas, and H. Wöhri  Eur. Phys. J. C 69 (2010) 657, doi:10.1140/epjc/s10052-010-1420-5, arXiv:1006.2738.
  * (5) P. Faccioli et al.  Phys. Lett. B 736 (2014) 98, doi:10.1016/j.physletb.2014.07.006, arXiv:1403.3970.
  * (6) P. Faccioli et al.  Phys. Lett. B 773 (2017) 476, doi:10.1016/j.physletb.2017.09.006, arXiv:1702.04208.
  * (7) ATLAS Coll.  Eur. Phys. J. C 76 (2016) 283, doi:10.1140/epjc/s10052-016-4050-8, arXiv:1512.03657.
  * (8) ATLAS Coll.  JHEP 09 (2014) 079, doi:10.1007/JHEP09(2014)079, arXiv:1407.5532.
  * (9) ATLAS Coll.  Phys. Rev. D 87 (2013) 052004, doi:10.1103/PhysRevD.87.052004, arXiv:1211.7255.
  * (10) CMS Coll.  Phys. Rev. Lett. 114 (2015) 191802, doi:10.1103/PhysRevLett.114.191802, arXiv:1502.04155.
  * (11) CMS Coll.  Phys. Lett. B 749 (2015) 14, doi:10.1016/j.physletb.2015.07.037, arXiv:1501.07750.
  * (12) CMS Coll.  Phys. Lett. B 780 (2018) 251, doi:10.1016/j.physletb.2018.02.033, arXiv:1710.11002.
  * (13) LHCb Coll.  Eur. Phys. J. C 71 (2011) 1645, doi:10.1140/epjc/s10052-011-1645-y, arXiv:1103.0423.
  * (14) LHCb Coll.  Eur. Phys. J. C 72 (2012) 2100, doi:10.1140/epjc/s10052-012-2100-4, arXiv:1204.1258.
  * (15) LHCb Coll.  JHEP 11 (2015) 103, doi:10.1007/JHEP11(2015)103, arXiv:1509.02372.
  * (16) LHCb Coll.  JHEP 06 (2013) 064, doi:10.1007/JHEP06(2013)064, arXiv:1304.6977.
  * (17) LHCb Coll.  JHEP 10 (2015) 172, doi:10.1007/JHEP10(2015)172, arXiv:1509.00771. [Erratum: JHEP 05, 063 (2017)].
  * (18) LHCb Coll.  Eur. Phys. J. C 80 (2020) 185, doi:10.1140/epjc/s10052-020-7638-y, arXiv:1908.03099.
  * (19) LHCb Coll.  JHEP 07 (2018) 134, doi:10.1007/JHEP07(2018)134, arXiv:1804.09214. [Erratum: JHEP 05, 076 (2019)].
  * (20) CMS Coll.  Phys. Lett. B 727 (2013) 381, doi:10.1016/j.physletb.2013.10.055, arXiv:1307.6070.
  * (21) CMS Coll.  Phys. Rev. Lett. 110 (2013) 081802, doi:10.1103/PhysRevLett.110.081802, arXiv:1209.2922.
  * (22) LHCb Coll.  Eur. Phys. J. C 73 (2013) 2631, doi:10.1140/epjc/s10052-013-2631-3, arXiv:1307.6379.
  * (23) LHCb Coll.  Eur. Phys. J. C 74 (2014) 2872, doi:10.1140/epjc/s10052-014-2872-9, arXiv:1403.1339.
  * (24) LHCb Coll.  JHEP 12 (2017) 110, doi:10.1007/JHEP12(2017)110, arXiv:1709.01301.
  * (25) ATLAS Coll.  JHEP 07 (2014) 154, doi:10.1007/JHEP07(2014)154, arXiv:1404.7035.
  * (26) CMS Coll.  Eur. Phys. J. C 72 (2012) 2251, doi:10.1140/epjc/s10052-012-2251-3, arXiv:1210.0875.
  * (27) CMS Coll.  Phys. Lett. B 743 (2015) 383, doi:10.1016/j.physletb.2015.02.048, arXiv:1409.5761.
  * (28) LHCb Coll.  Phys. Lett. B 714 (2012) 215, doi:10.1016/j.physletb.2012.06.077, arXiv:1202.1080.
  * (29) LHCb Coll.  JHEP 10 (2013) 115, doi:10.1007/JHEP10(2013)115, arXiv:1307.4285.
  * (30) CMS Coll.  Phys. Rev. Lett. 124 (2020) 162002, doi:10.1103/PhysRevLett.124.162002, arXiv:1912.07706.
  * (31) P. Faccioli and C. Lourenço  Eur. Phys. J. C 78 (2018) 731, doi:10.1140/epjc/s10052-018-6216-z, arXiv:1809.10488.
  * (32) S. Digal, P. Petreczky, and H. Satz  Phys. Rev. D 64 (2001) 094015, doi:10.1103/PhysRevD.64.094015, arXiv:hep-ph/0106017.
  * (33) F. Karsch, D. Kharzeev, and H. Satz  Phys. Lett. B 637 (2006) 75, doi:10.1016/j.physletb.2006.03.078, arXiv:hep-ph/0512239.
  * (34) Particle Data Group, M. Tanabashi et al.  Phys. Rev. D 98 (2018) 030001, doi:10.1103/PhysRevD.98.030001.
  * (35) HERA-B Coll.  Eur. Phys. J. C 49 (2007) 545, doi:10.1140/epjc/s10052-006-0139-9, arXiv:hep-ex/0607046.
  * (36) P. Faccioli, C. Lourenço, M. Araújo, and J. Seixas  Eur. Phys. J. C 78 (2018) 118, doi:10.1140/epjc/s10052-018-5610-x, arXiv:1802.01102.
  * (37) P. Faccioli, C. Lourenço, J. Seixas, and H. Wöhri  Phys. Rev. D 83 (2011) 096001, doi:10.1103/PhysRevD.83.096001, arXiv:1103.4882.
  * (38) P. Faccioli et al.  Eur. Phys. J. C 78 (2018) 268, doi:10.1140/epjc/s10052-018-5755-7, arXiv:1802.01106.
  * (39) J. E. Cavanaugh and A. A. Neath  WIREs Computational Statistics 11 (2019) e1460, doi:10.1002/wics.1460.
  * (40) C. S. Lam and W.-K. Tung  Phys. Rev. D 18 (1978) 2447, doi:10.1103/PhysRevD.18.2447.
  * (41) P. Faccioli, C. Lourenço, and J. Seixas  Phys. Rev. Lett. 105 (2010) 061601, doi:10.1103/PhysRevLett.105.061601, arXiv:1005.2601.
  * (42) NA10 Coll.  Z. Phys. C 37 (1988) 545, doi:10.1007/BF01549713.
  * (43) J. S. Conway et al.  Phys. Rev. D 39 (1989) 92, doi:10.1103/PhysRevD.39.92.
  * (44) NuSea Coll.  Phys. Rev. Lett. 99 (2007) 082301, doi:10.1103/PhysRevLett.99.082301, arXiv:hep-ex/0609005.
  * (45) NuSea Coll.  Phys. Rev. Lett. 102 (2009) 182001, doi:10.1103/PhysRevLett.102.182001, arXiv:0811.4589.
  * (46) E. Mirkes and J. Ohnemus  Phys. Rev. D 50 (1994) 5692, doi:10.1103/PhysRevD.50.5692, arXiv:hep-ph/9406381.
  * (47) CDF Coll.  Phys. Rev. Lett. 106 (2011) 241801, doi:10.1103/PhysRevLett.106.241801, arXiv:1103.5699.
  * (48) CMS Coll.  Phys. Lett. B 750 (2015) 154, doi:10.1016/j.physletb.2015.08.061, arXiv:1504.03512.
  * (49) D0 Coll.  Phys. Rev. D 63 (2001) 072001, doi:10.1103/PhysRevD.63.072001, arXiv:hep-ex/0009034.
  * (50) CDF Coll.  Phys. Rev. D 70 (2004) 032004, doi:10.1103/PhysRevD.70.032004, arXiv:hep-ex/0311050.
  * (51) CDF Coll.  Phys. Rev. D 73 (2006) 052002, doi:10.1103/PhysRevD.73.052002, arXiv:hep-ex/0504020.
  * (52) CMS Coll.  Phys. Rev. Lett. 107 (2011) 021802, doi:10.1103/PhysRevLett.107.021802, arXiv:1104.3829.
  * (53) ATLAS Coll.  Eur. Phys. J. C 72 (2012) 2001, doi:10.1140/epjc/s10052-012-2001-6, arXiv:1203.2165.
  * (54) HERA-B Coll.  Eur. Phys. J. C 60 (2009) 517, doi:10.1140/epjc/s10052-009-0957-7, arXiv:0901.1015.
  * (55) NuSea Coll.  Phys. Rev. Lett. 86 (2001) 2529, doi:10.1103/PhysRevLett.86.2529, arXiv:hep-ex/0011030.
  * (56) P. Faccioli and C. Lourenço  Eur. Phys. J. C 79 (2019) 457, doi:10.1140/epjc/s10052-019-6968-0, arXiv:1905.09553.
  * (57) P. Faccioli, C. Lourenço, J. Seixas, and H. Wöhri  JHEP 10 (2008) 004, doi:10.1088/1126-6708/2008/10/004, arXiv:0809.2153.