Centrality and transverse momentum dependence of inclusive J/$\psi$ production at midrapidity in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV

The inclusive J/$\psi$ meson production in Pb-Pb collisions at a center-of-mass energy per nucleon-nucleon collision of $\sqrt{s_{\rm NN}}$ = 5.02 TeV at midrapidity ($|y|$<0.9) is reported by the ALICE Collaboration. The measurements are performed in the dielectron decay channel, as a function of event centrality and J/$\psi$ transverse momentum $p_{\rm T}$, down to $p_{\rm T}$ = 0 GeV/$c$. The J/$\psi$ mean transverse momentum $\langle p_{\rm T} \rangle$ and $r_{\rm AA}$ ratio, defined as $\langle p^{\rm 2}_{\rm T} \rangle_{\rm PbPb}/\langle p^{\rm 2}_{\rm T} \rangle_{\rm pp}$, are evaluated. Both observables show a centrality dependence decreasing towards central (head-on) collisions. The J/$\psi$ nuclear modification factor $R_{\rm AA}$ exhibits a strong $p_{\rm T}$ dependence with a large suppression at high $p_{\rm T}$ and an increase to unity for decreasing $p_{\rm T}$. When integrating over the measured momentum range $p_{\rm T}$<10 GeV/$c$, the J/$\psi$ $R_{\rm AA}$ shows a weak centrality dependence. Each measurement is compared with results at lower center-of-mass energies and with ALICE measurements at forward rapidity, as well as to theory calculations. All reported features of the J/$\psi$ production at low $p_{\rm T}$ are consistent with a dominant contribution to the J/$\psi$ yield originating from charm quark (re)combination.


Introduction
The Quark-Gluon Plasma (QGP) is a state of strongly-interacting matter characterized by quark and gluon degrees of freedom predicted by Quantum Chromodynamics (QCD) to exist at high temperature and energy density [1,2]. Such conditions are realized during the initial hot and dense stages of ultrarelativistic heavy-ion collisions. The medium produced in these collisions has a short lifetime, which is of the order of 10 fm/c at the energies reached at the CERN Large Hadron Collider (LHC), see e.g. [3].
Due to their large masses, charm and beauty quarks are produced in hard partonic scatterings occurring during the early stage of the collision and therefore experience the full evolution of the medium. Charmonia, i.e. the bound states of charm and anti-charm quarks, are of particular interest for the understanding of the QGP, see e.g. [4,5]. In the framework of color-screening models, the suppression of the charmonium state J/ψ is an unambiguous signature of the QGP [6,7]. The high density of color charges prevents charm and anti-charm quarks from forming bound states. Therefore, the J/ψ yield is expected to be suppressed compared to probes unaffected by the hot and dense medium or from expectations of the incoherent superposition of nucleon-nucleon collisions at the same energy. This was experimentally observed in the most central heavy-ion collisions at SPS [8-10] and RHIC [11][12][13] energies.
At the significantly higher collision energies of the LHC, the suppression pattern of J/ψ mesons in heavy-ion collisions is fundamentally changed. In central Pb-Pb collisions at √ s NN = 2.76 TeV, where √ s NN is the center-of-mass collision energy per nucleon-nucleon pair, the suppression was found to be weaker [14][15][16] in comparison with the earlier measurements at lower energies mentioned above. The effect was measured by the ALICE Collaboration at both mid-and forward rapidity, dominantly for J/ψ mesons at a low transverse momentum (p T ). This phenomenon is understood as the result of the charmonium (re)generation due to copiously produced charm quarks, made possible by the deconfined nature of the QGP. In addition to the weaker nuclear suppression of charmonia, recent observations of non-zero elliptic flow of D [17,18] and J/ψ [19] mesons, suggest that charm quarks may thermalize and flow with the bulk particles during the QGP phase.
There are different phenomenological scenarios available for the description of charmonium production in heavy-ion collisions. In the framework of statistical hadronization, all charmonium states are created at chemical equilibrium at the phase boundary and their abundances are determined by thermal weights [20,21]. The transport approach considers a continuous production and dissociation of charmonium states already during the QGP phase governed by a set of rate equations [22]. Another approach includes charmonium dissociation by the scattering of comoving partons and hadrons with a (re)generation component at LHC collision energies [23]. All current models implementing statistical hadronization, microscopic transport approaches [24,25] or comover interactions take into account both the hot medium and the cold nuclear matter (CNM) [26] effects mainly originating from the modification of the gluon distribution function in the nucleus compared to the corresponding function of the free nucleon.
In this paper, we present the ALICE measurement of the inclusive J/ψ production at midrapidity in Pb-Pb collisions at a center-of-mass energy per nucleon pair of 5.02 TeV. The J/ψ mesons are reconstructed in the central barrel within the rapidity range |y| < 0.9 via the e + e − decay channel down to p T = 0 GeV/c. The J/ψ p T spectrum is measured in three centrality intervals. The J/ψ average transverse momentum p T and p 2 T are evaluated as a function of collision centrality: the latter is shown in comparison with the J/ψ p 2 T measured in pp collisions, via the ratio r AA = p 2 T PbPb / p 2 T pp . The nuclear modification factor R AA , which is defined by the ratio of the production yield in Pb-Pb collisions and the production cross section in pp collisions normalized by the nuclear overlap function T AA , as a function of event centrality and J/ψ p T , is obtained using the recent ALICE measurement of the inclusive J/ψ cross section in pp collisions at √ s = 5.02 TeV [27]. The new pp reference and the larger Pb-Pb data set allow for a significant reduction of the uncertainties compared to our previous measure-J/ψ nuclear modification factor in Pb-Pb collisions at √ s NN = 5.  [24,25], and comover [23] model calculations.
The measurements presented in this publication provide a significant complement to results in Pb-Pb collisions at the same collision energy by the ALICE Collaboration at forward rapidity [29], the measurements on J/ψ suppression at high p T by the ATLAS [30] and CMS [31] Collaborations around midrapidity, as well as to results at forward rapidity in Xe-Xe collisions at √ s NN = 5.44 TeV [32].

