Differential studies of inclusive J/$\psi$ and $\psi$(2S) production at forward rapidity in Pb-Pb collisions at $\mathbf{\sqrt{{\textit s}_{_{NN}}}}$ = 2.76 TeV

The production of J/$\psi$ and $\psi(2S)$ was measured with the ALICE detector in Pb-Pb collisions at the LHC. The measurement was performed at forward rapidity ($2.5<y<4 $) down to zero transverse momentum ($p_{\rm T}$) in the dimuon decay channel. Inclusive J/$\psi$ yields were extracted in different centrality classes and the centrality dependence of the average $p_{\rm T}$ is presented. The J/$\psi$ suppression, quantified with the nuclear modification factor ($R_{\rm AA}$), was studied as a function of centrality, transverse momentum and rapidity. Comparisons with similar measurements at lower collision energy and theoretical models indicate that the J/$\psi$ production is the result of an interplay between color screening and recombination mechanisms in a deconfined partonic medium, or at its hadronization. Results on the $\psi(2S)$ suppression are provided via the ratio of $\psi(2S)$ over J/$\psi$ measured in pp and Pb-Pb collisions.


Introduction
At high temperature, lattice quantum chromodynamics predicts the existence of a deconfined phase of quarks and gluons where chiral symmetry is restored [1].This state of matter is known as the Quark Gluon Plasma (QGP) [2], and its characterization is the goal of ultra-relativistic heavy-ion collision studies.
Among the probes used to investigate the QGP and quantify its properties, quarkonium states are one of the most prominent and have generated a large amount of results both theoretical and experimental.According to the color-screening model [3,4], measurement of the in-medium dissociation probability of the different quarkonium states could provide an estimate of the system temperature.Dissociation is expected to take place when the medium reaches or exceeds the critical temperature for the phase transition (T c ), depending on the binding energy of the quarkonium state.In the charmonium (cc) family, the strongly bound J/ψ could survive significantly above T c (1.5-2 T c ) whereas χ c and ψ(2S) melting should occur near T c (1.1-1.2T c ) [5,6].The determination of the in-medium quarkonium properties remains a challenging theoretical task.Intense and persistent investigations on the theory side are ongoing [7].Shortly after quarkonium suppression was suggested as a strong evidence of QGP formation, the first ideas of charmonium enhancement via recombination of c and c appeared [8,9].Since then, the J/ψ enhancement mechanism has been more formalized and quantitative predictions [10][11][12][13][14] were made.Since the charm quark density produced in hadronic collisions increases with energy [15], recombination mechanisms are predicted to give rise to a sizable J/ψ production at LHC energies, which is likely to partially compensate or exceed the J/ψ suppression due to color-screening in the QGP.The observation of J/ψ enhancement in nucleus-nucleus collisions via recombination would constitute an evidence for deconfinement and hence for QGP formation.In addition, information for the characterization of the QGP can come from the study of the ψ(2S) meson, a state which is less strongly bound and not affected by higher mass charmonium decays with respect to the J/ψ.In the pure melting scenario, the relative production of ψ(2S) with respect to J/ψ is expected to be very small at the LHC [4], which is not the case if recombination occurs [16,17].J/ψ suppression was observed experimentally in the most central heavy-nucleus collisions at the SPS [18,19], RHIC [20][21][22][23] and LHC [24][25][26][27][28], ranging from a center-of-mass energy per nucleon pair ( √ s NN ) of about 17 GeV to 2.76 TeV.The ψ(2S) suppression was measured at the SPS [29] and the LHC [30].The interpretation of these results is not straightforward as they are also subject to other effects, not all related to the presence of a QGP.A fraction of J/ψ originates from the strong and electromagnetic feed-down of the χ c and ψ(2S).Therefore, a melting of these higher mass states before they can decay into the J/ψ will lead to an effective suppression of the J/ψ yield already for a medium that does not reach the J/ψ dissociation temperature.Assuming charmonium states are initially produced with the same relative abundancies in Pb-Pb collisions as in pp collisions, the χ c and ψ(2S) melting would result in a reduction of the J/ψ yield of about 40% [31].In addition, a non-prompt J/ψ and ψ(2S) component from the weak decay of beauty hadrons also contributes to the inclusive measurements.Since the beauty hadrons decay outside the QGP volume, this contribution is not sensitive to the color-screening of charmonia.Finally, a fraction of the J/ψ and ψ(2S) suppression can be ascribed to cold nuclear matter (CNM) effects, also present in protonnucleus collisions [32,33].The CNM effects group together the nuclear absorption of the charmonia, the modification of the parton distribution functions (PDF) in the nuclei that leads to a reduction (shadowing) or an enhancement (anti-shadowing) of the cc pair production, and the energy loss of charm quarks in the nucleus.
Numerous studies of J/ψ production in different collision systems at different energies are now available.
Comparisons between experiments and to theoretical models can be made over wide kinematic ranges in rapidity and transverse momentum.We already published the centrality, transverse momentum (p T ) and rapidity (y) dependence of the J/ψ nuclear modification factor in Pb-Pb collisions at √ s NN = 2.76 TeV [26,

The ALICE detector
The ALICE detector is described in detail in [34].At forward rapidity (2.5 < y < 4) the production of quarkonium states is studied in the muon spectrometer via their µ + µ − decay channels down to zero p T .
In the ALICE reference frame, the positive z direction is along the counter-clockwise beam direction.The muon spectrometer covers a negative pseudo-rapidity (η) range and consequently a negative y range.However, due to the symmetry of the Pb-Pb system, the results are presented with a positive y notation, while keeping the negative sign for η.
The muon spectrometer consists of a ten-interaction-lengths (4.1 m) thick absorber, which filters the muons, in front of five tracking stations comprising two planes of cathode pad chambers each.The third station is located inside a dipole magnet with a 3 Tm field integral.The tracking apparatus is completed by a Muon Trigger system (MTR) composed of four planes of resistive plate chambers downstream from a seven-interaction-lengths (1.2 m) thick iron wall, which absorbs secondary hadrons escaping from the front absorber and low-momentum muons coming mainly from charged pion and kaon decays.A small-angle conical absorber protects the tracking and trigger chambers against secondary particles produced by the interaction of large rapidity primary particles with the beam pipe.Finally, a rear absorber protects the trigger chambers from the background generated by beam-gas interactions downstream from the spectrometer.
In addition, the Silicon Pixel Detector (SPD), scintillator arrays (V0) and Zero Degree Calorimeters (ZDC) were used in this analysis.The SPD consists of two cylindrical layers covering |η| < 2.0 and |η| < 1.4 for the inner and outer ones, respectively, and provides the coordinates of the primary vertex of the collision.The V0 counters, two arrays of 32 scintillator tiles each, are located on both sides of the nominal interaction point and cover 2.8 < η < 5.1 (V0-A) and −3.7 < η < −1.7 (V0-C).The ZDC are located on either side of the interaction point at z ≈ ±114 m and detect spectator nucleons at zero degree with respect to the LHC beam axis.The V0 and ZDC detectors provide triggering information and event characterization.

Data sample
The data sample analysed in this paper corresponds to Pb-Pb collisions at √ s NN = 2.76 TeV.These collisions were delivered by the LHC during 190 hours of stable beam operations spread over three weeks in November and December 2011.
The Level-0 (L0) minimum bias (MB) trigger was defined as the coincidence of signals in V0-A and V0-C detectors synchronized with the passage of two crossing lead bunches.This choice for the MB condition provides a high triggering efficiency (> 95%) for hadronic interactions.To improve the trigger purity, a threshold on the energy deposited in the neutron ZDC rejects the contribution from electromagnetic dissociation processes at the Level-1 (L1) trigger level.Beam induced background is further reduced at the offline level by timing cuts on the signals from the V0 and the ZDC.
The charmonium analysis was carried out on a data sample, where in addition to the MB prerequisite, a trigger condition of at least one or two reconstructed muon candidate tracks in the MTR (trigger tracks) was required in each event.The MTR logic allows for programming several L0 trigger decisions based on (i) the detection of one or two muon trigger tracks, (ii) the presence of opposite-sign or like-sign trigger track pairs and (iii) a lower threshold on the approximate transverse momentum (p trig T ) of the muon candidates.The latter selection is performed by applying a cut on the maximum deviation of the trigger track from an infinite momentum track originating at the nominal interaction point.Due to the finite spatial resolution of the trigger chambers, this does not lead to a sharp cut in p T , and the corresponding p A data sample of 17.3 • 10 6 Pb-Pb collisions was collected with the µ µ-MB trigger, defined as the coincidence of the MB and MUL conditions.A scaling factor F norm is computed for each run -corresponding to a few hours maximum of continuous data taking -in order to normalize the number of µ µ-MB triggers to the number of equivalent MB triggers.It is defined as the ratio, in a MB data sample, between the total number of events and the number of events fulfilling the µ µ-MB trigger condition.It should be noted that the MB sample used in this calculation was recorded in parallel to the µ µ-MB triggers.The F norm value, 30.56 ± 0.01(stat.)± 1.10(syst.),is given by the average over all runs weighted by the statistical uncertainties.A small fraction of opposite-sign dimuons were misidentified by the trigger algorithm as like-sign pairs.Although for the J/ψ it amounts to less than 1% when considering the full sample, it increases up to 4% at high p T in peripheral collisions.In this analysis, the missing fraction of opposite-sign dimuons was recovered by extracting the number of produced J/ψ and ψ(2S) from the union of the MUL and MLL data sample (MUL∪MLL).This is different from the selection applied in the former paper [27], where only the MUL data sample was used.On the other hand, the efficiency of the trigger algorithm to determine the sign of the muon pairs does not impact the normalization of the collected data sample to the number of equivalent MB events described above.This was cross-checked by computing the normalization factor of the MUL∪MLL data sample, resulting in less than 1% difference in the extracted number of equivalent MB events.

Definition of observables
The centrality determination is based on a fit of the V0 signal amplitude distribution as described in [36].
Variables characterizing the collision such as the average number of participant nucleons ( N part ) and the average nuclear overlap function ( T AA ) for each centrality class are given in Tab. 1.In this analysis a cut corresponding to the most central 90% of the inelastic nuclear cross section was applied as for these events the MB trigger is fully efficient and the residual contamination from electromagnetic processes is negligible.
J/ψ and ψ(2S) The average number of participant nucleons N part and the average value of the nuclear overlap function T AA with their associated systematic uncertainties for the centrality classes, expressed in percentages of the nuclear cross section [36], used in these analyses.
For each centrality class i, the measured number of J/ψ (N i J/ψ ) is normalized to the equivalent number of minimum bias events (N i events ).To obtain N i events , one simply multiplies the number of µ µ-MB triggered events by the F norm factor scaled by the width of the centrality class.Corrections for the branching ratio of the dimuon decay channel (BR J/ψ→µ + µ − = 5.93 ± 0.06%) and for the acceptance times efficiency (A × ε i ) of the detector are then applied.The J/ψ yield (Y i J/ψ ) in a centrality class i is given by It is then combined with the inclusive J/ψ cross section measured in pp collisions at the same energy to form the nuclear modification factor R AA defined as The p T and y integrated J/ψ cross section is σ pp J/ψ (p T < 8 GeV/c, 2.5 < y < 4) = 3.34 ± 0.13(stat.)± 0.24(syst.)± 0.12(luminosity) +0. 53 −1.07 (polarization)µb [37].The ALICE measurements reported here refer to inclusive J/ψ yields, i.e. include prompt J/ψ (direct J/ψ and feed-down from ψ(2S) and χ c ) and non-prompt J/ψ (decay of B-mesons).Contrary to prompt J/ψ, J/ψ from B-meson decays do not directly probe the hot and dense medium created in the Pb-Pb collisions.Beauty hadron decays occur outside the QGP, so the non-prompt J/ψ R AA is instead related to the energy loss of the beauty quarks in the medium.Although the prompt J/ψ R AA cannot be directly measured with the ALICE muon spectrometer, it can be evaluated via where F B is the fraction of non-prompt to inclusive J/ψ measured in pp collisions, and R non-prompt AA is the nuclear modification factor of J/ψ from B-meson decays in Pb-Pb collisions.The non-prompt and prompt J/ψ differential cross sections as a function of p T and y were measured by LHCb in pp collisions at √ s = 2.76 and 7 TeV [38,39] in a kinematic range overlapping with that of the ALICE muon spectrometer.Therefore, one can extract the p T and y dependence of F B from these data and use it in Eq. 3. A reliable determination of R non-prompt AA presents further complications.We have thus chosen two extreme hypotheses, independent of centrality, corresponding to the absence of medium effects on beauty hadrons (R non-prompt AA = 1) or to a complete suppression (R non-prompt AA = 0), to evaluate conservative limits on R prompt AA .An excess of J/ψ compared to the yield expected assuming a smooth evolution of the J/ψ hadro-production and nuclear modification factor was observed in peripheral Pb-Pb collisions at very low p T [40].This excess might originate from the photo-production of J/ψ.This contribution is negligible in pp collisions -from LHCb measurement at √ s = 7 TeV [41], it is O(10 −3 )% -but it is enhanced by a factor O(10 4 ) in Pb-Pb collisions, thus reaching the order of magnitude of the observed excess.The J/ψ coherent photo-production has been measured in ultra-peripheral Pb-Pb collisions [42].It is centered at very low p T , with ∼ 98% of these J/ψ below 0.3 GeV/c.An incoherent photo-production component is also observed in ultra-peripheral Pb-Pb collisions.About 30% of this contribution has a p T < 0.3 GeV/c, the rest being mainly located in the p T range 0.3-1 GeV/c.The influence of possible photo-production mechanisms on the inclusive J/ψ R AA presented in this paper has been evaluated by repeating the analysis placing a low p T threshold on the J/ψ of 0.3 GeV/c.Assuming that the observed excess in peripheral Pb-Pb collisions is indeed due to the photo-production of J/ψ, and that the relative contribution of the incoherent over coherent components is the same as the one estimated in ultra-peripheral collisions, then this selection would remove about 75% of the full photo-production contribution.Numerical values of R AA with the low p T threshold at 0.3 GeV/c are given in the Appendix A. All the figures and values presented in the paper refer to the inclusive J/ψ R AA but estimates of the difference between the inclusive and hadronic (without J/ψ photo-production) J/ψ R AA , are indicated where appropriate.
The results for the ψ(2S) analysis are given in terms of the ratio of their production cross sections (or, equivalently, of their production yields), expressed as When forming such a ratio the normalization factor N i events cancels out, as do most of the systematic uncertainties on A×ε corrections.The double ratio [ψ(2S)/J/ψ] Pb-Pb / [ψ(2S)/J/ψ] pp is used in order to directly compare the relative abundances of ψ(2S) and J/ψ in nucleus-nucleus and pp collisions.

Signal extraction
After a description of the muon selection procedure, we present here the two methods used to extract the J/ψ and ψ(2S) signals.The first one is directly based on fits of the µ + µ − invariant mass distribution while the second one makes use of the event mixing technique to subtract the combinatorial background.

Muon reconstruction
The muon reconstruction starts with the exclusion of parts of the detector that show problems during data taking such as high voltage trips, large electronic noise, pedestal determination issues.This selection is performed on a run-by-run basis to account for the time evolution of the apparatus.After pedestal subtraction, the adjacent well-functioning pads of both cathodes of each tracking chamber having collected a charge are grouped to form pre-clusters.These pre-clusters might be the superposition of several clusters of charges deposited by several particles crossing the detector close to each others.The number of clusters of charges contributing to the pre-cluster and their approximate location are determined with a Maximum √ s NN = 2.76 TeV ALICE Collaboration Likelihood -Expectation Maximization (MLEM) algorithm.It assumes that the charge distribution of each single cluster follows a two-dimensional integral of the Mathieson function [43].If the estimated number of clusters is larger than 3, the pre-cluster is split into several groups of 1, 2 or 3 clusters selected with the minimum total coupling to all the other clusters into the pre-cluster.Each group of clusters is then fitted using a sum of Mathieson functions, taking the MLEM results as a seed, to extract the precise location of where the particles crossed the detector.The overall spatial resolution is around 200 (550) µm in average in the (non-)bending direction.
The track reconstruction starts from the most downstream stations, where the multiplicity of secondary particles is smallest, by forming pairs of clusters in the two chambers of station 5(4), and deriving the parameters and associated errors of the resulting muon track candidates.The candidates are then extrapolated to the station 4 (5), validated if at least one compatible cluster is found in the station and duplicate tracks are removed.The procedure continues extrapolating the tracks to stations 3, 2 and 1, validating them by the inclusion of at least one cluster per station.The selection of compatible clusters is based on a 5σ cut on a χ 2 computed from the cluster and track local positions and errors.If several compatible clusters are found in the same chamber, the track is duplicated to consider all the possibilities and for each of them the track parameters and associated errors are recomputed using a Kalman filter.At each of the tracking steps, the track candidates, whose parameters indicate that they will exit the geometrical acceptance of the spectrometer in the next steps are removed.At the end of the procedure, the quality of the track is improved by adding/removing clusters based on a 4σ cut on the local χ 2 and fake tracks sharing clusters with others in the three outermost stations with respect to the interaction point are removed.The choice of the χ 2 cuts is a compromise between maximizing the tracking efficiency (< 1-2% muon rejection) and minimizing the amount of fake tracks (negligible background for this analysis).Finally, muon track candidates are extrapolated to the interaction vertex measured by the SPD taking into account the energy loss and the multiple Coulomb scattering in the front absorber.
An accurate measurement of the tracking chamber alignment is essential to reconstruct the tracks with enough precision to identify resonances in the µ + µ − invariant mass spectrum, especially the ψ(2S) for which the signal-to-background ratio is low.The absolute position of the chambers was first measured using photogrammetry before the data taking.Their relative position was then precisely determined using a modified version of the MILLEPEDE package [44], combining several samples of tracks taken with and without magnetic field.The small displacement of the chambers when switching on the dipole was measured by the Geometry Monitoring System (an array of optical sensors fixed on the chambers) and taken into account.The resulting alignment precision is ∼ 100 µm, leading to a reconstructed J/ψ invariant mass resolution of about 70 MeV/c 2 , and about 10% higher for the ψ(2S).The resolution is dominated by the energy loss fluctuation and multiple Coulomb scattering of the muons in the front absorber.More details on the muon spectrometer performances are given in [45].
In this analysis, the muon track candidates also have to fulfill the following requirements.First, the reconstructed track must match a trigger track with a p trig T above the threshold set in the MTR for triggering the event (1 GeV/c in this analysis).The trigger track is reconstructed from the average position of the fired strips on the two trigger stations, as computed by the trigger algorithm.The matching is based on a 4σ cut on a χ 2 computed from the tracker and trigger track parameters and errors including the angular dispersion due to the multiple Coulomb scattering of the muon in the iron wall.Second, the transverse radius coordinate of the track at the end of the front absorber must be in the range 17.6 < R abs < 89.5 cm.Muons exiting the absorber at small and large angles, thus outside the R abs cut range, have crossed a different amount of material, either the beam shield or the envelope of the absorber, affecting the precision of the energy loss and multiple Coulomb scattering corrections.Third, in order to remove muon candidates close to the edge of the spectrometer acceptance, a cut on the track pseudo-rapidity −4 < η < −2.5 is applied.
J/ψ and ψ(2S) production in Pb-Pb collisions at √ s NN = 2.76 TeV ALICE Collaboration 5.2 J/ψ signal J/ψ candidates are formed by combining pairs of opposite-sign tracks reconstructed within the geometrical acceptance of the muon spectrometer.The aforementioned cuts at the single muon track level remove most of the hadrons escaping from or produced in the front absorber, as well as a large fraction of low p T muons from pion and kaon decays, secondary muons produced in the front absorber, and fake tracks.The J/ψ peak becomes visible in the µ + µ − invariant mass spectra even before any background subtraction.At the dimuon level only cuts on rapidity (2.5 < y < 4) and transverse momentum (p T < 8 GeV/c) are applied.The J/ψ raw yields are extracted by using two different methods.
In the first method, the opposite-sign dimuon invariant mass distribution is fitted with a sum of two functions.The signal is described by a double-sided Crystal Ball function (CB2).This function is an extension of the Crystal Ball function [46], i.e. a Gaussian with a power-law tail in the low mass range, with an additional independent power-law tail in the high mass range.The CB2 function reproduces very well the J/ψ line shape in the Monte Carlo (MC) simulations.The underlying continuum is described by a variable width Gaussian function.This function is built on a Gaussian form, whose width is dependent on the invariant mass of the dimuon.It was checked that including or excluding a ψ(2S) contribution in the fitting procedure has a negligible effect on the extracted J/ψ yield within the present statistical and signal-extraction-related systematic uncertainties.Since the significance of the ψ(2S) signal in the centrality, p T and y intervals used for the J/ψ analysis is too small to extract its contribution, we do not include it in the fit for this analysis.
During the fitting procedure, the width of the J/ψ peak is kept as a free parameter as it cannot be reproduced perfectly in simulations, and its value varies from 65 to 76 MeV/c 2 (larger than those from MC by about 5-10%).The pole mass is also kept free although the differences observed between data and simulation are at the per mille level.The tail parameters cannot be constrained by the fit.Therefore they are fixed to values obtained from an embedding simulation (described in section 6) and adjusted for each p T and y interval under study in order to account for the observed dependence on the J/ψ kinematics.On the contrary, the J/ψ shape does not show a dependence on centrality, hence the CB2 tail parameters tuned on a centrality integrated MC sample are used in all the bins.Figure 1 presents fits of the opposite-sign dimuon invariant mass (m µ µ ) distributions for different p T ranges in central (top row) and peripheral (bottom row) collisions.
The signal-to-background ratio (S/B) and the significance (S/ √ S + B) of the signal are evaluated within 3 standard deviations with respect to the J/ψ pole mass.The S/B varies from 0.2 to 6.5 when going from the most central collisions to the most peripheral ones.Integrated over centrality and y (p T ), the S/B ranges from 0.2 (0.2) to 1.2 (0.6) with increasing p T (y).In all the centrality, p T or y intervals considered in this analysis, the significance is always larger than 8.
In the second method, the combinatorial background is subtracted using an event-mixing technique.The opposite-sign muon pairs from mixed-events are formed by combining muons from single muon low p T (MSL) triggered events.In order to limit the effect of efficiency fluctuations between runs and to take into account the dependence of muon multiplicity and kinematic distributions on the collision centrality, events in the same run and in the same centrality class are mixed together.The mixed-event spectra are normalized to the data using the combination of the measured like-sign pairs such as where pairs.The R factor in Eq. 5 is defined by and is introduced in order to correct for differences in acceptance between like-sign and opposite-sign muon pairs.Above a dimuon invariant mass of 1.8 GeV/c 2 , the R factor is equal to unity with deviations smaller than 1%.The accuracy of the mixed-event technique was assessed by comparing the distributions of likesign muon pairs from mixed-events to the same-event ones, which agree within 1% over the mass, p T and y ranges under study.This agreement justifies the use of the normalization given by Eq. 5, which implies that the correlated signal in the like-sign dimuon spectra is negligible with respect to the combinatorial background.The mass spectra of the opposite-sign mixed-event pairs are then subtracted from the data.The resulting background-subtracted spectra are fitted following the same procedure as in the first method, except that the variable width Gaussian function is replaced by an exponential function to account for residual background.Figure 2 shows fits of the background-subtracted opposite-sign dimuon invariant mass distributions for different p T ranges in central (top row) and peripheral (bottom row) collisions.

ψ(2S) signal
The invariant mass spectra used to extract the [ψ(2S)/J/ψ] ratio are obtained in the same way as described in the previous section, implementing the same cuts applied at the muon and dimuon levels.In order to improve the significance of the ψ(2S) signal, a wider centrality and p T binning than the one used for the J/ψ analysis was adopted, and the analysis is performed integrated over the full rapidity domain 2.5 < y < 4. The fits to the invariant mass spectra are performed by modeling the ψ(2S) signal with a CB2 function.Given the very low S/B ratio, the normalization is chosen as the only free parameter for ψ(2S).mass is fixed to the one of the J/ψ, shifted by the corresponding ∆m = m ψ(2S) − m J/ψ value taken from the PDG [47].The width of the ψ(2S) is fixed to the one of the J/ψ scaled by the ratio σ ψ(2S) /σ J/ψ estimated from MC simulations.
Fits of the invariant mass spectra showing the ψ(2S) are visible in Fig. 3 for the p T < 3 GeV/c interval in the centrality classes 20-40%, 40-60% and 60-90%.For the other intervals in centrality and p T the ψ(2S) signal could not be extracted, i.e. the ratio [ψ(2S)/J/ψ] is consistent with zero.In these cases, only the 95% confidence level upper limit is computed.

Acceptance and efficiency correction
In the J/ψ analysis, embedding simulations are used to compute the centrality, p T and y dependences of the acceptance times efficiency factor (A×ε).The Monte Carlo embedding technique consists of adding the detector response from a simulated signal event (charmonium in this case) to a real data event, and then performing the reconstruction as for real events.This has the advantage of providing the most realistic background conditions, which is necessary for Pb-Pb collisions where high multiplicities are reached: at η = 3.25, dN ch /dη ≈ 1450 for the 0-5% most central events [48].This leads to a large detector occupancy, which can reach about 3% in the most central collisions and alter the track reconstruction efficiency.
Monte Carlo J/ψ were embedded in MB triggered events recorded in parallel to the opposite-sign dimuon triggered events.Only one J/ψ was simulated per event at the position of the real event primary vertex reconstructed by the SPD.The shapes of the input MC p T and y distributions were tuned to match the measured distribution in Pb-Pb collisions (see discussion in section 7.2).The muons from the J/ψ decay were then transported through a simulation of the ALICE detector using GEANT3 [49].The detector simulated response was then merged with that of a real Pb-Pb event and the result was processed by the normal reconstruction chain.Embedding simulations were performed on a run-by-run basis to account for the time-dependent status of the tracking chambers.The residual misalignment of the detection elements, Fig. 3: Opposite-sign dimuon invariant mass distribution for the 20-40%, 40-60% and 60-90% centrality classes, for 2.5 < y < 4 and p T < 3 GeV/c, before (top row) and after background subtraction (bottom row) via event mixing.In these intervals the ψ(2S) signal is extracted whereas in all other centrality and p T intervals, only the 95% confidence level upper limits are provided.
whose amplitude is evaluated by analyzing the residual distance between the clusters and the tracks in data, was also taken into account.For the trigger chambers, the efficiency maps measured in data were used in the simulations.
The left panel of Fig. 4 shows the J/ψ A×ε as a function of collision centrality in the rapidity domain 2.5 < y < 4 and in the p T range p T < 8 GeV/c.We observe a relative decrease of 8% of the J/ψ reconstruction Centrality (%)  efficiency from the 80-90% centrality class to the 0-10% centrality class.This decrease is mostly due to a drop of about 3% of the single muon trigger efficiency in the most central collisions whereas the decrease of the single muon tracking efficiency is only on the order of 1%.When considering specific p T or y intervals, a maximum relative variation of ∼ 30% of the A×ε decrease with centrality is observed.The right panel of Fig. 4 shows the p T versus y dependence of A×ε.The rapidity dependence of A×ε reflects the geometrical acceptance of the muon pairs with a maximum centered at the middle of the rapidity interval and a decrease towards the edges of the acceptance.The p T dependence of A×ε is non-monotonic, with a minimum at p T ≈ 1.8 GeV/c corresponding to J/ψ kinematics for which one of the decay muons does not fall into the muon spectrometer acceptance.
For the ψ(2S) resonance, the embedding technique was not used.Since, in this case, only the ratio [ψ(2S)/J/ψ] is extracted, the A×ε correction factors for both resonances were evaluated through pure signal MC simulations, assuming that the dependence of the efficiency as a function of the centrality is the same for J/ψ and ψ(2S), and therefore cancels out in the ratio.The effect of possible differences in the centrality dependence of A×ε was studied and included as a source of systematic uncertainty.

Systematic uncertainties
In the following, each source of systematic uncertainty is detailed.Most of them affect the J/ψ and ψ(2S) results identically and vanish in the [ψ(2S)/J/ψ] ratio.Systematic uncertainties specific to the ψ(2S) analysis are explicitly mentioned.

Signal extraction
The systematic uncertainty on the signal extraction results from several fits of the invariant mass spectra, where signal line shape parameters, background description and fit range are varied as detailed below.In each centrality, p T and y intervals, the raw yield and the statistical uncertainty are given by the average of the results obtained from the different fits.The corresponding systematic uncertainty is defined as the RMS of these results.It was also checked that every individual result differs from the mean value by less than three RMS.
The J/ψ line shape is well described by the CB2 function, whose pole mass and width are constrained by the data while the tail parameters have to be fixed to values extracted from the embedding simulation.Alternatively, another set of tails was extracted from pp data, where a large statistics and a better S/B were available.In this case, the p T and y dependence of the tail parameters could not be determined with sufficient precision, so the same values were used for all p T and y intervals.In the event mixing approach, the influence of different normalizations of the opposite-sign mixed-event spectrum to the opposite-sign sameevent spectrum was investigated.We have tested a normalization performed on a run-by-run basis or after merging of all the runs, and a normalization based on the integral of the invariant mass spectrum in the intermediate mass region (1.5 < m µ µ < 2.5 GeV/c 2 ).None of these tests showed deviations larger than 1% in the number of extracted J/ψ, and thus were not included in the tests used to extract the systematic uncertainty on the signal extraction.The fit range of the invariant mass spectra was also varied considering a narrow (2.3 < m µ µ < 4.7 GeV/c 2 ) and a wide (2 < m µ µ < 5 GeV/c 2 ) interval.Finally, all the combinations of signal line shape, background description (with or without using the event-mixing technique) and fit range are performed to account for possible correlations.
The same procedure as above was applied when the ψ(2S) signal was included in the fit function for the specific centrality and p T intervals presented in this analysis.To account for the fact that the ψ(2S) width was fixed to the one of the J/ψ scaled by the ratio σ ψ(2S) /σ J/ψ estimated from MC simulations, all the fits were repeated varying the scaling factor by ±10%.This variation accounts for the fluctuations observed in pp data when fitting the invariant mass spectra leaving the width of the ψ(2S) free or fixing it as described above.

Monte Carlo input parametrization
The estimation of A×ε factors depends on the charmonium p T and y shapes used as input distributions in the MC simulation.In order to evaluate the sensitivity of the results on this choice, several MC simulations were performed, each one including modified p T and y distributions.For the J/ψ, the modification of the shapes was done in order to take into account the possible correlation between p T and y (as observed by LHCb in pp collisions [39]) and the correlation between p T (y) and the centrality of the collision (as reported in this paper).A systematic uncertainty of 3% is found for A×ε integrated over p T and y and is taken as correlated as a function of the centrality.The p T (y) dependence of this uncertainty varies in the range 0-1% (3-8%).The larger effect seen in the y dependence occurs at the low and high limits, where the acceptance falls steeply.
The same procedure was followed for the ψ(2S), assuming that the correlations between p T and y and with the centrality are of the same magnitude as those observed for the J/ψ.A systematic uncertainty of 2% is evaluated for the [ψ(2S)/J/ψ] ratio in the p T < 3 GeV/c interval.

Centrality dependence of the [ψ(2S)/J/ψ] A×ε
The embedding technique was not used for the ψ(2S) MC simulations as we have assumed the same A×ε dependence as a function of the centrality for the ψ(2S) and the J/ψ.In order to evaluate the systematic uncertainty introduced by this assumption, a conservative ±30% variation of the A×ε loss as a function of centrality was applied to the ψ(2S).This corresponds to the maximum variation of the A×ε loss between peripheral and central collisions observed for the J/ψ in different p T and y intervals.The effect on the (A×ε) J/ψ / (A×ε) ψ(2S) ratio is 1% or lower in all the centrality classes considered.Since this effect is much smaller than the systematic uncertainty on the signal extraction, it is neglected.

Tracking efficiency
The tracking algorithm, as described in section 5.1, does not require all the chambers to have fired to reconstruct a track.This redundancy of the tracking chambers can be used to measure their individual efficiencies from data, and since such efficiencies are independent from each other, we can combine them to assess the overall tracking efficiency.This evaluation of the tracking efficiency is not precise enough to be used to directly correct the data, because only the mean efficiency per chamber can be computed with the statistics available in each run.However, by comparing the result obtained from data with the same measurement performed in simulations, we can control the accuracy of these simulations and assess the corresponding systematic uncertainty on the A×ε corrections.
A 9% relative systematic uncertainty is obtained for the J/ψ by comparing the measured tracking efficiency in simulations and in peripheral Pb-Pb collisions.This uncertainty is constant and fully correlated as a function of centrality.From low to high p T (y), the systematic uncertainty varies from 9% to 7% (7% to 6% with a maximum of 12% at y ≃ 3.25).On top of that, a small difference was observed in the centrality dependence of this measurement between data and embedding simulations.This results in an additional 1% systematic uncertainty in the 0-10% centrality class and 0.5% in 10-20%.
Another systematic uncertainty can arise from correlated dead areas located in front of each other in the same station, which cannot be detected with the method detailed above.A dedicated study has shown that this effect introduces a 2% systematic uncertainty, fully correlated as a function of centrality and predominantly uncorrelated as a function of p T and y.
In the [ψ(2S)/J/ψ] ratio the systematic uncertainty on the tracking efficiency largely cancels out because the ψ(2S) and J/ψ decay muons have similar p T and y distributions and, therefore, cross about the same regions of the detector.Since the possible remaining systematic uncertainty is much smaller than that on the signal extraction, it is neglected in this analysis.

Trigger efficiency
The systematic uncertainty on the J/ψ A×ε corrections related to the trigger efficiency has two origins: the intrinsic efficiency of the trigger chambers and the response of the trigger algorithm.The first part was determined from the uncertainties on the trigger chamber efficiencies measured from data and applied to simulations.Propagating these efficiencies in J/ψ simulations results in a 2% systematic uncertainty on the A×ε corrections, fully correlated as a function of centrality and mainly uncorrelated as a function of p T and y.The effect of the systematic uncertainty on the shape of the trigger response as a function of the muon p T was determined by weighting MC J/ψ decay muons with different trigger response functions obtained in data and simulations.These functions were defined as the fraction, versus p T , of the single muons passing a 0.5 GeV/c p trig T threshold that also satisfy the 1 GeV/c p trig T threshold used in this analysis.The resulting systematic uncertainty on the J/ψ A×ε correction integrated over p T and y is 1%.As a function of p T , it amounts to 3% for p T < 1 GeV/c and 1% elsewhere.As a function of y, a 1% uncorrelated systematic uncertainty was obtained.These uncertainties are fully correlated as a function of centrality.
The systematic uncertainty on the modification of the trigger response as a function of centrality, i.e. for increasing multiplicity, was assessed by changing the detector response (space size of the deposited charge) to the passage of particles in embedding simulations.The corresponding uncertainties on the J/ψ A×ε corrections are 1% in the 0-10% and 10-20% centrality classes, and 0.5% in 20-30% and 30-40%.
As for the case of tracking efficiency, this source of systematic uncertainty largely cancels out in the [ψ(2S)/J/ψ] ratio and is neglected.

Matching efficiency
The systematic uncertainty on the matching efficiency between the tracking and the trigger tracks is 1%.It is given by the differences observed between data and simulations when applying different χ 2 cuts on the matching between the track reconstructed in the tracking chambers and the one reconstructed in the trigger chambers.This uncertainty is fully correlated as a function of the centrality and largely uncorrelated as a function of p T and y.
Also in this case, the effect on the [ψ(2S)/J/ψ] ratio is negligible.

pp reference
The statistical and systematic uncertainties on the measurement of the J/ψ differential cross section in pp collisions at √ s = 2.76 TeV are available in [37].The statistical uncertainty is combined with that of the Pb-Pb measurement when calculating the R AA as a function of p T and y, but is considered as a fully correlated systematic uncertainty as a function of the centrality.The correlated and uncorrelated part of the systematic uncertainty on the pp reference as a function of p T and y are both fully correlated as a function of the centrality.
The ψ(2S) statistics in the √ s = 2.76 TeV pp data sample are too low to be used for the normalization of the [ψ(2S)/J/ψ] Pb-Pb ratio.For this reason, pp results obtained at higher energy ( √ s = 7 TeV) [50] were used, thus introducing an additional source of systematic uncertainty.An interpolation procedure, as the one described in [33], was applied in order to extract the [ψ(2S)/J/ψ] pp ratio at √ s = 2.76 TeV.The discrepancy between the result of this interpolation in the kinematic range p T < 3 GeV/c 2.5 < y < 4 and the value obtained at √ s = 7 TeV is 10%: this relative difference is included in the systematic uncertainty on the pp reference.

Normalization
The systematic uncertainty on the normalization is the one attached to the scaling factor F norm and amounts to 4%.This value corresponds to one standard deviation of the distribution of the F norm computed for each run used in the analysis.This systematic uncertainty is fully correlated as a function of the centrality, p T and y.

Others
Systematic uncertainties on the nuclear overlap function T AA are available in Tab. 1.Another systematic uncertainty on the definition of the centrality classes arises from the V0 amplitude cut, which corresponds to 90% of the hadronic cross section [36].A maximum uncertainty of 5% is obtained in the centrality class (80-90%) vanishing with increasing centrality or in wider centrality classes.Systematic uncertainties due to the unknown polarization of the J/ψ are not propagated and we assume that J/ψ production is unpolarized both in pp and in Pb-Pb collisions.In pp collisions at √ s = 7 TeV, J/ψ polarization measurements at mid-rapidity (p T > 10 GeV/c) and forward-rapidity (p T > 2 GeV/c) are compatible with zero [51][52][53].In Pb-Pb collisions, J/ψ mesons produced from initial parton-parton hard scattering are expected to have the same polarization as in pp collisions and those produced from charm quarks recombination in the medium are expected to be unpolarized.

Summary
The systematic uncertainties related to the J/ψ analysis are summarized in Tab. 2. Concerning the ψ(2S) analysis, most of the systematic uncertainties cancel out in the [ψ(2S)/J/ψ] ratio and the main contributors are the signal extraction (14-45%) and the pp reference (10%).

Inclusive J/ψ mean transverse momentum
The p T dependence of the J/ψ yields per MB collision, defined by Eq. 1, was studied for three centrality classes (0-20%, 20-40% and 40-90%) and is displayed in Fig. 5.The statistical uncertainties appear as vertical lines.The systematic uncertainties uncorrelated as a function of p T are shown as open boxes, while the ones fully correlated as a function of p T but uncorrelated as a function of centrality are shown as shaded areas (mostly hidden by the points).The global systematic uncertainty, fully correlated as a function of centrality and p T , is quoted directly in the figure.Numerical values for the J/ψ yields can be found in Appendix A. The inclusive J/ψ mean transverse momentum was computed by fitting the p T distribution of  inclusive J/ψ yields with the function where C, p 0 and n are free parameters.This function is commonly used to reproduce the J/ψ p T distribution in hadronic collisions, see for instance [54][55][56].Fit results for the three centrality classes are displayed as full lines in the figure.An excess over this function is revealed in the lowest p T interval (corresponding to 0 < p T < 500 MeV/c) for peripheral Pb-Pb collisions.It could be caused by a residual contribution from J/ψ coherent photo-production, which was measured in ultra-peripheral collisions [42].A quantitative measurement of this contribution in hadronic collisions is reported in [40].Thus, in the most peripheral centrality class (40-90%) the fit was performed for p T > 500 MeV/c and extrapolated down to zero (dotted line).In the 0-20% and 20-40% centrality classes, no J/ψ excess was observed and fits were performed down to zero p T .As a cross-check, the same procedure as for the peripheral centrality class was tested and the obtained results are fully compatible within uncertainties.
Values of the mean transverse momentum ( p T ) and mean squared transverse momentum ( p 2 T ) obtained from the fits are given in Tab. 3 as a function of centrality.The statistical (systematic) uncertainty is extracted by fitting the p T distribution considering only the statistical (p T -uncorrelated systematic) uncertainty of the measurement.For comparisons, the p T and p 2 T results from PHENIX were recomputed with the function defined by Eq. 7, adjusted in the measured p T range and extrapolated to p T = 8 GeV/c to match our p T range.These results are also given in Tab. 3 along with the measurement in pp collisions at √ s = 2.76 TeV with updated uncertainties [57].In order to compare the evolution of p 2 T A-A at different energies, one can form the variable r AA defined as This variable was measured over the wide range of energies and colliding systems covered by NA50 and PHENIX experiments.The comparison with the ALICE results is done in Fig. 6 (right side).A very different N part dependence is seen, especially when comparing Pb-Pb collisions at the SPS and the LHC.At the SPS energy of √ s NN = 0.017 TeV [59], the increase of the J/ψ p 2 T with the centrality of the collision was attributed to the Cronin effect [61], interpreted as an extra p T kick due to multiple scatterings of the initial partons producing the J/ψ.At the LHC, a clear decrease of r AA is observed as a function of N part .This behavior could be related to the onset of recombination phenomena and to the thermalization of charm quarks.Theoretical calculations [13,60], based on transport models (described in the next section) are able to reproduce the r AA at SPS, RHIC and LHC energies.They correlate the specific dependence of r AA on collision centrality with the increased importance of recombination effects in the J/ψ production mechanism at the LHC.T at various energies and colliding systems.The statistical and systematic uncertainties are quoted separately, except for PHENIX measurements in Au-Au collisions where the quadratic sum is given.If the measurement is not available or not used in the range 0 < p T < 8 GeV/c, the fit function is extrapolated down to 0 and up to 8 GeV/c to compute p T and p 2 T .

Nuclear Modification Factor
Some of the R AA results presented here were already published in [27] and are shown again in this section, where they are compared with model calculations and with results from previous experiments.They include the centrality dependence of R AA (Fig. 7), the p T dependence of R AA for the full centrality range 0-90% and for the centrality class 0-20% (Fig. 9 top row) and the rapidity dependence of R AA (Fig. 10).The new results shown in this section include the centrality dependence of R AA for three p T intervals (Fig. 8) and the p T dependence of R AA for the centrality classes 20-40% and 40-90% (Fig. 9 bottom row).These new results were obtained using a slightly different trigger selection, as explained in section 3. The consistency of the results obtained with the two selections was verified.Fig. 6: Mean transverse momentum p T measured by ALICE [37,57] and PHENIX [21,55,58] as a function of the number of participant nucleons (left).r AA measured by NA50 [59], PHENIX and ALICE and compared to model calculations [13,60], as a function of the number of participant nucleons (right).
Tab. 1) is R 0-90% AA = 0.58 ± 0.01(stat.)± 0.09(syst.),indicating a clear J/ψ suppression.This suppression is significantly less pronounced than that observed at lower energy in PHENIX in a similar kinematic range, as previously discussed in [26,27].For N part larger than 70, corresponding to the 50% most central Pb-Pb collisions, the J/ψ R AA is consistent with being constant, within uncertainties.Such behavior was not observed in heavy ion collisions at lower energies (SPS, RHIC), where R AA is continuously decreasing as a function of centrality.
The impact of non-prompt J/ψ on the inclusive R AA analysis was studied.The R AA of prompt J/ψ is estimated (see Eq. 3) to be about 7% larger than the inclusive J/ψ R AA if the beauty component is fully suppressed.In the other extreme case, where the B-meson production is not affected by the medium and scales with the number of binary collisions, i.e.R non-prompt AA = 1, the R AA of prompt J/ψ would be about 6% smaller in central collisions and about 1% smaller in peripheral collisions.The excess in the inclusive J/ψ yield observed at very low p T [40] also influences the R AA in the most peripheral collisions.A large fraction of this contribution (about 75% as explained in section 4) can be removed by selecting J/ψ with a p T higher than 0.3 GeV/c.Assuming that the hadronic J/ψ R AA in the ranges 0 < p T < 0.3 GeV/c and 0.3 < p T < 8 GeV/c are the same, it becomes possible to estimate the impact of the J/ψ photo-production on the inclusive R AA .In the centrality classes 60-70%, 70-80% and 80-90%, the hadronic J/ψ R AA would be about 5%, 11% and 25% lower, respectively.Extreme hypotheses were made to define upper and lower limits, represented with brackets on the Figs. 7, 8 and 9.The upper limit calculation assumes no J/ψ from photo-production thus the inclusive measurement only contains hadronic production.The lower limit assumes that i) all J/ψ produced with a p T smaller than 0.3 GeV/c originate from photo-production and ii) the efficiency of the 0.3 GeV/c p T selection is reduced from 75% to 60% (corresponding to an increase by a factor two of the J/ψ photo-production above 0.3 GeV/c).
The comparison with theoretical models, shown on the right-hand side of Fig. 7, helps in the interpretation of the large difference observed between the PHENIX and the ALICE results.
The Statistical Hadronization Model (SHM) [62] assumes deconfinement and thermal equilibration of the bulk of the cc pairs.Charmonium production occurs at the phase boundary via the statistical hadronization of charm quarks.The prediction is given for two values of the charm cross section dσ cc /dy = 0.15 and 0.25  [13,60,62,63], which all include a J/ψ regeneration component (right).The brackets shown in the three most peripheral centrality classes on the right figure quantify the possible range of variation of the hadronic J/ψ R AA for two extreme hypotheses on the photo-production contamination in the inclusive measurement, see text for details.mb at forward rapidity.These values are derived from the measured charm cross section in pp collisions at √ s = 2.76 and 7 TeV [15] bracketing the expectation for gluon shadowing in the Pb-nucleus between 0.6 and 1.0.Production of non-prompt J/ψ from decays of B-mesons is not considered.
The two transport models from Zhao (TM1) [13] and Zhou (TM2) [60] mainly differ in the rate equation controlling the J/ψ dissociation and regeneration.In TM1, shadowing is implemented via a simple parametrization, leading to a 30% suppression in the most central Pb-Pb collisions.The charm cross section is assumed to be dσ cc /dy ≈ 0.5 mb at forward rapidity, the fraction of J/ψ from beauty hadrons to be 10% and no b-quenching is introduced in the calculation.This model is presented as a band connecting the results obtained with (lower limit) and without (upper limit) shadowing and is interpreted by the authors as the uncertainty of the prediction.In TM2, the shadowing is given by the EKS98 parametrization [64].The charm cross section is taken in the range dσ cc /dy ≈ 0.4 − 0.5 mb at forward rapidity; the calculations for these two values provide the lower and upper limits of the band displayed in the figure.The fraction of J/ψ from beauty hadrons is assumed to be 10% with a b-quenching of 0.8, increased to 0.4 for p T above 5 GeV/c.The Comover Interaction Model (CIM) [63] implements shadowing, interaction with a co-moving dense partonic medium and recombination effects.The shadowing is calculated within the Glauber-Gribov theory making use of the generalized Schwimmer model of multiple scattering.The J/ψ dissociation cross section due to comover interaction is taken as σ co = 0.65 mb from low-energy data.Recombination effects are included by adding a gain term proportional to σ co and to the number of c and c quarks, thus no additional parameter is added to the model.The charm cross section dσ cc /dy at forward rapidity is taken in the range 0.4 to 0.6 mb, which gives respectively the lower and upper limits of the calculation.Production of nonprompt J/ψ is not considered.
To match our J/ψ R AA results, all models above need to include in their calculation a sizeable J/ψ production from deconfined c and c quarks. ) and comparisons of the lowest and highest p T range to the transport and to the comover interaction models [13,60,63].The brackets quantify the possible range of variation of the hadronic J/ψ R AA for two extreme hypotheses on the photo-production contamination in the inclusive measurement.
A different test of these models was carried out by studying the J/ψ R AA centrality dependence in p T intervals.Figure 8 displays the measurement of the inclusive J/ψ R AA as a function of the number of participant nucleons measured in Pb-Pb collisions at √ s NN = 2.76 TeV for the three p T ranges 0-2, 2-5 and 5-8 GeV/c.The uncorrelated systematic uncertainties shown at each point were separated into uncorrelated as a function of centrality (open boxes) and fully correlated as a function of centrality but uncorrelated as a function of p T (shaded areas).For N part 150, the low p T J/ψ R AA is significantly larger than the mid and high p T ones.In the most central bin, the R AA values corresponding to the lowest and the highest p T are separated by 3.9σ .For N part 150, the centrality dependence exhibits similar trends for the 2-5 and 5-8 GeV/c ranges, while the most peripheral ( N part ∼ 20) R AA measurement in the low p T (0-2 GeV/c) range appears to deviate from the others.However, the J/ψ yield excess observed at very low p T may have a sizable effect in the 0-2 GeV/c interval.In the centrality classes 40-50%, 50-60% and 60-90%, based on the same assumptions made for the 0 < p T < 8 GeV/c case, the hadronic J/ψ R AA would be about 5%, 6% and 18% lower, respectively.Due to the increase of the non-prompt J/ψ component at large p T , the difference between the measured inclusive J/ψ R AA and the prompt J/ψ R AA increases with p T .If the beauty contribution is fully (not) suppressed, R AA of prompt J/ψ is estimated to be 6%, 8% and 11% larger (0-3%, 3-10% and 7-30% smaller, depending on centrality) for the p T ranges 0-2, 2-5 and 5-8 GeV/c, respectively.
Calculations from the transport models and the comover interaction model are plotted on top of the results shown in Fig. 8.For the most peripheral collisions ( N part 100), the models cannot correctly reproduce the R AA centrality dependence for both the low and high p T ranges.For the most central collisions ( N part 100), the R AA centrality dependence for high p T J/ψ is reasonably reproduced by all models.Concerning the low p T range in the most central events, the measurement is compatible with the upper side of the theoretical uncertainty band from the CIM and TM2 models.For these models, it corresponds to the highest value for dσ cc /dy, 0.6 and 0.5 mb respectively.

