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 [1,2,3]. 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 [4], 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. [5, 6].

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_\mathrm{T}\), and rapidity, y, have been measured in pp collisions at \(\sqrt{s} = 7\), 8 and 13 TeV for the \(\mathrm{J}/\psi \), \(\psi \mathrm{(2S)}\) and \(\varUpsilon \)(nS) vector states, by ATLAS [7,8,9], CMS [10,11,12] and LHCb [13,14,15,16,17,18,19]. Also the polarizations have been measured for these states, at mid-rapidity by CMS [20, 21] and at forward rapidity by LHCb [22,23,24]. 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 [25,26,27,28,29]. First experimental measurements on the polarizations of the \(\chi _{{c}1} \) and \(\chi _{{c}2} \) states have recently been reported by CMS [30], 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 [7,8,9,10,11,12,13,14,15,16,17,18,19, 25,26,27,28,29] 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 \(\mathrm{J}/\psi \) production. In fact, while the prompt \(\mathrm{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 \(\mathrm{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 [31]. 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 [32, 33].

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-\(\mathrm{J}/\psi \) feed-down fractions and \(\chi _{c}\) polarizations, and also of the properties of direct \(\mathrm{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 \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) data impose on the \(\chi _{c}\) cross sections and polarizations, through the feed-down contributions present in the \(\mathrm{J}/\psi \) case and absent in the \(\psi \mathrm{(2S)}\) case. The transition probabilities from heavier to lighter states needed in this work are all well known and listed in the PDG tables [34]. Our results are, therefore, fully data driven.

Besides reporting the \(\chi _{c}\) polarizations and feed-down contributions to \(\mathrm{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 \(\mathrm{J}/\psi \), \(\psi \mathrm{(2S)}\), \(\chi _{{c}1} \) and \(\chi _{{c}2} \) differential cross sections measured by ATLAS and CMS at 7 TeV [8, 10, 25], as well as the \(\chi _{{c}2} \) over \(\chi _{{c}1} \) cross section ratio [26] and the \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) polarizations [20] measured by CMS at the same energy. All of these measurements have been reported as functions of \(p_\mathrm{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 \(\mathrm{J}/\psi \)) in three \(\mathrm{J}/\psi \) \(p_\mathrm{T}\) bins. The total number of independent data points is 108 and they are all shown in Figs. 1, 2 and 3.

Fig. 1
figure 1

Mid-rapidity prompt quarkonium cross sections measured in pp collisions at \(\sqrt{s} = 7\) TeV by ATLAS (red markers) [8, 25] and CMS (blue markers) [10]. The inset shows the \(\chi _{{c}2} \) to \(\chi _{{c}1} \) cross section ratio [26]. The curves represent the result of the fit described in the text

Fig. 2
figure 2

Polar anisotropy parameter \(\lambda _\vartheta \), in the helicity frame, measured by CMS in pp collisions at \(\sqrt{s} = 7\) TeV, for prompt \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) dimuon decays [20]. Values corresponding to two (\(\mathrm{J}/\psi \)) or three (\(\psi \mathrm{(2S)}\)) rapidity bins were averaged. The curves represent the result of the fit described in the text

Fig. 3
figure 3

The \(\chi _{{c}2} / \chi _{{c}1} \) yield ratio vs. \(|\cos \vartheta |\) (in the helicity frame), for three \(\mathrm{J}/\psi \) \(p_\mathrm{T}\) bins (8–12, 12–18 and 18–30 GeV), as measured by CMS in pp collisions at \(\sqrt{s} = 8\) TeV [30]. The curves represent the result of the fit described in the text

The \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) \(p_\mathrm{T}\)-differential cross sections measured by ATLAS in the dimuon decay channel [7] 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_\mathrm{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_\mathrm{T}/M\), the ratio between the \(p_\mathrm{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_\mathrm{T}/M \) to simplify the notation.

The \(p_\mathrm{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

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

This functional shape describes very well the quarkonium transverse momentum distributions in different kinematic domains [6, 35, 36]. The parameter \(\gamma \) (having the meaning of the average \(p_\mathrm{T}/M\) squared) defines the function in the low-\(p_\mathrm{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_\mathrm{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 \mathrm{(2S)} \rightarrow \chi _{c1,2} \; \gamma \); \(\psi \mathrm{(2S)} \rightarrow \mathrm{J}/\psi \; X\); \(\chi _{c1,2} \rightarrow \mathrm{J}/\psi \; \gamma \). The rule for the momentum propagation from mother to daughter is, on average, \(p_\mathrm{T}/m = P_\mathrm{T}/M\), where M (m) and \(P_\mathrm{T}\) (\(p_\mathrm{T}\)) are, respectively, the mass and laboratory transverse momentum of the mother (daughter) particle [6]. 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 [37]. 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-\(\mathrm{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 [37]. Consequently, the terms longitudinal and transverse will here always refer to the yields of events where the daughter \(\mathrm{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_\mathrm{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 (\(\mathrm{J}/\psi \), \(\chi _{{c}1} \), \(\chi _{{c}2} \) and \(\psi \mathrm{(2S)}\)) would require eight \(\beta \) parameters, describing the \(p_\mathrm{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_\mathrm{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. [6, 36], the production cross sections of the \({^3\mathrm{S}_1}\) and \({^3\mathrm{P}_J}\) quarkonium states measured by ATLAS and CMS, at both 7 and 13 TeV, follow remarkably uniform patterns as a function of \(p_\mathrm{T}/M\). Such scaling is expected, from dimensional analysis considerations [36], 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_\mathrm{T}/M\) distribution shapes of the P-wave states and those of the S-wave states [38]. 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_\mathrm{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 [36]. 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 \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) have identical \(p_\mathrm{T}/M\)-dependent kinematic patterns, while keeping, obviously, two independent yield normalizations at \((p_\mathrm{T}/M)^{*}\). This assumption, suggested and supported by the observed universality of the \(p_\mathrm{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 [1]. Therefore, in the fit, a single parameter, \(\beta _{\mathrm {total}}^{\psi }\), describes the \(p_\mathrm{T}/M\) dependences of the total \(\mathrm{J}/\psi \) and \(\psi \mathrm{(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 \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) \(p_\mathrm{T}/M\) distributions and polarizations is attributed to the \(\chi _{{c}1} \) and \(\chi _{{c}2} \) feed-down contributions: \(\mathrm{J}/\psi \) and \(\psi \mathrm{(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 \mathrm{(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_\mathrm{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_\mathrm{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 \(\mathrm{J}/\psi \) and \(\psi \mathrm{(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 \(\mathrm{d} \sigma /\mathrm{d} p_\mathrm{T} \) values (in nb/GeV) at \((p_\mathrm{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. [34].

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_\mathrm{T}\) trends different from each other and from that of the \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) mesons, so that the results will provide a test of the importance of the \(p_\mathrm{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_\mathrm{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_\mathrm{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_\mathrm{T}/M)^{*}\) remaining free parameters in all options.

While the \(\mathrm{J}/\psi \) and \(\psi \mathrm{(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 \(\mathrm{J}/\psi \) transverse momentum. Those data points are parametrized with the expression

$$\begin{aligned} 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 }, \end{aligned}$$
(2)

where \(\xi _i\) (\(i=1,2,3\)) are the three average-\(p_\mathrm{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 \(\mathrm{J}/\psi \) from \(\chi _{{c}1} \) or \(\chi _{{c}2} \), calculated at the three \(p_\mathrm{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 \mathrm{(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 \mathrm{(2S)}\), \(\chi _{{c}1} \) and \(\chi _{{c}2} \) cross sections measured by ATLAS and between the \(\mathrm{J}/\psi \) and \(\psi \mathrm{(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. [34]. 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 \mathrm{(2S)} \rightarrow \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_\mathrm{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.

Table 1 Fit quality information. In all three fit variants there are 108 data points and 10 nuisance parameters

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, 2 and 3 for the baseline \(4\beta \) case. The smallness of the \(\chi ^2\) per degree of freedom, \(\chi ^2/\mathrm{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_\mathrm{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/\mathrm{ndf}\) value of the baseline analysis with those of the two variants mentioned before: \(\chi ^2/\mathrm{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) [39], the likelihood of the \(6\beta \) model is much smaller than that of the \(4\beta \) model.

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_\mathrm{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.

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-\(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\)-to-\(\mathrm{J}/\psi \) feed-down fractions are shown in Fig. 4. The results of the baseline (\(4\beta \)) fit are represented by the \(p_\mathrm{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_\mathrm{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.

Table 2 Values of the fitted parameters of interest, in the three considered scenarios
Fig. 4
figure 4

The fractions of the total prompt \(\mathrm{J}/\psi \) production rate due to feed-down decays from the \(\chi _{{c}1} \), \(\chi _{{c}2} \) and \(\psi \mathrm{(2S)}\) mesons, as a function of \(p_\mathrm{T}/M\)

Table 3 \(p_\mathrm{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_\mathrm{T}/M = 5\). The derived direct \(\mathrm{J}/\psi \) fraction is \((67.2 \pm 1.9)\)%
Fig. 5
figure 5

The polarization parameter \(\lambda _\vartheta \) of the \(\chi _{{c}1} \) and \(\chi _{{c}2} \) mesons (left), as well as of the \(\mathrm{J}/\psi \) mesons directly produced (same as of the \(\psi \mathrm{(2S)}\)) and produced in \(\chi _{{c}1} \) plus \(\chi _{{c}2} \) decays (right), as a function of \(p_\mathrm{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 \(\mathrm{J}/\psi \) mesons. The right panel also shows the polarization of the \(\mathrm{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_\mathrm{T}/M\)-independent results obtained in the \(1\beta \) fit, effectively representing averages over \(p_\mathrm{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 \(\mathrm{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_\mathrm{T}/M = 5\)) and in the \(1\beta \) option (averaged over \(p_\mathrm{T}/M\))
Fig. 6
figure 6

Ratio between the \(\chi _{{c}2} \) and \(\chi _{{c}1} \) cross sections, times the corresponding \(\chi _{{c}J} \rightarrow \mathrm{J}/\psi \, \gamma \) branching fractions, as a function of \(p_\mathrm{T}/M\), for two extreme polarization scenarios (left) and for the polarizations determined in our global-fit analysis (right)

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 \mathrm{(2S)}\) mesons. Finally, the ratio between the \(\chi _{{c}2} \) and the \(\chi _{{c}1} \) cross sections (times the corresponding \(\chi _{{c}J} \rightarrow \mathrm{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 6 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.

Fig. 7
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 \(\mathrm{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 \mathrm{(2S)}\) and of the directly produced \(\mathrm{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 [40,41,42,43,44,45] and vector boson [46,47,48,49,50,51,52,53] polarizations are known to be significantly non-zero and momentum dependent, as are those of low-\(p_\mathrm{T}\) quarkonia [54, 55]. 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. [41], 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 [56]. The other possibility is that the \(\mathrm{J}/\psi \) originates from the decay of a \(J = 0\) state, as expected to happen in the production from the \(^1\mathrm{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\mathrm{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_\mathrm{T}\)-independent unpolarized scenario has strong and rather remarkable physical implications.

Fig. 8
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 \(\mathrm{J}/\psi \) (left), the \(\psi \mathrm{(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 \(\mathrm{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 \mathrm{(2S)}\) polarization and, above all, of the difference between the \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) polarizations. This latter measurement can be performed, for example, by determining the ratio between the \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) angular distributions in a given \(p_\mathrm{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 \(\mathrm{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.

Fig. 9
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 \(\mathrm{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 \(\mathrm{J}/\psi \) production, and to extract the kinematic properties of the directly produced \(\mathrm{J}/\psi \). The analysis is fully data driven, not relying on any theoretical inputs. It uses LHC data on \(\mathrm{J}/\psi \), \(\psi \mathrm{(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 and 8 TeV.

A first result of the analysis is that all polarizations and cross section ratios are found to be perfectly compatible with being \(p_\mathrm{T}/M\)-independent. When the \(p_\mathrm{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 \(\mathrm{J}/\psi \) feed-down fractions are perfectly compatible with values obtained in a global analysis of low-\(p_\mathrm{T}\) data from fixed-target experiments [57], an observation that confirms the independence of such ratios on \(p_\mathrm{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. [30]) and all the indirect experimental information (red line).

Fig. 10
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 [30] (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 \(\mathrm{J}/\psi \) and \(\psi \mathrm{(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 \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) polarizations, potentially very precise given the cancellation of most systematic uncertainties. New measurements of the \(\psi \mathrm{(2S)}\) polarization, especially towards higher \(p_\mathrm{T}\), are also a top priority. In fact, the polarization of directly produced \(\mathrm{J}/\psi \) and \(\psi \mathrm{(2S)}\) states, accessible for the first time as a result of our global fit analysis, is found to be very small (\(\lambda _\vartheta ^{\mathrm{J}/\psi } = 0.04 \pm 0.06\)) and \(p_\mathrm{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_\mathrm{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_\mathrm{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.