$\Upsilon$ suppression at forward rapidity in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV

Inclusive $\Upsilon$(1S) and $\Upsilon$(2S) production have been measured in Pb-Pb collisions at the centre-of-mass energy per nucleon-nucleon pair $\sqrt{s_{_{\rm NN}}}=5.02$ TeV, using the ALICE detector at the CERN LHC. The $\Upsilon$ mesons are reconstructed in the centre-of-mass rapidity interval $2.5<y<4$ and in the transverse-momentum range $p_{\rm T}<15$ GeV/$c$, via their decays to muon pairs. In this Letter, we present results on the inclusive $\Upsilon$(1S) nuclear modification factor $R_{\rm AA}$ as a function of collision centrality, transverse momentum and rapidity. The $\Upsilon$(1S) and $\Upsilon$(2S) $R_{\rm AA}$, integrated over the centrality range 0-90%, are $0.37 \pm 0.02 {\rm{(stat)}}\pm 0.03 {\rm{(syst)}}$ and $0.10 \pm 0.04 {\rm{(stat)}}\pm 0.02 {\rm{(syst)}}$, respectively, leading to a ratio $R_{\rm{AA}}^{\Upsilon(\rm2S)}/R_{\rm{AA}}^{\Upsilon(\rm1S)}$ of $0.28\pm0.12\text{(stat)}\pm0.06\text{(syst)}$. The observed $\Upsilon$(1S) suppression increases with the centrality of the collision and no significant variation is observed as a function of transverse momentum and rapidity.


Introduction
A detailed study of the properties of the Quark-Gluon Plasma (QGP) [1] is the main goal of heavyion experiments at ultra-relativistic energies [2][3][4][5][6]. Quarkonia, i.e. bound states of charm or bottom quark-antiquark pairs, are sensitive probes of color deconfinement, due to the Quantum-Chromo Dynamics Debye screening mechanism [7][8][9] leading to quarkonium suppression. Also, since heavy quarks are produced in the initial parton-parton interactions, they experience the full evolution of the medium. Moreover, the various quarkonium states have different binding energies and therefore different dissociation temperatures in a QGP, leading to sequential suppression [7,10]. Theory estimates [11] indicate that bottomonium formation may occur before QGP thermalization [12] because of the large bottom quark mass. In this situation, a quantitative description of the influence of the medium on the bound states becomes challenging. While the dissociation temperatures vary significantly between different models [8,9], it is commonly accepted that the widths of the spectral functions of the bottomonium states increase compared to the widths in vacuum, due to the high temperature of the surrounding medium [13]. Finally, taking into account that feed-down processes from higher-mass resonances (around 40% for the ϒ(1S) and 30% for the ϒ(2S) [9]) are not negligible, the evaluation of the medium temperature via bottomonium measurements remains a complex endeavour.
The first studies of quarkonium production in heavy-ion collisions were devoted to charmonium states, and a suppression of their yields was observed at the SPS [14-16], at RHIC [17,18] and at the LHC [19][20][21][22]. The weaker J/ψ suppression observed at LHC energies, where the centre-of-mass energy per nucleon-nucleon pair ( √ s NN ) is one order of magnitude larger than at RHIC, is now explained by means of a competitive (re)generation mechanism, which occurs during the deconfined phase and/or at the hadronization stage [23][24][25][26]. This production mechanism strongly depends on the (re)combination probability of deconfined quarks present in the medium and thus on the initial number of produced cc pairs. The effect has been found to be more important at low p T and in the most central collisions [20,22,27].
The high-energy collisions delivered by the LHC allow for a detailed study of bottomonium states. For bottomonium production, perturbative calculations of production rates in elementary nucleon-nucleon collisions are more reliable than for charmonium yields due to the higher mass of the bottom quark with respect to charm. Since the the number of produced bb pairs in central heavy-ion collisions amount to a few pairs per event at the LHC, the probability for (re)generation of bottomonia through (re)combination is much smaller than in the case of charmonia.
A strong suppression of the ϒ(1S) state in Pb-Pb collisions with respect to properly scaled measurements in pp collisions has been observed at √ s NN = 2.76 TeV by ALICE [28] and CMS [29,30] in the rapidity ranges 2.5 < y < 4 and |y| < 2.4, respectively. The ϒ(1S) nuclear modification factor R AA is quantified as the ratio of the ϒ(1S) yield in nucleus-nucleus collisions to the production cross section in pp collisions times the nuclear overlap function T AA obtained via the Glauber model [31,32]. The suppression increases with the centrality of the collision, reaching about 60% and 80% for the most central collisions at mid-[30] and forward rapidity [28], respectively. Moreover, the ϒ(2S) suppression reaches about 90% and for ϒ(3S) data are compatible with a complete suppression [30]. As a function of p T the ϒ(1S) R AA , measured for p T < 20 GeV/c by CMS [30], is compatible with a constant value. When considering the y-dependence resulting from the comparison of ALICE and CMS results, there is an indication for a stronger suppression at forward y. Transport models [26,33] as well as an anisotropic hydro-dynamical model [34] fairly reproduce the experimental observations of CMS, while they tend to overestimate the R AA values measured by ALICE.
The bottomonium suppression due to the QGP should be disentangled from the suppression due to Cold Nuclear Matter (CNM) effects, such as the nuclear modification of the parton distribution functions due to shadowing [35,36], as well as parton energy loss [37]. These effects on the bottomonium production were studied in p-Pb collisions by ALICE [38] and LHCb [39], who reported for the ϒ(1S) a ϒ suppression at forward rapidity in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration nuclear modification factor slightly lower than unity at forward rapidity and compatible with unity at backward rapidity, although with significant uncertainties. Recently, ATLAS results indicate a significant suppression of the ϒ(1S) for p T < 40 GeV/c around mid-rapidity [40]. Additional measurements at forward/backward rapidity with higher statistics, are needed to fully constrain the models and perform a meaningful extrapolation of CNM effects to Pb-Pb collisions.
In this Letter we present the first results on the ϒ(1S) and ϒ(2S) R AA measured by the ALICE Collaboration in Pb-Pb collisions at √ s NN = 5.02 TeV. The pp reference cross sections used in the R AA calculations have been determined by an interpolation procedure based on various ALICE [41,42] and LHCb [43,44] results at different energies. The nuclear modification factor for the ϒ(1S) is presented as a function of the centrality of the collision and also differentially in p T and rapidity. For the ϒ(2S), an R AA value integrated over the centrality of the collision is quoted. Finally, the results are compared to theoretical calculations.

