$\Upsilon$ production and nuclear modification at forward rapidity in Pb-Pb collisions at $\mathbf{\sqrt{\textit{s}_{\textbf{NN}}}=5.02}$ TeV

The production of $\Upsilon$ mesons in Pb-Pb collisions at a centre-of-mass energy per nucleon pair $\sqrt{s_{\rm NN}}$ = 5 TeV is measured with the muon spectrometer of the ALICE detector at the LHC. The yields as well as the nuclear modification factors are determined in the forward rapidity region $2.5<y<4.0$, as a function of rapidity, transverse momentum and collision centrality. The results show that the production of the $\Upsilon$(1S) meson is suppressed by a factor of about three with respect to the production in proton-proton collisions. For the first time, a significant signal for the $\Upsilon$(2S) meson is observed at forward rapidity, indicating a suppression stronger by about a factor 2-3 with respect to the ground state. The measurements are compared with transport, hydrodynamic, comover and statistical hadronisation model calculations.


Introduction
The collisions of ultra-relativistic heavy ions are investigated to explore deconfined and chirally restored matter at high temperatures, the quark-gluon plasma (QGP) [1]. The characterisation of the degrees of freedom of the QGP as well as the transition from this new state of matter to ordinary hadrons are central questions to this field of research. The production rate of bound states of two heavy quark-antiquark pair, quarkonia, was proposed as a key observable for deconfinement [2]. Quarkonium states are the best approximations in nature of two static colour charges, hence representing a unique test system of the strong interaction. Since bottom and charm quarks have masses m b/c well above the temperature reached in heavy-ion collisions, they are produced dominantly within a short time scale of the order of 1/(2 · m b/c ) at the beginning of the collision by hard partonic interactions [3,4]. Therefore, heavy quarks experience the whole evolution of the thermodynamic system. As of today, the thermal properties of quarkonium are the subject of intense studies based on lattice quantum chromodynamics (QCD) and effective field theories of QCD [5]. These calculations observe strong modifications of the real and imaginary part of the potential between the heavy quark and its antiquark extracted with increasing precision [5,6]. Beyond the thermal properties of quarkonia in QCD matter, the investigation of the full quantum dynamical treatment of the time evolution involving the interaction between heavy-quark pairs and the medium has started [7], with the ultimate goal of deriving the relevant observables from QCD first principles [5].
Experimentally, charmonium production in nucleus-nucleus collisions was investigated at the SPS, RHIC and LHC [8]. The deviation of the production in nucleus-nucleus collisions with respect to the expectation from the proton-proton binary collision scaling is quantified via the nuclear modification factor. At the LHC, the suppression of J/ψ production in heavy-ion collisions is weaker than measurements at lower energies. This behaviour is commonly interpreted as a sign of (re)generation of charmonium either solely at the phase boundary [9] or during the deconfined stage of the medium evolution [10]. Both scenarios are advertised as signatures of deconfinement. In the bottomonium sector, the CMS collaboration at the LHC pioneered the measurements with the observation of a strong suppression of the ϒ(1S) state in Pb-Pb collisions [11][12][13][14]. Recently, the ALICE and CMS collaborations published the first measurement of the second Fourier coefficient of the azimuthal anisotropy of the ϒ(1S) production that indicates a weaker elliptic flow than the one measured for J/ψ [15,16]. The latter measurement provides a new experimental handle to which extent the emerging bound state participates in the collective motion or is sensitive to path-length dependence of the nuclear modification. The statistical precision of the present measurements is not yet sufficient for strong model constraints. Furthermore, the production of the excited ϒ(2S) state is much more strongly suppressed whereas the ϒ(3S) state has not been observed yet and is at least similarly suppressed as the ϒ(2S) state [11][12][13][14]17]. The interpretation of the data requires to account for the feed-down from the decay of excited states as well as the consideration of effects not related to the occurrence of the QGP, for example, differences between nuclear parton distribution functions (PDFs) and free nucleon PDFs. Nuclear modification factors were also measured in proton-nucleus collisions by the LHC collaborations [18][19][20][21][22] where no formation of a plasma phase was expected prior to the LHC and RHIC measurements [23][24][25][26][27][28][29]. The results indicate a significant modification of the ϒ(1S) production at midrapidity and forward rapidity, in line with expectations from nuclear PDFs [30][31][32] or energy loss considerations [33]. However, the experimental data hint of a stronger suppression than present in these approaches at backward rapidity. In addition, stronger suppressions of the excited states compared to the ϒ(1S) state were observed pointing towards other possible mechanisms, potentially analogue to nucleus-nucleus collisions [34].
Currently, the status of phenomenology in nucleus-nucleus collisions comprises frameworks implementing transport or rate equations [34][35][36][37][38] using a wide range of approaches for the created thermodynamic system and its interaction with the bb quark pair. The statistical hadronisation model, usually applied to the charm sector, was also proposed as a possible scenario for the bottom sector including the production of bottomonia [39].
The ALICE collaboration reported the suppression of ϒ production at forward rapidity 2.5 < y < 4.0 in Pb-Pb collisions at √ s NN = 2.76 TeV [40] and 5.02 TeV [41] based on the 2011 and 2015 data samples, respectively. This article presents the combined analysis of the 2015 and 2018 data sets recorded at √ s NN = 5.02 TeV, corresponding to a three times larger integrated luminosity than the previously published measurements at the same energy. This increase enables to perform a detailed measurement of the ϒ(1S) production as a function of the centrality, transverse momentum and rapidity. For the first time, a significant ϒ(2S) signal is observed in the forward rapidity region in heavy-ion collisions.

