Nuclear modification of $\Upsilon$ states in pPb collisions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV

Production cross sections of $\Upsilon$(1S), $\Upsilon$(2S), and $\Upsilon$(3S) states decaying into $\mu^+\mu^-$ in proton-lead (pPb) collisions are reported using data collected by the CMS experiment at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV. A comparison is made with corresponding cross sections obtained with pp data measured at the same collision energy and scaled by the Pb nucleus mass number. The nuclear modification factor for $\Upsilon$(1S) is found to be $R_\mathrm{pPb}(\Upsilon(1S))$ = 0.806 $\pm$ 0.024 (stat) $\pm$ 0.059 (syst). Similar results for the excited states indicate a sequential suppression pattern, such that $R_\mathrm{pPb}(\Upsilon(1S))$ $\gt$ $R_\mathrm{pPb}(\Upsilon(2S))$ $\gt$ $R_\mathrm{pPb}(\Upsilon(3S))$. The suppression is much less pronounced in pPb than in PbPb collisions, and independent of transverse momentum $p_\mathrm{T}^\Upsilon$ and center-of-mass rapidity $y_\mathrm{CM}^\Upsilon$ of the individual $\Upsilon$ state in the studied range $p_\mathrm{T}^\Upsilon$ $\lt$ 30 GeV$/c$ and $\vert y_\mathrm{CM}^\Upsilon\vert$ $\lt$ 1.93. Models that incorporate sequential suppression of bottomonia in pPb collisions are in better agreement with the data than those which only assume initial-state modifications.


