Universal kinematic scaling as a probe of factorized long-distance effects in high-energy quarkonium production

Dimensional analysis reveals general kinematic scaling rules for the momentum, mass, and energy dependence of Drell-Yan and quarkonium cross sections. Their application to mid-rapidity LHC data provides strong experimental evidence supporting the validity of the factorization ansatz, a cornerstone of non-relativistic QCD, still lacking theoretical demonstration. Moreover, data-driven patterns emerge for the factorizable long-distance bound-state formation effects, including a remarkable correlation between the S-wave quarkonium cross sections and their binding energies. Assuming that this scaling can be extended to the P-wave case, we obtain precise predictions for the not yet measured feed-down fractions, thereby providing a complete picture of the charmonium and bottomonium feed-down structure. This is crucial information for quantitative interpretations of quarkonium production data, including studies of the suppression patterns measured in nucleus-nucleus collisions.


Introduction
Quarkonium production is a paradigm case study for the understanding of hadron formation, but the theoretical description of its basic mechanisms remains a challenge [1]. Both the attractiveness and the complexity of the problem reside in the quark-antiquark (QQ) binding process, whose degrees of freedom manifest themselves in the observed spectrum of bound states of different masses and spin-angular-momentum properties, remaining unclear how they are interrelated with the momentum distributions measured for the different states. Non-relativistic QCD (NRQCD) [2], as well as other approaches to quarkonium production (e.g., colour-singlet and colour-evaporation models [3]) are based on a pivotal postulate: the long-distance effects related to boundstate formation can be factorized out in the calculation of the cross sections. In this way, short-distance QQ cross sections (SDCs), containing the dependence on transverse momentum (p T ), rapidity (y), and pp centre-ofmass collision energy ( √ s), can be calculated perturbatively [4][5][6][7][8], independently of the long-distance matrix elements (LDMEs), which are proportional to the probability of the transition from the 2s+1 l c j pre-resonance QQ, with c = 1 (colour singlet) or 8 (colour octet), to the finalstate observable 2S+1 L J quarkonium. Here, s, l, and j (S, L, and J) denote the spin, orbital angular momentum, and total angular momentum of the primordial QQ state (quarkonium), respectively.
The LDMEs (in practice, free theory parameters) are considered constant: they depend on the initial QQ quantum numbers and on the bound-state properties (quarkonium mass M , quark mass m Q , velocity of the heavy quark in the quarkonium rest frame v, etc.), but not on the "laboratory" variables p T , y, and √ s, nor on the collision system. The factorization postulate, believed to be valid in the limit of heavy quark mass, for which the time of creation of the QQ pair, ∆t ∼ 1/m Q , is distinctively smaller than the bound-state formation time, ∼ 1 fm, has not yet been demonstrated from QCD principles [1]. Also because of inconsistencies between the NRQCD predictions and the measurements, especially concerning polarization [9,10], it is often argued that the validity of factorization might be restricted to specific kinematic domains (high p T ) or to the (heavier) bottomonium family.
The analysis reported in this paper addresses the problem in a data-driven way, comparing the patterns of momentum, mass, and collision-energy dependences measured for the "reference" case of inclusive Drell-Yan production with the corresponding quarkonium patterns. It uses mid-rapidity pp LHC data and only relies on general kinematic scaling properties of the production cross sections.
2 Dimensional scaling in kinematics of particle production The p T and y double-differential cross section for the inclusive production of a given particle in pp collisions can be written as whereσ is the cross section of the parton-initiated process, P (x) are the parton density functions of the proton, x 1 and x 2 are the proton momentum fractions carried by the colliding partons, andŷ is the particle rapidity in the parton-parton centre-of-mass frame,ŷ = y − 1 2 ln(x 1 /x 2 ). An implicit delta function sets the particle's p T to the measured one. The cross section depends on the pp collision energy, √ s, through the relation √ŝ is the parton-parton collision energy, function of p T ,ŷ and M : We start by discussing inclusive Drell-Yan lepton-pair production. After integrating over the lepton emission angles in the dilepton rest frame (and over the degrees of freedom of the recoil system), the kinematics in the rest frame of the colliding partons is fully described by p T ,ŷ, and M . Since the dimensionality of the differential partonic cross section is [energy] −3 , we can write f being a dimensionless function of the dimensionless variablesŷ and ξ, where we use the definition ξ ≡ p T /M to have more compact equations. When the mixture of production mechanisms does not change with M , the partonic cross section, calculated at any arbitrary reference kinematic point (ξ ,ŷ ), scales as M −3 . This property can be tested using Drell-Yan cross sections, measured in a range of masses sufficiently smaller than the Z boson mass to stay unaffected by the strongly mass-dependent interference of γ-and Z-exchange mechanisms. The M scaling of the observable cross section depends on the pp collision energy via the parton densities, an effect that can be singled out in a purely data-driven way at LHC energies and mid rapidity, given the small average x 1,2 . For illustration, we parametrize P (x) as explicitly showing the factor describing the low-x behaviour and globally representing the remaining factors by the function q, typical examples of which being q(x) = , all tending to unity in the small-x limit. Therefore, Moreover, is the rapidity of the system of colliding partons. Using the previous relations and releasing the ξ integration in Eq. 1, we obtain the pp ("dressed") version of the partonic cross section, Eq. 3, where b = 2B − 2 and F ξ, y; The function g absorbs all the factors depending on ξ andŷ, including those contained in √ŝ /M ,Ê/M , andp/M (Eq. 2). F is a fully experimentally determinable dimensionless function of the dimensionless scaling variables ξ and y. It also depends on √ s and M , but always through the √ s/M ratio, via the terms q(x 1,2 ) = q( √ŝ / √ s e ±ŷ0 ) and the integration extremes ± ln( √ s/ √ŝ min ) = ± ln( √ s/M ) × (1 − ln(ξ + 1 + ξ 2 ) / ln( √ s/M )). The sensitivity of F on √ s/M decreases with increasing √ s and towards mid rapidity. We normalize F over the (ξ, y) range of the analysis and assimilate its integral I( √ s/M ) in a numerical redefinition of the b exponent (always possible to a very good approximation over a small range of √ s values, as in our analysis of √ s = 7 and 8 TeV data), with F dξ dy = 1. The (p T /M, y)-integrated cross section becomes σ = M −2 ( √ s/M ) β and its measurement at two or more √ s values provides a determination of β. These considerations lead to the following proposition. In a domain where Drell-Yan production mechanisms do not interfere in varying proportions, the joint distribution of the variables p T /M , y, and √ s/M is universal : its shape does not depend on M . At any given kinematic point ((p T /M ) , y ) and √ s, the dσ/dp T cross section scales like M −(3+β) , where β expresses the increase of the integrated cross section with √ s. The scaling of the cross section with M , integrated over p T , was previously discussed in Mid-rapidity Drell-Yan differential cross sections as a function of the dilepton mass, as measured by ATLAS [12] and CMS [13,14].
Ref. [11], but using process-specific arguments rather than completely general dimensional-analysis considerations.
We will now consider the Drell-Yan dσ/dM differential cross sections measured by ATLAS [12] and CMS [13,14], in pp collisions at √ s = 7 and 8 TeV, shown in Fig. 1 for the M < 50 GeV range. It can be recognized from the previous discussion that dσ/dM has the same M −3 ( √ s/M ) β scaling behaviour as dσ/dp T at fixed p T /M and y. Simultaneously fitting all four data sets to a single M −αDY function, only considering pointto-point uncorrelated uncertainties in order to determine the shape parameter with maximum significance, gives a remarkably precise result: α DY = 3.63 ± 0.03. Fitting each data set individually gives: 3.60 ± 0.05 (CMS 7 TeV), 3.63±0.05 (CMS 8 TeV), 3.75±0.07 (ATLAS 7 TeV, 2010) and 3.60 ± 0.07 (ATLAS 7 TeV, 2011).
The β exponent, reflecting the √ s dependence, can be derived from the ratio between the cross sections measured at 8 and 7 TeV, directly reported by CMS [14], leading to β = 0.73 ± 0.15, where the uncertainty is obtained assuming a relative uncorrelated luminosity uncertainty of 2%. The resulting α − β difference, 2.90 ± 0.15, is perfectly consistent with 3, as expected.
The Drell-Yan cross sections so far reported by the LHC experiments are integrated over p T . The future availability of p T -and p T /M -differential measurements will allow the realisation of more detailed and accurate tests. It is worth noting that the existence of these simple and general kinematic scaling rules has been ignored in all the published analyses of LHC data, an unfortunate situation because of the resulting loss in physics value and also because they can be very useful experimental tools. In fact, detector and trigger acceptances drastically sculpt the reconstructed dilepton distributions, especially in the low momentum and low mass regions. Verifying that the mea-sured p T /M and M spectra, after efficiency corrections, etc., satisfy the expected dimensional scaling constitutes a powerful cross-check of the analysis procedure and a validation of the reported systematic uncertainties.