Transverse momentum dependence of R AA
The p T dependence of the inclusive J/ψ R AA in the rapidity range 2.5 < y < 4 is shown in Fig. 9 for the full centrality range 0-90% and for three centrality classes 0-20% [27], 20-40% and 40-90%.In Fig. 9 top left corner, the inclusive J/ψ R AA in the centrality class 0-90% shows a decrease of about 50% from low to high p T .At low p T , the measurement is close to 0.8 showing very little suppression.At high p T , our R AA value is similar to that of CMS [25].They measured, in the different rapidity range 1.6 < |y| < 2.4, an inclusive J/ψ R AA = 0.41 ± 0.05 ± 0.04 for 3 < p T < 30 GeV/c.The corresponding mean p T is 6.27 GeV/c.When beauty contribution is fully (not) suppressed the prompt J/ψ R AA is expected to be 5% larger (2% smaller) for p T < 1 GeV/c and 17% larger (30% smaller) for 6 < p T < 8 GeV/c.The transport model calculations TM1 [13] and TM2 [60] are also shown in Fig. 9.Both models reproduce reasonably well the 0-90% centrality measurement at high p T .At low p T , TM1 reproduces rather well our measurement, while the data points sit on the upper limit of the TM2 calculation.One can also appreciate the relative contributions of the primordial (from the initial hard parton scattering) and regenerated (from coalescence of c and c quarks in the deconfined medium) components in these two calculations.The contribution of regenerated J/ψ is concentrated at low p T and its relative fraction with respect to the initial production differs between the models.In TM1, it is of the same order of the primordial J/ψ production, which is about constant over the full p T range.In TM2, the regenerated J/ψ contribution is almost three times larger than the primordial one in the lowest p T interval.For p T > 5 GeV/c, only the primordial production remains.
In the other panels of Fig. 9, the ALICE measurements are compared to those from PHENIX in Au-Au collisions at √ s NN = 0.2 TeV for the 0-20%, 20-40% and 40-60% centrality classes [21].For p T < 1 GeV/c, for all centrality ranges, the prompt J/ψ R AA is expected to be 5% larger (2% smaller) when the beauty contribution is fully (not) suppressed.For 6 < p T < 8 GeV/c the effect is much larger: if the beauty contribution is fully suppressed, the prompt J/ψ R AA would be 17% larger in all centrality ranges.If the beauty contribution is not suppressed, the prompt J/ψ R AA would be 44%, 15% and 8% smaller in the centrality ranges 0-20%, 20-40% and 40-90%, respectively.The very low p T excess in the inclusive J/ψ yield mentioned before has a non-negligible impact in the 0 < p T < 1 GeV/c range in the most peripheral centrality class 40-90%: following the same assumptions made for the 0 < p T < 8 GeV/c case, the estimated hadronic J/ψ R AA would be about 22% lower.In the most central collisions (0-20%), the inclusive J/ψ R AA at low p T is almost four times larger in Pb-Pb collisions at √ s NN = 2.76 TeV than in Au-Au collisions at √ s NN = 0.2 TeV.This difference cannot be explained only by the possible change in the size of the CNM effects that can be expected due to the different rapidity coverage and collision energy between the two measurements.Such a behavior, on the other hand, is expected by all the recombination models described in the previous section.The same trend is observed in the centrality class 20-40%, where the large difference between the PHENIX and ALICE results observed at low p T vanishes at high p T .Concerning the Fig. 9: Inclusive J/ψ R AA as a function of the J/ψ p T for 2.5 < y < 4 in the centrality class 0-90% [27] compared to transport models [13,60] (top left).The comparison is done with PHENIX results [21] and transport models in the 0-20% [27] (top right), 20-40% (bottom left) and 40-90% (bottom right) centrality classes.The brackets shown in the lowest p T interval for the centrality class 40-90% quantify the possible range of variation of the hadronic J/ψ R AA for two extreme hypotheses on the photo-production contamination in the inclusive measurement.Upper limits from PHENIX at high p T are not represented.most peripheral collisions, the inclusive J/ψ R AA is still slightly larger for ALICE results at low p T .However, here the comparison between the two experiments is done with different centrality classes, 40-90% (ALICE) and 40-60% (PHENIX), so that a firm conclusion, also because of the uncertainty size, cannot be drawn.Transport model calculations for Pb-Pb collisions at √ s NN = 2.76 TeV are also presented for the 0-20%, 20-40% and 40-90% centrality classes.TM1 shows a good agreement with the measurements in the 0-20% and 20-40% centrality classes, while TM2 tends to underestimate the data for p T < 5 GeV/c.In the most peripheral centrality class (40-90%), the two models follow significantly different trends, but the uncertainties from the measurement are too large to discriminate them.However, if the very low p T excess is taken into account, a rather flat p T dependence of the R AA is expected, pushing our measurement aside from TM1 calculations in this specific range.For the high p T region, both models reproduce well the experimental results in all the centrality classes.