Experimental apparatus and data sample
An extensive description of the ALICE apparatus can be found in [45,46]. The analysis presented in this Letter is based on muons detected at forward rapidity (2.5 < y < 4) 1 with the muon spectrometer [47]. The detectors relevant for ϒ measurements in Pb-Pb collisions are described below.
The Silicon Pixel Detector, corresponding to the two innermost layers of the Inner Tracking System [48], is used for the primary vertex determination. The inner and outer layer cover the pseudo-rapidity ranges |η| < 2 and |η| < 1.4, respectively.
The V0 scintillator hodoscopes [49] provide the centrality estimate. They are made of two arrays of scintillators placed in the pseudo-rapidity ranges 2.8 < η < 5.1 and −3.7 < η < −1.7. The logical AND of the signals from the two hodoscopes constitutes the Minimum Bias (MB) trigger. The MB trigger is fully efficient for the studied 0-90% most central collisions.
The Zero Degree Calorimeters (ZDC) are installed at ±112.5 m from the nominal interaction point along the beam line. Each of the two ZDCs is composed of two sampling calorimeters designed for detecting spectator protons, neutrons and nuclear fragments. The evaluation of the signal amplitude of the ZDCs allows for the rejection of events corresponding to an electromagnetic interaction of the colliding Pb nuclei [50].
The muon spectrometer covers the pseudorapidity range −4 < η < −2.5. It is composed of a front absorber, which filters muons upstream of the muon tracker, consisting of five tracking stations with two planes of cathode-pad chambers each, and of a dipole magnet providing a 3 T·m integrated magnetic field. Downstream of the tracking system, a 1.2 m thick iron wall stops efficiently the punch-through hadrons. The muon trigger system is located downstream of the iron wall and consists of two stations, each one equipped with two planes of Resistive Plate Chambers (RPC), with an efficiency higher than 95% [51]. The muon-trigger system is able to deliver single and dimuon triggers selecting muons with p T larger than a programmable threshold, via an algorithm based on the RPC spatial information [52]. Throughout its entire length, a conical absorber shields the muon spectrometer against secondary particles produced by the interaction of primary particles in the beam pipe.
The trigger condition used for data taking is a dimuon-Minimum Bias (µ µ-MB) trigger formed by the logical AND of the MB trigger and an unlike-sign dimuon trigger with a p T threshold of 1 GeV/c for each of the two muons.
The centrality estimation is performed using a Glauber fit to the sum of the signal amplitudes of the 1 In the ALICE reference frame, the muon spectrometer covers a negative η range and consequently a negative y range. We have chosen to present our results with a positive y notation ϒ suppression at forward rapidity in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration V0 scintillators [53,54]. Centrality ranges are given as percentages of the total hadronic Pb-Pb cross section. In addition to the centrality, the Glauber model allows an estimate of the average number of participant nucleons N part , of the average number of binary collisions N coll and of the nuclear overlap function T AA , for each centrality interval [55]. In the present analysis, the data sample corresponds to an integrated luminosity L int ≈ 225 µb −1 in the centrality interval 0-90% that has been divided into four centrality classes: 0-10%, 10-30%, 30-50% and 50-90%.

