Measurement of quarkonium production cross sections in pp collisions at $\sqrt{s}=$ 13 TeV

Differential production cross sections of J/$\psi$ and $\psi$(2S) charmonium and $\Upsilon$(nS) (n = 1, 2, 3) bottomonium states are measured in proton-proton collisions at $\sqrt{s} =$ 13 TeV, with data collected by the CMS detector at the LHC, corresponding to an integrated luminosity of 2.3 fb$^{-1}$ for the J/$\psi$ and 2.7 fb$^{-1}$ for the other mesons. The five quarkonium states are reconstructed in the dimuon decay channel, for dimuon rapidity $|y|<$ 1.2. The double-differential cross sections for each state are measured as a function of $y$ and transverse momentum, and compared to theoretical expectations. In addition, ratios are presented of cross sections for prompt $\psi$(2S) to J/$\psi$, $\Upsilon$(2S) to $\Upsilon$(1S), and $\Upsilon$(3S) to $\Upsilon$(1S) production.


Introduction
Since the discovery of heavy-quark bound states, quarkonium production in hadronic collisions has been the subject of many theoretical and experimental studies. A well established theoretical framework to describe quarkonium production is nonrelativistic quantum chromodynamics (NRQCD) [1][2][3], an effective theory that assumes that the mechanism can be factorized in two steps. In the first step, a heavy quark-antiquark pair is produced in a given spin and orbital angular momentum state, either in a color-singlet or color-octet configuration. The corresponding parton-level cross sections, usually called short-distance coefficients (SDCs), are functions of the kinematics of the state and can be calculated perturbatively, presently up to next-to-leading order (NLO) [4][5][6][7]. In the second step, the quark-antiquark pairs bind into the final quarkonium states through a nonperturbative hadronization process, with transition probabilities determined by process-independent long-distance matrix elements (LDMEs). Unlike the SDCs, the LDMEs are presently not calculable and must be obtained through fits to experimental data [4][5][6][7][8][9]. Until recently, for directly produced S-wave quarkonia, the color-octet 3 S 1 term was thought to dominate, which would result in a strong transverse polarization of the mesons relative to their direction of motion (helicity frame) at large transverse momentum, p T .
We report the measurement of double-differential cross sections of the five S-wave quarkonium states J/ψ, ψ(2S), and Υ(nS) in pp collisions at √ s = 13 TeV by the CMS detector at the LHC. The increased center-of-mass energy and production cross sections provide an extended reach in p T and improved statistical precision relative to similar measurements at 7 TeV [23-26, 34]. The measurements performed at 13 TeV also provide the opportunity to test the √ s dependence of the cross sections and to check the validity of the factorization hypothesis and LDME universality implied in NRQCD.
The product of the branching fraction of quarkonia to muon pairs, B(Q → µ + µ − ), and the double-differential production cross section, d 2 σ/(dp T dy), in bins of p T and rapidity, y, is given by where N(p T , y) is the number of prompt signal events in the bin, L is the integrated luminosity, ∆y and ∆p T are the bin widths, and 1/( (p T , y)A(p T , y)) represents the average of the product of the inverse acceptance and efficiency for all the events in the bin. Only prompt signal events are considered. The nonprompt components of the J/ψ and ψ(2S) mesons, i.e. originating from decays of b hadrons, are separated using the decay length defined as = L xy · m/p T , where L xy is the distance measured in the transverse plane between the average location of the luminous region and the fitted position of the dimuon vertex, m is the mass of the J/ψ (ψ(2S)) from Ref. [35], and p T the transverse momentum of the dimuon candidate. For the prompt signal events, we do not distinguish between feed-down decays of heavier quarkonium states and directly produced quarkonia.

