Quarkonium production at the LHC: A data-driven analysis of remarkably simple experimental patterns

Physics Letters B 773 (2017) 476–486

Contents lists available at ScienceDirect

Physics Letters B www.elsevier.com/locate/physletb

Quarkonium production at the LHC: A data-driven analysis of remarkably simple experimental patterns Pietro Faccioli a,b,∗ , Carlos Lourenço c,∗ , Mariana Araújo a,c , Valentin Knünz c , Ilse Krätschmer d , João Seixas a,b a

Physics Department, Instituto Superior Técnico (IST), Lisbon, Portugal Laboratório de Instrumentação e Física Experimental de Partículas (LIP), Lisbon, Portugal c European Organization for Nuclear Research (CERN), Geneva, Switzerland d Institute of High Energy Physics (HEPHY), Vienna, Austria b

a r t i c l e

i n f o

Article history: Received 29 November 2016 Received in revised form 28 July 2017 Accepted 4 September 2017 Available online 7 September 2017 Editor: L. Rolandi Keywords: Quarkonium Polarization NRQCD QCD Hadron formation

a b s t r a c t The LHC quarkonium production data reveal a startling observation: the J/ψ , ψ(2S), χc1 , χc2 and ϒ(nS) p T -differential cross sections in the central rapidity region are compatible with one universal momentum scaling pattern. Considering also the absence of strong polarizations of directly and indirectly produced S-wave mesons, we conclude 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. © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3 .