Data analysis
The evaluation of R AA is performed through the following expression: where N ϒ is the number of detected resonance decays to muon pairs, while BR ϒ→µ + µ − = (2.48 ± 0.05)% and (1.93 ± 0.17)% are the branching ratios for the dimuon decay of ϒ(1S) and ϒ(2S), respectively [56]. The (A × ε) ϒ→µ + µ − factor is the product of acceptance and detection efficiency for the ϒ state under study. The normalization factor N µ µ-MB · F norm is the product of the number of analyzed µ µ-MB events and the inverse of the probability to obtain an unlike-sign dimuon trigger in a MB-triggered event [22]. Finally, σ ϒ pp is the reference pp cross section and T AA represents the nuclear overlap function. The signal yields are evaluated by performing fits to the µ + µ − invariant mass distributions. In order to improve the purity of the dimuon sample a set of selection criteria [28] has been applied on the muon tracks, including the request of the matching between the tracks reconstructed in the trigger and tracking detectors of the muon spectrometer and a cut on the track transverse momentum (p T > 2 GeV/c). The latter cut has a small effect on the number of detected resonances. The raw ϒ yields are extracted using the sum of three extended Crystal Ball (CB) functions [57], one for each of ϒ(1S), ϒ(2S) and ϒ(3S). The extended CB function consists of a Gaussian core with non-Gaussian tails to take into account the radiative contributions of the ϒ production and the absorber effects, such as multiple Coulomb scattering and muon energy loss. The background is fitted with the sum of two exponential functions (see left panel of Fig. 1). Since the signal-to-background (S/B) ratio is low in the tail regions of the extended CB functions, the tail parameters are fixed to values obtained from the Monte Carlo (MC) simulation. The mass position and the width parameters of the ϒ(1S) are left free for the integrated spectrum (i.e. ϒ suppression at forward rapidity in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration centrality class 0-90%, p T < 15 GeV/c and 2.5 < y < 4). Whereas for the signal extraction as a function of centrality, the mass position and width of the ϒ(1S) are fixed to the values obtained in the fit to the centrality-integrated (0-90%) mass spectrum. Finally, for studies as a function of p T and y, the mass position and the width obtained for the centrality-integrated mass spectrum are scaled according to their evolution observed in the MC. Due to the poor S/B ratio for the higher mass states, the values of the mass of the ϒ(2S) and ϒ(3S) are fixed to the PDG mass differences with respect to the ϒ(1S), and the ratio of ϒ(2S) (ϒ(3S)) to ϒ(1S) widths is fixed to values from the MC simulation, i.e. 1.03 (1.06). In the fit shown in Fig. 1 only signals corresponding to the ϒ(1S) and ϒ(2S) are visible, since the ϒ(3S) contribution is compatible with zero events. Alternatively, the combinatorial background is modeled with the event-mixing method. In this approach, an invariant mass dimuon spectrum is constructed by pairing muons from different events with similar multiplicities as described in [22]. The combinatorial background is then subtracted from the raw dimuon spectrum (right panel of Fig. 1) and the resulting distribution is fitted with the sum of three extended CB and an exponential function to account for the residual background. Finally, the number of detected ϒ resonances, N ϒ , is obtained as the average [57] of the fitting methods described above (and also below in the discussion on signal systematics), leading to N ϒ(1S) = 1126 ± 53(stat) ± 47(syst) and N ϒ(2S) = 77 ± 33(stat) ± 17(syst).
The measured ϒ yields, N ϒ , are corrected for the detector acceptance and efficiency using MC simulations. Since the occupancy of the detector varies with the centrality of the collisions, the generated ϒ decays are embedded into real MB events to simulate the various particle multiplicity scenarios as in data.
The p T and y distributions of the generated ϒ are obtained from existing pp measurements [58-60] using the interpolation procedure described in [61]. The EKS98 nuclear shadowing parameterization [35] is used to include an estimate of CNM effects. Since available data favor a small or null polarization for ϒ(1S) [62-64], an unpolarized production is assumed. The variations of the performance of the tracking and triggering systems throughout the data-taking period as well as the residual misalignment of the tracking chambers are taken into account in the simulation. The A × ε values, for the range p T < 15 GeV/c, 2.5 < y < 4 and the 0-90% centrality class are 0.263 and 0.264 for the ϒ(1S) and ϒ(2S), respectively, with a negligible statistical uncertainty. A decrease of 2% is observed in A × ε for the 0-10% central collisions with respect to the 50-90% sample due to the higher occupancy in the most central events. The A × ε is higher by 20% in 3 < y < 3.5 compared to the values at 2.5 < y < 3 and 3.5 < y < 4, whereas it has no variation as a function of p T . The systematic uncertainty on A × ε is discussed below.
The systematic uncertainty on the signal extraction is evaluated using various functions for modelling the background shape, as well as adopting two fitting ranges, i.e. (7-14) GeV/c 2 and (7.5-14.5) GeV/c 2 . The tail parameters of the signal functions have been varied using estimates provided by two MC particle transport models: GEANT4 [65] and GEANT3 [66]. In the centrality, p T or y differential studies, the mass position and width are also varied by amounts, which correspond to the uncertainties on the mass position and the width returned by the fit to the centrality-integrated invariant mass spectrum. The ratio of ϒ(2S) (ϒ(3S)) to ϒ(1S) widths is varied from 1 (1) to 1.06 (1.12). The values of N ϒ and their statistical uncertainties are obtained by taking the average of N ϒ and of the corresponding statistical uncertainties from the various fits. This procedure is applied to both fits of the raw and combinatorial-background subtracted spectra. The systematic uncertainties are estimated as the root mean square of the distribution of N ϒ obtained from the various fits. The effect induced by the p T > 2 GeV/c cut on single muons on the A × ε-corrected ϒ yields was estimated by varying that cut by ±0.2 GeV/c. A ±2% maximum variation on N ϒ /(A × ε) was observed and included in the systematic uncertainties.
Various sources contribute to the systematic uncertainties of A × ε, such as the p T and y shapes of the input distributions for the MC simulations, the trigger efficiency, the track reconstruction efficiency and finally the matching efficiency between tracks in the muon tracking and triggering chambers. Various sets of simulations are produced with different ϒ input p T and y distributions, obtained from empirical parameterizations and/or extrapolations of available data sets at different energies. The maximum relative 6.6-11.3%(II) 5.5-11.5%(II) 6.3% 7.5% difference of A × ε for the various shapes is taken as the systematic uncertainty due to the input MC. In order to calculate the systematic uncertainty on trigger efficiency, the trigger response function for single muons is evaluated using either MC or data. The two response functions are then separately applied to simulations of an ϒ sample and the difference obtained for the ϒ reconstruction efficiency is taken as systematic uncertainty. The systematic uncertainty on the tracking efficiency is obtained starting from an evaluation of the single muon tracking efficiency in MC and data. This evaluation is performed via a procedure, detailed in [22], based on the redundancy of the tracking chamber information. The dimuon tracking efficiency is then obtained by combining the single muon efficiencies and the systematic uncertainty is taken as the difference of the values obtained with the procedure based on MC and data. The muon tracks for data analysis are chosen based on a selection on the χ 2 of the matching between a track segment in the trigger system with a track in the tracking chambers. The matching systematics are obtained by varying the χ 2 selection cut in data and MC and comparing the effects on the muon reconstruction efficiency [22].
The systematic uncertainty on the centrality measurement is evaluated by varying the V0 signal amplitude by ±0.5% corresponding to 90% of the hadronic cross section in Pb-Pb collisions, used as anchor point to define the centrality classes. The systematic uncertainty on the evaluation of σ pp ϒ is detailed in the next section. Finally, the systematic uncertainty evaluation of F norm and T AA are described in [22] and [53], respectively. The different systematic uncertainty sources on the R AA calculation are summarized in Table 1. If the above mentioned systematic uncertainty is correlated as a function of centrality, p T or y, it is quoted as correlated (type I) systematic uncertainty, otherwise it is treated as uncorrelated (type II).