3 Acceptance and efficiencies The CMS detector, data set, and event selection
The analysis uses dimuon events collected in pp collisions at √ s = 13 TeV with the CMS detector. The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid [36]. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [37].
The data were collected using a multilevel trigger system [38]. The first level (L1), made of custom hardware processors providing coarse momentum information, requires two muons within the range |η| < 1.6. Second (L2) and third (L3) levels, collectively known as the HLT (High-Level Trigger), are implemented in software. At these levels, the muon selection is refined, then opposite-charge muon candidates are paired and required to have an invariant mass in the regions 2.9-3.3, 3.35-4.05, or 8.5-11 GeV for the J/ψ, ψ(2S), and Υ(nS), respectively. The dimuon p T is required to be above 9.9 GeV for the J/ψ and above 7.9 GeV for the remaining states. For all five states, the dimuon rapidity is restricted to |y| < 1.25. A fit of the positions and momenta of the two muon candidates to a common vertex is performed, and the fit χ 2 probability is required to be above 0.5%. The sample collected with these triggers has a total integrated luminosity of 2.3 fb −1 for the J/ψ and 2.7 fb −1 for the other mesons. The lower value for the J/ψ is the consequence of the trigger prescaling that was applied to limit the rate during part of the data taking.
When reconstructing the five states offline, further requirements are applied: only muons with p µ T > 4.5 GeV in the range |η µ | < 0.3, or p µ T > 4.0 GeV in the range 0.3 < |η µ | < 1.4 are selected. The muons have to match the triggered pair and be identified as reconstructed tracks with at least five measurements in the silicon tracker and at least one in the pixel detector. The track is required to match at least one muon segment identified by a muon detector plane. Loose criteria are applied on the longitudinal and transverse impact parameters to reject cosmic rays and in-flight hadron decays. The dimuon vertex χ 2 probability is required to be greater than 1%. In the CMS magnetic field, the two muons can bend towards or away from each other; only the second type of event is considered in this analysis since the first type exhibits high trigger inefficiencies. The dimuon rapidity is restricted to |y| < 1.2.
The double-differential cross sections are presented in four (two) rapidity bins for the prompt J/ψ and ψ(2S) (Υ(nS)), and in several bins of p T , covering a p T range between 20 and 120 (100) GeV for J/ψ (ψ(2S), Υ(nS)), extending up to 150 (130) GeV for measurements integrated in rapidity.

Acceptance and efficiencies
The acceptance is calculated using simulated events produced with a single-particle event generator. The quarkonium states are generated with a flat y distribution and a realistic p T distribution derived from data [24,25], covering the analysis phase space. The PYTHIA 8.205 [39] Monte Carlo event generator is used to produce an unpolarized dimuon decay (corresponding to a flat dimuon angular distribution), also accounting for final-state photon radiation. The simulated events include multiple proton-proton interactions in the same or nearby beam crossings (pileup), with the distribution matching that observed in data, with an average of about 11 col-lisions per bunch crossing. The acceptance for events in a given (p T ,|y|) range is defined as the ratio of the number of generated events that pass the kinematic selection criteria described above to the total number of simulated events in that p T and |y| range. The acceptance depends on the quarkonium polarization. It is derived for the unpolarized scenario, which is compatible with experimental measurements within uncertainties. We also calculate multiplicative correction factors that allow, from the unpolarized case, to infer the acceptance that corresponds to three different values of the polar anisotropy parameter, λ HX θ , in the helicity frame: −1 (fully longitudinal), +1 (fully transverse), and k, with k reflecting the CMS measured value of λ HX θ for each quarkonium state [30,31], also used in Refs. [25,26]. The multiplicative factors to convert the cross sections calculated using the unpolarized scenario to the ones calculated employing one of the polarization scenarios described above are provided.
The single-muon trigger, reconstruction, and identification efficiencies are measured individually from data as a function of muon p T and |η|, applying a tag-and-probe [23, 34] technique on J/ψ and Υ(1S) candidates acquired with triggers that are independent from those used for the measurements of the yields. The individual efficiencies are multiplied and then parameterized using a sigmoid function. The dimuon efficiency is obtained as the product of the efficiencies of the two muons, multiplied by a correction factor, ρ, that takes into account the correlation between the two muons. The ρ factor is derived from data, using a trigger, independent from the ones used for the measurement of the yield, requiring a single muon at L1. ρ becomes increasingly important with higher dimuon p T , when the two muons are close to each other in space, causing the efficiency to decrease. Dimuon efficiencies are around 85% for the J/ψ and ψ(2S) up to a dimuon p T of 50 GeV and decrease slowly for higher p T due to the ρ factor. In the case of the Υ(nS) states, the dimuon efficiencies are nearly constant around 90%. The acceptance and efficiency term in Eq. (1) is obtained by averaging the values of the inverse of the acceptance times efficiency for all the individual dimuon candidates in each p T and |y| range.