Dimensional scaling in quarkonium production
Moving now to our main case study, quarkonium production, we write the differential cross section for the inclusive production of a given S-wave state as where the overall factor m −3 Q matches the global dimensionality of the observable. The only specification of the functions L i is that they have dimension [energy] 3 , formally compensated by the m 3 Q denominators, while the F i are dimensionless shape functions (defined with F dξ dy = 1, as in the Drell-Yan case). The √ s/M power-law factor represents the modification from partonic to observable level, as from Eq. 3 to Eq. 9. The coefficient β is the same as for Drell-Yan production. From the previous discussion, we can obtain a precise evaluation of its value at √ s 7 or 8 TeV as the difference between the average experimental mass scaling exponent α DY and the expected bare-cross-section scaling exponent: The expression above, Eq. 10, is built in analogy to the NRQCD factorized expansion. In this limit, F i and L i encode the information on, respectively, the hard scattering process producing a QQ pair, of mass 2 m Q , and its ensuing long-distance transition to the final bound state Q, of mass M . In NRQCD the LDMEs L i are independent of laboratory kinematics and only depend on degrees of freedom naturally defined in the Q rest frame (M , m Q , v, quantum numbers, spectroscopic energy levels, etc.), while the kinematic dependence is contained in the SDCs, In this formulation, a consequence of the factorization hypothesis is that each F i function should be independent of M , being, in particular, the same for all directly produced 3 S 1 charmonia and bottomonia. In fact, it is important to note that the mass difference between the Q and QQ states plays no role in the p T /M,ŷ dependence of F i , given that, similarly to the physical transitions between states of the same quarkonium family [15], the QQ → Q transition preserves the ratio between laboratory-vectormomentum and mass (and, therefore, both p T /M andŷ).
While these notations have been chosen to accommodate the NRQCD factorization expansion as a limit case, Eq. 10, with generic L i and F i both redundantly depending of the relevant variables, represents a fully general tem-   Fig. 2. Mid-rapidity prompt quarkonium cross sections in pp collisions at √ s = 7 TeV (left) and 13 TeV (right), as measured by ATLAS (red markers) [16][17][18] and CMS (blue and green markers) [19][20][21] (top panels), and after scaling up all the normalizations to match those of the J/ψ cases (middle panels). For each collision energy, the curves represent a single universal empirical function, of shape determined by a simultaneous fit to all data points of pT/M > 2 and normalizations specific to each quarkonium state. The bottom panels show the pulls between each data point and the fitted function. plate of the cross section, with no prejudice on its physical scaling properties.
Very interesting physical indications transpire from the unforeseen simplicity of the trends measured by ATLAS and CMS, as discussed in the next paragraphs.
The first experimental input to our considerations are the quarkonium cross sections measured at mid-rapidity (|y| 2) by the ATLAS and CMS experiments, shown in the top panels of Fig. 2 for √ s = 7 TeV (left) and 13 TeV (right). All the measurements, from the J/ψ to the Υ (3S), scale with p T /M in a state-independent way, at least for not very low p T (p T /M 2). This "universality" is well illustrated by the middle panels, which clearly show that, after rescaling the state-specific normalizations to those of the J/ψ, the cross section shapes become indistinguishable from each other. The goodness of this universal scaling can be better appreciated by looking at the bottom panels, which show (in a linear scale) the pull distributions, i.e. the differences between each data point and the universal fitted function, divided by the measurement uncertainty. No systematic trends are seen in the pull distributions and the observed deviations are very well compatible with statistical fluctuations. The cross section fits consider correlations between the luminosity uncertainties in each data set and the pulls, evaluated to check the consistency with one common shape, are calculated excluding such uncertainties.
A slightly broader p T /M distribution is observed at the higher energy, as shown in Fig. 3-top. The fact that the distinct P-and S-wave states show a compatible kinematic scaling is discussed in Ref. [22]. Here we focus on the precisely measured cross sections of the closely related 3 S 1 states. Their universal p T /M scaling, at a given energy, indicates that the same mixture of processes (or one common dominating process) describes the production of all these states.
The second experimental indication for our following discussion is the mass scaling of the cross section from charmonium to bottomonium. By exploiting the p T /M scaling we can determine it without relying on modeldependent extrapolations to low p T : it is sufficient to consider the fitted p T /M distributions at an arbitrary p T /M = (p T /M ) . We consider for each family the meson closest to the ground state, J/ψ and Υ (1S), at the two energies ( Fig. 3-bottom).
To obtain the yield of directly produced J/ψ mesons, we subtract from the fitted prompt-J/ψ normalization those of the χ c1 , χ c2 and ψ(2S) states (always considered at the same p T /M value: as mentioned before, p T /M is transferred unchanged from mother to daughter), weighted by the branching fractions [23] of the respective feed-down processes. The result is a J/ψ directproduction fraction of 31.9%. For the Υ (1S) we assume a direct-production fraction of (50 ± 10)%. Later in this paper we will need the corresponding fractions for the Υ (2S) and Υ (3S), which we assume to be 60% and 70%, respectively, also with ±10% uncertainties. These values were estimated from Υ (nS) yield ratios measured by ATLAS [17] and CMS [20], complemented by LHCb data on the (large) χ b (mP) → Υ (nS) feed-down fractions [24].
These experimental facts constrain and specify the functions L and F appearing in Eq. 10. The independence of the p T /M scaling on either M or m Q , at a given energy, indicates that the p T /M and (M, m Q ) dependences do not "mix": L × F = L(m Q , M, √ s/M ) × F(ξ, y, √ s/M ) (further studies considering forward-rapidity, low-p T data will assess the combined p T /M, y dependence of F, here effectively a function of only p T /M ). In other words, there is experimental evidence that, at mid rapidity and not too low p T , it is possible to describe quarkonium production "factorizing" short-and long-distance effects, defined as the dependences on, respectively, the laboratory momentum of the detected meson (F) and the specific boundstate properties (L).
The measured m −6 Q scaling of the differential cross section at a fixed p T /M (and y) further specifies L. By precisely equating the explicit m Q dependence of the expression in Eq. 10 (coming from the overall factor and the denominators of the L terms, with F now taken independent of m Q ), such result leaves no margin for a dependence of L on m Q , if not counterbalanced by a dependence on M and/or √ s: L = L(m Q /M, m Q / √ s). However, given that no significant difference in mass scaling is observed at 7 and 13 TeV, a m Q / √ s dependence would be experimentally indistinguishable from the analogous global energy scaling of all quarkonium cross sections already accounted for by the β power law in Eq. 10. It is therefore always possible to define the long-distance factors as √ sindependent, thereby reducing them to functions of the kind L = L(m Q /M ), as is illustrated by the discussion starting in the next paragraph.
It should be clear from these considerations that such L functions do not coincide with the LDMEs of NRQCD: instead of being defined by setting an energy scale within the theory, they are built on the basis of dimensionalanalysis and data patterns, becoming universal, experimentally definable quantities. Apart from differences in operative definitions, we note that the remarkably simple picture of S-wave quarkonium production emerging from the data mirrors well the primary concepts of NRQCD factorization and of universality of the long-distance boundstate formation effects. A future verification of the corresponding experimental patterns in different collision systems can further probe these fundamental concepts.
By adopting now the guiding idea of factorization of short-and long-distance effects, as supported by data, we can analyse in more detail the mass scaling, using all S-wave states. ture": within each quarkonium family, the cross section decreases faster than from one family to the other. In the factorization perspective, the global drop represents the short-distance scaling, reflecting the change from M (QQ) 2 m c to 2 m b . The steeper change within each family reflects the M -dependent QQ → Q transition probability L(m Q , M ). The factorized interpretation is represented by the fit curves shown in the figure. First (red and blue curves) the charmonium and bottomonium cross sections are extrapolated, respectively, to 2 m c and 2 m b (as defined above). From the resulting σ(2m b ) / σ(2m c ) ratio, the short-distance charmonium-to-bottomonium mass scaling (black curve) is determined as ∝ m −α QQ Q , with α QQ = 6.63 ± 0.08, practically identical to the previously obtained result, considering only the J/ψ and the Υ (1S), leading to the m −6 Q partonic-level scaling discussed above. Here we focus on the long-distance mass dependence. The bottom panel of Fig. 4 shows the dependence of the Interestingly, as shown in Fig. 5, this monotonic dependence can be seen as a tight correlation of the longdistance factors to the quarkonium binding energy, defined as 2M (D 0 ) − M (ψ(nS)) or 2M (B 0 ) − M (Υ (nS)) and calculated with mass values from Ref. [23]. As a function of this variable we see another "universal" trend: the data points are well described by the parametrization σ(Q)/σ(2 m Q ) ∝ E δ binding , common to charmonium and bottomonium and identical at 7 TeV (δ = 0.63 ± 0.02) and 13 TeV (δ = 0.63 ± 0.04). The central values of the δ and β parameters are identical by coincidence. This result gives further support to the idea that the dependence of the cross sections on the bound-state masses is a "factorizable" long-distance effect, independent of the laboratorymomentum.
We have until now considered only S-wave states. Given the absence of detailed χ b cross section data at mid rapidity, we can put to test the assumption that the "universal" σ(Q)/σ(2 m Q ) ∝ E δ binding dependence of the long-distance factors, with δ 0.63, can be extended to the P-wave states, even if with a different normalization constant multiplying E δ binding to account for a dependence on the orbital angular momentum. In this way, χ c data come to constrain the χ b (nP) cross sections, through the relation σ(χ bJ (nP)) = [E binding (χ bJ (nP))/E binding (χ cJ )] δ σ(χ cJ ), corresponding to the assumption that there is no dependence on J, n, or quark flavour.
Using branching-ratio measurements from Ref. [23], the feed-down structure of quarkonium production can be fully predicted. The result is presented in Table 1 (see also Appendix A). It should be kept in mind that these results correspond to the p T /M > 2 region. The predicted feeddown fractions of Υ (nS) production from χ b (nP) states are in reasonable agreement with forward-rapidity LHCb measurements merging the J = 1 and J = 2 signals [24], considered for p T /M > 2.