Detector, data samples and observables
A detailed description of the ALICE setup and its performance can be found in the references [42,43]. The analysis is based on the detection of muons within the pseudorapidity range −4.0 < η < −2.5 1 reconstructed and identified with the muon spectrometer [44]. In the following, we briefly introduce the detection systems relevant for the ϒ measurements in Pb-Pb collisions.
The primary vertex is determined with the silicon pixel detector, the two innermost layers of the inner tracking system of the central barrel of ALICE [45]. These two cylinders surrounding the beam pipe cover the pseudorapidity range |η| < 2 (first layer) and |η| < 1.4 (second layer) assuming the nominal interaction point (IP). The V0 detector [46] provides the centrality determination. It is made of two arrays of scintillators covering the pseudorapidity intervals −3.7 < η < −1.7 and 2.8 < η < 5. The logical AND of the signals from both subdetectors constitutes the minimum bias (MB) interaction trigger. This trigger decision is fully efficient for the 0-90% most central collisions. Zero-degree calorimeters [47], installed at ±112.5 m from the IP along the beam direction, are used for the rejection of events corresponding to electromagnetic interactions of the colliding lead nuclei.
The muon spectrometer of ALICE consists of a front absorber to filter muons, five tracking stations, a dipole magnet with a field integral of 3 T×m surrounding the third tracking station, an iron wall to reject further punch-through hadrons and low momentum muons from pion and kaon decays, and two trigger stations. These elements are traversed subsequently by the muons originating from the IP region. Each tracking station is composed of two planes of cathode-pad chambers. The two trigger stations are equipped with two planes of resistive plate chambers [48].
The trigger used for this analysis requires a coincidence of the MB signal and a dimuon trigger provided by the trigger stations. The dimuon condition consists of a positively and a negatively charged track candidate above a transverse momentum threshold of 1 GeV/c each. The analysed data set was recorded in 2015 and 2018 and corresponds to a total integrated luminosity of about 760 µb −1 .
The centrality determination relies on a Glauber fit to the sum of the signal amplitudes of the V0 detectors [49][50][51]. The centrality ranges are quoted as quantiles in percent of the total inelastic Pb-Pb cross section. The fit allows each centrality interval to be attributed an average number of participants N part , of binary nucleon-nucleon collisions N coll and the average nuclear overlap function T AA . The analysis comprises measurements integrated over the centrality interval 0-90% and differential in centrality. The Glauber fit quantities used in this analysis are tabulated in the note [51]. For the centrality intervals used in the present analysis and not reported in the note, the relevant values, including uncertainties, were obtained by averaging the corresponding values over the narrower ranges.
The nuclear modification factor is the main observable quantifying the difference of production between Pb-Pb and proton-proton (pp) collisions. It is defined as where N raw ϒ→µ + µ − denotes the raw number of ϒ candidates reconstructed via the dimuon decay channel within a given rapidity, centrality and transverse momentum interval, (A × ε) ϒ→µ + µ − the correction for the geometrical acceptance of the decay muons and the trigger, tracking and track quality selection inefficiencies. The quantity N MB represents the number of equivalent inelastic nucleus-nucleus collisions, T AA the average nuclear overlap function, ∆y and ∆p T the widths of the rapidity and transverse momentum intervals. The term d 2 σ pp ϒ→µ + µ − /dydp T corresponds to the product of the ϒ branching ratio for the dimuon decay channel BR(ϒ → µ + µ − ) and the ϒ production cross section in pp collisions at the same collision energy and in the same kinematic regime of the Pb-Pb measurement. The yields N ϒ→µ + µ − normalised by T AA and the yield ratio between the ϒ(2S) and ϒ(1S) states are also presented as complementary observables. The latter reads In this equation, the yields are corrected for the branching ratios for comparison with model calculations. The relative nuclear modification factor is written simply as In the following, the expression "integrated" stands for "integrated over the 0-90% centrality interval, in 2.5 < y < 4.0 and for p T < 15 GeV/c" unless otherwise specified.
3 Analysis description