Apparatus and data sample
A detailed description of the ALICE detector and its performance can be found in Refs. [33,34]. The ALICE central barrel detector allows for high resolution tracking and particle identification over the full azimuthal angle in the pseudorapidity range |η| < 0.9. The entire setup is placed inside a solenoidal magnet, which creates a uniform axial magnetic field of B = 0.5 T along the beam direction.
The main detectors used for the J/ψ meson reconstruction in the e + e − decay channel are the Inner Tracking System (ITS) [35] and the Time Projection Chamber (TPC) [36]. The ITS consists of 6 cylindrical layers of silicon detectors placed at radial distances to the beam line from 3.9 cm to 43 cm and provides high-precision tracking close to the interaction point as well as the determination of the primary vertex of the event. The two innermost layers form the Silicon Pixel Detector (SPD), the intermediate two layers are the Silicon Drift Detector (SDD), and the outermost layers are the Silicon Strip Detector (SSD).
Placed around the ITS, the TPC detector is a large cylindrical drift chamber extending radially from 85 cm to 250 cm from the nominal interaction point (x = y = z = 0 cm) and longitudinally between −250 cm and +250 cm. In addition to being the main tracking detector, the TPC also provides particle identification via the measurement of the specific energy loss (dE/dx) of charged particles in the detector gas.
The V0 detectors [37] consist of two scintillator arrays, V0A and V0C, which are located on both sides of the nominal interaction point at z = 329 cm and z = −90 cm and cover the pseudorapidity interval 2.8 ≤ η ≤ 5.1 and −3.7 ≤ η ≤ −1.7. The centrality of the events, expressed in fractions of the total inelastic hadronic cross section, is determined via a Glauber fit to the V0 amplitude as described in [38][39][40]. In addition, the V0 detectors are used to provide a minimum-bias trigger (MB), defined as the coincidence of signals in both V0 arrays and the beam crossing.
The results presented in this paper are obtained using the MB trigger data collected during the 2015 LHC Pb-Pb run at a center-of-mass energy per nucleon pair of 5.02 TeV. Beam-gas events are rejected using timing selections on the signals from the V0 and Zero Degree Calorimeters [41]. Pileup events are rejected online based on V0, but also in the offline analysis based on the correlation between the V0 multiplicity and the number of tracks reconstructed in the TPC. All events must have a reconstructed primary vertex with a longitudinal position within ±10 cm around the nominal interaction point. Only the events corresponding to the most central 90% of the Pb-Pb inelastic cross section (0-90%) are used in this analysis. For these events the MB trigger is fully efficient and the contamination by electromagnetic interactions is negligible. After all selections, a sample of 70 million events is available for analysis, corresponding to an integrated luminosity of L int ≈ 10 µb −1 .

