Measurement of nuclear modification factors of $\Upsilon$(1S), $\Upsilon$(2S), and $\Upsilon$(3S) mesons in PbPb collisions at $\sqrt{s_{_\mathrm{NN}}} =$ 5.02 TeV

The cross sections for $\Upsilon$(1S), $\Upsilon$(2S), and $\Upsilon$(3S) production in lead-lead (PbPb) and proton-proton (pp) collisions at $\sqrt{s_{_\mathrm{NN}}} =$ 5.02 TeV have been measured using the CMS detector at the LHC. The nuclear modification factors, R$_\mathrm{AA}$, derived from the PbPb-to-pp ratio of yields for each state, are studied as functions of meson rapidity and transverse momentum, as well as PbPb collision centrality. The yields of all three states are found to be significantly suppressed, and compatible with a sequential ordering of the suppression, R$_\mathrm{AA}$($\Upsilon$(1S)) $>$ R$_\mathrm{AA}$($\Upsilon$(2S)) $>$ R$_\mathrm{AA}$($\Upsilon$(3S)) . The suppression of $\Upsilon$(1S) is larger than that seen at $\sqrt{s_{_\mathrm{NN}}} =$ 2.76 TeV, although the two are compatible within uncertainties. The upper limit on the R$_\mathrm{AA}$ of $\Upsilon$(3S) integrated over $p_\mathrm{T}$ and rapidity is 0.094 at 95% confidence level, which is the strongest suppression observed for any hadron species in heavy ion collisions to date.