Signal extraction
The number of ϒ mesons is extracted from an unbinned log-likelihood fit to the dimuon invariant mass spectra. Opposite-charge pairs are formed with muon tracks meeting the selection criteria listed in the publication [40]. An extended Crystal Ball distribution [52] Fig. 1 for the two background descriptions. After subtraction of the mixed-event background from the raw distribution, the visible residual background is also fitted with an empirical function.
In total, N raw ϒ(1S)→µ + µ − = 3581 ± 119 (stat.) ± 156 (syst.) and N raw ϒ(2S)→µ + µ − = 325 ± 61 (stat.) ± 60 (syst.) are found. The quoted uncertainties are discussed in Section 3.4. No significant ϒ(3S) signal is observed in any of the kinematic regions. Assuming that the event-mixing technique provides the correct background estimation, this method improves the effective signal-to-background ratio of the ϒ(1S) by more than 2.5 and of the ϒ(2S) by almost a factor 2 with respect to the direct extraction of the raw yields.  However, as the description of the background shape is a delicate part of this analysis, the two methods are applied for the signal extraction.

Acceptance and efficiency corrections
The raw yields are corrected for the acceptance and the reconstruction efficiency with a Monte Carlo simulation. The ϒ signals are generated according to p T and y distributions extrapolated from collider data [55, 56], assuming an unpolarised production [57, 58]. These input shapes are adjusted by a nuclear PDF (nPDF) parametrisation to emulate the modification of the initial distributions. The generated resonances are then decayed into muon pairs using the EvtGen package [59], together with PHOTOS [60] to account for final-state radiation. Muons are transported through a realistic modelling of the ALICE apparatus via the GEANT3 code [61], and are injected in recorded MB events from data. This approach allows the experimental data-taking conditions and the occupancy of the detection elements to be accounted for.
The integrated acceptance and efficiency (A × ε), averaged over the 2015 and 2018 data sets, is about 25.7% for all ϒ states. This value varies with p T by less than 1%. Figure 2 shows the A × ε as a function of the collision centrality (left) and rapidity (right). An observed relative decrease of about 10% from peripheral to central events is due to the rise of the occupancy in the muon chambers. The variations as a function of rapidity are a direct consequence of the spectrometer geometry.