Analysis methods
The J/ψ candidates are reconstructed using the e + e − decay channel. The selected electron candidates are tracks reconstructed using both the ITS and TPC detectors. They must have a minimum transverse momentum of 1 GeV/c and pseudorapidity in the range |η| < 0.9. Primary electrons are selected using a maximum distance-of-closest-approach to the event vertex of at most 1 cm and 3 cm in the transverse and J/ψ nuclear modification factor in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration  longitudinal directions, respectively. Additionally, kink-daughters, i.e. secondary tracks from long-lived weak decays of charged particles, are removed from the analysis. In order to improve the resolution of track reconstruction and to reject secondary electrons from photon conversions in the detector material, the tracks are selected to have at least one hit in either of the two SPD layers. Electrons and positrons from photon conversions are further rejected using a prefilter method [27] in which track candidates forming pairs with an electron-positron invariant mass lower than 50 MeV/c 2 are removed from any further pairing. In the TPC, the electron candidates are required to have at least 70 out of 159 possible space points attached to the track, which ensures good tracking and particle identification resolution.
Electrons are identified by requiring that the measured dE/dx in the TPC lies within a ±3σ e band around the expected value for electrons estimated from a parameterization of the Bethe formula, where σ e is the particle identification (PID) resolution in the TPC for electrons. The hadron contamination is reduced by excluding tracks compatible with the proton or pion hypothesis within ±3.5σ p,π .
The number of observed J/ψ is obtained by constructing the invariant mass distribution of all combinations of opposite-sign (OS) electron pairs from the same event. The top panels of Figure 1 show the invariant mass distributions obtained in central (0-10%, left) and peripheral (60-90%, right) collisions together with the estimated background. The background is obtained using the distribution of OS pairs constructed by pairing electrons and positrons from different events, so-called mixed events (ME), which is scaled to match the same-event OS invariant mass distribution in two mass intervals on either side of the J/ψ signal region: 2.0 < m ee < 2.5 GeV/c 2 and 3.2 < m ee < 3.7 GeV/c 2 , where the J/ψ contribution is expected to be negligible. The raw J/ψ signal is then obtained by integrating the background-subtracted distribution in the mass window 2.92-3.16 GeV/c 2 . The lower panels of Figure 1 show the OS invariant mass distribution after background subtraction. Good agreement with the J/ψ invariant mass distribution from Monte Carlo (MC) simulation, normalized to the integral of the raw signal, is observed. The potentially remaining correlated background from semi-leptonic decays of cc and bb pairs is included in the systematic uncertainty. J/ψ nuclear modification factor in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration The corrected J/ψ p T -differential production yield is obtained in a given centrality class as where N J/ψ is the number of reconstructed J/ψ in the considered centrality class and p T and y intervals, N ev is the corresponding number of events, A × ε is the acceptance and efficiency correction factor and BR J/ψ→ee = (5.971±0.032)% is the branching ratio of the J/ψ decaying into the dielectron channel [42].
MC simulations of Pb-Pb collisions with embedded unpolarized J/ψ mesons are used to obtain the A × ε factors. The Pb-Pb collisions are generated using HIJING [43]. For the J/ψ, the prompt component is generated using a p T distribution tuned to the existing Pb-Pb measurements at forward rapidity while the non-prompt component is obtained from bb pairs generated with PYTHIA forced to decay into channels with J/ψ in the final state. The J/ψ decays into the e + e − channel are handled using PHOTOS [44]. The transport of the simulated particles in the detector material is performed using a GEANT3 [45] model of the ALICE apparatus and the same algorithm as for the real data is used to reconstruct the simulated tracks. The acceptance times efficiency correction factors include the kinematic acceptance, the reconstruction and PID efficiencies, and the fraction of signal in the integrated invariant mass window. The acceptance correction factor amounts to 33% and the fraction of the signal in the mass counting window is approximately 65%. Reconstruction and PID efficiencies are centrality dependent and together amount to approximately 24% in the most central collisions growing monotonically to approximately 32% in the most peripheral collisions. The correction factors are also p T dependent which, for large p T intervals, induces a dependence on the Monte Carlo p T distribution of the embedded J/ψ. This is taken as a source of systematic uncertainty and is discussed in the following section. The inclusive p T -integrated J/ψ production is measured in 5 different centrality classes: 0-10%, 10-20%, 20-40% , 40-60% and 60-90% while the p T -differential cross sections are obtained in larger centrality classes to ensure sufficient statistical significance: 0-20%, 20-40% and 40-90%.
The average transverse momentum of J/ψ, p T , is extracted using a binned log-likelihood fit of the p ee T distribution of all electron pairs as a function of the invariant mass. Each pair contribution is weighted by the (A × ε) −1 factor corresponding to its centrality and p T . The OS p ee T distribution is fitted with the function: where N bkg (m ee ) and N J/ψ (m ee ) are the mass-dependent distributions of background and signal pairs determined via the signal extraction procedure described above. The background mean transverse momentum, p bkg T , depends on the invariant mass and its shape is obtained from the ME technique, while its overall normalization can vary in the fit. Figure 2 illustrates the p T extraction procedure for the most central and most peripheral centrality intervals using dielectron pairs in the transverse momentum interval 0.15 < p T < 10 GeV/c. A similar procedure is employed also for the second moment of the transverse momentum distribution p 2 T . A low-p T cut-off is applied due to the observation of a J/ψ excess for reconstructed p T < 0.3 GeV/c at forward rapidity in peripheral Pb-Pb collisions at √ s NN = 2.76 TeV [46], which is found to originate from coherent photo-production. Since this production mechanism is not normally included in hadro-production models, the low-p T interval is excluded for enabling comparisons with theoretical calculations. At midrapidity, mainly due to a better momentum resolution, nearly all of the coherent yield is contained in the range of reconstructed p T < 0.15 GeV/c, as shown by the ALICE measurements of J/ψ photo-production in ultra-peripheral collisions [47]. A small component of incoherently photo-produced J/ψ is still present in the range p T < 1 GeV/c, but for the centrality intervals considered in this work it is negligible. Thus, in the following, unless otherwise specified, all the results refer to J/ψ with p T larger than 0.15 GeV/c.   17.9 ± 0.3 0.31 ± 0.01 0-90 125.9 ± 1.0 6.28 ± 0.12 The inclusive J/ψ nuclear modification factor is computed for a given centrality class as where d 2 N/dydp T is the inclusive J/ψ yield defined in Equation 1, the T AA is the average nuclear overlap function corresponding to the considered centrality class and d 2 σ pp /dydp T is the inclusive J/ψ cross section measured by ALICE in pp collisions at √ s = 5.02 TeV [27]. The values used for the nuclear overlap function are shown in Table 1 and are obtained from [40].