1. Introduction, motivation and conceptual remarks A large fraction of the recent publications based on LHC data are devoted to theory-driven analyses, aimed at seeing if the measurements agree with hypotheses derived in the context of specific theory models. Such analyses give the primary role to the theoretical predictions, which are based on calculations with a certain level of accuracy (often incorporating approximations) and precision (e.g., 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 theory commonly-considered in “global fits” of experimental data. While being the most sophisticated, complex, and conceptually profound theory presently available in this chapter of physics, it remains “work in progress”, with new improvements frequently taking place. The significant changes seen in some calculations from leading order (LO) to nextto-leading order (NLO) are a good example of the transient nature of the theoretical predictions.

*

Corresponding authors. E-mail addresses: [email protected] (P. Faccioli), [email protected] (C. Lourenço).

It has been previously shown [2] that global fits of quarkonium production data in the framework of NRQCD can lead to puzzling results, if one assumes that a certain superposition of the presently available NLO perturbative QCD calculations of shortdistance coefficients (SDCs) is able to describe the measured differential cross sections down to very low quarkonium transverse momentum (p T ). The highest (statistical) precision of the lowest-p T data drives the result of the fit and leads to what has been historically considered an inescapable prediction: quarkonium production must be transversely polarized, already at not-so-high p T values. The puzzling nature of the measured (absence of) quarkonium polarizations, in clear contradiction with the predictions, disappears by considering that the existing calculations are not sufficiently accurate at low p T . Indeed, Ref. [2] shows that the NLO 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 T region. This example clearly shows that fitting (a suited word) experimental measurements within a given theoretical mould might lead to results determined by the theoretical corset, that forces the data into predetermined shapes, preventing us from exploring a richer spectrum of options and, worse, potentially blinding us from simpler and more natural interpretations.

http://dx.doi.org/10.1016/j.physletb.2017.09.006 0370-2693/© 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3 .

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

477

√

Fig. 1. Mid-rapidity prompt quarkonium cross sections measured in pp collisions at s = 7 TeV by ATLAS (red markers) [3–5] and CMS (blue markers) [6,7], as a function of p T / M. The curves represent a single empirical function, with shape parameters determined by a simultaneous fit to all data (of p T / M > 2) and normalizations specific to each state (left panel) or adjusted to the J/ψ points (right panel) to directly illustrate the universality of the kinematic dependences. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The study reported in the present Letter goes one step further in our data-driven approach: perturbative calculations of the production kinematics are not used as ingredients anywhere in our analysis, the outcome of the fit being exclusively determined by the measurements, which are already sufficiently precise to provide stable and significant results; the theory calculations are only used a posteriori, in comparisons with the distributions determined by the fit. Remarkably, the (high-p T ) quarkonium measurements made at the LHC provide strong modelindependent indications regarding the mechanisms of prompt quarkonium production, allowing us to explore a broader landscape than in previous phenomenological analyses, 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 should 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 remarkably simple scenario, where all quarkonium states, despite their different quantum numbers, follow an almost universal production trend. This observation raises the central physics question of this Letter: can we derive significant statements, from the experimental observations, regarding possible differences in the production mechanisms behind 3 S1 and 3 P J quarkonium production? In Section 3 we 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 T distributions and relative normalizations are parametrized by an empirical function. The cross-section measurements are complemented by the corresponding polarizations and all the relevant feed-down decays are properly accounted for, using general kinematic relations for the modelling of the momentum distributions and polarizations of the indirectly produced states. Section 4 shows the results of the

analysis, restricted to the charmonium family given the lack of experimental information on bottomonium feed-down fractions. We then discuss, in Section 5, how the seemingly universal data patterns compare to the NRQCD theory framework, with its relatively complex structure of hierarchies and constraints, and differences in the calculated kinematic behaviours of the participating processes. 2. The universality of the measured patterns Fig. 1 shows the first of the two interesting observations we want to discuss: a seemingly universal pattern in the shapes of the p T distributions of all prompt quarkonia. Indeed, when presented as p T / M distributions, where M is the mass of the quarkonium state, the prompt J/ψ , χc1 , χc2 , ψ(2S) and ϒ(nS) production cross sections are all compatible with a single kinematic dependence, at least for mid-rapidity (| y |  1) and for not-too-small p T / M. This observation confirms, with higher-p T data, a trend first noticed in Ref. [2]. Fig. 2 offers a linear representation of the p T / M scaling, in the form of pulls, i.e. the differences between each data point and the function shown in Fig. 1-right, normalized to the measurement uncertainty. Only a small fraction of all the data points depart from the universal function by more than two standard deviations, fluctuating without any systematic trends. The good fit quality, expressed by the normalized χ 2 of 1.11 for 193 degrees of freedom (ndf), provides a quantitative measure of the p T / M universality. Fig. 3 presents a clearer view of the χc measurements, also adding χb data: the χc2 /χc1 and χb2 /χb1 yield ratios, as well as the ratio of J/ψ from χc decays to prompt J/ψ yield, show a flat dependence vs. p T / M. Since the prompt ψ(2S) mesons are fully directly produced while the J/ψ and ϒ(nS) states are significantly affected by feed-down contributions from χc and χb decays ( 25% [5] and  40% [12], respectively), the perfect compatibility of their p T / M shapes is another observation supporting the similarity of P- and S-wave quarkonium production, even in p T / M ranges uncovered by the existing χ data. The universality of the production kinematics could indicate that a single parton-level process (or mixture of processes) describes the production of all states, irrespectively of their masses and quantum numbers. Indeed, 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

478

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

Fig. 2. ATLAS (red markers) and CMS (blue markers) prompt quarkonium cross sections, represented in the form of pulls with respect to the single empirical function shown in Fig. 1-right. The horizontal bars denote the measurement bins. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

quarkonium states of different masses and angular momentum properties, at least in the mid-rapidity region. Such a “universal production” scenario is an unforeseen experimental observation: in principle, partonic production cross sections should be different for states of different quantum numbers, as will be discussed in Section 5 in the context of NRQCD. 3. Fit methodology and feed-down decays

Fig. 3. Mid-rapidity χc2 /χc1 and χb2 /χb1 prompt yield ratios, as well as the fraction of the prompt J/ψ yield coming from χc decays, measured in pp collisions by ATLAS [5] and CMS [8,9], as a function of p T / M.

invariance of the p T / M distribution, for each given rapidity range; a different scaling might apply to the forward-rapidity region of LHCb (Refs. [13–15] and references therein), as compared to the one in the mid-rapidity region covered by the ATLAS and CMS data we are using here. 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 T [16]. The most precise measurements of the λϑ polar anisotropy parameter of the J/ψ , ψ(2S) and ϒ(nS) dilepton decay distributions, shown in Fig. 4 (for the helicity frame) as a function of p T / M and averaged over rapidity, suggest a small (almost negligible) and constant polarization. While no direct χ polarization data exist, we remind that the three states shown in Fig. 4 have very different χ feed-down components, suggesting that weak polarizations are also to be expected for the contribution from the χc and χb states. We conclude that, today, there is no experimental evidence of differences in production cross sections and polarizations between

In this section we describe a more quantitative and detailed analysis of the data, performing a fit of the cross sections and polarization measurements, properly considering all the feed-down connections between the different states. Inspired by the pattern of slightly transverse polarizations seen in Fig. 4, we consider that the production cross section can be decomposed into an “unpolarized” term and a “transversely polarized” term, letting the data determine, in a model-independent way, the two contributions. This procedure could easily be generalized to scenarios including observed longitudinal polarizations by replacing the unpolarized contribution with a longitudinally polarized one. The cross section for the direct yield of a given 3 S1 state is parametrized as ∗ σdir ( p T / M ) = σdir [(1 − f p∗ ) g u ( p T / M ) + f p∗ g p ( p T / M )] ,

(1)

∗ and f ∗ are, respectively, the total direct-production where σdir p

cross section and the fractional contribution of the polarized process, both calculated at a fixed reference value, which we choose to be ( p T / M )∗ = 2. The p T / M dependences of the unpolarized and polarized yields are described by the shape functions g u ( p T / M ) and g p ( p T / M ), respectively, both normalized to unity at ( p T / M )∗ , g ( p T / M ) = h( p T / M )/h(( p T / M )∗ ), with

h( p T / M ) =

pT M

 · 1+

1

β −2

·

( p T / M )2

γ

−β .

(2)

The parameter γ (having the meaning of the average p T / M squared) defines the function in the low-p T turn-on region and is only mildly sensitive to the data we are considering here; hence, in the fit we consider γ as a common free parameter. The β power-law exponent, instead,  characterizes the high-p T shape: h ∝ ( p T / M )1−2β for p T / M  γ (β − 2). Therefore, we distinguish the unpolarized and polarized cross sections with two different powers, βu and βp , respectively.

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

479

√

Fig. 4. Polar anisotropy parameter λϑ measured by CMS in pp collisions at s = 7 TeV, for prompt J/ψ , ψ(2S) and ϒ(1S) dilepton decays [10,11]. For improved clarity, values corresponding to two or three rapidity bins were averaged, assuming uncorrelated systematic uncertainties, and the very uncertain ϒ(2S) and ϒ(3S) data are not shown.

The discrimination between the two physical contributions relies on the experimental polarization constraint, the unpolarized and polarized processes for vector quarkonia being characterized by λϑ = 0 and λϑ = 1, respectively. The polarized yield fraction can be expressed, as a function of p T / M, as f p = 3λϑ /(4 − λϑ ). As will be discussed later, 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 T dependences 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 T / M, the exclusive role of constraining both the relative normalizations and the differences between the momentum dependences 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 affected by theoretical uncertainties and will improve as the experimental panorama evolves. Our previous analysis [2] addressed the production of the ψ(2S) and ϒ(3S) quarkonia, the ones less affected by feed-down decays from P-wave quarkonia. We now concentrate on the charmonium system, complementing the ψ(2S) measurements with J/ψ , χc1 and χc2 data. The bottomonium system, with its more complex feed-down structure, not yet sufficiently constrained by LHC data (no mid-rapidity χb measurements exist), is left for future consideration. We also profit from much improved ψ(2S) and J/ψ cross section measurements recently made available by ATLAS and CMS, extending to significantly higher p T . The similarity of the p T / M dependences shown by these precise data, discussed in Section 2, justifies assuming that the direct production mechanism is the same for the J/ψ and ψ(2S) mesons. This translates, in our fit, in the use of a single “polarized fraction”, common to the two states. No significant experimental indications exist, so far, regarding χc polarization. This observable is indirectly constrained by the difference between the measured J/ψ and ψ(2S) polarizations, the former including about 25% of χc contribution while the latter is fully directly produced. Future measurements of the χc polarization and/or more precise J/ψ and ψ(2S) polarizations will allow us to perform a fully data-driven fit, without conjectures on the χc1 and χc2 polarizations. At the present moment, however, and given the precision of the existing data, it is convenient to ensure

a sufficiently constrained fit by limiting the number of free parameters, which we do by adopting plausible scenarios for the P-wave polarizations. It is only at this point, and only for the production of P-wave quarkonia, that our fit acquires an ingredient not exclusively driven by the data. As a baseline scenario, we compensate the missing χc polarization information by a calculation [17,18] based on the NRQCD framework, assuming that χc1 and χc2 production is dominated by [8] [1] the 3 S 1 octet and the 3 P 1,2 singlet terms (Section 5). Given that, in this procedure, the χc1 and χc2 polarizations are imposed, their direct cross sections are not split in two terms, contrary to what is done for the S-wave states. In short, we consider four contributions to direct quarkonium production, the unpolarized and polarized ψ terms plus the χc1 and χc2 direct cross sections, in all characterized by a common γ and four independent β parameters. In our analysis we model the charmonium feed-down decay chains ψ(2S) → χc1,c2 γ , ψ(2S) → J/ψ + X and χc1,c2 → J/ψ γ taking into account momentum and polarization transformations. Calculating the decay kinematics gives the relation connecting the mother’s and daughter’s laboratory momenta,

( M 2 − m2 )2 ( M 2 − m2 )2  p 2 /m2 = 1+ + × (1 + cos2 ) , (3) 2 2 P /M 4m2 P 2 4m2 M 2 where M (m) and P (p) are, respectively, the mass and laboratory momentum of the mother (daughter) particle, and is the emission angle of the daughter in the mother’s rest frame. The average symbols denote an integration over cos , making linear terms in this variable disappear. The term cos2 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 M 2 + m2 , 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  ( M − m), almost indistinguishable directions in the laboratory, the relation is also valid vectorially and we can assume p T /m = P T / M. In these conditions, we can model the kinematics of an indirectly produced particle using the p T / M distribution of its mother particle. Therefore, we account for the momentum

480

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

conversion along any feed-down chain by simply considering all observables as functions of p T / M. The total prompt cross sections of the states affected by feeddown are obtained as χc(1,2 χc(1,2)

( p T / M ) = σdir

σ

( pT /M )

+ B (ψ(2S) → χc(1,2) γ ) σ ψ(2S) ( p T / M ) , J/ψ σ J/ψ ( p T / M ) = σdir ( pT /M )

(4)

+ B (χc1 → J/ψ γ ) σ χc1 ( p T / M ) + B (χc2 → J/ψ γ ) σ χc2 ( p T / M ) + B (ψ(2S) → J/ψ X ) σ ψ(2S) ( p T / M ) ,

while σ = σdir in the ψ(2S) case. We now address the relevant polarization transfer rules. Measurements [19–21] of the ψ(2S) → J/ψ ππ decay angular distribution reveal that, in very good approximation, the J/ψ retains the angular momentum alignment of the mother ψ(2S), being the ππ system dominated by its J = 0 component. It can be assumed, therefore, that in the ψ(2S) to J/ψ feed-down decay the polarization is transmitted without modification, and that the shape parameters of the dilepton decay distribution (in particular, λϑ ) do not change. Concerning the J/ψ mesons produced by χc feed-down decays, it has been shown [22] that, while the angular momentum projection changes from mother to daughter, the χc → J/ψ + γ decay and the J/ψ 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 λϑ remains the same (even if obviously referring to two different decay channels) in the χc to J/ψ transition. Finally, the χc mesons produced by radiative decays of the ψ(2S) (as well as the J/ψ mesons from the subsequent χc radiative decays) have the following polarization parameters: χ

ψ(2S)

λϑ 1 = λϑ

ψ(2S)

/(4 + λϑ

),

(5)

χ2

ψ(2S) ψ(2S) /(60 + 13λϑ ). λϑ = 21λϑ

To calculate the observable polarizations of the J/ψ and χc states, combining their direct and indirect contributions (Eq. (4)), we use the general addition rule from Ref. [23], which, in the example of two processes contributing with fractions f (1) and f (2) of the cross section, reads as

λϑ(1) ⊕ λϑ(2) =

(1 ) (1 ) (2 ) (2 ) f (1) λϑ /(3 + λϑ ) + f (2) λϑ /(3 + λϑ ) (1 ) (2 ) f (1) /(3 + λϑ ) + f (2) /(3 + λϑ )

.

(6)

4. Fit results The fit has six free shape parameters: γ , βu (ψ), βp (ψ), β(χc1 ), β(χc2 ) and f p∗ (ψ). The other parameters are the normalizations of the direct production cross sections of the four charmonia. Additional nuisance parameters are used to model the correlations induced by the luminosity uncertainties. The branching ratios used to extract the cross sections from the yields measured in specific decay channels, as well as those weighting the corresponding feed-down terms of the prompt cross sections, are also treated as constrained parameters, with central values and uncertainties taken from Ref. [24]. 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 T -rapidity bin of the cross section measurement, we recalculate the average λϑ parameter for the considered state and apply the corresponding polarization-dependent acceptance correction. It is worth noting that analyses of measured cross sections and polarizations that do not follow this correction procedure are effectively using inconsistent data and are likely to obtain biased results. The fit includes the charmonium data shown in Figs. 1, 3 and 4. The cross sections are only considered for p T / M > 2. The χc -to-J/ψ feed-down fraction and χc2 -to-χc1 yield ratio of ATLAS are not considered, to avoid correlations with the χc cross section data. The results are shown in Fig. 5. The quality of the fit is exceptionally good, with a total χ 2 of 28 for 80 degrees of freedom. The best values of the shape parameters are γ = 0.73 ± 0.19, βu (ψ) = 3.42 ± 0.05, βp (ψ) = 2.67 ± 0.29, β(χc1 ) = 3.46 ± 0.08 and β(χc2 ) = 3.49 ± 0.10. The polarized fraction is quite small at low p T / M, reflecting the polarization data: f p∗ = (1 ± 1)% (at ( p T / M )∗ = 2). It increases almost linearly with p T / M, becoming (10 ± 7)% at p T / M = 15. To evaluate the sensitivity of the results to the χc1 and χc2 polarizations assumed in the fit, we have redone it considering an alternative scenario, inspired by the observation that the P-wave and S-wave p T / M-differential cross sections have seemingly identical shapes. In this option, P-wave states are produced with the same polarized and unpolarized cross section components that characterize S-wave production: g u , p (χc1,2 ) ≡ g u , p (ψ) and f p (χc1,2 ) ≡ f p (ψ). This new fit has the same quality (normalized χ 2 of 29/82) and essentially identical shape parameters: γ = 0.80 ± 0.17, βu = 3.47 ± 0.05 and βp = 2.97 ± 0.15. The similarity of the two sets of values indicates that the χc1 and χc2 polarizations do not have a strong influence on the results of the fit. In particular, it is clear that the χc1 and χc2 p T / M distributions are constrained to be very similar to those of the S-wave states. This follows the observation that the measured p T / M-differential J/ψ cross section, 25% of which arises from χc feed-down decays, is very similar to the ψ(2S) one, fully directly produced. In fact, the results remain unchanged even when the χc1 and χc2 spectra are replaced with single (average) points, thereby completely removing the experimental constraints on the χc cross section shapes, which are then freely and independently determined by the global charmonium fit. 5. Comparison with theory 5.1. The NRQCD framework Current studies of quarkonium production phenomenology are based on NRQCD, a non-relativistic effective field theory whose pillar is the hypothesis of factorization of the long- and shortdistance 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 (S ) to produce heavy quark-antiquark pairs ( Q Q ) of different colours (singlet or octet) and spin-angular momentum configurations:





σ A+B→H+X =



  S A + B → Q Q [2S +1 L CJ ] + X

S , L ,C



×L Q Q

[2S +1 L CJ ]



(7)

→H .

The coefficients L are the so-called long-distance matrix elements (LDMEs), representing the probabilities that the different Q Q states, with spin, angular-momentum and colour configuration described by the indices S, L, J and C , evolve into the

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

Fig. 5. Comparison between the data points and the best fit curves: J/ψ (top left) and ψ(2S) (top right) cross sections; λϑ for the J/ψ and ψ(2S) (bottom right).

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 T and rapidity), while the SDCs are functions of the collision type and energy, as well as of the Q 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 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 Q preresonance 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 T and rapidity regions. Assuming the limit of heavy component quarks strongly reduces and constrains, for each quarkonium, the wide range of processes entering Eq. (7). 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

481

χc cross sections and χc2 /χc1 ratio (bottom left); and

the linear combination. Moreover, some of these terms can be neglected because their SDCs lead to negligible yields in comparison to the observed prompt cross sections, a well known case being [1] the 3 S1 singlet contribution to J/ψ and ψ(2S) production [25]. As a result, the prompt production cross sections of the 3 S1 quarkonia, J/ψ , ψ(2S) and ϒ(nS), are expected to be dominated by the [8] 1 [8] 3 [8] S0 , S1 and 3 P J octet contributions, with comparable magni[8]

tudes (the 3 S1

being suppressed at the short-distance level), the

χc and χb production should mainly proceed through the 3 S[18] and [1] 3 [1] P1,2 channels, and the ηc production should reflect the 1 S0 sin[8]

[8]

[8]

glet and the 3 S1 octet channels (the 1 S0 and 1 P1 terms being suppressed [26]). 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 χc2,b2 and χc1,b1 production is expected to be equal to the corresponding ratio of spin states,

 [8 ]   [8 ]  L 3 S1 → χc2,b2 / L 3 S1 → χc1,b1 = 5/3 .

(8)

5.2. The complex SDC patterns In the present study, we are mainly interested in the characteristic kinematic dependences of the different contributing processes and on how they compare to the simplicity of the experimental scenario. Fig. 6 shows the individual contributions of

482

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

Fig. 6. SDCs (top) and corresponding helicity-frame λϑ polarization parameters (bottom), for the main octet and singlet components of directly-produced ψ mesons (left), √ [8] and for J/ψ mesons from χc1 (middle) and χc2 (right) decays, calculated at NLO [17,27,28] for pp collisions at s = 7 TeV and a charmonium mass of 3 GeV. The 3 P J SDC [8]

[8]

[8]

[8]

[1]

is defined as S (3 P0 ) + 3S (3 P1 ) + 5S (3 P2 ). The 3 P J and 3 P1,2 SDC functions are multiplied by mc2 , the squared charm-quark mass; they become negative at high p T / M and are plotted with flipped signs. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the dominating Q Q configurations, calculated at NLO [17,27,28]. The top panels show the p T / M dependence of each SDC; the P-wave SDCs are shown (here and throughout the paper) multiplied by mc2 . The bottom panels show the polarization parameter λϑ = (1 − 3 β0 )/(1 + β0 ) for S-wave dilepton decays, where β0 is the longitudinal cross section fraction (angular momentum projection J z = 0) in the helicity frame, for directly-produced J/ψ or ψ(2S) mesons (left) and for mesons resulting from χc1 (middle) or χc2 (right) feed-down decays. The P-wave contributions have rather peculiar kinematic behaviours, with cross sections becoming negative above certain p T / M values and unphysical polarization parameters (|λϑ | > 1). In fact, the calculations of the singlet and octet P-wave SDCs present singularities that match (and can be cancelled by) the one appear[8] [1,8] + 3 S[18] sum ing in the 3 S1 term [1,29]. Therefore, only the 3 P J is required to be well-defined and physical, while the two “frag[1,8] [8] ments” 3 P J and 3 S1 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, [1,8] on their own, unobservable. In particular, the 3 P J components at NLO contain large and negative transverse and/or longitudinal contributions, leading to the kinematic peculiarities mentioned above. While such characteristics can be seen as technical details of a theoretical calculation involving an expansion over purely math-

ematical objects, they are, nevertheless, in seeming contradiction with the universality of the p T / M dependences of the cross sections and polarizations measured for the different quarkonium states. The non-trivial p T / M dependences and sign-changing normalizations of the P-wave contributions require suitable cancellations to yield results that, besides being physical, must reproduce the comparatively trivial experimental trends. The different final states come from different pre-resonance [8] mixtures, of characteristic kinematic trends: the 3 S1 octet term, containing the contributions of processes where the Q Q comes from a single gluon, is fully transversely polarized (at high-enough [8] p T ), while the 1 S0 octet is unpolarized (by definition) and has a much steeper p T / M shape. At high p T / M, all three P-wave components seem to asymptotically tend to the cross section and po[8] larization shapes of their complementary-term 3 S1 . Instead, at low p T / M (including the region of the current χ data) the trends [1] [1] strongly deviate and, for example, the 3 P1 and 3 P2 terms sharply differentiate the J = 1 and 2 states, possibly leading to strong polarizations. Comparing the SDCs of Fig. 6 with the experimental trends of Fig. 3, we see, in particular, that the different turnover points of [1] [1] [8] the 3 P1 , 3 P2 and 3 P J SDCs may prevent a smooth description of the flat χ2 /χ1 and χc /J/ψ 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

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

data. In NRQCD, however, v-scaling rules predict an approximate [8] [8] [8] equality in magnitude between the 1 S0 , 3 S1 and 3 P J LDMEs for [1]

[8]

the production of 3 S1 quarkonia, and between the 3 P1,2 and 3 S1 ones for the production of χ mesons. Moreover, the HQSS relation [8] in Eq. (8) would lead to a χ2 /χ1 ratio close to 5/3 if the 3 S1 octet term dominated χ production. The measured ratio (Fig. 3) violates that expectation by a factor of 2, indicating that the singlets must be very important. On the other hand, large singlet terms would lead to a significant difference between the J = 1 and 2 p T / M dependences, probably excluded by the measured ratio. 5.3. Drawbacks of traditional data-to-theory comparisons The variety of pre-resonance contributions and their associated production and decay kinematics, as implied by v-scaling hierarchies, HQSS and current perturbative calculations, seems to be redundant with respect to the observed universal p T / M scaling and lack of polarization. In other words, the current state-of-the-art theory does not seem to naturally reflect the experimental scenario, considering its simplicity as an accidental coincidence. On the other hand, the responsibility of achieving the necessary cancellations is left to the LDMEs (Eq. (7)), which are obtained from fits to experimental data. This procedure may lead to unstable results, given that the fits are driven by the necessity of cancelling kinematic SDC dependences not visible in the data. Seemingly successful cancellations may be a transient result, allowed by the experimental uncertainties on the observables most sensitive to the critical P-wave contributions: the available χ cross sections are much more uncertain than the J/ψ and ψ(2S) cross sections, and no χ polarization measurements exist so far. The obtained LDMEs may, therefore, be fortuitous and unphysical, a situation not immediately seen as a problem, in the absence of alternative approaches. Furthermore, it is difficult to estimate the precision of the current calculations, at least in the case of the P-wave terms, which change from positive SDCs and reasonable polarizations at LO to mostly negative SDCs and unphysical polarizations at NLO [17,27, 28,30–32]. Moreover, leading-power fragmentation corrections [33, 34], accounting for a subset of next-to-next-to-leading order processes, still bring startling shape and normalization modifications [8] to the 3 P J SDC, indicating a slow convergence of the perturbative expansion. It is worth emphasising the potential fragility of the fit procedures adopted in some NRQCD analyses of quarkonium production. They rely on free LDMEs, representing unphysical quantities not calculable and not directly comparable to observations, to obtain precise cancellations between the individually-unobservable terms 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 χc and χb measurements and, to a lesser extent, in the polarization measurements, still lacking precision and completeness. Equally worrisome, it may well be that the currently available SDCs have kinematic dependences (shapes) that differ from the ultimate result by more than what is covered by the usual uncertainties, evaluated with the standard variations of the theory inputs (changes of factorization and renormalization scales, etc.). The dramatic shape changes from LO to NLO of some SDCs suggest that we might be far from the ideal calculations, in which case the present term-to-term cancellations are a temporary solution, a situation that seemingly invalidates, in practice if not conceptually, the pivotal postulate that the coefficients of the factorization expansion are universal constants.

483

5.4. A new data-to-theory comparison approach Our analysis uses a data-driven method, faithfully reflecting the current experimental inputs and ensuring that the fit outcome is stable and amenable to rigorous statistical interpretations, even if at the cost of larger uncertainties in the results. In this new fitting framework, only observable parameters and cross section contributions are considered, the perturbative kinematic calculations being replaced by empirical functions, determined by the data. The results, essentially independent of theoretical inputs, are then compared to the calculations, thereby providing an unbiased test of the theory. The shapes of the cross section terms fitted in our baseline scenario (unpolarized ψ , polarized ψ , χc1 and χc2 ) are shown in Fig. 7 as coloured bands and compared with relevant SDC combinations. The unpolarized and polarized ψ bands are compared [8] [8] with, respectively, the 1 S0 and 3 S1 SDC shapes, calculated at NLO and including fragmentation corrections [33]. The dotted curve il[8] lustrates how the addition of a small 3 P J contribution leads to steeper shapes, disfavoured by the measured polarized ψ band [8] with respect to the 3 S1 term alone. The widths of the χc1 and χc2 SDC bands reflect the uncertainty in the relative contribution [8] [1] of the 3 S1 octet and 3 P1,2 singlet terms, as evaluated [18] from the χc2 over χc1 cross section ratios measured by ATLAS and CMS. In general, we observe a good agreement between the datadriven bands and the calculated SDC shapes. The unpolarized ψ band measured in scenario 2, with a barely visible uncertainty, [8] shows an astonishing agreement with the 1 S0 octet SDC shape. In principle, theory uncertainty bands should also be shown in Fig. 7, surrounding the lines representing the SDCs (or SDC combinations). However, only the uncertainties affecting the shape should be considered here; variations of the factorization and renormalization scales, for instance, lead to visible changes in the normalization of the curves but have a small effect on their shapes. Naturally, adding such theory uncertainties would increase the overlap between the data-driven bands and the calculated SDCs. The top panels of Fig. 7 support the hypothesis that the universal unpolarized and polarized ψ cross sections have a physical [8] [8] affinity with, respectively, the 1 S0 and 3 S1 terms of the NRQCD factorization expansion. It is worth noting that the shapes of the two SDCs are significantly different from each other, their ratio changing by a factor 10 in the range 3 < p T / M < 13. 5.5. The LDME determination Making specific hypotheses about the physical nature of the unpolarized and polarized contributions allows us to extract LDME values from the comparison with the calculated SDCs. For exam[8] ple, identifying the unpolarized cross section with the 1 S0 term, the corresponding LDME for a given state is

 [8 ]     [8 ]  L 1 S0 = 1 − f p σdir / S 1 S0 .

(9)

Using the SDCs calculated at NLO, complemented with fragmentation contributions, for the production of a Q Q pair with 3 GeV of rest energy, the fit results for the direct J/ψ cross section in [8] the baseline scenario lead to the 1 S0 LDME represented by the grey band in Fig. 8. The width of the band only represents the experimental uncertainty; it does not reflect possible changes of the calculations due to scale dependences (mostly a normalization shift) or to potential higher-order improvements. As evident in the figure and implied by Eq. (9), where f p , σdir [8]

and S (1 S0 ) are functions of p T / M, the (supposedly constant) co[8]

efficient L(1 S0 ) becomes, technically, a function of p T / M. This illustrates how the universality of the LDMEs, and hence the validity

484

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

Fig. 7. Direct production cross sections resulting from the fit of the data, in the two considered χc polarization scenarios, with uncertainty bands reflecting correlated variations in the fit parameters. Suitable NLO SDC combinations are also shown, described in the text. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

L(1 S[08] ) results obtained in three different analyses of charmonium

[8]

Fig. 8. 1 S0 LDME for direct J/ψ production in the baseline scenario (grey band), compared to results obtained in other analyses: BK [35], CMSWZ [27] and BCKL [33].

of the factorization hypothesis itself, are in reality only approximations relying on the goodness of the perturbative SDC calculations at a given order. The extracted LDME increases by a factor of 2 in the low-p T / M region, presumably more affected by calculation uncertainties, and then stabilizes. Also for this reason, the LDME, when expressed by a constant number, necessarily depends on the method used for theory-data comparison. Fig. 8 also shows the

production, CMSWZ [27], BCKL [33] and BK [35], represented as uncertainty bands illustrating the p T / M ranges of the used data. Since these results were extracted from prompt-J/ψ cross sections, neglecting the ψ(2S) and χc feed-down, we rescaled them by 2/3, the approximate ratio of direct-over-prompt J/ψ yields (a practically p T / M-independent ratio, given the almost universal p T / M scaling of the cross sections). Remarkably, the three results follow a trend seemingly determined by the lower p T / M limit chosen for the input data, which also reflects the highest event population and most constraining data points. That trend is well described by our curve, implicitly showing that all results are conceptually compatible, while clearly revealing that the specific numerical results are somewhat arbitrary and mostly determined by the choice of the starting p T / M value (or, in our case, by the normalization point in Eq. (9)). Our analysis method keeps data and theory disentangled, by using the data to obtain the best possible experimental determinations of global observables, which are consecutively compared to theory. Such comparisons can, thus, be done as a function of kinematics, which is not possible when the theory calculations are intrinsically embedded in the fits. 5.6. Comments on an alternative polarization hypothesis We conclude this section by discussing a new hypothesis on long-distance polarization effects in quarkonium production and how it relates to the scenarios we have been considering. A cru-

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

[8]

cial ingredient of NRQCD is the “standard” conjecture that the 3 S1 Q Q transfers its polarization (fully transverse at not-too-low p T ) to the observable J/ψ , ψ(2S) or ϒ(nS) without attenuation, if no intermediate feed-down transitions are involved. The process is modelled as a transition of the Q Q pair with the simultaneous emission of two gluons having total angular momentum J = 0, in analogy with the observable process ψ(2S) → J/ψ ππ , where the J/ψ inherits the ψ(2S) polarization. It has been recently observed [36] that another kind of process may be at play, consisting of two successive electric dipole transi[8] [8] [8] tions: 3 S1 → 3 P J + g, followed, for example, by 3 P J → J/ψ + g. [8]

In this case, the resulting J/ψ and ϒ polarization from 3 S1 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 ψ(2S) → χc + g and χc → J/ψ + g. In the first step the polarization is attenuated, according to the relations in Eq. (5), while the second one preserves the value of λϑ . [8] According to the calculations of Ref. [36], the 3 S1 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 unknown proportions, bringing the NRQCD predictions closer to an unpolarized scenario of J/ψ , ψ(2S) and ϒ(nS) production. The new hypothesis does not affect χc po[8] larizations, since the transition of 3 S1 to χc is already modelled according to the minimal assumption of a single electric dipole transition. [8] A further possible justification of an unpolarized 3 S1 contribution has been provided in the kT -factorization approach [37], which takes into account the off-shellness of the initial gluons. Our general decomposition of J/ψ and ψ(2S) cross sections into “unpolarized” and (fully transversely) “polarized” contributions can [8] fit both scenarios, of transverse or unpolarized 3 S1 , since neither involves possible longitudinal contributions. In the new conjecture [8] the 3 S1 term contributes to both the unpolarized (σu ) and polarized (σp ) cross sections, depending on the effective polarization [8]

[8]

(0 < λϑ (3 S1 ) < 1) transmitted by the 3 S1 octet to the quarkonium state. The comparison in Fig. 7, for example, should be made with a different association of the unpolarized and polarized terms to the two octet channels,

1

[8 ] 

σu → σ S 0 σp →

3 λϑ

3

4 − λϑ



+

4 − 4 λϑ



3 S[ 8 ] 1

[8] 



[8 ] 

S1

[8 ] 4 − λ ϑ 3 S1

[8 ] 

S1

3

3

[8 ] 

 σ S1 

3 [ 8 ]   σ S1 , (10)

, [8] 



[8] 

where σ 2S+1 S J = S 2S+1 S J × L 2S+1 S J . The second equation shows that the polarized term should now be compared to [8] a rescaled 3 S1 function, leaving the pure shape comparison unaltered. The comparison term for the unpolarized function, how[8] [8] ever, becomes a combination of the 1 S0 shape with a 3 S1 contribution, depending on the assumed

λϑ (3 S[18] ).

3 [8] 1 [8] S1 / S0

yield ratio and

Independently of these two assumptions, such a combination would be a less steep function of p T / M. It is interesting to notice that, given the very good agreement between the un[8] polarized shape and the 1 S0 term alone (shown in Fig. 7), the present measurements already disfavour a non-negligible contribution from an additional unpolarized contribution having the shape [8] of the 3 S1 SDC. In other words, the data do not indicate that the [8]

[8]

3 S1 polarization, always considered to be λϑ (3 S1 ) ≈ +1, should instead be closer to 0.

485

6. Discussion and conclusions In this Letter we described a data-driven approach to the study of quarkonium production, replacing perturbative calculations of differential cross sections by empirical functions determined by the data. The responsibility of discriminating between the different elementary processes behind quarkonium production is assigned to the polarizations, which are mainly determined by basic physics considerations and, hence, are less dependent on fixed-order QCD calculations. This new method reflects a change of perspective with respect to traditional analyses, where the process discrimination is exclusively based on the p T dependences of the calculated individual cross-section contributions (the SDCs). 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 of multidimensional decay-angle analyses with respect to the much simpler yield counting. The fit results only reflect the data trends and uncertainties, independently of theoretical ingredients, such as the existence of unobservable or even unphysical partial contributions, or the factorization into short- and long-distance factors. Comparisons with SDC calculations and LDME evaluations are only made a posteriori. This is a safer procedure for data-theory comparisons, given that only the data values are known with uncertainties that have a rigorous statistical interpretation. At this moment, given the lack of χc polarization data, we must resort to some model assumptions regarding χc production. Once χc polarization measurements will become available, the procedure will no longer need any model input. The disentangling of theory calculations from the fit to data allows us to more easily face problems specific to one of the two heterogeneous worlds. In particular, tensions originating from mutual experimental inconsistencies (such as those seen in Ref. [38]) can be addressed before and independently of the comparison to theory, without affecting the judgement of how well any given model describes the data. And since the fit results only represent the data, with their uncertainties, they can be confronted simultaneously to different models or multiple variants of a model, the best way to appreciate the diversity and uncertainty of the theory world. Instead, analyses that embed model calculations in the fit usually neglect their uncertainties, or assign somewhat arbitrarily values, so that it is difficult to evaluate how well a model describes the data, given that the resulting χ 2 does not reliably express the goodness of the fit. In such cases, interesting indications can only come from dramatic divergences, 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 blurred the crucial evidence that the perturbative SDC calculations are steeper, at low p T , than the experimental distributions. This feature, studied in detail in Ref. [2] by performing fits as a function of a selection cut progressively removing low-p T data, is immediately visible when our data-driven results are compared to calculations, as done in Fig. 7. Also the LDME determination shown in Fig. 8 illustrates how the new approach brings to light the limitations of the current calculations. It is rewarding to see that the present experimental knowledge on quarkonium production is significantly more satisfactory than in the pre-LHC years, when mutual inconsistencies shed doubts on the reliability of some data. With forthcoming LHC measurements filling the few remaining gaps and further improving the constraining power of the experimental picture, quarkonium pro-

486

P. Faccioli et al. / Physics Letters B 773 (2017) 476–486

duction could soon become an area of precision QCD studies, a prospect that calls for the use of rigorous theory-data comparison methods, also accounting for the correlations between the cross section and polarization measurements, common uncertainties among data sets, feed-down decays, etc. A remarkable observation can already be made today: there is no evidence at all, even in the latest ATLAS and CMS cross-section measurements, which extend up to p T around 100 GeV, of differences among the short-distance mechanisms behind the production of the several quarkonium states. Indeed, all the seven measured states, from the J/ψ to the ϒ(3S), exhibit identical p T / M cross-section shapes, and all the states for which the polarizations have been measured (the five S-wave quarkonia) show perfectly compatible trends (almost unpolarized production). While more data are needed for the χ states (in particular, first polarization measurements), the comparison of the kinematic distributions of the several S-wave particles, which receive very different contributions from χ decays, strongly suggests that the P-wave and S-wave quarkonia follow the same p T / M scaling. These simple patterns in the data should stimulate theoretical considerations towards a natural explanation. Acknowledgements We are indebted to H.-S. Shao, who kindly gave us the NLO calculations of the short distance coefficients. The work of I.K. is supported by FWF, Austria, through the grant P 28411-N36. References [1] G.T. Bodwin, E. Braaten, P. Lepage, Phys. Rev. D 51 (1995) 1125, http:// dx.doi.org/10.1103/PhysRevD.51.1125, arXiv:hep-ph/9407339; Erratum: Phys. Rev. D 55 (1997) 5853, http://dx.doi.org/10.1103/PhysRevD.55.5853. [2] P. Faccioli, et al., Phys. Lett. B 736 (2014) 98, http://dx.doi.org/10.1016/ j.physletb.2014.07.006, arXiv:1403.3970. [3] ATLAS Collaboration, J. High Energy Phys. 09 (2014) 079, http://dx.doi.org/ 10.1007/JHEP09(2014)079, arXiv:1407.5532. [4] ATLAS Collaboration, Phys. Rev. D 87 (2013) 052004, http://dx.doi.org/10.1103/ PhysRevD.87.052004, arXiv:1211.7255. [5] ATLAS Collaboration, J. High Energy Phys. 07 (2014) 154, http://dx.doi.org/ 10.1007/JHEP07(2014)154, arXiv:1404.7035. [6] CMS Collaboration, Phys. Rev. Lett. 114 (2015) 191802, http://dx.doi.org/ 10.1103/PhysRevLett.114.191802, arXiv:1502.04155. [7] CMS Collaboration, Phys. Lett. B 749 (2015) 14, http://dx.doi.org/10.1016/ j.physletb.2015.07.037, arXiv:1501.07750. [8] CMS Collaboration, Eur. Phys. J. C 72 (2012) 2251, http://dx.doi.org/10.1140/ epjc/s10052-012-2251-3, arXiv:1210.0875. [9] CMS Collaboration, Phys. Lett. B 743 (2015) 383, http://dx.doi.org/10.1016/ j.physletb.2015.02.048, arXiv:1409.5761. [10] CMS Collaboration, Phys. Lett. B 727 (2013) 381, http://dx.doi.org/10.1016/ j.physletb.2013.10.055, arXiv:1307.6070.

[11] CMS Collaboration, Phys. Rev. Lett. 110 (2013) 081802, http://dx.doi.org/ 10.1103/PhysRevLett.110.081802, arXiv:1209.2922. [12] LHCb Collaboration, Eur. Phys. J. C 74 (2014) 3092, http://dx.doi.org/10.1140/ epjc/s10052-014-3092-z, arXiv:1407.7734. [13] LHCb Collaboration, Phys. Lett. B 714 (2012) 215, http://dx.doi.org/10.1016/ j.physletb.2012.06.077, arXiv:1202.1080. [14] LHCb Collaboration, Phys. Lett. B 718 (2012) 431, http://dx.doi.org/10.1016/ j.physletb.2012.10.068, arXiv:1204.1462. [15] LHCb Collaboration, Eur. Phys. J. C 74 (2014) 2872, http://dx.doi.org/10.1140/ epjc/s10052-014-2872-9, arXiv:1403.1339. [16] P. Faccioli, C. Lourenço, J. Seixas, H. Wöhri, Eur. Phys. J. C 69 (2010) 657, http://dx.doi.org/10.1140/epjc/s10052-010-1420-5, arXiv:1006.2738. [17] H.-S. Shao, Y.-Q. Ma, K. Wang, K.-T. Chao, Phys. Rev. Lett. 112 (2014) 182003, http://dx.doi.org/10.1103/PhysRevLett.112.182003, arXiv:1402.2913. [18] P. Faccioli, C. Lourenço, M. Araújo, J. Seixas, I. Krätschmer, V. Knünz, χc1 and χc2 polarizations: a distinctive experimental probe of NRQCD, submitted to Phys. Rev. Lett. (2017). [19] BES Collaboration, Phys. Rev. D 62 (2000) 032002, http://dx.doi.org/10.1103/ PhysRevD.62.032002, arXiv:hep-ex/9909038. [20] CLEO Collaboration, Phys. Rev. D 30 (1984) 1433, http://dx.doi.org/10.1103/ PhysRevD.30.1433. [21] CLEO Collaboration, Phys. Rev. D 58 (1998) 052004, http://dx.doi.org/10.1103/ PhysRevD.58.052004, arXiv:hep-ex/9802024. [22] P. Faccioli, C. Lourenço, J. Seixas, H. Wöhri, Phys. Rev. D 83 (2011) 096001, http://dx.doi.org/10.1103/PhysRevD.83.096001, arXiv:1103.4882. [23] P. Faccioli, C. Lourenço, J. Seixas, Phys. Rev. D 81 (2010) 111502, http:// dx.doi.org/10.1103/PhysRevD.81.111502, arXiv:1005.2855. [24] K.A. Olive, Chin. Phys. C 40 (2016) 100001, http://dx.doi.org/10.1088/ 1674-1137/40/10/100001. [25] CDF Collaboration, Phys. Rev. Lett. 79 (1997) 572, http://dx.doi.org/10.1103/ PhysRevLett.79.572. [26] M. Butenschön, Z.-G. He, B. Kniehl, Phys. Rev. Lett. 114 (2015) 092004, http://dx.doi.org/10.1103/PhysRevLett.114.092004, arXiv:1411.5287. [27] K.-T. Chao, et al., Phys. Rev. Lett. 108 (2012) 242004, http://dx.doi.org/10.1103/ PhysRevLett.108.242004, arXiv:1201.2675. [28] H.-S. Shao, Comput. Phys. Commun. 198 (2016) 238, http://dx.doi.org/10.1016/ j.cpc.2015.09.011, arXiv:1507.03435. [29] A. Petrelli, et al., Nucl. Phys. B 514 (1998) 245, http://dx.doi.org/10.1016/ S0550-3213(97)00801-8, arXiv:hep-ph/9707223. [30] M. Butenschön, B. Kniehl, Nucl. Phys. B, Proc. Suppl. 222 (2012) 151, http://dx.doi.org/10.1016/j.nuclphysbps.2012.03.016, arXiv:1201.3862. [31] M. Butenschön, B. Kniehl, Mod. Phys. Lett. A 28 (2013) 1350027, http:// dx.doi.org/10.1142/S0217732313500272, arXiv:1212.2037. [32] M. Butenschön, B. Kniehl, Phys. Rev. Lett. 108 (2012) 172002, http:// dx.doi.org/10.1103/PhysRevLett.108.172002, arXiv:1201.1872. [33] G.T. Bodwin, H.S. Chung, U.-R. Kim, J. Lee, Phys. Rev. Lett. 113 (2014) 022001, http://dx.doi.org/10.1103/PhysRevLett.113.022001, arXiv:1403.3612. [34] G.T. Bodwin, et al., Phys. Rev. D 93 (2016) 034041, http://dx.doi.org/10.1103/ PhysRevD.93.034041, arXiv:1509.07904. [35] M. Butenschön, B.A. Kniehl, Phys. Rev. D 84 (2011) 051501, http://dx.doi.org/ 10.1103/PhysRevD.84.051501, arXiv:1105.0820. [36] S. Baranov, Phys. Rev. D 93 (2016) 054037, http://dx.doi.org/10.1103/ PhysRevD.93.054037. [37] S. Baranov, A. Lipatov, N. Zotov, Eur. Phys. J. C 75 (2015) 455, http://dx.doi.org/ 10.1140/epjc/s10052-015-3689-x, arXiv:1508.05480. [38] ATLAS Collaboration, Eur. Phys. J. C 76 (2016) 283, http://dx.doi.org/10.1140/ epjc/s10052-016-4050-8, arXiv:1512.03657.