Determination of the yields
The signal and background yields are obtained through an extended unbinned maximumlikelihood fit to the dimuon invariant mass distribution in the case of the Υ(nS) states, and to the dimuon invariant mass and decay length distributions for the J/ψ and ψ(2S) mesons. In both cases, the number of signal and background candidates are free parameters in the fit.
The three Υ(nS) signal peaks are modeled with Crystal Ball (CB) functions [40], composed of a Gaussian core, characterized by a mean m, a width σ m , and a tail characterized by two parameters, n and α. The CB function is used to account for the energy loss due to the final-state radiation of the muons. The mean mass values are fixed to those of the Particle Data Group [35], multiplied by a common factor that calibrates the mass scale, left as a free parameter in the fit. The width of the CB function is a free parameter only in the case of the Υ(1S), while the width of the CB functions describing the Υ(2S) and Υ(3S) peaks are fixed to the width of the Υ(1S), scaled by the ratio of their masses to the mass of the Υ(1S). The Υ(nS) dimuon mass resolution σ m is a function of rapidity and spans the range 60 to 90 MeV for |y| < 1.2 in the case of the Υ(1S). The tail parameters n and α are the same for all three CB functions; n is fixed and α is constrained to a Gaussian probability distribution. Both constraints are derived from a fit of the Υ(1S) dimuon invariant mass shape, using the p T -integrated distribution to reduce the statistical fluctuations. The background is modeled using an exponential function.
For the J/ψ and ψ(2S) mesons, an additional nonprompt component originating from the decay of b hadrons must be taken into account. The prompt and nonprompt yields are measured by fitting the dimuon invariant mass and decay length distributions. The J/ψ dimuon invariant 4 5 Systematic uncertainties mass distribution is modeled by the sum of a CB and a Gaussian function with common mean, while the corresponding ψ(2S) distribution is described using only a CB function. The widths of the CB and Gaussian functions, as well as the α of the CB functions, are free parameters. The σ m varies as a function of rapidity between 20 and 50 (40) MeV for the J/ψ (ψ(2S)) state. The m and n parameters are fixed to values derived from fits to the invariant mass distribution of the p T -integrated data. An exponential is used to describe the dimuon mass background. The decay length distribution is modeled by a prompt signal component represented by a resolution function, a nonprompt term given by an exponential function convolved with the resolution function, and a background term represented by the sum of a resolution function plus an exponential decay function to take into account prompt and nonprompt background components. The resolution function is modeled by the sum of two Gaussian functions whose widths are taken as the event-by-event decay length uncertainty, multiplied by global scale factors. The two scale factors are free parameters in the fit and are constrained with Gaussian probability distributions that are derived from fits to the p T -integrated data, less affected by statistical fluctuations. The effective width of the two Gaussian functions is approximately 25 µm.
To verify that the fits to the quarkonium states are unbiased and the uncertainties are correctly modeled, 1000 pseudo-experiments were produced from simulation. Similarly, simulated events were used to test the hypotheses made on the constraints of the parameters. Differences in the event-by-event uncertainty information between signal and background candidates could introduce biases in the fitting of the decay length using the simplified model described above, but we verified that these effects are negligible. Examples of fits to the invariant mass and decay length distributions are provided in Figs. A.1-A.2 of Appendix A.