Systematic uncertainties
The systematic uncertainties on the measured J/ψ yields, p T , and p 2 T originate from uncertainties on tracking, electron identification, signal extraction procedure, the kinematics used in the MC simulation for estimating the A × ε corrections, and the J/ψ decay branching ratio. For the R AA and the r AA , the uncertainties on the J/ψ cross section measurement in pp collisions [27] and (only in the case of the R AA ) on the nuclear overlap function T AA [40] need to be considered in addition. A summary of the uncertainties on the p T -integrated and p T -differential yields is given in Table 2.
The systematic uncertainty on the tracking of the candidate electrons is mainly due to uncertainties on the ITS-TPC track matching and on the track reconstruction efficiency in both the ITS and the TPC. These uncertainties, mainly due to differences in the reconstruction efficiency between data and MC, are √ s NN = 5.02 TeV ALICE Collaboration Table 2: Systematic uncertainties on the p T -integrated and on p T -differential J/ψ yields for different centrality intervals. Only the ranges of uncertainty are quoted over the considered centrality intervals. The individual contributions and the total uncertainties are given as percent values. Source n.a. n.a. pp reference 7 9-12 3 5 estimated by varying the main track selection criteria and repeating the whole analysis chain. All variations which provide a corrected yield that deviates from the yield obtained with the standard selection criteria by more than one standard deviation are considered [48]. The tracking uncertainty is then obtained as the root-mean-square of the distribution of all the valid variations, while the distribution mean is used as the central value. For the J/ψ yields, this uncertainty ranges between 2% and 7% as a function of centrality (integrated over p T ) and between 4% and 9% as a function of transverse momentum. The tracking systematic uncertainty on the p T and p 2 T are smaller than those for the corrected yields and detailed in Table 2.
Uncertainties on the electron identification are due to the TPC electron PID response and the hadron rejection. A data-driven procedure is used to improve the matching between data and simulation for the electron selection by employing a pure sample of electrons from tagged photon conversions in the detector material. The residual mismatches are estimated by varying all the PID selection criteria following a similar procedure as for the tracking systematic uncertainty. The extracted uncertainty on the J/ψ yields ranges between 1% and 6%, depending on the centrality and transverse momentum interval.
The uncertainty from the signal extraction procedure includes two components, one due to the J/ψ signal shape and one due to the residual correlated dielectron background in the J/ψ mass region. In order to estimate the uncertainty on the signal shape, the corrected J/ψ yields are estimated using variations of the mass region where the number of raw J/ψ are counted. The correlated dielectron background in the invariant mass range used to extract the signal has generally a different shape compared to the combinatorial background. So matching the ME background in the sidebands of the same-event OS distribution may lead to a bias in the estimation of the raw yields. This is taken as a systematic uncertainty and is estimated by varying the mass ranges of the sidebands where the ME background is matched. The total uncertainty on the signal extraction ranges between 2% to 7% as a function of centrality and between 5% and 7% as a function of p T .
The acceptance and efficiency correction is p T dependent which makes correction factors averaged over large p T intervals sensitive to the J/ψ p T spectrum used in the simulation. Since precise measurements of the J/ψ transverse momentum spectra at midrapidity down to p T = 0 GeV/c are not available, the simulations used for corrections rely on the ALICE measurement at forward rapidity (2.5 < y < 4.0) in Pb-Pb collisions at √ s NN = 5.02 TeV [29]. The measured spectrum including the statistical and systematic uncertainties is fitted using a power law function and the fit parameters are varied randomly within their allowed uncertainties taking into account their correlation matrix. The resulting uncertainty amounts to 2% for the p T -integrated corrected yields and ranges between 1% to 2% in the considered p T intervals.
Systematic uncertainties on the extraction of p T and p 2 T are obtained by repeating the fit procedure with similar variations of tracking and PID selections as for the yield estimation. Since for this measurement the A × ε correction is applied for each dielectron pair using a fine-binned distribution of the J/ψ nuclear modification factor in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration correction factors, the systematic uncertainty due to the kinematics of the J/ψ used in the MC simulation is negligible. In addition, the p T and p 2 T are also extracted by directly fitting the corrected J/ψ spectrum with a power law function and are found to be compatible to the values obtained from the fit with Equation 2.
Systematic uncertainties on the tracking, PID and MC kinematics are considered to be partly correlated over both the centrality and the transverse momentum. The systematic uncertainties on signal extraction are considered as uncorrelated. The uncertainties on the nuclear overlap function are taken as uncorrelated over centrality and fully correlated over p T within a given centrality interval. The uncertainty on the p T -integrated pp reference is considered to be fully correlated over centrality, while the uncertainties on the p T -differential values are fully correlated over centrality and highly correlated over p T .