Introduction
Properties of the color-deconfined quark-gluon plasma (QGP) created in high-energy collisions of heavy nuclei can be studied using heavy-quark resonances produced by initial hard scatterings [1][2][3][4][5][6]. Yields of various quarkonium states, which have a short formation time in their rest frames and can typically escape the QGP before they decay, encode information on the evolution of the plasma starting from its early stages [1,2,4,[6][7][8][9]. Debye screening and gluodissociation [10][11][12][13][14] in the QGP produced in lead-lead (PbPb) collisions are understood to modify yields of quarkonium states hierarchically, according to their binding energies. Each state dissociates when a high enough temperature is reached in the QGP [4-6, 9, 15]. To interpret the quarkonium-state suppression patterns observed in heavy ion collisions as signals of color deconfinement in the hot plasma, it is essential to understand initial state and final state "cold nuclear matter" (CNM) effects. In this context, initial state refers to the partons in the relevant quantum chromodynamics process that stem from the colliding proton (p) or nucleus and scatter to produce a heavy quark pair but before it hadronizes into a quarkonium state. Examples of CNM effects that have been discussed in pA collisions include shadowing of the parton distribution functions in the nucleus (initial state) [16], energy loss in the nucleus (initial and final states) [17], and interactions with hadronic comovers (final state) [18]. For a recent review, see Ref. [7]. Traditionally, all modifications observed in pPb collisions were assumed to be due to CNM effects. However, it is worth noting that this assumption has been questioned given recent evidence of collective behavior in pp and pPb collisions with the highest amount of emitted particles, further referred to as high-activity [19][20][21][22][23][24][25]. This might be explained by assuming the formation of a QGP-like medium [26].
Bottomonia serve as particularly powerful probes for studying the QGP, since their high masses require that their production be dominated by initial hard scattering of partons in the collisions [7,[27][28][29][30]. When compared to charmonia, their yields are considerably less modified by regeneration or recombination in the QGP [30][31][32]. Measurements by the CMS experiment showing sequential modification of Υ(nS) (where n = 1, 2, 3) decaying via the dimuon channel in PbPb compared with pp collisions at a nucleon-nucleon center-of-mass energy of √ s NN = 2.76 TeV [33, 34] and 5.02 TeV [35,36] were used to infer model-dependent [32, 37] QGP temperatures. This effect is consistent with models that incorporate sequential suppression due to color screening [5,8,30]. Similar measurements of Υ(nS) production in pPb collisions can help to disentangle hot and cold nuclear matter effects and to investigate various CNM mechanisms.
Nuclear modification factors R pPb are ratios of particle production cross sections in pPb collisions over the corresponding cross sections in pp collisions scaled to account for the number of nucleons in the Pb nucleus. The R pPb values quantify the modification of hard probe production in pPb collisions due to the nuclear environment created by a single lead nucleus in the initial state. In this analysis, these factors are determined for Υ(nS) under the assumption that the cross sections scale as σ pPb = Aσ pp , where A is the mass number of Pb. With this assumption, also known as the A-scaling hypothesis, values of R pPb different from unity indicate modifications that go beyond simple superposition of binary nucleon-nucleon collisions. These R pPb values, together with measurements of the nuclear modification factors R AA in PbPb collisions [36], can be used to investigate the relative contributions of hot and cold nuclear matter effects.
Since pPb collisions create an imbalance of nuclear matter in the proton-going (forward rapidity) and lead-going (backward rapidity) directions, they can be used to investigate differences in CNM effects in these regions of varying nuclear matter density within the same collision sys-tem. In the charmonium sector, CMS has found hints of differences in the level of suppression between the excited and ground state in the lead-going region [38,39]. One CNM modification mechanism that relies on the abundance of nuclear matter is dissociation by interaction with comoving particles, where the cross section of interaction increases with particle multiplicity in the rapidity region of the produced Υ meson [18,40]. This is quantified by measuring the forward-backward production ratios R FB of Υ states in pPb collisions.
The LHCb [41] and ALICE [42] Collaborations reported measurements of the Υ(nS)/Υ(1S) yield ratios (LHCb for n = 2 and 3; ALICE for n = 2), along with R pPb and R FB for Υ(1S) in pPb collisions at √ s NN = 5.02 TeV using Υ mesons detected in the forward rapidity region. In those studies, the proton reference was obtained by interpolating results from event samples collected at other collision energies, i.e., 2.76, 7, and 8 TeV. In the midrapidity region, the ATLAS Collaboration studied bottomonia in pPb collisions using same-energy pp reference data [43], reporting Υ(nS)/Υ(1S) (for n = 2 and 3), as well as Υ(1S) yields self-normalized to their activity-integrated values, and R pPb (Υ(1S)). The CMS Collaboration previously reported the Υ(nS)/Υ(1S) (for n = 2 and 3) yield ratios versus event activity in  16 TeV for direct comparison. These bottomonium measurements in pPb have focused on the ground state and indicate that the level of suppression is consistent with that expected from shadowing calculations, but they provide little information on the behavior of the excited states.
In this Letter, we analyze pPb and pp collision data from the CERN LHC collected with the CMS detector at the same nucleon-nucleon center-of-mass (CM) energy of √ s NN = 5.02 TeV. The yields of Υ(nS) mesons are measured using their decay to two muons. By comparing the yields measured in the two colliding systems, the R pPb and R FB factors are determined including all bottomonium states for the first time. Because models that incorporate final-state CNM effects are the only ones to predict different modifications for the excited states, these measurements for the ordering of excited state R pPb values may reveal these types of final-state mechanisms. Ordered suppression could arise from various causes -e.g., the size of the states, their cross section with potential comovers, or their binding energy. These results are compared with measurements of the Υ(nS) nuclear modification factors R AA in PbPb collisions [36] using PbPb data also collected at 5.02 TeV with the CMS detector, allowing a model-dependent comparison of bottomonia in hot and cold nuclear matter.

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. Muons are detected in the range |η lab | < 2.4 in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. In the barrel region |η lab | < 1.2 muon detection planes are based on drift tube technology, while the endcap region 0.9 < |η lab | < 2.4 uses cathode strip chambers. Resistive plate chambers provide additional muon detection capability in the range |η lab | < 1.6. Matching muons to tracks measured in the silicon tracker leads to a relative transverse momentum p T resolution on the order of 1% for a typical muon used in this analysis [48]. In addition, two steel and quartz-fiber hadron forward calorimeters cover the range 2.9 < |η lab | < 5.2. A 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. [49].
A two-tiered system is used to select collision events of interest from the detector. The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 µs [50]. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage [51].

Data selection and simulated samples
The events used for this analysis are selected using the trigger systems described above, requiring two muon candidates in the muon detectors with no explicit cuts in muon transverse momentum, p µ T , or muon pseudorapidity measured in the laboratory, η µ lab . The event samples used in this analysis correspond to integrated luminosities of 28.0 ± 0.6 pb −1 and 34.6 ± 1.2 nb −1 for pp [52] and pPb [53] collisions, respectively. The uncertainties in the integrated luminosity determination are considered as a global uncertainty in all results. All recorded pPb events are required to have an energy deposit above 3 GeV in the hadron forward calorimeters on each side of the interaction point in order to suppress background from ultra-peripheral collisions and beam-gas events, while having a high efficiency for the selection of beam-beam hadronic collisions.
In the case of pPb collisions, the value of the integrated luminosity represents the combined luminosity of collisions with proton and lead beams traveling in either direction. While in the symmetric pp and PbPb collision systems the CM and laboratory (lab) reference frames coincide, in the case of pPb collisions the difference between the energy-per-nucleon of the two beams induces a shift between the two frames. For pPb collisions at √ s NN = 5.02 TeV, the rapidity y is shifted in the CM frame by δy = 0.465 compared to the lab frame. The rapidity range of the reconstructed dimuons in the lab frame |y µµ lab | < 2.4 corresponds to a CM frame rapidity range of either −2.87 < y µµ CM < 1.93 (Pbp) or −1.93 < y µµ CM < 2.87 (pPb), depending on the direction of the proton beam. In order to minimize the influence of asymmetric detector conditions, data are taken with both beam directions and then combined by inverting the rapidity of one of the datasets.
For both pp and pPb data, we select events with muon candidates in the kinematic range p µ T > 4 GeV/c, |η µ lab | < 2.4. The muon tracks are required to have at least 6 hits in the silicon tracker, at least one hit in the silicon pixel detector, and match with at least one segment in any detection plane of the muon system. The distance of the track from the closest primary vertex [54] must be less than 20 cm in the longitudinal direction and 0.3 cm in the transverse direction. When forming a muon pair, each of the two muons is required to match the hardware trigger that prompted recording of the event and to originate from a common vertex with a χ 2 probability larger than 1%, as obtained by a Kalman vertex filter algorithm [55]. For pPb data, an additional filter is used to remove events that contain multiple interactions per bunch crossing (pileup) [38]. This filter reduces the fraction of pileup events from 3% to less than 0.2%, and reduces the effective luminosity of pPb collisions by 4.1% compared to the numbers noted above.
Dedicated Monte Carlo (MC) simulations of collision data are used to validate fitting techniques and to correct the extracted Υ(nS) yields for losses due to finite detector acceptance and efficiency. Simulated samples are independently generated for the Υ(1S), Υ(2S), and Υ(3S) mesons, in pp collisions using PYTHIA 8.209 [56], assuming no polarization based on measurements at the LHC [57, 58]. To simulate pPb collisions, the rapidities of all particles in the generated pp events are boosted by δy = 0.465 in the Pb-going direction to mimic the y CM shift in data. The CMS detector response is simulated using GEANT4 [59]. The reconstructed p Υ T distributions of the simulated Υ states are weighted using a fit to the ratio of the p Υ T spectra in data and simulation. The rapidity distributions in simulation are consistent with those in data. Figure 1 shows the invariant mass distributions of opposite-sign muon pairs for pp (left) and pPb (right) collisions. The dimuon data are integrated in the dimuon range p µµ T < 30 GeV/c and |y µµ CM | < 1.93. The yields of the Υ states, uncorrected for detector acceptance and efficiency, are obtained via unbinned maximum-likelihood fits to the invariant mass spectra, shown as solid blue lines. A dashed red line is used in Fig. 1 (right) to depict the expected Υ(1S), Υ(2S), and Υ(3S) yields under the R pPb = 1 hypothesis, obtained by scaling the signal shape of each state by the inverse of its finally measured R pPb value (including the ratio of the efficiencies corresponding to pp and pPb collisions). This comparison illustrates that the Υ(nS) yields are suppressed in pPb relative to pp collisions in the integrated kinematic region. We bin the data in the dimuon kinematic variables p µµ T and y µµ CM , as well as in event activity variables which we discuss below. Quarkonium peaks can be modeled by a Crystal Ball (CB) function [60], whose low-mass power-law tail accounts for dimuons that undergo bremsstrahlung radiation in the detector material as well as final-state radiation. We model the shape of each Υ state with a sum of two CB functions. A parameter representing the relative CB widths is left free in the fit, to accommodate muons with different momentum resolutions (depending on their η µ lab ). The relative contributions of the two CBs are also allowed to vary, both in the kinematic and event activity variables.

Signal extraction
To eliminate unnecessary degrees-of-freedom in the fits, the relative widths and relative contributions of the two CB functions are constrained to be the same for all three Υ states, consistent with fits to simulated samples. Furthermore, parameters governing the shape of the radiative tail are constrained to be the same for all six CB functions in the fit. The mass parameter of the Υ(1S) is left free to account for possible systematic shifts in the momentum scale of the reconstructed tracks. The final value of this parameter is consistent between fits to pp and pPb data. Since any changes in the momentum scale should affect all measured Υ(nS) similarly, we constrain the masses of the excited states such that their ratio matches the Particle Data Group (PDG) world-average values [61] as follows: (m(nS)/m(1S)) fit = (m(nS)/m(1S)) PDG . Similarly, the CB widths are also scaled by the ratio of the PDG mass values.
The parameters governing the tail shapes and ratio of CB widths are found to be correlated across kinematic bins. Rather than allowing the parameters to be completely unconstrained, we instead allow them to vary around their mean values within an interval estimated from a set of preliminary fits. The deviation of each parameter is translated into a Gaussian probability that is multiplied with the fit likelihood. The width of each Gaussian function is set to the RMS value of the corresponding parameter in the preliminary fits. In the case of pPb collisions, the central value of the parameter determining the relative contributions of the two CB functions is constrained in this manner as well.
As a result of the large number of free parameters in the preliminary fits, it is possible for parameters to converge to different values in repeated fits. By fitting the data across all the analysis bins, we find certain parameter values to be normally distributed. Each normally distributed parameter is first restricted to its mean value across the preliminary fits, in order to enable the rest of the parameters to converge consistently across the bins. We take an iterative approach to this constraining technique to avoid biasing the final parameter values. The mean values of parameters from preliminary fits are obtained separately in different rapidity regions to allow for differences in Υ meson reconstruction resolution in the barrel and end-cap regions of the detector, where muons pass through different amounts of material and are detected using different technologies.
The background is modeled with a shifted and scaled error function multiplied by an exponential. The exponential function models the dominant combinatorial background, which falls with increasing dimuon invariant mass according to a statistical phase space factor. The use of an error function is motivated by the effect of the p µ T > 4 GeV/c selection applied to single muons, which produces a hump-like feature in the combinatorial background at low invariant masses and at low dimuon p µµ T . For dimuon p µµ T > 6 GeV/c, this feature moves to lower invariant masses outside the fit region and we model the background solely with an exponential function.

Acceptance and efficiency corrections
The Υ(nS) yields that are extracted using fits to the invariant mass spectra are corrected to account for geometric limitations of the detector and inefficiencies of the online and offline selection algorithms. The dedicated MC simulations of Υ(nS) decays are used to determine the acceptance, which is the fraction of generated Υ mesons in a given kinematic region that decay to muons satisfying the kinematic requirements applied in this analysis.
The efficiency of dimuon reconstruction, event triggering, and muon identification are studied using dedicated MC simulations of Υ(nS) decays, after they have undergone full detector response simulation. The dimuon efficiency is determined as the fraction of generated Υ mesons in simulation that are identified as such, having satisfied all the same conditions that are required of muon pairs in collision data. Since pure PYTHIA-based MC samples are used for pPb collisions, we verify that the efficiency correction does not exhibit any dependence on multiplicity. This was also found in the related study of charmonium states reconstructed via muon pairs in pPb collisions [38].
Additional corrections are estimated to compensate for possible discrepancies between simulation and data efficiencies. To estimate such discrepancies, muon triggering, track reconstruction, and identification efficiencies are measured using single muons from prompt J/ψ meson decays in both simulation and data, as described in Ref. [48]. The ratios of the single-muon detection efficiencies between J/ψ data and simulation are estimated. These ratios differ significantly from unity in the case of muon triggering and identification in pPb collisions and muon triggering and track reconstruction in pp collisions. These ratios are used to correct the simulation-based efficiencies. For the bulk of muons in this analysis this correction to the efficiencies is small (∼1%, but for some regions of phase space it can grow to at most 10%).

Systematic uncertainties
Typical ranges of total and individual sources of systematic uncertainties in R pPb and R FB for all three Υ(nS) states are tabulated in Table 1. Two important sources of systematic uncertainty in the Υ(nS) yields originate from an incomplete knowledge of the signal and background shapes. The signal shape systematic uncertainty is estimated using an alternative fit model of the signal, consisting of a single CB function in combination with a Gaussian function. This alternative fit model provides a comparable goodness-of-fit to the nominal one. For the background uncertainty estimation, a similar method of recalculating yields using an alternative fit model for the background is used. Because the background shape evolves with p µµ T , the model was varied in different kinematic regions. In higher p µµ T regions, a power law is used as the alternative background fit model. In lower p µµ T regions, the background model is constructed from a linear combination of four invariant-mass fits to four p µµ T subintervals of a MC simulation of dimuon decays.
When estimating systematic uncertainties in the Υ(nS) yields using nominal and alternative models for signal and background distributions, we employ a method that helps to reduce the contribution of statistical fluctuations. We perform pseudo-experiments where we generate a set of invariant mass distributions by MC sampling the shape fitted to the dimuon invariant mass spectra in each analysis bin. Each generated invariant-mass distribution is fitted separately with the nominal and alternative signal models, using the nominal background model in both cases. The systematic uncertainty is evaluated as the mean of absolute values of relative differences between yields extracted using the nominal and alternative models. Similarly, additional pseudo-experiments are performed to estimate the uncertainty associated to the choice of the background model.
The procedure for constraining the parameters of the signal model introduces another source of systematic uncertainty. In order to estimate the systematic uncertainty on the yields we perform a set of preliminary fits to the dedicated Υ(nS) MC simulations, in which the parameter phase space is iteratively reduced in the same way as for data. For each parameter, we use the mean value obtained from the last iteration of MC fits as an alternative value for the mean of the Gaussian probability function used to constrain the parameters when fitting to data. In the case of pPb, the value of the parameter determining the relative contributions of the two CB functions in the free-parameter fit to MC is used as its alternative value. We compare the yields from the nominal fits with those extracted using the alternative mean of the constraining Gaussians and calculate the deviations of the yields for each constrained parameter. The largest of these deviations is assigned as the systematic uncertainty.
Systematic uncertainties in the acceptance corrections are estimated by varying the parameters of the fit used to weight MC p µµ T spectra within uncertainties, and recording the largest produced deviation. For the R pPb measurement, the Y states are assumed to be unpolarized in both pp and pPb collisions. We also assume that if the production mechanism of Υ(nS) leads to a different polarization, then it remains the same for both pp and pPb collisions. As a consequence, the acceptance ratio under a different polarization scenario will remain equal to unity. Therefore the uncertainty in the polarization would not affect the nuclear modification factor.
Systematic uncertainties in efficiency corrections are estimated by combining two sources in quadrature. The first is the uncertainty in weighting MC p µµ T spectra, which is estimated by determining the efficiency using the MC samples with and without weighting. The second source of uncertainty arises from the data-to-simulation ratio of single-muon detection efficiencies that are used to correct the purely simulation-based dimuon efficiencies. The systematic uncertainty in the efficiency of each stage of muon detection (triggering, tracking, muon identification) is studied by varying the selection criterion for that stage, while the statistical uncertainty is determined by repeating such variations one hundred times and estimating the standard deviation.
The total systematic uncertainty from uncorrelated sources is obtained by combining the uncertainties in quadrature in the signal and background extractions, as well as in the acceptance and efficiency corrections. The combined systematic uncertainty in the results increases slightly with increasing |y Υ CM | and with decreasing p

Results
The product of the branching fraction of Υ(nS) to muon pairs, B(Υ(nS) → µ + µ − ), and the double-differential production cross section, d 2 σ/dp Υ T dy Υ CM , is obtained as where N

Υ(nS)
Fit is the yield of Υ(nS) mesons extracted from the fit in a given (p µµ T , y µµ CM ) bin, a is the dimuon acceptance correction, ε is the dimuon efficiency correction, and L int is the integrated luminosity. Figure 2 Fig. 2 (lower row).
The quantity R pPb is calculated as where A = 208 is the mass number of the Pb nucleus. Results for R pPb are shown in Fig. 3 as functions of p Υ T and y Υ CM . We observe that all three Υ(nS) states are suppressed in pPb relative to pp collisions throughout the kinematic region explored, suggesting modification by CNM effects in pPb collisions. Similar to the PbPb case [36], the level of suppression for each Υ state in pPb collisions is consistent with a constant value in the kinematic region studied, although the level of suppression seen in PbPb is much stronger. The ATLAS Collaboration reported an increasing R pPb with p Υ T for Υ(1S) [43] in a similar midrapidity region as in CMS. The CMS data is consistent with no dependence, but the overall p Υ T dependence of the R pPb (Υ(1S)) in the two experiments is consistent within uncertainties. Moreover, our data shows no y Υ CM dependence, which is consistent with the ATLAS result.
In the charmonium sector, the CMS Collaboration found hints of an ordered suppression pattern. The R pPb of ψ(2S) was found to be smaller than that of J/ψ [38] in pPb collisions at √ s NN = 5.02 TeV for backward rapidity and p J/ψ T < 10 GeV/c [39]. The results presented here suggest a similar ordered suppression of the Υ states in the backward rapidity region as well as across the entire p Υ T region studied. The measured R pPb (Υ(1S)) is systematically larger than that of Υ(2S), which in turn is systematically larger than the R pPb (Υ(3S)), suggesting different levels of modification to the three states by final-state effects in these regions. In the forward rapidity region, the measured R pPb of the three states appear more mutually consistent.
We further compare the y Υ CM dependence of the measured R pPb to predictions from three CNM models: shadowing, energy loss, and comover interaction. The shadowing calculations incorporate next-to-leading order nuclear modifications of the PDFs (nPDFs) [16], according to  EPS09 [62]. Predictions using coherent energy loss [17] are made with and without using EPS09. Since they affect the quarkonium system before hadronization, both shadowing and energy loss are initial-state effects. Finally, predictions using the comover interaction model (CIM) [18] are provided with two different leading-order nPDF calculations: EPS09 and nCTEQ15 [63]. Since the comovers in the CIM interact with the quarkonium system after hadronization, it is deemed to be a final-state CNM effect.
In Fig. 4, the measured R pPb (Υ(1S)) is compared to predictions from shadowing [16] (left) and predictions using energy loss only and energy loss with shadowing [17] (right). The uncertainty in the models comes from the nPDFs. The combined energy loss with shadowing model is in better agreement with our data, but given the current uncertainties in theory and experiment, the models using only shadowing or energy loss cannot be ruled out. The shadowing model, which includes only initial-state effects, predicts equal modification of all bottomonium states and therefore is incompatible with the Υ(2S) and Υ(3S) data.
In contrast to shadowing and energy-loss models, the CIM predicts different degrees of modification for the Υ(1S), Υ(2S), and Υ(3S) states [18, 40], since higher excited states have a larger size and hence increased comover interactions. In addition, comover modification of quarkonium states is expected to be stronger in regions where the comover densities are larger, such as in the nucleus-going direction in asymmetric proton-nucleus collisions and in regions of higher event activity. Figure 5 shows comparisons of predicted R pPb in the CIM [18], including shadowing corrections from both nCTEQ15 and EPS09 nPDFs, with the measured R pPb for Υ(1S) (upper left), Υ(2S) (upper right), and Υ(3S) (lower). The CIM R pPb predictions show similar ordered suppression to that found in the data, an effect missing in models with only initial-state effects.
By comparing the R pPb (Υ(nS)) in the forward (proton-going) and backward (lead-going) directions, we can investigate the dependence of bottomonium suppression on the amount of nuclear matter present. Figure 6 shows the R pPb of Υ(nS) states for −1.93 < y Υ CM < 0 and We study the forward-backward production ratio of Υ mesons in pPb collisions defined as follows, where y Υ CM is positive. We measure event activity near the measured Υ meson using the number of reconstructed tracks, N tracks , in the region |η lab | < 2.4 (a detailed discussion of the eventactivity variables can be found in Ref. [44]). To measure event activity further from the Υ meson, we use the sum of deposited transverse energy E T in 4 < |η lab | < 5.2. Figure 7 shows the R FB as a function N tracks (left), and E T (right). The uncorrected mean values of the event activity variables in minimum bias pPb collisions are N tracks = 41 and E T = 14.7 GeV. The measured R FB remains consistent with unity at all levels of event activity for all three Υ states. This observation is independent of the η region used to measure event activity. The ALICE Collaboration determined a value of R FB consistent with unity for Υ(1S) for integrated event activity for Υ mesons in the forward (2.03 < y    Figure 8 shows the integrated R pPb of Υ states as well as the R AA observed in PbPb collisions [36] at √ s NN = 5.02 TeV. The 95% confidence level upper limit on the Υ(3S) R AA is depicted using an arrow. The data indicate an ordering of nuclear modification for the Υ family with R pPb (1S) > R pPb (2S) > R pPb (3S) : R pPb (Υ(1S)) = 0.806 ± 0.024 (stat) ± 0.059 (syst), R pPb (Υ(2S)) = 0.702 ± 0.041 (stat) ± 0.058 (syst), R pPb (Υ(3S)) = 0.536 ± 0.058 (stat) ± 0.050 (syst).
We determine the p-value of the observed suppression of Υ states in pPb relative to pp collisions against the hypothesis of A-scaling, which predicts no nuclear modification. Given the uncertainties in the measured R pPb values for Υ(1S), Υ(2S), and Υ(3S), we determine these p-values to be 1.16 × 10 -3 , 1.36 × 10 -5 , and 6.84 × 10 -10 , respectively.
Given that initial-state CNM models predict equal nuclear modification to all three Υ states in contrast to final-state CNM, which result in different levels of nuclear modification, we can determine the p-value of the observed additional suppression of each excited state compared to the ground state. This can be done under the hypothesis that no final-state CNM effects are evident, and the R pPb of the excited and ground states are equal. The p-values of the measured lower R pPb values of the excited states relative to Υ(1S) are 1.24 × 10 -1 for Υ(2S) and 1.02 × 10 -3 for Υ(3S), corresponding to significances of 1.2 and 3.1 standard deviations, respectively. Figure 8 illustrates that the measured modifications in Υ(nS) production in pPb collisions are considerably smaller than those seen in PbPb collisions [36]. A direct comparison of the R AA to the R pPb requires model-dependent scaling of the R pPb to reflect modification by two lead nuclei in PbPb collisions instead of one. Such a comparison of the observed modification effect of CNM on bottomonia in pPb to the nucleus-nucleus collision environment is needed to determine if hot nuclear matter effects in the QGP result in additional suppression of bottomonia in PbPb. Additional modification in PbPb compared to pPb collisions is expected from the presence of color deconfinement as predicted by Refs. [2,3,5,7,64], and by larger comover interaction effects in the dense medium [18].
Tabulated results are provided in the HEPData record for this analysis [65].

Summary
The Υ(nS) (where n = 1, 2, 3) family is studied in proton-lead (pPb) collisions at √ s NN = 5.02 TeV and the production cross sections are presented. Using pp collision data obtained at the same collision energy, the nuclear modification factors R pPb in pPb collisions for the three Υ states are measured. Compared to the hypothesis of scaling by the number of nucleons A, we find the Υ(nS) yields to be suppressed. This suppression is observed over the entire kinematic range that is studied, i.e., transverse momentum p values are consistent with unity for all states, independent of the region used to measure the event activity.
The integrated nuclear modification factors for Υ(nS) in pPb collisions are compared with those measured in PbPb collisions. The nuclear modification factors R AA in PbPb collisions are much smaller than the corresponding R pPb value for each state. However, a similar ordering of the measured R pPb (Υ(nS)) is observed, with Υ(1S) the least suppressed. This suggests the presence of final-state effects in pPb collisions, consistent with predictions from models that break up the bound quarkonium states via interactions with comoving particles from the underlying event. These results will help us to understand how bottomonia are modified in heavy-ion collisions.      [45] CMS Collaboration, "Investigation into the event-activity dependence of Υ(nS) relative production in proton-proton collisions at √ s = 7 TeV", JHEP 11 (2020) 001, doi:10.1007/JHEP11(2020)001, arXiv:2007.04277.