Systematic uncertainties
Systematic uncertainties are due to the measurement of the integrated luminosity (2.3%) [41], the determination of the signal yields, and the dimuon efficiencies and acceptances. Uncertainties in the estimation of the yields are evaluated by changing the signal and background models used in the maximum-likelihood fits. To assess the systematic uncertainty in the modeling of the signal invariant mass distribution of each state, the n and α parameters of the CB function are varied by up to ±5 standard deviations, one at a time, while the mean, which is constrained in the nominal fit, is allowed to float. The half-differences between the largest resulting deviations of the signal yields measured in the fit from the nominal yields are added in quadrature to obtain an uncertainty in the modeling of the signal. The systematic uncertainty originating from a possibly imperfect description of the background is evaluated by changing the background model from an exponential to a linear function. The observed differences from the nominal signal yields are taken as a systematic uncertainty. The total uncertainty in the determination of the yields is obtained as the sum in quadrature of the uncertainties in signal and background, and is about 2.0% for all quarkonium states.
Uncertainties in the discrimination between charmonia that are promptly produced rather than originating from b hadron decays arise from the determination of the primary vertex position (the production point of the mesons, which enters in the calculation of the decay length) and from the modeling of the signal and background in the decay length distributions. We assess the uncertainty originating from the choice of the primary vertex by using an alternative to the average position of the luminous region, the position of the collision vertex closest to the dimuon vertex extrapolated towards the beam line. The systematic uncertainty related to the description of the background is evaluated by measuring the difference between the prompt fractions using the nominal fit and a fit modeling the background by the sum of four expo-nential functions and a simplified resolution function composed of only one single Gaussian function. To study the impact of imperfect modeling of the resolution function, the scale parameters of the Gaussian functions that had Gaussian constraints in the nominal fit are varied by ±1 standard deviation. Similarly, we assess the impact of modeling the nonprompt signal by fixing the parameterization of the exponential decay function. The systematic uncertainty stemming from the choice of the primary vertex is added in quadrature with the uncertainty derived from the fit strategy. The latter is calculated as half of the difference between the maximum deviations observed from the nominal fit when the above variations are applied one by one. The total systematic uncertainty in the determination of the nonprompt yield is less than 3% in almost all the (p T , y) bins for the J/ψ meson, without a dominant contribution from any one of the sources described above. The largest systematic uncertainty in the ψ(2S) measurements can reach a maximum of 16%, mostly owing to the uncertainty in the modeling of the background decay length distribution. The effect of pileup on the analysis results has been studied using both data and simulation, and found to be negligible.
Uncertainties in the single-muon efficiencies, reflecting their statistical precision as well as possible imperfections of the parametrization, are evaluated by varying the three parameters of the sigmoid function used to parameterize the single-muon efficiencies within their uncertainties. The resulting systematic uncertainties are nearly constant as a function of p T and are around 2.5% for the J/ψ and ψ(2S), and 1.8% for the Υ(nS) in the central regions |y| < 0.6 and around 1% in the remaining rapidity regions. The L3 single-muon efficiencies are calculated from simulations because of the low number of collected events useful for their measurements. The corresponding uncertainty is estimated to be 3%.
Systematic uncertainties related to the ρ factor are of three kinds. The first originates from the number of events available in the control sample collected with the independent trigger used to evaluate the ρ factor. The relative uncertainty is about 1% from 20 to 50 GeV and increases to about 5% near 100 GeV, with no dependence on rapidity. The measurement of the ρ factor also requires the evaluation of an additional single-muon efficiency using the tag-and-probe method, which introduces an uncertainty of about 1% at low p T (below 50 GeV) and up to 4% at high p T . Moreover, we assign the fractional difference in the ρ factor obtained from data and simulation as a systematic uncertainty. The difference is in the range 2-5% up to 60 GeV and increases slowly for higher p T , reaching a value of up to 15%, in the worst case. This is the dominant uncertainty for all the quarkonium states except the ψ(2S).
The finite number of events generated for the acceptance calculation imposes a systematic uncertainty of 0.5% at low p T and up to 6% at high p T . Other sources of systematic uncertainties, like the kinematic modeling of simulated events, are found to have a negligible influence on the acceptance calculation. The effect of the quarkonium polarization on the acceptance is not treated as a systematic uncertainty; instead correction factors are provided in Appendix Ato recalculate the cross sections according to different polarization scenarios.
For the cross sections measured in the rapidity-integrated range |y| < 1.2, we conservatively assign the total systematic uncertainties of the most-forward rapidity range, which are larger than the uncertainties for central rapidities. Taking advantage of the larger yields in the integratedrapidity range, an additional p T bin was added for each state. The systematic uncertainty in the yields for this bin was evaluated as described above for the other bins, while for other uncertainties the same value as in the neighboring lower-p T bin was used. It was verified that systematic uncertainties extrapolated to the additional p T bin have either negligible p T dependence in that region or are negligibly small compared to other systematic or statistical uncertainties.

Results
For the measurement of the ratios of the cross sections of the prompt ψ(2S), Υ(2S), and Υ(3S) states relative to their ground states, the systematic uncertainties in the yields, the ρ factor, the single-muon efficiencies, and the acceptance are the only ones considered. Uncertainties in the yields for the ratio of ψ(2S) and J/ψ cross sections are treated as uncorrelated, because their corresponding yields are determined from independent fits. In contrast, yield uncertainties are treated as correlated for the ratio of the Υ(nS) to Υ(1S) cross sections, as they are extracted from a combined fit to the three states, as shown in Fig. A.2 of Appendix A. The correlation factors are found to be negligible, and are on the order of 0.05. The same single-muon efficiencies are used for all the measured cross sections, therefore their uncertainties are treated as correlated in all the ratios. The systematic uncertainties in the ratios are determined by consistently varying the efficiencies in the numerator and the denominator by their uncertainties and recalculating the ratios. The resulting effect is less than 0.4%. The uncertainty in the integrated luminosity is fully correlated, and is not included in the ratios. Uncertainties in the ρ correction factor are treated as uncorrelated.
The statistical uncertainty in the ψ(2S) to J/ψ cross section ratio is more important than any systematic uncertainty except for the high-p T region, where the ρ factor uncertainty is the dominant one, reaching 28%. For the Υ(2S) to Υ(1S) and Υ(3S) to Υ(1S) cross section ratios, the uncertainty in the ρ factor dominates across the entire p T region, ranging from 3% to 12%.

Results
The measured double-differential cross sections times the dimuon branching fractions are presented in Fig. 1 as a function of p T , for four rapidity ranges in the case of the prompt J/ψ and ψ(2S) states, and two rapidity ranges for the Υ(nS). The top panels of Fig. 2 show the measured cross sections times branching fractions for the rapidity-integrated range |y| < 1.2. The presented results are obtained under the assumption of unpolarized production, which is very close to the polarization that was measured by CMS [30,31]. If the quarkonium states are fully polarized, the cross sections can change by up to 25%. The numerical values of the cross sections for all five quarkonium states in the chosen bins of p T and |y| in the unpolarized scenario are reported in Tables A.1-A.5 of Appendix A. Tables A.6-A.10 list the multiplicative scale factors needed to recalculate the cross sections in the three different polarization scenarios described in Section 3. The conversion to a new polarization scenario is achieved by multiplying the unpolarized cross section result in each (p T , |y|) bin by the corresponding scale factor.
The NLO NRQCD predictions [42,43] are in agreement with the measured cross sections times branching fractions within uncertainties, as shown in the top panels of Fig. 2. The ratios of the measured to predicted values are plotted in the middle panels of Fig. 2, where the vertical bars represent the experimental uncertainties. The shaded bands show the theoretical uncertainties stemming from the extraction of the LDMEs, renormalization scales, and the choice of c and b quark masses, added in quadrature with the uncertainties in the dimuon branching fractions [35]. The theory tends to underestimate (overestimate) the cross section for the J/ψ (ψ(2S)), while staying within the one-standard-deviation uncertainty band. The bottom panels of Fig. 2 show the ratios of the p T differential cross sections times branching fractions measured at √ s = 13 TeV and 7 TeV [25, 26] for |y| < 1.2. The 13 TeV cross sections of all five quarkonium states are factors of 1.5 to 3 larger than the corresponding 7 TeV cross sections, changing slowly as a function of dimuon p T . An increase of this order is expected from the evolution of the parton distribution functions.  Table A.11 of Appendix A.