Production cross sections in proton-proton collisions
The inclusive production of ϒ mesons in pp collisions at √ s = 5.02 TeV has been recently measured at forward rapidity, 2.5 < y < 4.0 [62]. The relative statistical uncertainties are of the order of the ones obtained in the present Pb-Pb analysis and limit the potential of the R AA study. To overcome this limitation, the ϒ production cross sections in pp collisions at √ s = 5.02 TeV are interpolated from measurements performed at various centre-of-mass energies at the LHC.
The first procedure exploits the ALICE [63, 64] and LHCb [65-67] data within the 2.5 < y < 4.0 acceptance. Empirical functions, such as an exponential or a power law, are fitted to the cross sections as a function of the centre-of-mass energy and then evaluated for √ s = 5.02 TeV. This energy interpolation procedure, described in Ref.
[68], is employed to estimate the integrated and the p T -dependent differential cross sections. For the rapidity-differential cross sections, a further interpolation is performed as in the note [69]. The rapidity dependence of dσ pp ϒ→µ + µ − /dy is fitted and integrated over the ranges matching the present Pb-Pb analysis. The inclusion of CMS measurements [14] constrains the curvature of the fit function down to y = 0. The corresponding uncertainties are calculated as the quadratic sums of the un-  certainties of the data propagated through the interpolation procedures and the spread of the results from the fitted functions. The interpolation results are in agreement with the measured cross sections [62] with systematic uncertainties smaller by more than a factor 2.

Systematic uncertainties
This section describes the different sources of systematic uncertainties to compute the yields and the nuclear modification factors. The raw yields as well as the ϒ(2S)-to-ϒ(1S) yield ratio are evaluated as the average of the results obtained from the following fit variations. The parameters of the Crystal Ball distributions are estimated from MC simulations using either the GEANT3 [61] or GEANT4 [70] transport package. A set of tail parameters from simulations of pp collisions at √ s = 5.02 TeV is also used for the signal extraction. The ratios of widths characterising the signal shapes of ϒ(2S) and ϒ(3S) are varied from 1 to 1.08 and from 1 to 1.14, respectively, in order to account for discrepancies between data and simulation. Several empirical functions are used to model the background shape, whether the event-mixing technique is applied or not. The raw spectrum is fitted with the sum of two decreasing exponentials as well as a pseudo-Gaussian function whose width varies linearly with mass. These functions are defined in the note [52]. Once the mixed-event background is subtracted from the raw distribution, a single exponential or a power-law function are employed. The fits are alternatively performed within two mass ranges, [6][7][8][9][10][11][12][13] and [7][8][9][10][11][12][13][14] GeV/c 2 , to cover different invariant mass regions. The final statistical uncertainty is calculated as the linear average over the uncertainty returned by the fit whereas the systematic uncertainty is estimated as one standard deviation of the distribution of the results. The systematic uncertainty for the signal extraction arises predominantly from the background description uncertainty. At maximum, the relative systematic uncertainty ranges from 3.6% to 10% for the ϒ(1S) and from 9.7% and 39% for the ϒ(2S) as a function of rapidity. The signal extraction uncertainties are uncorrelated to any kinematic variable.
The uncertainties of the calculations of the A × ε corrections have multiple origins. The p T and y distributions in the initial MC conditions are modified by nuclear PDFs or by not considering any shadowing effect. The uncertainty associated to the choice of nPDF set is estimated as the maximum relative difference across the A × ε values obtained when switching the input shapes. The uncertainties due to the tracking, trigger and matching efficiencies are treated as in Ref. [71], the dominant source being the dimuon reconstruction efficiency. Discrepancies between data and MC studies are propagated from single muon tracks to ϒ mesons by multiplying the uncertainties with a factor two. The uncertainties of the MC input shapes as well as of the tracking, trigger and matching efficiencies are uncorrelated with p T and rapidity and fully correlated with the centrality. In addition, the tracking systematic uncertainty con-tains an uncorrelated component as a function of the centrality. The latter is negligible in the peripheral collisions and increases up to 1% for the most central events.
Furthermore, the yields can be sensitive to the uncertainty associated to the definition of the centrality intervals. The corresponding uncertainties are taken from a previous J/ψ analysis [71] as the statistical fluctuations are too large to make a sensitive estimate with the present ϒ measurement. This uncertainty is negligible for the 0-90% interval. The uncertainty assigned to the T AA calculation is detailed in the references [50,51]. The number of equivalent minimum bias events is obtained from the number of analysed dimuon-triggered events normalised by the probability of having a dimuon trigger in a MB event. Its uncertainty is estimated in Ref. [72] and is correlated with all kinematic variables. The uncertainties of the production cross sections in pp collisions are fully correlated with the centrality and uncorrelated as a function of p T and rapidity.
All the sources of systematic uncertainties entering in the computation of the nuclear modification factor are summarised in Tables 1 and 2 for the ϒ(1S) and ϒ(2S), respectively. The ranges indicate the minimum and maximum relative uncertainties as a function of the given kinematic variable. The uncertainties are common to the normalised yields. For the ratio of yields, only the uncertainty attributed to the signal extraction as well as the branching ratio uncertainties [53] are assumed not to cancel.