Proton-proton reference cross sections
The pp reference cross section for ϒ(1S) and ϒ(2S) production are computed by means of an interpolation procedure as described for ϒ(1S) in [67]. The energy interpolation for the ϒ cross section, as a function of rapidity and for the p T and y integrated result, uses the measurements of ϒ production cross sections in pp collisions at √ s = 7 and 8 TeV by ALICE [41, 42] and at √ s = 2.76, 7 and 8 TeV by LHCb [43, 44]. The interpolation is performed by using various empirical functions and, in addition, the shape of the energy dependence of the bottomonium cross sections calculated using two theoretical models, i.e. the Leading Order Colour Evaporation Model (LO-CEM) [68] and the Fixed Order Next-to-Leading Logs (FONLL) model [69]. The latter gives cross sections for open beauty, which is here used as a proxy to study the evolution of the bottomonium cross section. The energy interpolation for the ϒ(1S) cross section as a function of p T is based on LHCb measurements only, since the p T coverage of the results of this analysis (p T < 15 GeV/c) is more extended than that of the corresponding ALICE pp data (p T < 12 ϒ suppression at forward rapidity in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration  GeV/c). The result of the interpolation procedure gives BR ϒ(1S)→µ = 302 ± 23(syst) pb assuming unpolarized quarkonia and integrating over the ranges 2.5 < y < 4 and p T < 15 GeV/c. The uncertainties correspond to the quadratic sum of two terms. The first term dominates the total uncertainty on the interpolated value and reflects the statistical and systematic uncertainties on the data points used in the interpolation procedure. The second term is related to the spread among the interpolated cross sections obtained by using either the empirical functions or the energy dependence estimated from the theoretical models mentioned above. The numerical values obtained from the interpolation procedure are summarized in Table 2 for the various kinematic ranges used in the analysis.

