From identical S- and P-wave pT/M spectra to maximally distinct polarizations: probing NRQCD with chi states

A global analysis of ATLAS and CMS measurements reveals that, at mid-rapidity, the directly-produced $\chi_{c1}$, $\chi_{c2}$ and J/$\psi$ mesons have differential cross sections of seemingly identical shapes, when presented as a function of the mass-rescaled transverse momentum, $p_{\rm T}/M$. This identity of kinematic behaviours among S- and P-wave quarkonia is certainly not a natural expectation of non-relativistic QCD (NRQCD), where each quarkonium state is supposed to reflect a specific family of elementary production processes, of significantly different $p_{\rm T}$-differential cross sections. Remarkably, accurate kinematic cancellations among the variegated NRQCD terms (colour singlets and octets) of its factorization expansion can lead to a surprisingly good description of the data. This peculiar tuning of the NRQCD mixtures leads to a clear prediction regarding the $\chi_{c1}$ and $\chi_{c2}$ polarizations, the only observables not yet measured: they should be almost maximally different from one another, and from the J/$\psi$ polarization, a striking exception in the global panorama of quarkonium production. Measurements of the difference between the $\chi_{c1}$, $\chi_{c2}$ and J/$\psi$ polarizations, complementing the observed identity of momentum dependences, represent a decisive probe of NRQCD.


Introduction
The mechanisms behind hadron production continue to challenge our understanding: analytical perturbative QCD calculations are insufficient to tackle all the aspects of the strong interactions driving the binding of quarks into observable particles. Studies of quarkonium production can provide crucial progress towards solving this problem [1]. According to non-relativistic QCD (NRQCD) [2], one of the theory approaches in this area of QCD phenomenology, S-and P-wave quarkonia are produced from the binding of quark-antiquark pairs created with a variety of quantum numbers, in color singlet or octet configurations. These terms are characterized by significantly different kinematic dependences and polarizations, determined by the short-distance cross sections (SDCs), presently calculated at next-to-leading order (NLO) [3][4][5][6]. They contribute with probabilities proportional to long distance matrix elements (LDMEs), extracted from fits to experimental data. While conceptually appealing and successful in several respects, it has been confusing to see that different groups performing global fits to experimental data extract significantly different matrix elements, despite using identical theory calculations, as a result of using different data fitting strategies [3][4][5][6]. These puzzles and a potential solution were discussed in Ref. [7], mostly devoted to the quarkonia least affected by feed-down contributions, the ψ(2S) and Υ (3S) states. A detailed data-driven analysis of the cross sections and polarizations of five S-wave and two P-wave states, complemented by an original comparison with theory calculations, was presented in Ref. [8]. That analysis is extended in this paper to address two main questions: how different and experimentally recognizable are the χ c production mechanisms with respect to those of the J/ψ and ψ(2S) mesons; and how important will be, for the understanding of quarkonium production, new or improved χ c measurements.

