Studies of J/ψ production at forward rapidity in Pb–Pb collisions at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV

The inclusive J/ψ production in Pb–Pb collisions at the center-of-mass energy per nucleon pair sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV, measured with the ALICE detector at the CERN LHC, is reported. The J/ψ 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/ψ yield in Pb–Pb collisions with respect to binary-scaled pp collisions is quantified by the nuclear modification factor (RAA). The RAA at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV is presented and compared with previous measurements at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 2.76 TeV as a function of the centrality of the collision, and of the J/ψ transverse momentum and rapidity. The inclusive J/ψ RAA shows a suppression increasing toward higher transverse momentum, with a steeper dependence for central collisions. The modification of the J/ψ average transverse momentum and average squared transverse momentum is also studied. Comparisons with the results of models based on a transport equation and on statistical hadronization are 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 socalled 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 (charm (c) and beauty (b)) 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 c and 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 P C = 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) collisions was observed at the SPS at √ s NN = 17 GeV [6,7], at RHIC up to -1 -
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 ) [20]. The observed suppression does not increase with increasing collision energy as expected in the colorscreening 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/ψ [21,22]. 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 [21,22]. The regeneration scenario was further supported by the measurement of a positive J/ψ elliptic flow (v 2 ) [23][24][25][26][27] which, at low p t , can be acquired via charm-quark recombination [28,29]. 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 [30,31]) 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 [32][33][34][35][36][37]. 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 b quark [38]. 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 nonprompt 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. The paper is organized as follows: section 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 section 3. Results are presented and compared to available measurements at √ s NN = 2.76 TeV and model calculations in section 4.

Apparatus and data sample
The ALICE detector and its performance are extensively described in refs. [39] and [40], 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 [41] 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 [42], 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 [43], 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 [44].
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 refs. [45,46]. 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 [47].
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 .
1 In the ALICE reference frame, the muon spectrometer covers a negative η interval and consequently negative y values. We have chosen to present our results with a positive y notation.

JHEP02(2020)041
3 Data analysis For a centrality class i, the double-differential 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 [48], ∆p t and ∆y are respectively the widths of the p t and y intervals, (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. The values of N i MB are 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 [20]. The values of T AA i in different centrality classes were obtained using a Glauber calculation [46,49,50]. 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 refs. [17,51] and are used here as a reference. In addition, 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 [52].
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. [53]. The selection criteria applied to both single muons and dimuons are identical to the ones used in refs. [14,17], requiring a match between tracks reconstructed in the tracking system and track segments in the muon trigger system.

JHEP02(2020)041
The signal extraction was performed over the US dimuon invariant mass ranges [2.2,4.5] and [2.4,4.7] GeV/c 2 using two methods. In the first one, the 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 [14], 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 a double-sided Crystal Ball function (CB2, where independent non-Gaussian tails are present on both sides of a Gaussian core) or a pseudo-Gaussian with a mass-dependent width [54]. For both functions, the position of the J/ψ pole 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 [55] and GEANT4 [56]) 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 [51], the sample with the largest significance of the J/ψ signal. The ψ(2S) signal was included in the fits using the same signal functions as for the J/ψ, with mass and width tied to those of the J/ψ [57].
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 and originating mainly from semi-muonic decays of pairs of charm hadrons, 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 figure 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. 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. [33]. Unpolarized J/ψ production was assumed, consistently with the measurements of inclusive J/ψ polarization in pp collisions [58,59]. 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. [51]. 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. [60].
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  Matching efficiency

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 refs. [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 √ s = 5.02 TeV by interpolating in energy the corresponding LHCb crosssection measurements in pp collisions at √ s = 2.76 and 7 TeV [61][62][63]. The limits on R prompt AA correspond to the two extreme hypotheses of total non-prompt J/ψ suppression (R non-prompt AA = 0) and absence of suppression (R non-prompt AA = 1). 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 [64] and related to J/ψ photo-production [65]. 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 to the region p t > 0.3 GeV/c 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 and lower limit brackets correspond to the extremest hypotheses on the contribution from photo-produced J/ψ and on the efficiency of the aforementioned p t selection as described in ref. [14].
4.1.1 Centrality-integrated R AA as a function of y and p t Figures 2 and 3 show the inclusive J/ψ nuclear modification factor as a function of transverse momentum and rapidity, integrated over the centrality class 0-90%. The results are compared with those obtained at √ s NN = 2.76 TeV [14] and with the results of the calculation of a transport model [28,66]. A significant increase of R AA is visible with decreasing p t , which was already observed for the most central events (0-20%) and reported in ref. [17]. Within uncertainties, the results are compatible with those obtained, in a more restricted p t interval, at the lower LHC energy, with a possible hint (1.2σ) of a weaker suppression in the region 2 < p t < 6 GeV/c. The transport model calculations are in qualitative agreement with the data. In this model, a competition between suppression and regeneration of charmonia is assumed, choosing a cc production cross section dσ cc /dy = 0.57 mb and dσ pp J/ψ /dy = 3.35 µb for 2.5 < y < 4. The latter value is ∼10% smaller than our measurement of the same quantity [17]. The model also 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. Figure 3 shows that, in the explored rapidity interval, there is no significant variation of the R AA values. The calculations of the transport model are in good agreement with the experimental results. The comparison of the results with those obtained at √ s NN = 2.76 TeV [14] hints (1.5σ) for a weaker suppression at √ s NN = 5.02 TeV at large y (3.75 < y < 4).
4.1.2 Centrality-differential R AA as a function of y and p t   noting that the results for 0-20% were already published in ref. [17]. In this paper, the corresponding values were updated with the improved T AA uncertainties reported in ref. [49]. In Figure 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) [67]. A good agreement with this calculation, which does not include contributions from non-prompt J/ψ production, 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 figure 5, the R AA values exhibit a very weak rapidity dependence in all the centrality classes, as also observed in 0-90% (figure 3). The calculation of the transport model is able to describe the data, in 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 figures 6 and 7 the 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 figure 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. [68,69], the origin may be from the bias introduced by the event selection 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 figure 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. The vertical error bars represent statistical uncertainties, the boxes around the points uncorrelated systematic uncertainties, while the correlated uncertainty is shown as a filled box around R AA = 1. The corresponding measurements in Pb-Pb collisions at √ s NN = 2.76 TeV [14] are also shown, as well as the ratio of the R AA values, which is depicted in the bottom panel of the figure. The R AA values at √ s NN = 5.02 TeV and the ratios to lower energy results are compared with transport model calculations [28] and, for 0-20% and 20-40% centrality, with the results of the SHM [67]. The brackets around R AA values for 40-90% centrality in the lowest p t interval represent an estimate of the maximum influence of J/ψ photo-production, as detailed in section 4.  The vertical error bars represent statistical uncertainties, the boxes around the points uncorrelated systematic uncertainties, while the correlated uncertainty is shown as a filled box around R AA = 1. The R AA values are compared with transport model calculations [28]. The brackets around R AA values for 40-90% centrality represent an estimate of the maximum influence of J/ψ photoproduction, as detailed in section 4.1. -

J/ψ average transverse momentum and r AA
A complementary insight into the modification of J/ψ transverse momentum distributions in Pb-Pb collisions can be obtained by the study of the J/ψ average transverse momentum p t and the average squared 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 where C, p 0 and n are free parameters. This function is widely used to reproduce the J/ψ  p t distribution in hadronic collisions (e.g refs. [70,71]). 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 figure 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.
In the left panel of figure 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 figure 4. The p t 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 figure 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 hypothe-  ses 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.  [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 [28] 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.

JHEP02(2020)041
The inclusive J/ψ R AA as a function of transverse momentum and rapidity for the centrality range 0-90%, is compatible with previously published results at √ s NN = 2.76 TeV [14]. A suppression of the 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 centralitydifferential 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 mea- 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.

Acknowledgments
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex.

JHEP02(2020)041
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.