Results and discussions
The inclusive J/ψ p T -differential yields evaluated using Equation 1 are shown in the left panel of Figure 3 (left) for the 0-20%, 20-40%, and 40-90% centrality intervals. The vertical error bars indicate statistical uncertainties while the systematic uncertainties, independently of their degree of correlation, are shown as boxes around the data points. The horizontal error bars show the evaluated p T -range with the data point placed in the center.
In the SHM, all heavy quarks are produced during the initial hard partonic interactions followed by their thermalization in the QGP and the subsequent formation of bound states at the phase boundary according to their thermal weights. The p T -integrated charmed-hadron yields depend only on the total cc cross section in heavy-ion collisions and on the chemical freeze-out parameters, which are determined by fitting measured light-flavored hadron yields. In addition to the high-density core part in the QGP, a corona contribution is added for the case that the nuclear density decreases below 10% of its maximal value, where no QGP is assumed and the number of J/ψ is calculated from yields in pp collisions scaled by the number of binary nucleon-nucleon collision. A recent update of the SHM [49] uses a MUSIC (3+1)D [51] hydrodynamical simulation to extract the transverse flow velocity and the radial velocity profile of the freeze-out hyper-surface, such that the J/ψ p T can be extracted from a blast-wave J/ψ nuclear modification factor in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration parameterization which follows a Hubble-like expansion [52].
The CIM [23] was developed specifically for the description of charmonium suppression in heavy-ion collisions via its interactions with a comoving medium, either hadronic or partonic. The hot medium effects are modeled using a rate equation which contains a loss term for charmonium dissociation, and a gain term for (re)generation. In this model, the charmonium dissociation rate depends on the density of comovers, obtained from experimental measurements and on the charmonium dissociation cross section which is an energy-independent parameter of the model, fixed from fits to low energy data. Charmonium dissociation is balanced by the (re)generation component which depends on the primordial charm-quark cross section.
Both microscopic transport models considered here, TM1 [24] and TM2 [25], solve the Boltzmann equation for charmonia (J/ψ, χ c and ψ ) with dissociation and recombination terms. Each model considers the fireball evolution using implementations of ideal hydrodynamics which include both the deconfined and the hadronic phase separated by a first order phase transition. The dissociation rate in both models depends on the medium density and on a lattice-QCD-inspired charmonium binding energy (in TM1) or squared radius of the bound state (in TM2), all of them being functions of temperature. The (re)generation component is implemented using different approaches. In the TM1 calculations, it is based on the assumption that the charm quarks reach statistical equilibration after a relaxation time of about a few fm/c, while in the TM2 calculations the charm quarks are recombined using the same cross section as for the dissociation process and a thermalized distribution of charm quarks.
The primordial cc production cross section in Pb-Pb collisions is a common input for all of the above mentioned models. There is so far no measurement of the cc cross section in Pb-Pb or pp collisions at √ s NN = 5.02 TeV at midrapidity, which lead dominantly to the uncertainty of the models. The cross section in Pb-Pb collisions is obtained from the total cc cross section in pp collisions dσ cc /dy scaled by the average number of nucleon-nucleon collisions N coll in a given centrality class of Pb-Pb collisions with additional CNM effects taken into account. For the rapidity interval used in this work, |y| < 0.9, the value of dσ cc /dy estimated for MB Pb-Pb collisions is 0.53 ± 0.10 mb for the SHM, 0.76 ± 0.13 mb for TM1, 0.78 ± 0.09 mb for TM2 and 0.56 ± 0.11 mb for CIM.
The right panel of Figure 3 shows a comparison of the inclusive J/ψ transverse momentum spectrum in the 20% most central Pb-Pb collisions to calculations from the SHM and TM2 models. The bands indicate model uncertainties mainly due to the assumptions on the dσ cc /dy. Good agreement between data and the SHM predictions is observed in the low-p T region, while for p T 5 GeV/c the calculations underestimate the data. The TM2 calculations underestimate the measured yields over the measured p T range.
In order to facilitate the comparison of the J/ψ p T spectra obtained in this work with other measurements or theory calculations, the J/ψ p T and p 2 T are extracted in several centrality intervals, using the method described in Sec. 3. The left panel of Figure 4 shows the J/ψ p T dependence on the mean number of participant nucleons N part . The p T in Pb-Pb collisions at √ s NN = 5.02 TeV shows a monotonic decrease from the most peripheral collisions, where it is compatible to the measurement in pp collisions at √ s = 5.02 TeV, to the most central collisions, which hints towards a strong contribution from (re)combination processes. This trend is not clearly visible for the measurement at √ s NN = 2.76 TeV, which suffered from large statistical and systematic uncertainties.
The r AA ratio, defined as p 2 T PbPb / p 2 T pp , which is shown in the right panel of Figure 4, is a measure of the broadness of the p T spectra in heavy-ion collisions relative to pp collisions at the same energy. A strong decrease of the r AA is observed in Pb-Pb collisions at √ s NN = 5.02 TeV between peripheral, where it is consistent with unity, and central collisions where r AA reaches a value of 0.6 at midrapidity and 0.75 at forward rapidity [53]. When comparing with measurements at lower energies from RHIC [11,55,56] and SPS [54], a very different picture emerges. While the RHIC measurements for both p T and r AA are J/ψ nuclear modification factor in Pb-Pb collisions at √ s NN = 5.   [28], Pb-Pb collisions at 5.02 at forward rapidity [53], and with those at lower collision energies at SPS [54] and RHIC [11,55,56]. The red box around unity at N part ≈ 400 in the right panel indicates the correlated uncertainty of the ALICE data points due to the p 2 T in pp collisions.  [27] and Pb-Pb collisions at √ s NN = 5.02 TeV at midrapidity as a function of the mean number of participating nucleons. The ALICE results are compared with calculations from the transport models [57, 58] and the SHM [49]. The colored bands represent model uncertainties. As in Figure 4, the red box around unity at N part ≈ 400 in the right panel indicates the correlated uncertainty of the ALICE data points due to the p 2 T in pp collisions.
compatible with a constant trend as a function of N part [25], the SPS results show a monotonic increase of both p T and r AA as a function of collision centrality which, at this energy, can be explained by a broadening of the p T distribution due to the Cronin effect [59].
The results for the p T and r AA in Pb-Pb collisions at √ s NN = 5.02 TeV are compared with model calculations in Figure 5. The statistical hadronization model agrees with the data only for the most central collisions but underestimates the measurements for more peripheral collisions. A good description of the centrality trend is obtained with the transport model TM1 calculation, which includes a detailed implementation of the fireball evolution, with the exception of most central collisions where the model overestimates both the J/ψ p T and r AA . Besides the most central collisions where there is a hint of an increase of the R AA with collision energy, the results at the two energies are compatible within uncertainties. A comparison of the experimental results at √ s NN = 5.02 TeV with calculations based on the models described before is shown in the right panel of Figure 6. The calculations are shown as bands that indicate model uncertainties, dominated by the uncertainties on the cc cross section and on the CNM effects. The SHM calculation shows a good agreement with the data over the entire centrality range. CIM, TM1 and TM2 calculations underestimate the experimental results towards the data points corresponding to the most central collisions despite the fact that the total cc cross section assumed in TM1 and TM2 is significantly larger compared to the SHM and the CIM. The large model uncertainties do not allow a conclusion to be made on the phenomenology of charmonium production in nuclear collisions. This emphasizes the importance of a precise measurement of the total cc cross section, but also the need of using consistent model inputs, including the total cc cross section, the pp reference J/ψ cross section and CNM effects.
The inclusive J/ψ nuclear modification factor in Pb-Pb collisions as a function of p T is shown in the left panel of Figure 7 for the centrality intervals 0-20%, 20-40% and 40-90%. The systematic uncertainties shown as boxes around the data points include the systematic uncertainties from the Pb-Pb analysis while the uncertainties from the pp reference, correlated over centrality, are shown as the gray band around unity. The colored boxes at high p T around unity indicate the correlated uncertainties due to the T AA values used for the R AA calculation. These results are compatible with binary scaling for p T < 3 GeV/c, with the exception of the data point around 2 GeV/c which shows a downward statistical fluctuation for 40-90% centrality, while the J/ψ production is suppressed at higher p T . With the current uncertainties it is difficult to extract a centrality trend except for the highest p T interval, 5-10 GeV/c, where a stronger suppression is observed in the most central collisions relative to the more peripheral centrality intervals at a significance level of about 3σ . The results for the 20% most central collisions are compared with model calculations and shown in the right panel of Figure 7. Both the SHM and TM1 models describe qualitatively the data. In these models, the increasing R AA towards low p T is a consequence of the dominant contribution of (re)generated J/ψ. At high p T , the contribution from recombination decreases, and the J/ψ production is suppressed due to color charges in the medium. The main J/ψ sources at high p T are primordial production and feed-down from beauty decays. The SHM, J/ψ nuclear modification factor in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration where the J/ψ at high p T are produced only in the corona, overestimates the degree of J/ψ suppression.
Since the charm quark density, i.e. the cc cross section, is expected to decrease towards larger rapidity, the comparison to the forward-rapidity measurements is a valuable source of information. In the left panel of Figure 8, the p T dependence of the J/ψ R AA in the 20% most central Pb-Pb collisions at midrapidity is compared with the ALICE results measured at forward rapidity (2.5 < y < 4) [29]. The boxes around the data points represent systematic uncertainties, while the boxes drawn around R AA = 1 show global uncertainties on the pp reference due to uncertainties on the beam luminosity and T AA . In the low-p T range (p T < 5 GeV/c) these data indicate larger R AA values at midrapidity compared to those at forward rapidity, with a combined statistical significance of nearly 4σ , compatible with expectations from a (re)generation scenario due to the larger primordial cc density at midrapidity. The rapidity dependence of the inclusive J/ψ suppression, integrated over p T , is shown in the right panel of Figure 8 for the 0-90% centrality interval. The value of the R AA at midrapidity is 0.97 ± 0.05(stat.) ± 0.1(syst.), and a monothonic decrease is observed towards forward rapidity [29,53].