Data-driven considerations
As shown in the top panel of Fig. 1, the 3 S 1 and 3 P J charmonium and bottomonium cross sections measured at the LHC at mid-rapidity show a remarkably uniform pattern as a function of p T /M , the ratio between the quarkonium transverse momentum and its mass (   [9][10][11] and CMS (blue markers) [12,13]. The normalizations were adjusted to the J/ψ points to directly illustrate the universality of the pT/M dependence. The curve represents a fit to all points of pT/M > 2 [8], with a normalized χ 2 of 215/193 = 1.11. The inset shows the NLO SDCs [5,14,15]. The 3 P 1,2 SDCs are multiplied by m 2 c , the mass of the charm quark squared; they are negative and plotted with flipped signs. Bottom: Polar anisotropy parameter λ ϑ , in the helicity frame, measured by CMS in pp collisions at √ s = 7 TeV, for prompt J/ψ, ψ(2S) and Υ (1S) dimuon decays [16,17]. For improved visibility, values corresponding to two or three rapidity bins were averaged. The curves represent the calculated λ ϑ = (ST −SL)/(ST +SL) values, where ST (SL) is the transverse (longitudinal) short distance cross section, in the helicity frame (HX).
Ref. [8] shows the seven independent distributions). Moreover, the corresponding measurements of the quarkonium decay distributions indicate similar polarizations for all Swave states, independently of their different P-wave feeddown contributions, as expressed by their polar anisotropy parameters, in the helicity frame, shown in the bottom panel.
This seemingly "universal" picture of quarkonium production is an unexpected result, when compared to the wide variety of kinematic shapes of the differential cross sections (NLO SDCs) contributing to the observable pat-terns within the NRQCD framework, as shown in the inset of Fig. 1. The most surprising aspect is that the χ c1 and χ c2 P-wave states have, at least at mid-rapidity, p T /M distribution shapes indistinguishable from those of the S-wave states. According to the SDCs calculated at NLO, on the other hand, the singlet and octet P-wave terms of the NRQCD expansion, which contribute differently to J/ψ (ψ(2S)), χ c1 and χ c2 production, have rather peculiar and differentiated kinematic behaviours, with cross section terms becoming negative above characteristic p T /M thresholds and having unphysical polarization parameters (|λ ϑ | > 1). In advance of any detailed numerical analysis, the qualitative comparison between data and theory illustrated in Fig. 1 indicates that the theory requires precise and seemingly unnatural cancellations between terms of the expansion, in order to reproduce observable cross sections and polarizations that are not only physical but also identical for states of different quantum numbers.
It should be noted that comparing the shapes of seven different quarkonium states, including five S-wave states affected by very different fractions of P-wave feed-down contributions, provides a stronger (more precise) statement regarding the overall equality between S-and P-wave quarkonium production than one might initially expect, given the uncertainties of the χ c1 and χ c2 measurements on their own. We will quantify this observation when presenting the results of the charmonium fit.
The χ c1 and χ c2 polarizations are the main missing element in the current experimental landscape, but two data-driven observations provide indirect indications. First, the ψ(2S), J/ψ and Υ (1S) polarizations are very similar ( Fig. 1-bottom), despite the diversity of χ feeddown fractions (0, ∼ 25% [11,18] and ∼ 40% [19], respectively). Assuming that the directly-produced S-wave mesons have very similar production mechanisms, as indicated by the seemingly identical shapes of the p T /Mdifferential cross sections ( Fig. 1-top), the χ c1 plus χ c2 summed feed-down contributions cannot have a large impact in the observed J/ψ polarization. The second observation derives from comparing χ c2 /χ c1 cross-section ratios measured in different experimental acceptances, profiting from their strong sensitivity to the polarization hypothesis used in the acceptance corrections. Interestingly, as seen in Fig. 2, the J z = 0 alignment hypothesis gives the best mutual agreement between the χ c2 /χ c1 ratios reported by ATLAS and CMS, as well as between the LHCb values obtained using photons detected in the calorimeter or with conversions to e + e − pairs. When both are polarized in the J z = 0 limit, the χ c1 and χ c2 decays produce strongly polarized J/ψ mesons, with, respectively, λ ϑ = +1 and λ ϑ = −3/5, leading to a weighted λ ϑ ∼ 0.3 when the feed-down fractions and the cross-section ratio itself are taken into account. These observations suggest that the χ c1 and χ c2 polarizations might be a striking exception in the global panorama of high-energy quarkonium production, at least at mid rapidity. The χc2/χc1 ratio measured in pp collisions at 7 TeV by ATLAS [11], CMS [20] and LHCb [21,22], with acceptance corrections calculated with two extreme polarization hypotheses: spin alignments Jz(χc1) = ±1, Jz(χc2) = ±2 (top) and Jz(χc1) = Jz(χc2) = 0 (bottom). The unpolarized hypothesis leads to intermediate values.

Analysis method
To quantify our previous data-driven considerations and compare the results with theory, we perform a simultaneous fit of the mid-rapidity differential cross sections and polarizations, including a detailed account of how the mother's momentum and polarization are transferred to the daughter in the relevant feed-down decays: ψ(2S) → χ c1,2 γ; ψ(2S) → J/ψ X; χ c1,2 → J/ψ γ. The analysis is restricted to the charmonium family, given the lack of experimental information on bottomonium feed-down fractions. The rule for the momentum propagation from mother to daughter is, approximately, p T /m = P T /M , where M (m) and P T (p T ) are, respectively, the mass and laboratory transverse momentum of the mother (daughter) particle [8]. 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 [23]. The fit is exclusively based on empirical parametrizations. 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 and, therefore, only affected by statistical and systematic experimental uncertainties.
Inspired by the pattern of slightly transverse polarizations seen in Fig. 1, we parametrize the directly-produced J/ψ and ψ(2S) cross section shapes as a superposition of unpolarized (λ ϑ = 0) and transversely polarized (λ ϑ = +1) processes, λ ϑ being the polar anisotropy parameter of the dilepton decay in the helicity frame [24]: where f p , identical for the two charmonia, is the fractional contribution of the polarized process considered at an arbitrary reference point (p T /M ) * . The shape functions g u (p T /M ) and g p (p T /M ) describe the p T /M dependences of, respectively, the unpolarized and polarized yields. Both are normalized to unity at the The parameter γ (having the meaning of the average p T /M squared) defines the function in the low-p T turnon 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: . Therefore, we distinguish the unpolarized and polarized cross sections with two different powers, β u and β p , respectively, identical for the J/ψ and ψ(2S). The relative contributions and shapes of the g u and g p functions are constrained by the polarization data. In fact, the polarized yield fraction, equal to f p at (p T /M ) * , can be expressed as a function of p T /M as For the χ c1 and χ c2 direct cross sections we use the same general p T /M shape parametrization, but without discriminating between polarized and unpolarized contributions, which, in the absence of χ polarization data, would not be individually constrained by the fit. In short, we consider four contributions to direct quarkonium production, the unpolarized and polarized ψ terms plus the χ c1 and χ c2 cross sections, altogether characterized by one γ and four β parameters, β u , β p , β(χ 1 ) and β(χ 2 ). Their theoretical counterparts are, respectively, 1 S 2 (where each term indicates the SDC function times the LDME constant), the four leading cross section components foreseen by NRQCD hierarchies for 3 S 1 and 3 P J quarkonium production. However, this parallelism is only a guidance in the parametrization of the fit, not a theoretical input. As discussed in more detail in Ref. [8], our approach is very different with respect to fits using the calculated SDC shapes, where the fit results are mostly determined by the p T -differential cross sections; the less precise polarization data are not included in the fits or have a negligible effect. In our fit, the polarization data, versus p T /M , have the exclusive role of constraining both the relative normalizations and the differences in momentum dependence of the polarized and unpolarized contributions. The precision of these data-driven results will evolve as new measurements become available, remaining insensitive to specific theoretical calculations and uncertainties.
Without χ polarization measurements, the χ c1 and χ c2 cross sections cannot help, today, discriminating the g u and g p contributions to J/ψ production and, therefore, relate the parametrized direct-J/ψ polarization to the mea-sured prompt one. However, the two data-driven observations mentioned above allow us to implement such a relation by adopting an approximate constraint on the total χ c polarization contribution to J/ψ production. Given that, on average, λ J/ψ ϑ λ ψ(2S) ϑ , we can infer that λ J/ψ←χc ϑ should be positive, under the assumption that the direct J/ψ and ψ(2S) polarizations are equal. On the other hand, the extreme hypothesis, discussed above, according to which both χ c1 and χ c2 are polarized in the J z = 0 limit leads to λ J/ψ←χc ϑ 0.3 (which, when weighted by the 25% feed-down fraction of J/ψ from χ c , is comparable to the average difference λ The ATLAS and CMS integrated-luminosity uncertainties are (independently) varied as nuisance parameters, following Gaussian functions centred at unity and of widths equal to the relative uncertainties of the published luminosities. These two nuisance parameters multiply all the data points (cross sections) of the respective experiment, thereby correlating the several datasets within each experiment. Moreover, the experiments measured products of cross sections times branching ratios, so that the uncertainties of the branching ratios have also been treated as nuisance parameters (with central values and uncertainties taken from Ref. [25]), multiplying all relevant data points and representing correlations between ATLAS and CMS.
Another source of correlation between all the points being fitted 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 T , using the shape-parametrization functions described above. The expected λ ϑ values can be immediately compared to the measured ones, for the determination of the corresponding χ 2 terms, while for the calculation of the cross-section χ 2 terms we first scale the measured cross sections by acceptance-correction factors calculated for the λ ϑ 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 polarization.
The fit has 100 experimental constraints and 20 parameters: 5 shape parameters, 4 normalizations and the fraction f p , plus 2 luminosity and 8 branching-ratio nuisance parameters.

Analysis results
As shown in Fig. 3, the charmonium cross sections and polarizations are described by the fit just presented, with a χ 2 per degree of freedom of 28/80. Figure 4 shows the fitted cross section terms as bands of widths reflecting the experimental uncertainties. A very interesting and non-trivial indication of this purely data-driven fit is that the χ c1 and χ c2 p T /M distributions are very similar to the unpolarized term dominating ψ production, as quantified by the compatibility of the β parameters: β u = 3.42 ± 0.05, β(χ 1 ) = 3.46 ± 0.08 and β(χ 2 ) = 3.49 ± 0.10. This very clear experimental observation is predominantly the result of the perfect compatibility of the (high precision) J/ψ and ψ(2S) p T /M shapes, even in p T /M ranges beyond those covered by the existing χ data, reflecting the fact that the prompt ψ(2S) mesons are fully directly produced while 25% [11] of the J/ψ yield comes from χ c decays. In fact, the χ c cross sections are measured at relatively low p T and with comparatively poor precision. To verify this conclusion, we repeated the fit keeping only one experimental point for each of the two χ c cross sections (chosen in the middle of the measured range), so that these measurements constrain the feed-down fractions at that point but not the p T /M shapes. As expected, the fit results for the χ c cross section shapes do not change significantly, with shape parameters remaining the same within the one-sigma range. The experimental bands for the four observable cross sections are compared with the corresponding NRQCD terms, the dashed/dotted lines corresponding to (combinations of) SDCs calculated at NLO. We emphasize that the two terms of comparison are completely independent, the first being the result of a model-independent fit of experimental data and the second a pure theoretical calculation. The unpolarized and polarized ψ bands are compared with, respectively, the 1 S [8] 0 and 3 S [8] 1 SDC shapes, calculated at NLO and also including fragmentation corrections representing a partial account of next-to-nextto-leading order processes [26,27]. As illustrated by the dot-dashed line, corresponding to ( 1 ) , adding a 3 P [8] J term leads to steeper shapes, departing from the polarized ψ band more than the 3 S [8] 1 term alone. The 1 S [8] 0 SDC shape is in remarkable agreement with the experimental "unpolarized" band. Moreover, adding the negative 3 P term results in shapes approximating the 1 S [8] 0 term, reproducing relatively well the observed similarity between the unpolarized-ψ and χ c1,2 patterns.
Before discussing in more detail this data-theory comparison, we will now describe the derivation of the NRQCD curves for the χ c distributions and corresponding polarization predictions. In NRQCD the χ c1,2 polarizations and cross sections are functions of one common parameter, equal for all χ c states, Top: χc2/χc1 ratio measured by ATLAS [11] and CMS [20], before (open markers) and after (filled markers) accounting for the dependence of the detection acceptances on the (simultaneously calculated) χc1 and χc2 polarizations. The grey bands reflect the theory fits to the data, with widths reflecting the Kχ uncertainty. Bottom: χc1,2 polarizations calculated adding the 3 S with O denoting the LDME. The χ c production cross sections σ J and the spin-density matrix elements σ ij J have the general form where S (ij) denotes the SDC or its spin projection. The λ ϑ are calculated as λ χ1 ϑ = (σ 00 1 − σ 11 1 )/(σ 00 1 + 3σ 11 1 ) and λ χ2 ϑ = (−3σ 00 2 − 3σ 11 2 + 6σ 22 2 )/(5σ 00 2 + 9σ 11 2 + 6σ 22 2 ), where the σ ij J depend on K χ through Eq. 3. The χ c1,2 λ ϑ parameters refer to the corresponding J/ψ dilepton decay distributions, which are the ones directly measured and fully reflect the χ polarization state, while being insensitive to the uncertain contributions of higher-order photon multipoles [23].
We determine K χ from the measured χ c2 /χ c1 ratio, taking into account that the published values strongly depend on the χ c1 and χ c2 polarizations assumed for the corrections of the detector's acceptance. For each K χ considered, we calculate the χ c1 and χ c2 polarizations us-ing NLO SDCs and correct the published ratio by the corresponding acceptance ratio. The fit χ 2 is then calculated comparing the corrected measurement (with statistical and systematic uncertainties, but no "polarization uncertainties") with the prediction for that K χ value. Figure 5-top shows how the χ c2 /χ c1 ratio changes, and the theoretical fit improves, when we use the NRQCD polarization conjecture instead of the unpolarized scenario that the experiments assume to report the measurements. The result of our fit is K χ = 4.60 ± 0.06, much more precise than the value 3.7 +1.0 −0.7 , derived [14] using the unpolarized ratios and including the entire spectrum of polarization hypotheses in the experimental uncertainty. The corresponding polarization predictions are shown in Fig. 5bottom. Interestingly, as p T /M decreases, λ ϑ tends to the extreme physical values +1 (χ c1 ) and −3/5 (χ c2 ), in agreement with the alignment scenario suggested by the measured χ c2 /χ c1 cross-section ratios (Fig. 2): these limit values correspond to two very different decay distribution shapes, but to the same pure J z = 0 angular momentum configuration of the χ c . The J/ψ λ ϑ from the weighted χ c1 and χ c2 feed-downs (blue band) is close to the values measured by CMS (squares) for the prompt sample, implying that the direct and feed-down terms have similar polarizations.
It is quite remarkable to observe that the difference ∆λ ϑ ≡ λ ϑ (χ c2 ) − λ ϑ (χ c1 ) is predicted with a rather high precision and, furthermore, reaches extreme values (around −1). In particular, in the p T ≈ 20 GeV region, where experimental measurements will be provided in the near future, the prediction is ∆λ ϑ = −0.80 ± 0.05, implying a strong deviation from the mild polarizations shown in Fig. 1-bottom. Comparing the discriminating power of this result to the corresponding predictions of Ref. [14] ( Fig. 4), λ ϑ (χ c1 ) = 0.25 +0.08 −0.05 and λ ϑ (χ c2 ) = 0.10 +0.15 −0.20 , one can see the crucial importance of a proper treatment of the uncertainties and correlations affecting the experimental data. It is also relevant to note that, thanks to the cancellation of most experimental systematic uncertainties, the difference λ ϑ (χ c2 ) − λ ϑ (χ c1 ) can be measured with maximal significance and accuracy.
We will now discuss in more detail the theory-data comparison. To discern shape differences more easily than in the logarithmic-scale plots of Fig. 4, we present in Fig. 6 some of the results in the form of ratios, in a linear scale. In the top panel, we can see that the ratio between the experimental "unpolarized" band and the state-of-the-art 1 S This effect might represent a residual limitation of current finite-order perturbative calculations, as suggested by the observation that the ratio shows a more pronounced nonflatness when the SDC is calculated at NLO without fragmentation contributions ("NLO" band), and is not flat at all when we use the LO SDC as reference ("LO" band). The differences between these three ratios provide a pedagogical illustration of the improvements made in the successive evolutions of the calculations. The bottom panel shows that the ratios between the measured χ c1 (blue ratio 0 2 4 6 8 10 LO [8] Fig. 6. Ratios of direct-production charmonium cross section shapes for different combinations of the measured and/or calculated terms already presented in Fig. 4. For visibility reasons, all ratios are normalized to unity at pT/M = 4. filled band) or χ c2 (pink open band) p T /M -differential cross sections and the corresponding unpolarized-ψ cross section are practically identical to each other, and essentially flat, offering an effective representation of the strong experimental observation mentioned above. It is interesting to compare these two bands, exclusively determined by the measurements, with the two corresponding (and completely independent) theory ratios, here represented by the (blue and pink) solid lines, calculated as the ratios between suitable combinations of the 3 S SDCs (analogous to the χ c1 or χ c2 ) and the 1 S [8] 0 SDC (analogous to the unpolarized-ψ). The dashed curves surrounding the solid ones reflect the 1.3% uncertainty on K χ , already shown in Fig. 4. Also this ratio deviates from a flat function in the lower part of the p T /M range, but this deviation is a relatively small effect, as can be judged by comparing it with the corresponding rate of increase of the individual components, 3 S  1,2 / 1 S [8] 0 (brown and green curves). It is actually quite remarkable to see how effective is the mutual cancellation of the individual (steep) variations, in the combinations pertinent to the χ c1 and χ c2 states. Also taking into consideration that the P-wave SDCs seem to be affected by a slower convergence of the perturbative series than the S-wave SDCs [27], the present level of agreement between the shapes of the χ c -to-1 S [8] 0 predicted ratios and the corresponding χ c -to-unpolarized-ψ measured bands can be considered very promising. As a matter of fact, and despite the initial impression of unnecessary complexity expressed by Fig. 1, we see that NRQCD provides predictions that are, already today, very close to reproducing the uniformity of the observed p T /M trends, as well as the small measured S-wave polarizations. This unexpected agreement is the result of a series of cancellations, which, given their fragile and unstable nature, must be tested with precise ingredients. Further improvements in the perturbative calculations, especially for the P-wave SDCs, are needed for more conclusive statements.

Summary
The χ c1 and χ c2 states have, both, p T /M distributions well compatible with being identical to the one of the J/ψ mesons. This conclusion results from the study of the full set of charmonium data and has a much higher significance than one would obtain if only considering the χ c cross section measurements, given their limited precision and p T coverage in comparison to the J/ψ and ψ(2S) measurements. This is a very specific and non-trivial experimental observation, seemingly in contradiction, at least a priori, with the expectations of NRQCD, given the significantly different shapes of the relevant SDCs. For example, the ratios between the 3 S 1,2 SDCs, dominant terms of the factorization expansion for χ production, to the 1 S [8] 0 SDC, very well describing J/ψ and ψ(2S) production, have an order-of-magnitude increase from low to high p T /M . Therefore, within NRQCD, one would a priori expect different p T /M dependences for the χ c1 , χ c2 and ψ states. Remarkably, thanks to mutual cancellations of the steep SDC shapes differences, NLO NRQCD calculations approximately reproduce the similarity betwen the χ c1 , χ c2 and ψ cross setions shapes, giving a satisfactory description of charmonium production as measured at mid-rapidity by the ATLAS and CMS experiments.
On the other hand, at the present state of the SDC calculations and within the limits of the currently adopted v 2 expansion, this agreement comes with a definite prediction of strong and opposite χ c1 and χ c2 polarizations. It is worth restating this conclusion with different words. The χ c1 and χ c2 differential cross sections have p T /M dependences compatible with being identical to the J/ψ distribution, an experimental observation that NRQCD (given the presently available NLO SDCs) can only reproduce in a very specific and non-trivial configuration, leading to a remarkable prediction: the χ c1 and χ c2 polarizations are as different from each other as physically possible.
If confirmed experimentally, through an accurate measurement of the variable λ ϑ (χ c2 ) − λ ϑ (χ c1 ), the existence of strong χ c1 and χ c2 polarizations (an exception among all quarkonia observed by high-p T experiments) would be a big step forward to confirm the existence of the diversified and polarized processes that are at the heart of NRQCD. If, instead, similar and weak χ c1 and χ c2 polarizations will be measured, it will be crucial to investigate if the predicted strong and opposite polarizations, experimentally falsified, are caused by approximations and inaccuracies of the presently available fixed-order perturbative calculations or from problems in the conceptual foundations of the theory. In that case, NRQCD would be facing a big challenge: even if future improvements of the P-wave SDC calculations would eventually make the χ c1 and χ c2 polarization predictions compatible with the measurements (e.g., building upon the recent progress on fragmentation corrections [27]) one would still think that the homogeneity of the observed kinematic patterns deserves a more natural theoretical explanation than a series of "coincidences" cancelling out the variegated complexity of NRQCD. In either case, accurate measurements of the χ c1 and χ c2 polarizations constitute a decisive test of NRQCD.