Results and discussion
The rapidity and p T -differential yields of inclusive ϒ production in Pb-Pb collisions, normalised by T AA = 6.28 ± 0.06 mb −1 for the 0-90% centrality interval [51], are presented in Fig. 3. The results are shown together with the CMS measurements performed at midrapidity [14]. The vertical error bars represent the statistical uncertainties whereas the boxes correspond to the uncorrelated systematic uncertainties described in the previous section. Henceforth, this convention is adopted for all figures. The rapidity dependence indicates that the forward rapidity measurement spans an interesting dynamic range where the signal falls off significantly with respect to the approximate plateau reached around midrapidity. For the first time, a significant ϒ(2S) signal is observed in Pb-Pb collisions at forward rapidity. The measured p T spectrum of ϒ(1S) within the ALICE forward acceptance is softer than at midrapidity as also observed in pp and p-Pb collisions [8].   The measurements are compared with calculations based on transport and rate equations. Within the comover picture, quarkonia are dissociated via the interaction with surrounding particles in the final state [34]. The revisited version of this model aims to explain the suppression of bottomonium production in both p-Pb and Pb-Pb collisions with the same assumptions. It takes into account the nuclear modification of parton distribution functions (nPDFs). Uncertainties from the nCTEQ15 shadowing [32] and the comover-ϒ interaction cross sections are depicted together in the figures as grids. Predictions are also derived from the thermal modification of a complex heavy-quark potential inside an anisotropic plasma [35]. The survival probability of bottomonia is evaluated based on the local energy density, integrating the rate equation over the proper time of each state. The background medium is described with viscous hydrodynamics for three values of the shear viscosity-to-entropy density ratio η/s. These calculations do not include any modification of nuclear PDFs or any regeneration phenomenon. The transport approaches describe an interplay of dissociation and regeneration mechanisms regulating the production of bottomonia at the QGP stage. For the transport model [36], the medium evolves as an expanding isotropic fireball. Results are provided with and without the presence of a regeneration com- The filled boxes at unity correspond to the relative uncertainties correlated with centrality. The results are compared with calculations from the comover and the hydrodynamic models [34,35] in the left panel and with the transport descriptions [36,37] in the right panel.
ponent. The width of the bands represents the modification of the PDF modelled by an effective scale factor on the initial number of bb pairs. This scale factor is neglected in peripheral collisions and is varied between no modification up to 30% suppression at forward rapidity in the most central collisions. In the framework of coupled Boltzmann equations [37], the regeneration is dominated by real-time recombinations of correlated heavy-quark pairs. The simulation of the collision system includes the EPPS16 nPDF parametrisation [31]. In the figures, the calculations are shown with a band due to the nPDF uncertainty and with three curves from the variation of the coupling constants.
Inclusive production can be decomposed into the direct production component and the production from feed-down of higher bottomonium resonances. Direct ϒ(1S) production constitutes approximately 70% to the total inclusive cross section in pp collisions at the LHC [8]. The feed-down contributions of P-wave states to excited ϒ production are not measured at low transverse momentum and are extrapolated. All models treat the excited states and their decay chains to ϒ mesons. They account for a feed-down contribution to ϒ(1S) production in pp collisions consistent with the measured values, but with varying assumptions for the feed-down of the excited states. An uncertainty based on the variation of the feed-down contributions is only considered for the comover interaction model [34].
The predictions from the comover and the hydrodynamic models are compared with the data in the left panel of Fig. 4 while the calculations from the transport approaches are shown in the right panel of the same figure. The various calculations reproduce the trend of the data within uncertainties. For the ϒ(1S), the measurement points lie on the lower limit of the comover interaction model [34] and of the coupled Boltzmann equations [37]. The sharp slope expected in all cases for the R AA of ϒ(2S) is not measurable because of statistical limitations. The current models do not account for the spatial dependence of the nPDF modification: stronger effects are expected for nuclei probed close to their centre. The impact on measurements has been discussed early on [73], a more recent extraction can be found in Ref. [74]. Future models should consider this phenomenon in particular for the discussion of centrality-dependent studies.
Before studying the different approaches via the relative suppression, the ϒ(2S)-to-ϒ(1S) yield ratio in Pb-Pb collisions shown in the left panel of Fig. 5 is compared with the statistical hadronisation model [39]. It assumes that the final state light-flavour hadron yields are calculable from hadron resonance gas densities with a common chemical freeze-out temperature and the baryochemical potential tracing closely the phase boundary from the QGP to hadrons. The approach has been extended to heavy-flavour production assuming a kinetic equilibration of the heavy quarks prior to freeze-out and total heavy quark conservation using the production from initial hard scatterings as an input. The model calculates the abundances of various heavy-flavour species assuming thermal weights. For non-central collisions, a contribution from pp-like scattering behaviour at the surface of the interaction zone is introduced. The calculation underestimates the measured ϒ(2S)-to-ϒ(1S) yield ratio for the 0-30% most central collisions. Taking into account all the uncertainties, the deviation is about one sigma. Comparisons with other measurements in the bottom sector are required to further test the applicability of the model.  [39]. The two curves represent the uncertainty of the pp-like contribution of the corona of the nuclear overlap. (Right) Relative nuclear modification factor along with the predictions from the comover interaction model [34], hydrodynamic calculations [35], from the transport model [36] and calculations based on the coupled Boltzmann equations [37]. The filled red box at unity denotes the uncertainty on the ϒ(2S)-to-ϒ(1S) cross section ratio in pp collisions.
The relative nuclear modification factor defined by Eq. 3 is an appropriate observable to confront the different approaches. Considerations of effects common to both states are thus expected to disappear as indicated by the smaller uncertainties in the right panel of Fig. 5 compared to Fig. 4. The relative uncorrelated systematic uncertainties are 4% smaller when accounting for correlations in the signal extraction, while the systematic uncertainty correlated with the centrality is reduced from 9.2% to 2.5%. Integrated over the 0-90% class, the ϒ(2S)-to-ϒ(1S) R AA is 0.360 ± 0.069 (stat.) ± 0.055 (syst.) i.e. 7.2σ from unity. The comover interaction model shows a deviation with respect to the measurement for the 0-30% most central collisions, also noticed in a comparison with CMS data [34]. The hydrodynamic calculations [35] describe well the data within the present experimental uncertainties. Interestingly, within the transport approaches [36,37], larger double ratio values are achieved by the large regeneration component of the ϒ(2S) production. With more precise measurements, the relative R AA could serve as a model discriminator thanks to the cancellation of sources of uncertainty and the different slopes between the models.
In the following, differential studies of the nuclear modification factor are carried out to scrutinise the suppression features. The dependence of the R AA on the transverse momentum is investigated in Fig. 6 for the 0-90% centrality interval. No significant variation is observed up to 15 GeV/c, in line with hydrodynamic and transport model expectations. The present ϒ(1S) measurement disfavours the hydrodynamic calculation for the highest shear viscosity-to-entropy density ratio. It is interesting to note that the nuclear modification factor in p-Pb collisions exhibits a significant p T dependence in both forward and backward regions [21,22]. The difference in spectral shape between proton-nucleus and nucleusnucleus collisions should receive a close attention from a phenomenological point of view. Figure 7 shows the rapidity dependence of the nuclear modification factors measured by ALICE and  Figure 6: Nuclear modification factor of ϒ(1S) as a function of the transverse momentum. The red box at unity corresponds to the global uncertainty correlated with p T . Predictions from the hydrodynamic [35] and transport [36] models are also shown.

CMS [14]
, and the results are compared with the hydrodynamic model as well as with the calculations based on the coupled Boltzmann transport equations. The experimental data indicate a plateau value of the ϒ(1S) nuclear modification factor of around 0.4 between midrapidity and y ≈ 3. The two most forward measurement points hint of a decrease of the nuclear modification factor down to a value close to 0.3: the R AA is lower by 2σ in the most forward rapidity interval with respect to the central range of the ALICE measurement. The hydrodynamic calculations indicate the opposite behaviour. In this model, the rapidity profile inherits from the initial conditions of the simulated medium [35]. The results from the coupled Boltzmann equations exhibit a structure induced by the rapidity-dependent impact of the used nPDF [37]. The curves cannot describe the CMS and ALICE measurement consistently, albeit the most forward data points lie on the limit of the nPDF uncertainty band. These discrepancies may point towards a physical mechanism not captured in the presently available models. This behaviour will need to be scrutinised further in future analyses improving the current experimental and theoretical uncertainties. Interestingly, the nuclear modification factor of ϒ(1S) calculated with the coherent energy loss model shows a decreasing trend towards forward rapidity at √ s NN = 2.76 TeV [75], even though the model does not reproduce the overall suppression at this centre-of-mass energy. This different behaviour is caused by the conjunction of a non-constant rapidity distribution and the rapidity shift due to the heavy-quark pair medium interaction. All other models do not consider the latter effect. The R AA of the ϒ(2S) is independent of the rapidity within uncertainties in the measured interval with values between 0.07 and 0.20. The rapidity dependence of the models for ϒ(2S) is similar to the one observed for the ground state and is compatible with the experimental measurements.