Conclusions
The measurements of the inclusive J/ψ yields and nuclear modification factors at midrapidity (|y| < 0.9) were performed in the dielectron decay channel in Pb-Pb collisions at a center-of-mass energy √ s NN = 5.02 TeV using an integrated luminosity of L int ≈ 10 µb −1 collected by the ALICE Collaboration. The results were presented as a function of transverse momentum in different collision centrality classes.
The J/ψ transverse momentum dependent yields in central Pb-Pb collisions are well reproduced in the low p T range by an updated SHM calculation [49] and underestimated for large p T . The TM2 [25] transport calculations underestimate the J/ψ yields over the entire measured p T range. The J/ψ p T and p 2 T show a decrease from peripheral collisions, where they are similar to the values observed in pp collisions, towards most central collisions. This centrality-dependent behavior is qualitatively different compared to the observations at lower energies from RHIC and SPS and can be explained through the interplay between the (re)generation process, dominant at low p T for central events at the LHC, color screening, and CNM effects like gluon shadowing. A good description of the observed trends is provided by the TM1 calculations, while the SHM calculations agree with the data for central collisions only.
The p T -integrated nuclear modification factor as a function of the number of participant nucleons shows a moderate level of suppression in the range 50 < N part < 300, and indicates an increase towards cen-J/ψ nuclear modification factor in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration TeV as a function of p T , at midrapidity and at forward rapidity [29,53]. Right: Rapidity dependence of the inclusive J/ψ R AA in the centrality interval 0-90%. The error bars represent statistical uncertainties, while the boxes around the data points represent systematic uncertainties. The boxes around unity represent global uncertainties on the pp reference due to normalization and T AA . In the right panel, the correlated uncertainty of the point at midrapidity is included in the box around the data point. tral collisions. In the most peripheral collisions, our results are compatible with binary scaling of the J/ψ production. The nuclear modification factor as a function of the transverse momentum shows a strong suppression, centrality dependent, for p T > 3 GeV/c but is compatible with unity or with a small enhancement at small p T , suggestive of the large contribution from the (re)generation process. Furthermore, from these measurements we observe significantly larger values for R AA compared to the results at forward rapidity [29] for both the p T -integrated values in the 0-90% centrality interval and for the p T -differential R AA in the low p T region (p T < 5 GeV/c) in the centrality interval 0-20%.
Consequently, these results strenghten the hypothesis that charmonium at low p T is produced predominantly via (re)generation in the late stages of the collision at the LHC. However, due to the remaining experimental and theoretical uncertainties, the exact phenomenology leading to these observations cannot be determined yet.        [20] P. Braun-Munzinger and J. Stachel, "(Non)thermal aspects of charmonium production and a new look at J/ψ suppression", Phys. Lett. B490 (2000)  [57] X. Zhao