Introduction
The measurement of quarkonium production in heavy ion collisions is one of the most promising ways to study the properties of strongly interacting matter at high energy density and temperature. It has been predicted that in such an environment, a strongly interacting medium of deconfined quarks and gluons (the quarkgluon plasma, QGP) is formed [1,2]. Bottomonium states have been the subject of studies in heavy ion collisions for several reasons. Bottomonia are produced during the early stages of collisions via hard parton scattering. Their spectral functions are modified as a consequence of Debye screening of the heavy-quark potential at finite temperatures [3,4], as well as by thermal broadening of their widths due to interactions with gluons [5,6]. These in-medium effects have been studied in numerical simulations of quantum chromodynamics (QCD) on a space-time lattice, and captured as real and imaginary components of the heavy-quark potential [7]. One of the most remarkable signatures of these interactions with the medium is the sequential suppression of quarkonium states in heavy ion collisions compared to the production in proton-proton (pp) collisions, both in the charmonium (J/ψ , ψ(2S), χ c , etc.) and the bottomonium (ϒ(1S), ϒ(2S), ϒ(3S), χ b , etc.) families [8].
This scenario follows from the expectation that the suppression of quarkonia is stronger for states with smaller binding energy. ⋆ E-mail address: cms -publication -committee -chair @cern .ch.
The quarkonium yield can also increase in the presence of QGP, from the recombination of uncorrelated quarks [9][10][11][12]. However, recombination-like processes for bottomonia are expected to be negligible compared to the charmonium family [13][14][15], because these processes are driven by the number of heavy-quark pairs present in a single event, which is much smaller for beauty than for charm. The dissociation temperatures for the ϒ states, above which suppression occurs, are expected to be correlated with their binding energies, and are predicted to be T dissoc ≈ 2T c , 1.2 T c and 1T c for the ϒ(1S), ϒ(2S), and ϒ(3S) states, respectively, where T c is the critical temperature for deconfinement [16]. Therefore, measurements of the yields of each ϒ state can provide information about the thermal properties of the medium during its hot early phase.
Modifications of particle production in nucleus-nucleus (AA) collisions are quantified using the nuclear modification factor, R AA , which is the ratio of the yield measured in AA to that in pp collisions, scaled by the mean number of binary NN collisions. Comparisons of the bottomonium data with dynamical models incorporating the heavy-quark potential effects found in high-temperature lattice QCD are thus expected to extend our understanding of the nature of colour deconfinement in heavy ion collisions. Measurements of both the charmonium (J/ψ and ψ(2S))[ [17][18][19][20] and bottomonium (ϒ(1S), ϒ(2S), and ϒ(3S)) [21,22] families have been carried out at a nucleon-nucleon (NN) center-of-mass energy of laboration show strong suppression of J/ψ and ψ(2S) mesons [26,27], as well as of both ϒ(2S) and ϒ(3S) relative to the ϒ(1S) ground state [28]. The suppression of the excited ϒ(2S) relative to the ϒ(1S) ground state persists at very forward rapidity, 2.5 < y < 4 [ 29]. These measurements provide new constraints for theoretical models of the medium [9,11].
In this Letter, we report measurements of the differential cross sections and nuclear modification factors for ϒ(1S), ϒ(2S), and ϒ(3S) mesons using their decay into two oppositely charged muons in lead-lead (PbPb) and pp collisions at √ s NN = 5.02 TeV.
Results are presented as functions of the ϒ transverse momentum (p T ) and rapidity ( y), as well as PbPb collision centrality (i.e., the degree of overlap of the two lead nuclei). The data were collected with the CMS detector at the CERN LHC in 2015.

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 tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Muons are detected in the pseudorapidity interval of |η| < 2.4u s i n g gas-ionization chambers made of three technologies: drift tubes, cathode strip chambers, and resistiveplate chambers. These are embedded in the steel flux-return yoke of the solenoid. The silicon tracker is composed of pixel detectors followed by microstrip detectors. The p T of muons matched to tracks reconstructed in the silicon detector is measured with a resolution between 1% and 2% for typical muons used in this analysis [30]. In addition, CMS has extensive forward calorimetry, including two steel and quartz-fiber Cherenkov hadron forward (HF) calorimeters that cover the range of 2.9 < |η| < 5.2.
The HF calorimeters are segmented into towers and the granularity is η × φ = 0.175 × 0.175 radians. These are used in the present analysis to select PbPb collision events and to define their centrality class. Centrality, defined as the fraction of the total inelastic hadronic cross section with 0% representing collisions with the largest overlap of the two nuclei, is determined experimentally using the total energy in both HF calorimeters [31]. A more detailed description of the CMS detector, together with a definition of the coordinate system and the kinematic variables, can be found in Ref. [32].

Data selection and simulation samples
The ϒ mesons are identified using their dimuon decay channel.
In both pp and PbPb collisions, the dimuon events are selected by a fast hardware-based trigger system, which requires two muon candidates in a given bunch crossing with no explicit requirement on the muon momentum beyond the intrinsic selection due to the acceptance coverage of the CMS muon detectors. In pp collisions, this trigger registered an integral luminosity of 28.0 pb −1 . The PbPb data were taken with two triggers based on the same algorithm used for pp data. The first mode, designed to enhance the event count for muon pairs from peripheral events, added an additional selection that the collision centrality be in the 30-100% range. This trigger sampled the full integrated luminosity of 464 µb −1 . The second mode, using just the pp trigger alone, was prescaled during part of the data taking and therefore sampled a smaller effective integrated luminosity of 368 µb −1 . Data taken with this latter trigger were used to analyze the yields in the 0-30% and 0-100% centrality bins.
In order to keep hadronic collisions and reject beam-related background processes (beam-gas collisions and beam scraping events), an offline event selection is applied. Events are required to have at least one reconstructed primary vertex. In pp collisions at least 25% of the tracks have to pass a tight track-quality selection [33]. A filter on the compatibility of the silicon pixel detector cluster width and the vertex position is also applied [34]. The PbPb collision events have an additional requirement of the presence of at least three towers in the HF on both sides of the interaction point with an energy above 3G e V . The combined efficiency for this event selection, and the remaining contamination due to non-hadronic ultra-peripheral events which can raise the efficiency above 100%, is (99 ± 2)% [ 35,36]. The minimum-bias trigger requirement removes a negligible fraction of the events with a hard collision needed to produce ϒ mesons. We also studied a possible contamination from photoproduction processes in the peripheral region and found it to be negligible. Multiplecollision events (pileup) have a negligible effect on the measurement, since the average number of additional collisions per bunch crossing is approximately 0.9 for pp and much smaller for PbPb data.
Muons are selected in the kinematic range of p μ T > 4 GeV and |η μ | < 2.4, and are also required to be reconstructed using the combined information of the tracker and muon detectors (so-called "global muons" defined in Ref. [30]). To remove cosmic ray muons, the distance of the muon track from the closest primary vertex must be less than 20 cm in the beam direction and 3m m in the transverse direction. Pairs of oppositely charged muons are fitted with a common vertex constraint and kept if the fit χ 2 probability is larger than 1%. The studied dimuon kinematic range is limited to p μ + μ − T < 30 GeV and |y μ + μ − | < 2.4. Dimuons in this p T range comprise 99% of those passing all of the analysis selection criteria.
Simulated Monte Carlo (MC) ϒ events are used to calculate correction factors for all of the results presented, including the geometrical acceptance and reconstruction efficiency, as well as the trigger and offline selection efficiency. The samples are generated using pythia 8.209 [37]f o r the pp collisions and pythia 8.209 embedded in hydjet 1.9 for the PbPb events [38]. The PbPb simulation is tuned to reproduce the observed charged-particle multiplicity and p T spectrum in PbPb data. The CMS detector response is simulated using Geant4 [39]. Since the simulated p T spectrum of ϒ is not identical to the spectrum observed in data, an event-by-event weight is applied to the simulations in order to match the two distributions. The weight is given by a function fit to the ratio of data over MC p T spectra.

Signal extraction
The yields of ϒ mesons are extracted using unbinned maximum-likelihood fits to the invariant mass spectra, following the same procedure for pp and PbPb data. The signal of each ϒ state is modeled by a double Crystal-Ball (CB) function which is the sum of two CB functions [40]. This choice together with leaving the width parameter for the first CB free in the fit, is made in order to account for the different mass resolution in the barrel compared to the endcap region of the detector. A parameter relates the widths of the two CB functions, the second one being constrained to be narrower. The mass and the two radiative-tail parameters of both CB functions for a given state are kept the same, as these are not affected by the detector resolution. The mass parameter of the ground state is left free to allow for possible shifts in the absolute momentum calibration of the reconstructed tracks. For the excited states (ϒ(2S) and ϒ(3S)), the yields can vary while all other fit parameters are fixed to be identical to those for the ground state except for the mean and width which are fixed to values found by multiplying those for ϒ(1S) by the ratio of the published masses of the states [41]. In the pp data fits, the two radiative-tail parameters and the parameter for the ratio of the two widths are allowed to vary within a Gaussian probability density function (PDF). The mean and the width of the constraining Gaussian function represent the average and its uncertainty, respectively, from the fits in all the rapidity bins of the analysis with no fixed parameters. In the PbPb fits, in addition, the parameter for the fraction of the two CB functions is also constrained. In this case, the mean and the width of the constrained parameters represent the corresponding parameter values and their uncertainties from the pp fits for each kinematic region. The background PDF is an error function multiplied by an exponential, with the yield, the error function's two parameters, and the decay parameter of the exponential all allowed to vary in the final fit. For bins with p T > 6G e V , an exponential without the error function provides the best fit, and was used for the nominal result. Fig. 1 shows the dimuon invariant mass distributions in pp and PbPb collisions along with the fits using the model described above, for the kinematic range p

Corrections
In order to obtain the normalized cross sections, the yields extracted from the fits to the dimuon invariant mass spectra are corrected for acceptance and efficiency, and scaled by the integrated luminosity. The acceptance corresponds to the fraction of dimuon events originating from ϒ mesons within the kinematic range of the analysis. The acceptance values for the considered kinematic region are 22.5% (ϒ(1S)), 27.8% (ϒ(2S)), and 31.0% (ϒ(3S)) for PbPb collisions and differ by <1% from the corresponding pp data values, with the small difference being due to a small residual difference in the kinematic spectra after weighting the MC to data.
The dimuon efficiency is defined as the probability that a muon pair within the acceptance is reconstructed offline, satisfies the trigger condition, and passes the analysis quality criteria described in Section 3. The dimuon efficiency is calculated using MC. The individual components of the efficiency (track reconstruction, muon identification and selection, and triggering) are also measured using single muons from J/ψ meson decays in both simulated and collision data, with the tag-and-probe (T&P) method [30]. For the muons used in this analysis, data and MC efficiencies are seen to differ only in the case of the trigger efficiency, and there only by 1%. For this case, scaling factors (SF), calculated as the ratio of data over simulated efficiencies as function of p μ T and η μ , are applied to each dimuon on an event-by-event basis. The other components of the T&P efficiency are used only for the estimation of systematic uncertainties. The average efficiencies integrated over the full kinematic range are 73.5% (ϒ(1S)), 74.4% (ϒ(2S)), and 75.0% (ϒ(3S)) in PbPb collisions, and they are 8-9% higher for pp collisions.
The integrated luminosity of 28.0 pb −1 with an uncertainty of 2.3% [42]i s used to normalize the yields for pp data. For PbPb collisions, the number of minimum bias collision events sampled by the trigger (N MB ), together with the average nuclear overlap function (T AA ), are used for the normalization. The overlap function T AA is given by the number of binary NN collisions divided by the inelastic NN cross section, and can be interpreted as the NN-equivalent integrated luminosity per heavy ion collision. Values of T AA are calculated with a Glauber model MC simulation [43,44], which is also used to obtain the average number of participating nucleons, N part . This latter number is highly correlated with the impact parameter of the collision, and is used as the abscissa when plotting results as a function of PbPb collision centrality.

Systematic uncertainties
Point-to-point systematic uncertainties arise from the choices of signal and background PDFs and of the central value in the fit constraints, as well as from acceptance and efficiency corrections. Larger relative uncertainties are obtained when the background level is higher (at lower p T or more forward y regions), and, in particular for the ϒ(3S), when the absolute yield is small.
The uncertainty from the choice of signal model is estimated by fitting the data using a single CB function in combination with a Gaussian function instead of a double CB function. The uncertainties are determined by calculating the difference between the yield obtained with the alternative model compared to the nominal one. For the PbPb (pp) yields, the differences are in the range of 1-7% (0.1-4.6%) for the ϒ(1S), 2-19% (0.1-1.3%) for the ϒ(2S), and 5-78% (0.7-7%) for the ϒ(3S) mesons.
The systematic uncertainty from the choice of the central value in the fit constraints is estimated by using instead of the average parameter values from the pp fits, the values in each pp analysis bin when all parameters were left floating. The differences in the PbPb (pp) signal yields, typically below 4% (4.5%) for the ϒ(1S), below 8% (3%) for the ϒ(2S), and 45% (2%) for the ϒ(3S), are quoted as a systematic uncertainty.
The systematic uncertainty due to the choice of background model is estimated using two alternative background functions. One is in the form of a fourth-order polynomial function and the other is an exponential plus an additional linear function. The maximal deviations of the PbPb (pp) yield between these two models compared to the nominal are quoted as the uncertainty and are typically in the range of 1-6% (1-5%) for the ϒ(1S), 2-23% (2-4%) for the ϒ(2S), and 5-200% (3-5%) for the ϒ(3S) mesons.
For the estimation of systematic uncertainties due to acceptance and efficiency corrections, the source of uncertainty is the imperfect knowledge of the simulated p T distribution shape. To take this source into account, the function used to weight the MC p T spectra event-by-event is modified within its fit uncertainty. The acceptance and efficiency obtained from the simulated p T distribution are compared with and without the variation of the function, with the difference between the two used as an estimate of the systematic uncertainty. In addition, there is a systematic uncertainty for the efficiency in the T&P correction arising from the uncertainty in the SFs of the single-muon efficiency. The systematic uncertainties of the SFs are taken into account for trigger, tracking, and muon identification. The uncertainties in the single muon efficiencies are propagated to the dimuon efficiency values to estimate the systematic uncertainty from this source. The statistical uncertainty inherent in the data set used for the T&P studies is also considered as an additional component of the systematic uncertainty in the corrected yields. The PbPb (pp) systematic uncertainties are in the range of 3.5-6.4% (2.6-3.9%) in the case of the total efficiency correction and in the range of 0.1-3.0% (0.1-0.8%) for the acceptance correction.
Finally, several sources of correlated uncertainties (i.e., global uncertainties common to all points) are considered: for the pp dataset from the pp integrated luminosity, and for the PbPb dataset from the T AA and the N MB estimations. The uncertainty on the integrated luminosity measurement for the pp dataset is 2.3% [42]. The uncertainty for N MB in PbPb collisions is 2%, which accounts for the inefficiency of trigger and event selection. For the R AA calculation, T AA uncertainties ( Table 2 in Appendix A) are estimated by varying the Glauber model parameters within their uncertainties ( Table 1 in Appendix A) [ 36]. The total combined uncertainty is calculated by adding the results from the various sources in quadrature. The global uncertainty for the differential cross section results arises from the integrated luminosity in pp collisions and N MB in PbPb collisions. For the R AA results, the global uncertainty combines the uncertainties from T AA , pp luminosity, and PbPb N MB for the bins integrated over centrality. For the centrality dependent R AA results, the uncertainty from T AA is included bin-by-bin, while the total uncertainty from the pp measurement is included in the global uncertainty. Using the updated uncertainties of the Glauber model parameters in Ref. [45], instead of those from Ref.
[36], would reduce the T AA uncertain-ties by 0.1-1.1% and the total systematic uncertainties for R AA by less than 0.7% (with the largest change for the 70-100% centrality bin). However, in order to allow direct comparisons to previous results [27,36,46,47], these updated parameters are not used in this analysis. The bin migration effect due to the momentum resolution is negligible for the kinematic range of this measurement.

Results
The ϒ cross sections and values of R AA are measured in sev-

Differential cross sections in pp and PbPb collisions
The differential production cross section of ϒ mesons decaying in the dimuon channel in pp collisions is given by The branching fraction for the decay ϒ → μ + μ − is denoted by B.
The quantity N corresponds to the extracted yield of ϒ mesons in a given (p T , y) bin, (A ε) represents the average acceptance and efficiency in the given bin, L int is the integrated luminosity, and p T and y are the widths of the given bin. For PbPb data, L int is replaced by (N MB T AA ), as explained in Section 4.2, to compare the pp and PbPb data under the hypothesis of binary-collision scaling.

Nuclear modification factor R AA
The nuclear modification factor is derived from the pp cross sections and PbPb normalized yields as where T AA is the average value of T AA computed in each centrality bin. The quantities N AA and σ pp refer to the normalized yield of ϒ mesons in PbPb collisions corrected by acceptance and efficiency, and the pp cross section for a given kinematic range, respectively.  range explored here. The kinematic dependence of R AA is useful to constrain models of ϒ meson suppression in a deconfined medium [9].
The dependence of R AA on PbPb collision centrality, as quantified using the average N part , is depicted in Fig. 5. The strong suppression of the ϒ(3S) meson is observed in both centrality bins studied, 0-30% and 30-100%. The R AA decreases with increasing centrality in the case of the ϒ(1S) and ϒ(2S) mesons. A hint of this centrality dependence of R AA for ϒ(2S) was first seen in data at √ s NN = 2.76 TeV [21] and is now confirmed using the larger data sample at 5.02 TeV. Fig. 6 shows a comparison between the measured R AA for ϒ(1S) and ϒ(2S) mesons and two models of bottomonium suppression from Krouppa and Strickland [9], and from Du, He, and Rapp [12]. Both models incorporate color-screening effects on the bottomonium family and feed-down contributions from decays of heavier quarkonia. No regeneration in QGP or cold nuclear matter effects are considered by the first model, but are included in the sec-     while no significant dependence on p T or y is found in the measured region. The suppression of ϒ(1S) is larger than that seen at √ s NN = 2.76 TeV, although the two are compatible within uncertainties. The R AA of the ϒ(3S) state is measured to be below 0.096 at 95% confidence level, making this the strongest suppression observed for a quarkonium state in heavy ion collisions to date. heavy ion nuclear physics. Without his work, our understanding of the modifications of particle production in heavy ion collisions, of which this paper is only one of very many, would not have been possible. We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses.

Appendix A. Glauber model values
Centrality variables computed using a Glauber model [44]a r e summarized in Tables 1 and 2