Rapidity dependence of R AA
The rapidity dependence of the inclusive J/ψ R AA in Pb-Pb collisions [27] is shown in Fig. 10.The inclusive J/ψ R AA measured in the rapidity range |y| < 0.8 is about 0.7, consistent with the value measured at y ∼ 3.
From y ∼ 3 to y ∼ 4, the J/ψ R AA shows a decreasing trend leading to a drop of about 40%.The influence of non-prompt J/ψ on this result is small, as the prompt J/ψ R AA is expected to be only 8% larger (5% smaller) for 2.5 < y < 2.75 and 6% larger (9% smaller) for 3.75 < y < 4 if the beauty contribution is fully (not) suppressed.The Pb-Pb measurements are compared to theoretical calculations, which only consider shadowing and coherent energy loss.The break-up of the cc pair and nuclear absorption are not taken into account in any of the models.Shadowing only predictions are made within the Color Singlet Model at Leading Order [65] and the Color Evaporation Model at Next to Leading Order [66], with the EKS98 [64] and the EPS09 [68] parametrizations of the nPDF, respectively.For EKS98 (EPS09) the upper and lower limits correspond to the uncertainty in the factorization scale (uncertainty of the nPDF).Finally, a theoretical prediction, which includes a contribution from coherent parton energy loss processes in addition to EPS09 shadowing [67] is also shown.All models show a fair agreement with our measurements over a wide rapidity range, |y| 3.If the amplitude of CNM effects is correctly given by the calculations shown in Fig. 10, the observed J/ψ suppression due to CNM effects could be as large as 40%.Moreover, if an additional J/ψ suppression occurs in the hot nuclear matter (as expected from lower energy measurement and observed at high p T by both CMS and ALICE), other mechanisms compensating this suppression are needed to explain the R AA measurements.Figure 10 supports this scenario, where suppression effects in hot matter are qualitatively counterbalanced by recombination.This is indeed what is expected from all models featuring recombination discussed in this paper.At higher rapidity, for |y| 3, the models implementing only CNM effects tend to deviate from the data, although the one combining shadowing with coherent energy loss seems to match the decreasing trend of the R AA better.Such a decrease of the R AA values can also be explained by recombination models, where a reduction of the recombination effects is expected with increasing rapidity, due to the decrease of dσ cc /dy.

[ψ(2S)/J/ψ] ratio
The ratio between inclusive ψ(2S) and J/ψ yields measured in Pb-Pb collisions at √ s NN = 2.76 TeV is shown in the left side of Fig. 11 as a function of N part .In the interval p T < 3 GeV/c, the ψ(2S) signal was extracted in three centrality classes (20-40%, 40-60% and 60-90%) while only the 95% confidence level upper limit was established for the centrality class 0-20%.At higher p T , in the interval 3 < p T < 8 GeV/c, the yield of ψ(2S) could not be extracted and the 95% confidence level upper limit is shown for the 0-20% and 20-60% most central collisions.The results do not allow a firm conclusion since statistical fluctuations inside one standard deviation allow our data points to range between very low double ratios (strong ψ(2S) suppression with respect to J/ψ) to values higher than unity (less ψ(2S) suppression with respect to J/ψ).Nevertheless, the limit set for the lowest p T bin for the 0-20% most central collisions points to a larger suppression of the ψ(2S) in that region.A transport model calculation [17] for inclusive ψ(2S) and J/ψ production is shown for the two p T intervals considered.The theoretical uncertainty band is due to different choices of the quenching factor for the b-quark.A qualitative agreement can be appreciated for both p T intervals.
It is worth underlying that our result is for inclusive ψ(2S) and J/ψ production, while SHM predictions and CMS results are for prompt charmonia production.The impact of the B-mesons feed-down on the ratio was extensively studied in [17], showing a very strong influence of the non-prompt ψ(2S) component on the final result.According to this study, removing this non-prompt contribution would lead to a significantly lower double ratio at high N part : in the 0 < p T < 3 GeV/c bin a 60% decrease is expected, while in the 3 < p T < 8 GeV/c bin the effect could be even stronger, leading to a 80% decrease.

Conclusions
We have presented a study of J/ψ and ψ(2S) production in Pb-Pb collisions at √ s NN = 2.76 TeV in the transverse momentum and rapidity ranges p T < 8 GeV/c and 2.5 < y < 4.This analysis was carried out in the muon spectrometer system, whose tracking and triggering capabilities were described in detail.
The [ψ(2S)/J/ψ] Pb-Pb ratio was measured in two p T ranges as a function of centrality.In some intervals, only the 95% confidence level upper limits could be obtained.The suppression pattern of the ψ(2S) is compatible with that of the J/ψ in most of the centrality and p T intervals studied.The large uncertainties leave open the possibility of strong enhancement or suppression factors.An accurate ψ(2S) measurement in Pb-Pb would require significantly more statistics than the one presented in this analysis.
The J/ψ signal was extracted as a function of p T , y and the collision centrality.We have computed the J/ψ p T and p 2 T .The J/ψ p T in Pb-Pb collisions decreases significantly (5σ effect) from peripheral to central collisions.In addition we have studied r AA defined as the ratio of the J/ψ p 2 T measured in Pb-Pb and pp collisions at the same energy.The r AA exhibits a clear decrease as a function of centrality for Pb-Pb collisions.
The nuclear modification factor, R AA , of inclusive J/ψ was measured as a function of centrality.A constant suppression of about 40% was observed for N part larger than 70 [27].New studies of the J/ψ suppression pattern as a function of centrality for three p T ranges were presented.Above N part ∼ 150, the low p T J/ψ R AA clearly differs from the high p T J/ψ R AA and is about three times larger for N part > 250, corresponding to a 3.9σ separation.Complementary to this, the p T dependence of the suppression pattern was analysed for the different centrality classes.An increase of the inclusive J/ψ R AA with decreasing p T is observed below 5 GeV/c in the most central Pb-Pb collisions (0-20%), while no significant p T dependence is seen in the most peripheral collisions (40-90%).As a function of rapidity, the results published in [27] show compatible R AA values for |y| < 0.9 and 2.5 < y < 3.For larger rapidity, a decreasing trend is visible.
Comparisons of the r AA and R AA measured in ALICE with lower energy experiments show significant differences.The decreasing trend of r AA observed as a function of centrality is opposite to NA50 and PHENIX measurements.The R AA in the most central collisions is three times larger than the one measured by PHENIX, and the difference reaches a factor four in the p T region below 1 GeV/c.If the suppression sources observed at lower energies, which were related to color screening in hot nuclear matter on top of CNM effects, are still present at the LHC, then other mechanisms compensating the J/ψ suppression are needed to explain the ALICE measurements.This conclusion is further substantiated, in the region |y| < 3, by the comparison of the inclusive J/ψ R AA measurements as a function of y to models implementing only CNM effects, which shows a qualitative agreement.
The inclusive J/ψ r AA and R AA measurements were also compared to various theoretical calculations including hot and cold nuclear matter effects.The hadronic part of the J/ψ R AA was estimated when needed to allow for a direct comparison to models, which do not implement the J/ψ production mechanism at the origin of the observed very low p T excess [40].All these models feature a full or partial J/ψ production from charm quarks recombination and are in fair agreement with the experimental results.The transport models considered in this paper are also able to generate an amount of J/ψ elliptic flow comparable to the one measured in ALICE [69].The double differential studies of the inclusive J/ψ R AA as a function of centrality and p T brings new constraints to the models.Reproducing the suppression pattern in peripheral collisions for both low and high p T J/ψ is challenging for all models.Some tensions also appear in describing the R AA evolution at low p T for all centrality classes.However, the uncertainties on the measurements on one side, and on the CNM and dσ cc /dy in the theoretical calculations on the other side, do not allow for drawing a firm conclusion.The large uncertainties on the model predictions also show the limit of using the R AA as an observable to measure the J/ψ suppression due to hot medium effects.Ideally one should, in Pb-Pb collisions, compare the J/ψ production to the charm production to cancel out the cold nuclear matter effects affecting the initial cc dynamics.However, the measurement of the charm cross section in Pb-Pb collisions is very ambitious and still remains to be done at the LHC.

= 1 = 1 = 1
in simulation as the p T value for which the muon trigger probability is 50%.The following muon-specific L0 triggers were used: -Single muon low p T (p trig T GeV/c): MSL -Opposite-sign dimuon low p T (p trig T GeV/c on each muon): MUL -Like-sign dimuon low p T (p trig T GeV/c on each muon): MLL

c
Entries

Fig. 4 :
Fig. 4: The J/ψ acceptance times efficiency, shown as a function of centrality (left) and as a function of p T and y for the centrality class 0-90% (right).The vertical error bars in the left panel represent the statistical uncertainties.

Fig. 5 :
Fig. 5: Differential yields of inclusive J/ψ in Pb-Pb collisions at √ s NN = 2.76 TeV as a function of p T for three centrality classes.Solid lines correspond to the results from the fit described in the text.
The p T of inclusive J/ψ measured in pp and Pb-Pb collisions at √ s NN = 2.76 TeV is shown in Fig. 6 (left side) as a function of N part .The error bars (open boxes) represent the statistical (systematic) uncertainties.A clear downward trend in p T is observed when going from pp to the most central Pb-Pb collisions.The p T decrease from peripheral (40-90%) to central (0-20%) collisions is significant, the two values being separated by more than 5σ .These results are compared to the ones obtained by PHENIX in pp, Cu-Cu and Au-Au collisions at √ s NN = 0.2 TeV.There is no evidence for a decreasing trend, contrary to what is observed in the ALICE measurement.

9. 1
Centrality dependence of R AA Our measurement of the inclusive J/ψ R AA at √ s NN = 2.76 TeV in the range 2.5 < y < 4 and p T < 8 GeV/c is shown in Fig. 7 as a function of N part .Statistical (uncorrelated systematic) uncertainties are represented by vertical error bars (open boxes).A global correlated systematic uncertainty affecting all the values by the same amount is quoted in the legend.The same convention is applied in the following figures, unless otherwise specified.The J/ψ R AA in the centrality class 0-90% (corresponding to N part ∼ 124, see 〉

Fig. 8 :
Fig.8: Inclusive J/ψ R AA as a function of the number of participant nucleons measured in Pb-Pb collisions at √ s NN = 2.76 TeV for three p T ranges (0-2, 2-5 and 5-8 GeV/c) and comparisons of the lowest and highest p T range to the transport and to the comover interaction models[13,60,63].The brackets quantify the possible range of variation of the hadronic J/ψ R AA for two extreme hypotheses on the photo-production contamination in the inclusive measurement.

Fig. 10 :
Fig. 10: Inclusive J/ψ R AA as a function of the J/ψ rapidity measured in Pb-Pb collisions at √ s NN = 2.76 TeV [27], compared to theoretical calculations of CNM effects due to shadowing and/or coherent energy loss [65-67].

Fig. 11 :
Fig. 11: Inclusive [ψ(2S)/J/ψ] ratio measured as a function of N part in Pb-Pb collisions at √ s NN = 2.76 TeV for two p T intervals, compared to NA50 results [29] and to a theoretical calculation [16] (left).Double ratio, as a function of N part , between the ψ(2S) and J/ψ measured in Pb-Pb at √ s NN = 2.76 TeV and pp collisions at √ s = 7 TeV, compared to theoretical calculations [17] (right).

Table 2 :
Summary of the systematic uncertainties (in %) entering the J/ψ yield and/or R AA calculation as a function of centrality, p T and y.Numbers with an asterisk correspond to the systematic uncertainties fully correlated as a function of the given variable.
p T range y range Centrality p T ± stat.± syst.

Table 3 :
Values of p T and p 2 Inclusive J/ψ R AA as a function of the number of participant nucleons measured in Pb-Pb collisions at √ s NN = 2.76 TeV [27], compared to the PHENIX measurement in Au-Au collisions at √ s NN = 0.2 TeV [21] (left) and to theoretical models

Table A . 2 :
Inclusive J/ψ R AA and Pb-Pb yields as a function of centrality, for p T < 8 GeV/c and 2.5 < y < 4.0.Statistical and systematic uncertainties are also reported.A global systematic uncertainty of 15% (12%) affects all the R AA (yields) values.

Table A . 4 :
Inclusive J/ψ R AA and Pb-Pb yields as a function of centrality, for 2 < p T < 5 GeV/c and 2.5 < y < 4.0.Statistical and systematic uncertainties are also reported.A global systematic uncertainty of 14% (11%) affects all the R AA (yields) values.

Table A . 12 :
Inclusive J/ψ R AA as a function of centrality, for 0.3 < p T < 8 GeV/c and 0.3 < p T < 2 GeV/c in the rapidity range 2.5 < y < 4.0.Statistical and systematic uncertainties are also reported.A global systematic uncertainty of 15% affects all the R AA values.Centrality R AA ± (stat.) ± (syst.) for 0.3 < p T < 1 GeV/c