Conclusions
The results presented in this article provide a detailed measurement of the ϒ(1S) production as well as the first significant measurement of the ϒ(2S) state at forward rapidity in Pb-Pb collisions at the LHC. For the 0-90% centrality class, the nuclear modification factor is 0.353 ± 0.012 (stat.) ± 0.029 (syst.) for ϒ(1S) and 0.128 ± 0.024 (stat.) ± 0.026 (syst.) for ϒ(2S). The corresponding excited-to-ground state ratios are in agreement with hydrodynamic calculations [35], a transport model with a regeneration component [36], and predictions from coupled Boltzmann equations [37]. The results are in tension with the comover interaction model [34] and with the statistical hadronisation model [39] for the 0-30% most central collisions. Taken together with the CMS data, these measurements constrain the rapidity dependence of ϒ suppression with respect to proton-proton collisions. The available models describing bottomonium production in heavy-ion collisions do not capture the rapidity dependence observed for the   Figure 7: Nuclear modification factor of ϒ(1S) and ϒ(2S) as a function of rapidity. The red and violet filled boxes at unity correspond to the global uncertainties common to both ϒ states from the ALICE and CMS measurements. The experimental data are compared with hydrodynamic [35] and coupled Boltzmann equations [37] calculations on the left and right panel, respectively.
R AA of ϒ(1S) in the ALICE acceptance.
Bottomonia are privileged observables to understand the formation and the interaction of bound states in strongly interacting matter and hence to learn about the degrees of freedom of the QGP. This measurement is one of the starting points for more differential studies of bottomonium production. Flow, ϒ(3S) measurements as well as better constraints on parton densities in nuclear collisions, feed-down chains and beauty production cross sections will become available in the upcoming years at the LHC [76].

Acknowledgements
We thank the phenomenology groups for their predictions. The authors would like to extend special thanks to Xiaojian Du, Elena Gonzalez Ferreiro and Xiaojun Yao for enlightening discussions.
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The         [59] D. Lange, "The EvtGen particle decay simulation package", Nucl. Instrum. Meth. A 462 (2001) 152-155.