Summary
The double-differential production cross sections of the J/ψ, ψ(2S), and Υ(nS) (n = 1, 2, 3) quarkonium states have been measured, using their dimuon decay mode, in pp collisions at √ s = 13 TeV with the CMS detector at the LHC. The production cross sections of all five S-wave states are presented in a single analysis. The measurement has been performed as a function of transverse momentum (p T ) in several bins of rapidity (y), covering a p T range 20-120 GeV for the J/ψ meson and 20-100 GeV for the remaining states. The cross sections integrated over |y| < 1.2 are also presented, and extend the p T reach to 150 and 130 GeV, respectively. Also presented are the ratios of cross sections measured at √ s = 13 (this analysis) and 7 TeV (from Refs. [25, 26]), as well as the cross sections of the prompt ψ(2S), Υ(2S), and Υ(3S) mesons relative to their ground states. These results will help in testing the underlying hypotheses of nonrelativistic quantum chromodynamics and in providing further input to constrain the theoretical parameters.   [11] LHCb Collaboration, "Measurement of the χ b (3P) mass and of the relative rate of χ b1 (1P) and χ b2 (1P) production", JHEP 10 (2014) 88, doi:10.1007/JHEP10(2014)088, arXiv:1409.1408.
[13] CMS Collaboration, "Measurement of the production cross section ratio σ(χ b2 (1P))/σ(χ b1 (1P)) in pp collisions at              Table A.1: Double-differential cross section times the dimuon branching fraction of the J/ψ meson for different ranges of p T , in bins of |y| and for the full |y| range, for the unpolarized decay hypothesis, with their statistical and systematic uncertainties in percent. The average p T value in each bin is also given. The global uncertainty in the integrated luminosity of 2.3% is not included in the systematic uncertainties.
B dσ 2 /dp T dy   Table A.3: Double-differential cross section times the dimuon branching fraction of the Υ(1S) meson for different ranges of p T , in bins of |y| and for the full |y| range, for the unpolarized decay hypothesis, with their statistical and systematic uncertainties in percent. The average p T value in each bin is also given. The global uncertainty in the integrated luminosity of 2.3% is not included in the systematic uncertainties.   Table A.5: Double-differential cross section times the dimuon branching fraction of the Υ(3S) meson for different ranges of p T , in bins of |y| and for the full |y| range, for the unpolarized decay hypothesis, with their statistical and systematic uncertainties in percent. The average p T value in each bin is also given. The global uncertainty in the integrated luminosity of 2.3% is not included in the systematic uncertainties.  Table A.6: Multiplicative scaling factors to obtain the J/ψ differential cross sections for different polarization scenarios (λ HX θ = +1, k, −1) from the unpolarized cross section measurements given in Table A Table A.9: Multiplicative scaling factors to obtain the Υ(2S) differential cross sections for different polarization scenarios (λ HX θ = +1, k, −1) from the unpolarized cross section measurements given in Table A.4. The parameter k corresponds to a linear interpolation of the CMS measured value of λ HX θ [31] as a function of p T for p T < 50 GeV. For p T > 50 GeV, where no measurements of λ HX θ exist, k is taken as the average of all the measured values of λ HX θ for p T < 50 GeV. p T |y| < 0.6 0.6 < |y| < 1.2 |y| < 1.2 [GeV] λ θ = +1 λ θ = k λ θ = −1 λ θ = +1 λ θ = k λ θ = −1 λ θ = +1 λ θ = k λ θ = −1 20-22  Table A.10: Multiplicative scaling factors to obtain the Υ(3S) differential cross sections for different polarization scenarios (λ HX θ = +1, k, −1) from the unpolarized cross section measurements given in Table A.5. The parameter k corresponds to a linear interpolation of the CMS measured value of λ HX θ [31] as a function of p T for p T < 50 GeV. For p T > 50 GeV, where no measurements of λ HX θ exist, k is taken as the average of all the measured values of λ HX θ for p T < 50 GeV, which are all consistent with a single value.