Summary and discussion
At the current level of experimental precision, midrapidity LHC proton-proton data for inclusive charmonium and bottomonium production are well described by a simple parametrization reflecting a universal (i.e. stateindependent) scaling with two variables: the shapes of the p T distributions of all states become one common shape as a function of p T /M and the normalization of this shape scales in a simple monotonic way with E binding , at least for the S-wave states. While the cross section shape as a function of p T /M depends on √ s, the normalization scaling with E binding is found to be identical at 7 and 13 TeV.
Having a simple empirical parametrization faithfully describing the proton-proton measurements is certainly very useful for model-independent studies of quarkonium production, especially when considering more complex collision systems such as proton-nucleus and nucleus-nucleus collisions. It is also very convenient for the tuning of Monte Carlo simulations. More importantly, the analysis reported in this paper, exclusively based on non-trivial (albeit potentially misleadingly simple) dimensional analysis arguments, reveals significant experimental evidence supporting that quarkonium production can be understood as being factorized between short-and long-distance phases. This data-driven result mirrors very well the general concept of NRQCD factorization into process-dependent SDCs and universal LDMEs. More generally, seeing evidence of factorization in the LHC data is an important step towards establishing that the QQ pairs are produced uncorrelated from the rest of the event, allowing us to experimentally probe how they evolve and bind into the observable spectra of quarkonium states.
The experimental verification of the factorization ansatz can be extended by including in the analysis the available forward-rapidity LHCb measurements. Such study, currently ongoing, has a higher order of complexity: since the p T /M distributions show a significant y dependence, the determination of scaling patterns is no longer an effectively one-dimensional problem as in the case where only mid-rapidity data are considered. Further tests, which may become possible in the future, include the study of the production of quarkonia associated to specific final state particles and the search for simplicity patterns in different kinds of collisions.
More data on the production of P-wave quarkonia are needed, especially in the bottomonium family: the χ b absolute cross sections have not yet been measured at mid rapidity, and the forward rapidity results, reported relative to Υ (nS) production, do not distinguish between the J = 0, 1, and 2 states. It is reassuring, however, to see that those forward rapidity measurements are in reasonable agreement with the pattern of χ bJ (nP) → Υ feeddown fractions that we obtained assuming that the simple E binding scaling pattern can be extended to the χ states, an assumption that follows very naturally from the work presented in this paper. The detailed feed-down structure of the charmonium and bottomonium families we have determined in this work needs to be tested by new measurements, filling gaps in the experimental picture. It cannot be excluded that future measurements of direct χ cJ and χ bJ production will show that the scaling with binding energy is not the same for the P-and S-wave states. Such a result would still be very interesting, naturally. For example, the difference in scaling trends, when known, could be used to see if new resonances, as the X(3872), align better with the S-or the P-wave mass trends, thereby contributing to the understanding of their nature.
A complete account of all the feed-down fractions is particularly important in quantitative analyses of the