Results
The nuclear modification factors for inclusive ϒ(1S) and ϒ(2S) production in Pb-Pb collisions at is 0.28 ± 0.12(stat) ± 0.06(syst). Since the decay kinematics of the two ϒ states is very similar, most of the systematic uncertainty sources entering the ratio cancel out except those on the signal extraction and on the pp cross section, which are the dominant contributions to the total systematic uncertainty. The measurements show a strong suppression for both bottomonium states with the more weakly bound state being significantly more suppressed. The ratio between the ϒ(1S) R AA at √ s NN = 5.02 TeV and 2.76 TeV is 1.23 ± 0.21(stat) ± 0.19(syst). The sources of systematic uncertainties entering the calculation of the ratio are considered uncorrelated, except for the T AA component, whose uncertainty cancels out. The ratio is compatible with unity within uncertainties.
The centrality, p T and y dependences of the ϒ(1S) R AA at forward rapidity at √ s NN = 5.02 TeV are shown in Fig.2. A decrease of R AA with increasing centrality is observed down to R ϒ(1S) AA = 0.33 ± 0.03(stat) ± 0.03(syst) for the 0-10% most central collisions. No significant p T -dependence is observed up to p T = 15 GeV/c within uncertainties. The nuclear modification factor shows no significant dependence on rapidity.
The inclusive ϒ(1S) R AA measurements are compared in Fig. 2 to several calculations: two transport models (TM) [33,70] and one hydro-dynamical model [34]. To describe the quarkonium motion in the medium, both transport codes use a rate-equation approach which accounts for both suppression and (re)generation mechanisms in the QGP. In the TM1 model [33] the evolution of the thermal medium is based on a thermal-fireball expansion while the TM2 model [70] uses a 2+1 dimensional version of the ideal hydrodynamic equations. The two models use different rate equations and both models include a feed-down contribution from higher-mass bottomonia to the ϒ(1S). In TM2, two sets of feeddown fractions are assumed. Finally, the ϒ(1S) production cross section in pp collisions at √ s = 5.02 TeV in the rapidity range 2.5 < y < 4 is taken as dσ  using the pp interpolation method reported in the previous section. TM1 predictions are shown as bands accounting for shadowing effects as calculated in [71]. The upper limit shown in Fig. 2 corresponds to the extreme case of the absence of shadowing while the lower limit reflects a reduction of 30% due to shadowing. The TM1 model implements the feed-down fractions reported in [9]. In the TM2 model, the shadowing parameterization is based on EKS98 [35] and the band edges correspond to two different sets of feed-down fractions (27% from χ b ; 11% from ϒ(2S + 3S) and 37% from χ b ; 12% from ϒ(2S + 3S)) adopted by the authors. In the third model [34], a thermal suppression of the bottomonium states is calculated using a complex-valued heavy-quark potential parametrized by means of lattice QCD and embedded in a medium evolving according to 3+1d anisotropic hydrodynamics. In this recent study, the R AA shows no sensitivity to the plasma shear viscosity-to-entropy density ratio (4πη/s) parameter of the hydro evolution, which is therefore set to 4πη/s = 2 consistent with particle spectra fits. The band of the model quantifies the heavy-quark potential uncertainty, which has been estimated by including a ±15% variation of the Debye mass of the QCD medium that is tuned by a fit to the real-part of the lattice in-medium heavy-quark potential. Furthermore, the predictions shown are referring to the initial momentum-space anisotropy parameter ξ 0 = 0, which corresponds to a perfectly isotropic QGP at the starting point of the hydro-dynamical evolution at τ 0 = 0.3 fm/c. Finally, this model accounts for feeddown contributions but it includes neither a (re)generation mechanism nor CNM effects. The centrality dependence of the ϒ(1S) R AA is fairly reproduced by the model calculations in the top panel of Fig. 2.
ϒ suppression at forward rapidity in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration The data are best described by TM1 when (re)generation is included and by TM2 when (re)generation is not taken into account. The hydro-dynamical model describes the trend of the data, the fact that the data lie on the upper edge of the uncertainty band for N part > 70 could indicate a smaller Debye mass and thus a stronger heavy-quark potential. The data as a function of p T (bottom left panel of Fig. 2) can be described with or without the (re)generation scenario of the TM1 model while showing agreement with the hydro-dynamical model for the upper edge of the uncertainty band. Finally, the y-dependence of the ϒ(1S) R AA is described, within uncertainties, by the hydro-dynamical model in the bottom right panel of Fig. 2 despite the possibly different trend between data and calculations.
The low ϒ(1S) R AA reported in this Letter raises the important question whether direct ϒ(1S) are suppressed at LHC energies or only the feed-down contribution from higher mass states. However, the large uncertainties of the current measurements of CNM effects [38-40] prevent a firm conclusion.

Summary
The nuclear modification factors of inclusive ϒ(1S) and ϒ(2S) production at forward rapidity (2.5 < y < 4) and for p T < 15 GeV/c in Pb-Pb collisions at √ s NN = 5.02 TeV have been measured using the ALICE detector. The observed ϒ(1S) suppression increases with the centrality of the collision and no significant variation is observed as a function of transverse momentum or rapidity. A larger suppression of the ϒ(2S) bound state compared to the ground state is also reported. Transport and dynamical model calculations reproduce qualitatively the centrality and kinematic dependence of the ϒ(1S) nuclear modification factor.            [55] D. G. d'Enterria, "Hard scattering cross-sections at LHC in the Glauber approach: From pp to pA and AA collisions," arXiv:nucl-ex/0302016 [nucl-ex].   ϒ suppression at forward rapidity in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration ϒ suppression at forward rapidity in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration ϒ suppression at forward rapidity in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration