Studies of J/$\psi$ production at forward rapidity in Pb-Pb collisions at $\sqrt{s_{\rm{NN}}}$ = 5.02 TeV

The inclusive J/$\psi$ production in Pb-Pb collisions at the center-of-mass energy per nucleon pair $\sqrt{s_{\rm{NN}}}$ = 5.02 TeV, measured with the ALICE detector at the CERN LHC, is reported. The J/$\psi$ meson is reconstructed via the dimuon decay channel at forward rapidity ($2.5<y<4$) down to zero transverse momentum. The suppression of the J/$\psi$ yield in Pb-Pb collisions with respect to binary-scaled pp collisions is quantified by the nuclear modification factor ($R_{\rm{AA}}$). The $R_{\rm{AA}}$ at $\sqrt{s_{\rm{NN}}}$ = 5.02 TeV is presented and compared with previous measurements at $\sqrt{s_{\rm{NN}}}$ = 2.76 TeV as a function of the centrality of the collision, and of the J/$\psi$ transverse momentum and rapidity. The inclusive J/$\psi$ $R_{\rm{AA}}$ shows a suppression increasing toward higher $p_{\rm{T}}$, with a steeper dependence for central collisions. The modification of the J/$\psi$ average $p_{\rm{T}}$ and $p_{\rm{T}}^{2}$ is also studied. Comparisons with the results of models based on a transport equation and on statistical hadronization are also carried out.


Introduction
The study of ultra-relativistic heavy-ion collisions aims to investigate the properties of strongly-interacting matter at high temperature and energy density. Lattice Quantum Chromodynamics calculations predict that a deconfined state of partonic matter, the so-called Quark-Gluon Plasma (QGP), can be created in such collisions [1][2][3]. Among the many possible probes to study this phase of matter, heavy quarks are particularly interesting as they are expected to be produced in the initial stage of the collisions, by hard partonic scatterings, and to experience the full evolution of the system. In particular, it was predicted that bound states of charm (c) and anti-charm (c) quarks (known as charmonia) should be suppressed due to the color-screening mechanism [4]. The suppression probabilities of the quarkonium (cc or bb) states in the QGP depend on their binding energies and the medium temperature. Therefore, the measurement of the relative production rates of the quarkonium states should give indications on the temperature of the system [5]. Among the different charmonium states, the study of the ground state with quantum numbers J PC = 1 −− (J/ψ) is comparatively more accessible due to its larger abundance and to the relatively large branching ratio to dileptons, and has led to several important results.
Over the past decades, the J/ψ production in heavy-ion collisions was measured at the SPS, RHIC and the LHC, covering a wide range of center-of-mass energies per nucleon pair ( √ s NN ) from about 17 GeV to 5.02 TeV. A suppression of the J/ψ production yield in nucleus-nucleus (AA) relative to that expected from measurements in proton-proton (pp) [17,18]. The suppression is evaluated through the calculation of the nuclear modification factor (R AA ), corresponding to the ratio of the production yields in AA and the cross section in pp collisions, normalised by the nuclear overlap function T AA [19]. The observed suppression does not increase with increasing collision energy as expected in the color-screening picture considering the increasing temperature of the formed QGP. This observation is naturally explained by a further production mechanism known as regeneration, in which abundantly produced cc pairs recombine into J/ψ [20,21]. The contribution of the regeneration to J/ψ production has to increase with the density of cc pairs and consequently with the collision energy. It is worth noting that the regeneration contribution should favour low transverse momentum (p T ) J/ψ, as the bulk of charm quarks are produced at small momenta [20,21]. The regeneration scenario was further supported by the measurement of a positive J/ψ elliptic flow (v 2 ) [22,23] which, at low p T , can be acquired via charm-quark recombination [24,25]. It is important to note that in addition to the effects discussed above, related to the production of a high energy-density medium, the so-called cold-nuclear-matter effects may also have a sizeable influence on the charmonium yields. In particular, the modification of the parton distribution functions in the nucleus (e.g. nuclear shadowing [26,27]) may modify the initial yields of charm quarks and has to be taken into account in the interpretation of the results. Quantitative estimates of these effects are carried out via the study of proton-nucleus collisions [28][29][30][31][32][33]. Finally, a quantitative interpretation of the results requires taking into account that the observed J/ψ are produced either promptly, i.e. as direct J/ψ or via decay of higher-mass charmonium states (χ c , ψ(2S)), or non-promptly through the weak decay of hadrons containing a beauty (b) quark [34].
For a better assessment of the suppression-regeneration scenario, extensive studies of the centrality, p T and rapidity dependence of the J/ψ nuclear modification factor have to be carried out. The first ALICE measurement of the inclusive (sum of prompt and non-prompt sources) J/ψ production at √ s NN = 5.02 TeV at forward rapidity [17] has shown a hint for an increase of R AA with respect to the √ s NN = 2.76 TeV results in the region 2 < p T < 6 GeV/c, while the results were consistent elsewhere.
In this paper, we complement the results obtained in Ref. [17]. The J/ψ R AA is simultaneously obtained in different collision centrality classes and p T or rapidity intervals. In addition, to further assess the kinematic region of influence of the J/ψ regeneration mechanism, results on the J/ψ average p T and p 2 T as a function of centrality are presented. J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration The paper is organized as follows: Sec. 2 is dedicated to the description of the ALICE detector systems used in this analysis. The analysis procedure is briefly explained and a summary of the systematic uncertainties is also given in Sec. 3. Results are presented and compared to available measurements at √ s NN = 2.76 TeV and model calculations in Sec. 4.

Apparatus and data samples
The ALICE detector and its performance are extensively described in Ref. [35] and [36], respectively. J/ψ mesons are reconstructed in the muon spectrometer (covering the pseudo-rapidity interval −4 < η < −2.5 1 ) via their dimuon decay channel down to zero p T . The muon spectrometer consists of a 4.1 m (10 interaction lengths, λ int ) thick front absorber which is used to filter out hadrons coming from the interaction point (IP), followed by tracking (MCH) and triggering (MTR) systems. Each of the five tracking stations is composed of two planes of Cathode Pad Chambers. The third tracking station is located inside a dipole magnet with a field integral of 3 Tm. A 1.2 m (7.2λ int ) thick iron wall, which absorbs secondary hadrons escaping from the front absorber and low-momentum muons produced predominantly from π and K decays, is located between the tracking system and the trigger stations. Each of the two trigger stations consists of two planes of Resistive Plate Chambers. Finally, a small-angle conical absorber around the beam-pipe protects the spectrometer from secondary particles produced by interactions of large-η primary particles with the beam-pipe.
The other detectors used in this analysis are the Silicon Pixel Detector (SPD), the V0 scintillator detectors, the Cherenkov detectors T0 and the Zero Degree Calorimeters (ZDC). The SPD [37] provides the coordinates of the primary vertex of the collision, and consists of two cylindrical layers covering |η| < 2 (inner layer) and |η| < 1.4 (outer layer). The V0 [38], composed of two arrays of 32 scintillator tiles each, and located on both sides of the IP, covers 2.8 < η < 5.1 (V0A) and −3.7 < η < −1.7 (V0C), and is used as a trigger detector, for the centrality determination and to remove beam-induced background. It is also used for the measurement of the luminosity along with the T0 detector [39], which consists of two quartz Cherenkov counters, located on each side of the IP and covering the pseudo-rapidity intervals −3.3 < η < −3 and 4.6 < η < 4.9. The ZDCs, located on either side of the IP at ± 114 m along the beam axis, detect spectator nucleons emitted at zero degrees with respect to the LHC beam axis, and are used to reject electromagnetic Pb-Pb interactions [40].
The centrality determination and the evaluation of the average number of participant nucleons in the collision N part for each centrality class is based on a Glauber model fit to the V0 signal amplitude distribution as described in Ref. [41,42]. The events are classified in centrality classes corresponding to percentiles of the nuclear hadronic cross section. In this analysis, events corresponding to the most central 90% of the inelastic cross section were selected. For these events the minimum bias (MB) trigger is fully efficient and the residual contamination from electromagnetic processes is negligible. The MB trigger is defined by a coincidence of the signals from both sides of the V0 detector.
The analysis presented here is based on dimuon-triggered events which require, in addition to the MB condition, the detection of two Unlike-Sign (US) tracks in the triggering system of the muon spectrometer. The muon trigger selects muon candidates having a transverse momentum larger than a given threshold which corresponds to the value for which the trigger efficiency reaches 50%. In Pb-Pb collisions the p T threshold is ≈ 1 GeV/c with the single-muon trigger efficiency reaching a plateau value of 98% at ≈ 2.5 GeV/c [43].
The current analysis exploits the data samples of Pb-Pb collisions at √ s NN = 5.02 TeV collected during 2015. This corresponds to an integrated luminosity L PbPb int ≈ 225 µb −1 .
J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration

Data analysis
Below, the definitions of the measured observables are given, followed by a presentation of the main analysis steps and a discussion of the systematic uncertainties. For a centrality class i, the doubledifferential J/ψ invariant yield Y i J/ψ is defined as where N i J/ψ (p T , y) is the number of J/ψ for a given p T and y interval, BR J/ψ→µ + µ − = (5.96 ± 0.03)% is the branching ratio of the dimuon decay channel [44], (A × ε) i (p T , y) is the product of the detector acceptance and the reconstruction efficiency for that p T and y interval, and N i MB is the equivalent number of minimum-bias events. It is obtained as the product of the number of dimuon-triggered events times the inverse of the probability of having a dimuon trigger in a MB event (F i ). The F i values correspond to those quoted in Ref. [17]. For the centrality integrated sample, the value of the normalization factor is F 0−90% = 11.84 ± 0.06. The quoted uncertainty is systematic and corresponds to the difference between the results obtained with two methods, either by calculating the ratio of the counting rates of the two triggers, or by applying the dimuon trigger condition in the analysis of MB events. The nuclear modification factor R AA is given by where T AA i is the average of the nuclear-overlap function [19]. The values of T AA i in different centrality classes were obtained using a Glauber calculation [42,45,46]. The systematic uncertainty on the T AA i calculation, which ranges from 1% in the most central class to 3% in the most peripheral one, was determined by varying the density parameters of the Pb nucleus and the nucleon-nucleon inelastic cross section within their uncertainties. The systematic uncertainty on the definition of the centrality intervals is evaluated by varying by ±0.5% the fraction (90%) of the hadronic cross section selected with the chosen minimal cut on the V0 signal amplitude, and redefining accordingly the centrality intervals, following the approach detailed in Ref. [17]. Values of the J/ψ cross section in pp collisions (d 2 σ pp J/ψ /dydp T ) at √ s = 5.02 TeV were already reported in Ref. [17,47]. Here, following the same analysis procedure as detailed in those papers, the cross section was evaluated in four p T intervals (0.3-2, 2-5, 5-8, and 8-12 GeV/c) for the interval 2.5 < y < 4 and three y intervals (2.5-3, 3-3.5, and 3.5-4) for the interval p T < 12 GeV/c. The integrated luminosity of the pp sample is L pp int = (106.3 ± 2.2) nb −1 [48]. The J/ψ candidates were formed by combining US muons reconstructed within the geometrical acceptance of the muon spectrometer using the tracking algorithm described in Ref. [49]. The selection criteria applied to both single muons and dimuons are similar to the ones used in Ref. [14,17], requiring a match between tracks reconstructed in the tracking system and track segments in the muon trigger system.
The signal extraction was performed over various centrality, p T and y intervals using two methods. In the first one, the US dimuon invariant mass distributions were fitted with the sum of a signal and a background function. In the second one, the combinatorial background (dominant in central Pb-Pb collisions) was first estimated using an event-mixing technique, and then subtracted from the raw invariant mass distribution. Finally, the resulting distributions were fitted with the sum of a signal and a residual background component.
The signal component of the fitting function is either an extended Crystal Ball (CB2) function or a pseudo-Gaussian with a mass-dependent width [50]. For both functions, the position of the J/ψ pole J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration mass, as well as the width of the resonance, are free parameters of the fits, while the non-Gaussian tail parameters were fixed. Two sets of tail parameters were obtained from Monte Carlo (MC) simulations using different particle transport codes (GEANT3 [51] and GEANT4 [52]) to account for the sensitivity of these parameters to the description of the detector materials. In addition, another set of CB2 tail parameters was extracted from the pp collisions at √ s = 13 TeV data sample [47], the sample with the largest significance of the J/ψ signal. Finally, the contribution of the ψ(2S) signal was also taken into account in the fits, using the same signal functions as for the J/ψ, with mass and width tied to those of the J/ψ [53].
The background was parametrized either as a pseudo-Gaussian with a width quadratically dependent on the mass or as a ratio of a 2 nd to 3 rd order polynomial. When using the event-mixing technique, the continuum component of the correlated background remaining in the US dimuon distributions after the background subtraction was parametrized using the sum of two exponential functions. Examples of fits to the US dimuon invariant mass distributions, without and with subtraction of mixed-event background, are shown in Fig. 1 for different centrality classes and p T intervals. For each centrality class, p T and y interval, the number of J/ψ and the statistical uncertainty are given by the average of the results from the considered fit configurations obtained by varying the signal and background functions, the tail parameters and the invariant mass fit range ([2.2,4.5] or [2.4,4.7] GeV/c 2 ). The systematic uncertainty is defined, for each centrality, p T and y interval, as the RMS of the various fit results. It varies between 1.5% and 3.6% as a function of centrality or p T and between 1.5% and 5% as a function of y. The J/ψ A × ε was obtained using MC simulations, where the p T and y distributions for the generated J/ψ were matched to the ones extracted from data using an iterative procedure as done in Ref. [29]. Unpolarized J/ψ production was assumed [54,55]. The misalignment of the detection elements as well as the time-dependent status of each electronic channel during the data taking period were taken into account in the simulation. Generated J/ψ → µ + µ − signals were embedded into real minimum bias events in order to properly reproduce the effect of detector occupancy and its variation from one centrality class to another, and reconstructed as for real events. A relative decrease by ∼14% of A × ε was observed in the most central Pb-Pb collisions with respect to the most peripheral ones.
The following sources of systematic uncertainty on A × ε were considered: (i) the parametrization of the input p T and y shapes, (ii) the uncertainty on the tracking efficiency in the muon tracking chambers, (iii) the uncertainty on the MTR efficiency and (iv) the matching between tracks reconstructed in the tracking and triggering systems.
For the parametrization of the MC input distributions, two sources of systematic uncertainty were considered: the effect of the finite data sample used to parametrize these distributions and the correlations between p T and y (more explicitly, the fact that the p T distribution of the J/ψ varies within the rapidity interval in which it is measured). The former turns out to be negligible. For the latter, different MC simulations were performed by varying the input p T and y distributions within limits that correspond to this effect and re-calculating the A × ε in each case as done in Ref. [47]. The uncertainties on the tracking efficiency in the MCH, trigger efficiency in the MTR, and on the matching efficiency between MTR and MCH tracks were evaluated by comparing the efficiencies obtained in data and MC at the single muon level and propagating the observed differences to the J/ψ candidates, as done in Ref. [56].
In each centrality, p T and y interval, the total systematic uncertainty on the yield and R AA is determined as the quadratic sum of the uncertainties from the different sources listed in Table 1. Correlations of various uncertainties vs centrality, p T or y are also reported. The values in the last row correspond to the sum of the statistical and systematic uncertainties on the pp reference.

Nuclear modification factor
This section summarizes the results for the inclusive J/ψ R AA at forward rapidity in Pb-Pb collisions at √ s NN = 5.02 TeV as a function of: -rapidity and transverse momentum, integrated over the centrality (class 0-90%); -rapidity and transverse momentum, for the centrality classes 0-20%, 20-40% and 40-90%; -centrality, in four transverse momentum intervals and in three rapidity intervals.
When possible, the ratio between the results of this analysis and the measurements in Pb-Pb collisions at √ s NN = 2.76 TeV [14], in the same kinematic interval, is computed. Only the uncertainty related to the T AA cancels out in the ratio, as discussed in Ref.
[17]. Following the same approach as in Ref. [14,17], an estimate of the R AA of prompt J/ψ was determined by making conservative assumptions on the size of the non-prompt R AA . The relation between the inclusive (R AA ), prompt (R prompt AA ) and non-prompt (R non-prompt AA ) nuclear modification factors can be expressed as: where F B is the fraction of non-prompt to inclusive J/ψ in pp collisions. This quantity is evaluated at The effect is small at moderate transverse momentum ( 10% for p T 5 GeV/c) and then increases at higher p T . Numerical values for the limits on R prompt AA can be found in the HepData record associated to this paper. Another effect which may influence the interpretation of the inclusive J/ψ results is the presence of an excess at very low p T , observed at √ s NN = 2.76 TeV [60] and related to J/ψ photo-production [61]. This source was shown to be significant with respect to hadronic production for peripheral Pb-Pb collisions and has a strong influence on the measured R AA values. For this reason, the region p T < 0.3 GeV/c was excluded when dealing with peripheral collisions. The remaining contribution of this source was evaluated following the procedure detailed in Ref. [14] and the maximum effect on R AA is explicitly shown in the following figures by use of bracket symbols. The upper limit bracket corresponds to the assumption of no contribution from photo-produced J/ψ, while the lower one corresponds to the hypothesis that all the J/ψ with p T < 300 MeV/c are photo-produced. includes contributions from both prompt and non-prompt J/ψ. The upper (lower) limit of this calculation corresponds to a 10% (25%) contribution of nuclear shadowing. 4.1.2 Centrality-differential R AA as a function of y and p T Figures 4 and 5 show, respectively, the p T and y dependence of the inclusive J/ψ R AA , for events corresponding to the centrality classes 0-20% 2 , 20-40% and 40-90%. In Fig. 4, moving from central to peripheral collisions, a weaker p T dependence of the R AA is observed, up to an almost constant nuclear modification factor for 40-90% centrality. When comparing results at √ s NN = 5.02 TeV and 2.76 TeV, a slight increase of the R AA is visible for the most central collisions and for 2 < p T < 6 GeV/c at the higher collision energy, while the results are compatible in the 20-40% and 40-90% samples. A fair agreement with the transport model calculations is observed. The results for the 0-20% and 20-40% centrality classes are also compared with a model based on statistical hadronization (SHM) [63]. A good agreement can be found up to p T ∼ 4 GeV/c, while at higher transverse momentum R AA is underestimated. This feature could partly be due to additional production mechanisms, not implemented in the model, such as J/ψ production from gluon fragmentation in jets.
In Fig. 5, the R AA values exhibit a very weak rapidity dependence in all the centrality classes, as also observed in 0-90% (Fig. 3). The calculation of the transport model is able to describe the data, in J/ψ production in Pb-Pb collisions at √ s NN = 5.  particular when a weak nuclear shadowing scenario (10%, corresponding to the lower limit chosen by the authors) is adopted.

Centrality dependence of R AA
In Figs. 6 and 7 the nuclear modification factor R AA as a function of the average number of participant nucleons N part is shown for various transverse momentum and rapidity intervals, respectively. The N part intervals correspond to the centrality selections 0-10%, 10-20%, 20-30%, 30-40%, 40-50%, 50-60%, and 60-90%, from larger to smaller N part values. The results of Fig. 6 clearly show that moving from low to high p T the centrality dependence of R AA becomes steeper, with R AA reaching a minimum value of 0.29 ± 0.02(stat) ± 0.01(syst) for the 0-10% centrality class and 8 < p T < 12 GeV/c. In the low-p T region (0.3 < p T < 2 GeV/c), the R AA has a weak N part dependence and is compatible with being constant (∼ 0.7) for N part > 150. In the most peripheral centrality class, a deviation from unity can be observed, in particular for p T > 2 GeV/c, not seen in the theoretical calculations. As discussed in Refs [64, 65], the origin may be from the bias introduced by the event select and collision geometry, which causes an apparent suppression. When comparing the results with those corresponding to Pb-Pb collisions at √ s NN = 2.76 TeV [14], systematically higher R AA values are found in the p T interval 2 < p T < 5 GeV/c, even if the maximum observed difference is only at 1.5σ level, for the centrality region 0-10%. In all other p T intervals where the comparison is possible, the results at the two energies are compatible. When comparing the results with the transport model calculations, the agreement is good at low p T (0.3 < p T < 2 GeV/c), while the data lie close to the upper edge of the calculation at higher p T .
In Fig. 7 the centrality dependence of the nuclear modification factor is shown for 3 rapidity intervals. No variation of the suppression pattern against rapidity is observed. The same weak dependence can also be observed with the transport model calculations.     momentum p 2 T as a function of the collision centrality. By normalizing p 2 T to the corresponding pp value, one obtains an adimensional quantity, r AA = p T 2 AA / p T 2 pp , useful for comparisons between various collision energies and/or theory calculations.
As a first step, the J/ψ invariant yields as a function of p T are fitted in various centrality classes with the following function This function is widely used to reproduce the J/ψ p T distribution in hadronic collisions (e.g Ref. [66,67]). The quantities to be determined, p T and p 2 T , are then computed as the first and second moment of f (p T ) respectively. In Fig. 8, the J/ψ invariant yields as a function of p T are shown for various centrality classes together with the fitted functions. In order to limit the influence of the J/ψ production excess at low p T , due to photo-production, the interval p T < 0.5 GeV/c was excluded from the fit. The statistical (systematic) uncertainties on p T and p 2 T were obtained from fits to the invariant yield distributions, considering only statistical (p T -uncorrelated systematic) uncertainties on the J/ψ yields. √ s NN = 5.02 TeV ALICE Collaboration  In the left panel of Fig. 9, the centrality dependence of p T is shown and compared with previous results at √ s NN = 2.76 TeV [14]. The centrality dependence of the √ s NN = 5.02 TeV results is weak up to N part ∼ 150, followed by a significant decrease towards central events. This softening of the J/ψ p T distributions is a direct consequence of the smaller suppression observed at low p T when considering the transverse-momentum dependence of the nuclear modification factors, shown in Fig. 4. The p T J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration values are systematically larger than those at √ s NN = 2.76 TeV, an effect due to the increase of the collision energy, but the decrease of p T with increasing centrality is similar at the two energies. A more direct comparison with lower energy results and theoretical calculations can be performed by studying the quantity r AA . The results are shown in the right panel of Fig. 9, and compared with those obtained in Pb-Pb collisions at √ s NN = 2.76 TeV and the transport model calculations. In peripheral collisions, and up to N part ∼ 150, the r AA value is compatible with unity within uncertainties. A maximum decrease of ∼25% is observed for central collisions. The brackets around the p T and r AA in peripheral collisions represent the possible variation of the hadronic J/ψ p T and r AA for two extreme hypotheses on the J/ψ photo-production contamination. The lower limit bracket corresponds to the assumption of no contribution from photo-produced J/ψ, while the upper one corresponds to the hypothesis that all the J/ψ with p T < 300 MeV/c are photo-produced. The results are compatible with those at √ s NN = 2.76 TeV, with a hint for larger r AA values at √ s NN = 5.02 TeV. A very different centrality dependence was observed at lower collision energies ( √ s NN = 200 GeV at RHIC [68] and √ s NN = 17 GeV at SPS [69]) as the r AA increases (especially at the SPS energy) towards more central collisions (the comparison was shown in Ref. [14]). The different behaviors of r AA at different energies can be explained by the increasing amount of J/ψ regeneration with collision energy. Finally, the comparison with the transport model calculation [24] shows good agreement for peripheral and central collisions, but an underestimation of the data points is observed in the intermediate centrality class, reaching a significance up to 2.5σ for N part ∼ 150.

Conclusions
This paper reports on ALICE measurements of the inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV in the kinematic range 2.5 < y < 4 up to p T < 12 GeV/c. Results on the nuclear modification factor R AA , the average transverse-momentum p T , and the ratio r AA were presented. A systematic comparison with the calculation of a transport model was carried out and, for the p T dependence of R AA , with the results of a statistical hadronization model.  J/ψ production is observed (R AA < 1), mild at low p T but increasing towards higher p T , and not strongly depending on rapidity. The centrality-differential studies show that the y dependence of R AA is weak and fairly independent of centrality, while the p T dependence of R AA grows steeper for more central events. All the R AA results are fairly reproduced by the calculation of a transport model, with a tendency to underestimate the observed R AA at intermediate p T . The statistical hadronization model reproduces, although with larger uncertainties, the p T dependence of R AA for various centrality classes, but shows a discrepancy in the high-p T region.
A complementary study was also carried out by measuring the centrality dependence of R AA for different p T and y intervals. A suppression strongly increasing with centrality is visible at high p T , while at low p T the suppression is relatively weak (R AA ∼ 0.7) and practically independent of centrality. On the contrary, the shape of the R AA as a function of centrality does not vary significantly in the studied rapidity ranges, showing a mild decrease until N part ∼ 100, followed by a plateau.
Finally, the r AA ratio decreases with increasing centrality, similarly to previous observations at √ s NN = 2.76 TeV. The transport model calculation underestimates the measurement at intermediate N part values.
The results shown in this paper confirm, with better accuracy, the observations carried out at √ s NN = 2.76 TeV and strengthen the evidence for the presence of a mechanism that leads to a significant increase of R AA at low p T . Recombination of charm-quark pairs during the deconfined QGP phase, as implemented in the transport model compared with our results, is a strong candidate for explaining the features of the data.                    J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration