Measurements of inclusive J/ ψ production at midrapidity and forward rapidity in Pb–Pb collisions at √ s NN = 5.02 TeV

The measurements of the inclusive J / ψ yield at midrapidity ( | y | < 0.9) and forward rapidity (2.5 < y < 4) in Pb–Pb collisions at √ s NN = 5 . 02 TeV with the ALICE detector at the LHC are re-ported. The inclusive J / ψ production yields and nuclear modification factors, R AA , are measured as a function of the collision centrality, J / ψ transverse momentum ( p T ), and rapidity. The J / ψ average transverse momentum and squared transverse momentum ( ⟨ p T ⟩ and ⟨ p 2T ⟩ ) are evaluated as a function of the centrality at midrapidity. Compared to the previous ALICE publications, here the entire Pb–Pb collisions dataset collected during the LHC Run 2 is used, which improves the precision of the measurements and extends the p T coverage. The p T -integrated R AA shows a hint of an increasing trend towards unity from semicentral to central collisions at midrapidity, while it is flat at forward rapidity. The p T -differential R AA shows a strong suppression at high p T with less suppression at low p T where it reaches a larger value at midrapidity compared to forward rapidity. The ratio of the p T -integrated yields of J / ψ to those of D 0 mesons is reported for the first time for the central and semicentral event classes at midrapidity. Model calculations implementing charmonium production via the coalescence of charm quarks and antiquarks during the fireball evolution (transport models) or in a statistical approach with thermal weights are in good agreement with the data at low p T . At higher p T , the data are well described by transport models and a model based on energy loss in the strongly-interacting medium produced in nuclear collisions at the LHC


Introduction
Quantum chromodynamics (QCD) is the theory describing the strong interaction.Lattice QCD, i.e. the discrete formulation of QCD, predicts the existence of a state of deconfined matter at high energy density that is characterised by quark and gluon degrees of freedom [1,2].This quark-gluon plasma (QGP) is created during the early hot and dense stage of heavy-ion collisions at ultra-relativistic energies.
The heavy quarks, charm (c) and beauty (b), are unique probes for this phase of matter [3,4].Due to their large masses, they are produced as quark-antiquark pairs in hard partonic scattering processes in the early stage of the collision, and they thus experience the full evolution of the system.While the majority of the produced heavy quarks and antiquarks hadronise independently into open heavy-flavour hadrons, bound quarkonium states can also be formed [5].However, it has been predicted that the formation of bound states should be suppressed due to the mechanism of colour screening where the large density of colour charges in the QGP hinders the production of bound quarkonia [6,7].The degree of suppression of the various quarkonium states depends on their binding energy along with the medium properties, such as its temperature.Consequently, the measurements of quarkonium production rates in heavy-ion collisions have been considered as a potential thermometer of the medium.
The production of J/ψ, the charmonium ground state with quantum numbers J PC = 1 −− , has been studied extensively in heavy-ion collisions over the last several decades.Suppression of the J/ψ yield was observed in nucleus-nucleus collisions with respect to the expectation from proton-proton collisions at the SPS (up to centre-of-mass energy per nucleon pair √ s NN = 17 GeV) [8-10], at RHIC ( √ s NN up to 200 GeV) [11][12][13][14] and at the LHC ( √ s NN = 2.76 TeV [15-18] and 5.02 TeV [19][20][21][22][23][24]).However, contrary to the prediction from the colour-screening scenario, the measured suppression does not increase with increasing collision energy from RHIC to LHC despite of an increased energy density of the produced QGP.At LHC energies, J/ψ production is found to be less suppressed than at the lower RHIC energies, in particular at low transverse momentum, p T , of the J/ψ [18,[25][26][27].In addition, a significant azimuthal anisotropy in the J/ψ production was observed via the elliptic and triangular flow measurements reported in Refs.[28,29].These observations are explained by an additional J/ψ production mechanism, referred to as (re)generation in the following, in which copiously produced uncorrelated charm quarks and antiquarks bind into J/ψ mesons [30,31].This process can only take place in a deconfined medium, and its contribution to the measured J/ψ yield increases with the density of cc pairs and, therefore, with increasing collision energy and decreasing p T [32][33][34].With increasing p T , (re)generation becomes less relevant for J/ψ production and, instead, charmonium dissociation and the fragmentation of high-energy partons into charmonia become dominant.In the latter case, the suppression of high-p T J/ψ yields should reflect the energy loss of partons [35], which is mostly of radiative nature in this kinematic regime.
The formation process of charmonia in heavy-ion collisions is complex and various phenomenological approaches are considered.In the statistical hadronization scenario, the relative abundances of charmomium states with respect to other charmed hadrons are determined by thermal weights [30,32] at the system chemical freeze-out.In microscopic transport and comover interaction models, charmonia are continuously produced and broken up during their propagation through the QGP [31,33,34,36].Furthermore, it is important to consider cold nuclear matter (CNM) effects.In particular, the modification of the parton distribution functions in nuclei with respect to nucleons [37] has to be taken into account for the interpretation of the results.These CNM effects were investigated in ALICE especially with proton-nucleus collisions [38][39][40][41][42][43][44].
For a better assessment of the production mechanisms, systematic measurements of the centrality, p T , and rapidity dependence of J/ψ production are pivotal.In this article, the ALICE results on inclusive J/ψ production at midrapidity (|y| < 0.9) for 0.15 < p T < 15 GeV/c and forward rapidity (2.5 < y < 4) for 0.3 < p T < 20 GeV/c, from the full Run 2 data sample at the LHC, are reported.Inclusive J/ψ At midrapidity, the main detectors employed in the analysis are the Time Projection Chamber (TPC) [47] and the Inner Tracking System (ITS) [48], both immersed in a uniform magnetic field of 0.5 T provided by a solenoid magnet.The TPC is used for tracking and particle identification.It covers the pseudorapidity range |η| < 0.9 for tracks with full radial length and has full coverage in azimuth.It provides excellent momentum resolution and electron-hadron separation in a wide range of track transverse momentum.The ITS is a cylindrical six-layer silicon detector, with the innermost layer located at 3.9 cm from the beam pipe, providing additional space points for tracking that enhance the spatial resolution in the reconstruction of primary and secondary vertices.
At forward rapidity, the muon spectometer [49,50] detects muons in the range −4 < η < −2.5.It consists of a 3 Tm dipole associated with a tracking and a trigger system.A front absorber with a thickness of 10 interaction lengths is placed before the tracking system in order to filter out hadrons produced in the interaction.The tracking system consists of five tracking stations, each one made of two planes of cathode pad chambers.An iron wall with 7.2 interaction length thickness is located between the tracking and trigger stations in order to stop secondary hadrons escaping the front absorber and low momentum muons produced predominantly from π and K decays.The trigger system consists of two stations, each one made of two planes of resistive plate chambers.Finally, a conical absorber around the beam pipe protects the spectrometer against secondary particles produced by the interaction of primary particles with large pseudorapidity at the beam pipe.In the analysis at forward rapidity, the determination of the primary vertex of the collision is provided by the Silicon Pixel Detector (SPD) that constitutes the two innermost layers of the ITS.
Both analyses, at midrapidity and forward rapidity, use the V0 [51] and Zero Degree Calorimeter (ZDC) [52] detectors.The V0 detector consists of two scintillator detector arrays and covers the full azimuth in the pseudorapidity regions −3.7 < η < −1.7 and 2.8 < η < 5.1, respectively.It is used for triggering, beam-gas background rejection, and characterisation of the event centrality.The ZDC detectors are located at a distance of 112.5 m on both sides of the interaction region along the beam direction, and they detect spectator nucleons emitted at zero degree with respect to the LHC beam axis.They are used to reject electromagnetic Pb-Pb interactions.
The trigger for minimum-bias (MB) events was provided by the coincidence of signals in the two scintillator arrays of the V0 detector.The dimuon analysis relies on a dimuon trigger which requires, in addition to the MB trigger, the detection of two opposite-sign tracks in the muon trigger system.The muon trigger Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration selects muon candidates with a p T larger than a threshold of ∼ 1 GeV/c.The trigger efficiency reaches 50% at this threshold value and a plateau value of 98% at p T ∼ 2.5 GeV/c [53].
The results presented in this article are based on the data sample collected by ALICE from Pb-Pb collisions at √ s NN = 5.02 TeV in 2015 and 2018 during Run 2 at the LHC.During the 2018 data taking, the Pb-Pb dataset for the central barrel was enhanced with central (0-10%) and semicentral (30-50%) events.In total, the integrated luminosity corresponding to the analysed data sample was about 105 µb −1 and 51 µb −1 for the central and semicentral events, respectively.For the other centrality intervals, the integrated luminosity of the data sample was 22 µb −1 .For the analysis at forward rapidity, the dimuon triggered sample corresponds to an integrated luminosity of 756 µb −1 .At midrapidity, triggered events containing collisions that overlap within a time window smaller than the readout time of the TPC were removed to preserve a uniform particle identification performance of the TPC, which is sensitive to the total charge produced by the ionising tracks in the sensitive volume.Only events with the primary vertex, reconstructed within ±10 cm from the nominal interaction point in the beam direction, were considered for further analysis at midrapidity.In the forward analysis, there was no selection on the primary vertex.

Analysis details
The primary observable is the p T -differential J/ψ yield per unit of rapidity d 2 N/(dy dp T ).For a given interval of centrality, rapidity (∆y), and transverse momentum (∆p T ), this is obtained as where N J/ψ is the number of reconstructed J/ψ mesons, N ev is the number of events corresponding to the analysed centrality interval, and (A × ε) is the acceptance times efficiency factor.The branching ratio (BR) corresponds to either the dielectron or the dimuon J/ψ decay channel.Since the analysis at forward rapidity was based on a sample of dimuon-triggered events, the equivalent N ev was 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 triggered event, F norm [40,54].The number of equivalent N ev was first obtained for the 0-90% centrality class and was then scaled to the centrality classes considered in the analysis.
The nuclear modification factor, R AA , is obtained as where ⟨T AA ⟩ is the average nuclear overlap function as described in Ref. [55] and given in Table 1 for the centrality intervals used for the analyses at midrapidity and forward rapidity.The R AA is evaluated as a function of the average number of participant nucleons, N part , corresponding to a given centrality class (as shown in Table 1), and as a function of the J/ψ p T .The differential J/ψ production cross section in pp collisions, d 2 σ pp /(dy dp T ), was measured at both midrapidity and forward rapidity as reported in Refs.[56] and [57], respectively.At midrapidity, the data sample size does not allow to obtain the cross section for p T > 10 GeV/c and an extrapolation is applied for the last p T interval (10 < p T < 15 GeV/c), the result of the extrapolation is 6.82 ± 0.91 nb.The details of this approach are described in Ref. [58,59], and the corresponding R AA is shown as the open points in Figs. 6, 7 and 9.
The p T -differential J/ψ yields at midrapidity were studied further by extracting the ⟨p T ⟩ and ⟨p 2 T ⟩ in fine centrality intervals, as described later in this section.For quantitative comparisons of the p T distributions, the ratio r AA of the J/ψ ⟨p 2 T ⟩ measured in Pb-Pb collisions to the one obtained in pp collisions at the same energy is calculated as Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration The J/ψ mesons are reconstructed employing the e + e − decay channel at midrapidity and the µ + µ − decay channel at forward rapidity.The analysis techniques are discussed in detail in Refs.[18,20,21].
Here, only a brief overview is given and differences with respect to previous analyses are highlighted.
Electron candidates for the analysis at midrapidity are tracks reconstructed in both the ITS and TPC in the pseudorapidity range |η| < 0.9 and with a p T > 1 GeV/c to suppress combinatorial background.All tracks are required to have at least one hit in the SPD layers and at least 70 out of a maximum of 159 clusters reconstructed in the TPC.These and other quality criteria that were applied in addition (see Ref. [20]) ensure good tracking resolution and particle identification.They reduce the background electrons from the conversion of photons in the detector material or from long lived weakly-decaying hadrons as well as tracks from pileup collisions, occuring within the readout time of the TPC.Electrons and positrons are identified using selections on the specific energy loss, dE/dx, in the TPC gaseous volume.The measured dE/dx is required to be within 3 standard deviations (σ ) relative to the expected electron specific energy loss corresponding to the track momentum, and more than 3.5σ different from either the π or proton specific energy loss hypotheses.Electrons from photon conversions surviving the track quality criteria are rejected using a technique where candidate electrons are paired with other electrons selected with less strict criteria to enhance the probability of finding the conversion partner, as described in detail in Ref. [21].
Muon candidates were selected such that the track pseudorapidity is within the geometrical acceptance of the muon spectrometer, −4 < η < −2.5, and that the reconstructed track matches a track segment reconstructed in the trigger system.The transverse position of the tracks at the end of the absorber is required to be 17.6 < R abs < 89.5 cm in order to reject tracks crossing the thickest part of the absorber.In addition, a selection was applied on the product of the track momentum and the transverse distance to the primary vertex in order to reduce the contamination produced by particles that do not originate from the interaction point.
The number of reconstructed J/ψ mesons, i.e. the raw J/ψ yield, was obtained by constructing the invariant-mass distribution of all the possible opposite-sign dileptons with rapidity selections of 2.5 < y < 4 for dimuons and |y| < 0.9 for dielectrons.At midrapidity, the signal extraction was performed in two steps.First, the combinatorial background was estimated using an event-mixing technique [21] and subtracted from the invariant-mass distribution.In the second step, the remaining distribution was fitted with a two-component function, one corresponding to the J/ψ signal and the other to the residual Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration background, which mainly arises from correlated semileptonic decays of heavy-flavour hadrons.The J/ψ signal line shape was obtained from Monte Carlo (MC) simulations of J/ψ mesons decaying in the dielectron channel embedded in simulated Pb-Pb collisions, as described below, while for the residual background a second-order polynomial function was employed.The raw J/ψ yield was obtained by first counting the dielectron pairs in the mass range 2.92 < m e + e − < 3.16 GeV/c 2 in the combinatorialbackground subtracted invariant-mass distribution, and then subtracting the residual background based on a two-component fit.Finally, the raw J/ψ-meson yield is corrected for the fraction of J/ψ reconstructed outside of the counting mass interval, as described in more detail in Section 3.3.This procedure is illustrated in the left panels of Fig. 1 for the collision centrality interval 0-5% and p T > 0.15 GeV/c.The upper left panel shows the invariant-mass distributions for the opposite-sign dielectrons constructed from the same event (black) and mixed events (red).The fitting procedure of the combinatorial background subtracted invariant-mass distribution discussed above is illustrated in the lower left panel.At forward rapidity, two different methods were used to extract the number of J/ψ counts.In the first method, the invariant-mass distributions were fitted with a sum of a signal and a background function, while in the second method the event-mixing technique was employed, as described in Ref. [18].The fit functions corresponding to the signal are either a double-sided Crystal Ball function (CB2) or a pseudo-Gaussian with a mass-dependent width [60].In both cases, the J/ψ pole mass and width were free parameters of the fit, while the non-Gaussian tail parameters were fixed.Two sets of tail parameters were obtained, one based on MC simulations and one extracted from a large data sample of pp collisions at √ s = 13 TeV [61].The MC simulations were embedded into real MB events in order to properly account for the effect of the detector occupancy.The ψ(2S) resonance was also included in the fit to the invariantmass spectrum, using the same signal function as for the J/ψ with mass and width bound to those of the J/ψ [61,62].The background functions employed in the first method were either a variable-width Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration Gaussian [60] or a ratio of a second order to a third order polynomial.The residual background in the second method was parameterised with a sum of two exponential functions.Finally, two invariant-mass ranges were considered for the fit procedure: 2.2 < m µ + µ − < 4.5 and 2.4 < m µ + µ − < 4.7 GeV/c 2 .An example of the signal extraction fit is shown in the right panel of Fig. 1 before (upper plot) and after (lower plot) the subtraction of the combinatorial background estimated with the event-mixing technique.For each p T and centrality interval, several fits were performed with the two different approaches, different combinations of signal and background functions, signal tail parameters, and fitting ranges.The number of J/ψ was obtained as the average of the results from the various fitting methods [57].These various fitting methods are used to determine the systematic uncertainties on the yield extraction as described in Section 3.4.
About 9.0×10 5 and 8.2×10 4 raw J/ψ counts are measured at forward rapidity and midrapidity, respectively, integrated over all available centrality and p T intervals.
3.2 J/ψ ⟨p T ⟩ and ⟨p 2 T ⟩ extraction At midrapidity, a quantitative study of the J/ψ p T spectrum in fine centrality intervals is conducted by extracting the J/ψ mean p T , ⟨p T ⟩ J/ψ , and mean p T squared, ⟨p 2 T ⟩ J/ψ .For a given centrality interval, these quantities are obtained based on a fit to the mass dependent ⟨p T ⟩ and ⟨p 2 T ⟩ distributions after efficiency correction, using a function defined as: where X stands for either the ⟨p T ⟩ or ⟨p 2 T ⟩, and f is the invariant-mass dependent fraction of J/ψ signal determined in the signal-extraction procedure explained above.The invariant-mass dependent background component X bkg(m e + e − ) is determined from the event-mixing procedure plus a second order polynomial function for the residual background.Examples of these fits for the ⟨p T ⟩ observable are shown in Fig.The data points correspond to opposite-sign e + e − pairs from the same event, the blue line to the e + e − pairs from mixed events, and the red line is the combined fit that includes the mixed events and residual background which is described by the polynomial function.

Acceptance and efficiency correction
The acceptance times reconstruction efficiency factor (A × ε), which enters into Eq. 1, was computed employing both MC and data-driven methods.At midrapidity, this factor includes the kinematic ac-Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration ceptance, track-reconstruction and particle-identification efficiencies, and the fraction of J/ψ with an invariant mass in the signal counting range.With the exception of the particle identification, obtained with a data-driven method, the corrections were obtained using a MC simulation of J/ψ embedded in simulated Pb-Pb collisions.The Pb-Pb collisions were generated using the HIJING 1.0 model [63], while the J/ψ were generated using a cocktail of prompt J/ψ with a kinematic distribution tuned to existing measurements and non-prompt J/ψ from beauty hadrons forced to decay into channels containing J/ψ, using PYTHIA 6.4 [64].The e + e − decay of the embedded J/ψ was handled using PHOTOS [65].
Both the prompt J/ψ and the beauty hadrons forced to decay into non-prompt J/ψ were assumed to be unpolarised, in agreement with existing measurements, which indicate small or no polarisation [66,67].All generated particles were transported through the ALICE detector setup using GEANT3 [68], taking into account the time dependence of detector conditions during the 2015 and 2018 data-taking periods.
For the determination of the particle-identification efficiency, a clean sample of electrons from photon conversions, passing similar quality selection criteria as primary electrons in the TPC, is used to compute differential maps in pseudorapidity, azimuthal angle ϕ and momentum p for the single-electron selection efficiency.These were then propagated to the J/ψ dielectron pairs using the phase-space distribution of the J/ψ decay simulated in the above mentioned MC simulations.The total (A × ε) for the p T -integrated J/ψ yields is about 6.5% in the 0-10% centrality interval, and it slightly increases towards more peripheral collisions.As a function of p T , the (A × ε) has a non-monotonic behaviour, with a minimum value of 5.6% around p T = 2 GeV/c, and a maximum of about 9% towards zero and high p T .
At forward rapidity, the acceptance and reconstruction efficiency values were determined using simulated J/ψ mesons forced to decay via the dimuon channel, embedded into real events.The J/ψ p T -and y-differential distributions used in the simulation were adjusted to measurements via an iterative procedure, and separately for all centrality intervals employed in this analysis.The J/ψ were assumed to be unpolarised, in agreement with the small polarisation, compatible with zero, measured in Pb-Pb collisions for 2 < p T < 10 GeV/c [67].As in the analysis at midrapidity, the simulations take into account the time dependence of the detector conditions, such as the status of the tracking chambers and the residual detector element misalignment.The trigger-chamber efficiency was determined from data and used as input in the simulations.The (A × ε) reaches a minimum of 11% at p T ≈ 2 GeV/c and increases up to 13% at low p T and up to 46% at high p T in the 0-20% centrality interval.It increases towards peripheral collisions by a few percent [46].

Systematic uncertainties
The considered sources of systematic uncertainty for the analysis at midrapidity include central-barrel tracking, electron identification, signal extraction, and the kinematics of the J/ψ injected in the MC simulations.For the analysis at forward rapidity, the main systematic uncertainties originate from the signal extraction, the muon tracking and trigger efficiencies, and the kinematics of the J/ψ used in the embedded MC simulations.In addition, for the R AA , uncertainties on the pp reference and the nuclear overlap function are included in both analyses [55].The uncertainty on the J/ψ decay branching ratio and the evaluated systematic uncertainties are summarised in Table 2 and Table 3 for the analysis at midrapidity and forward rapidity, respectively.
At midrapidity, the tracking uncertainty is the largest source of systematic uncertainty and it is dominated by the ITS-TPC track matching.This was determined based on the difference in the matching efficiency observed for single tracks between data and MC simulations.The systematic uncertainty due to the electron identification takes into account the residual miscalibration of the TPC particle identification (PID) response and also the statistical uncertainty of the clean electron sample used to compute the identification efficiency.For its estimation, the PID selection criteria (both electron inclusion and hadron rejection) are varied, each time obtaining a new set of raw yields and corresponding PID efficiencies.The assigned systematic uncertainty is taken as the standard deviation of the distribution of the corrected results obtained in this procedure.The systematic uncertainty of the signal extraction includes Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration one component from the J/ψ signal shape obtained from MC simulations and one component related to the description of the dielectron background.The former was determined by varying the mass range in which the signal is counted and recomputing each time the corresponding signal fraction correction, while the latter was determined by repeating the fit to the invariant-mass distribution in different mass ranges.The standard deviation of the corrected yield distribution was then taken as the systematic uncertainty.Since the J/ψ efficiency depends on p T , the average efficiency computed over wide p T intervals depends in turn on the underlying J/ψ p T distribution used in the MC simulations.The corresponding uncertainty was minimised by iteratively tuning the injected p T spectrum to match the corrected spectrum measured in this analysis.A systematic uncertainty which takes into account possible variations of the p T spectrum, statistically compatible with the finally measured p T spectrum, was assigned and it is typically below 1%.The total systematic uncertainty on the J/ψ corrected yield varies in the range 6-10% in different p T and centrality intervals.The systematic uncertainties, which are dominated by the tracking uncertainties, are correlated over centrality and p T to a very large extent.The systematic uncertainties for ⟨p T ⟩ and ⟨p 2 T ⟩ are evaluated via similar procedures as for the p T integrated yields.The uncertainties from signal extraction and track selection criteria range, respectively, from 0.2 to 1.2% and from 0.5 to 1.3%.The electron identification and ITS-TPC matching systematic uncertainties are calculated by propagating the p T -differential systematic uncertainties to the ⟨p T ⟩ based on the measured p T -differential spectrum.
In the analysis at forward rapidity, the systematic uncertainty corresponding to the signal extraction was determined using several variations of the fit to the invariant-mass spectra, including the fit method, the signal and background functions, and the fitted mass range.This uncertainty varies in the range 1.5-10.7%depending on the p T interval and centrality class.The systematic uncertainty on (A × ε) depends on the uncertainty on the p T and y distributions of the simulated J/ψ, and on the tracking, trigger, and matching efficiency.The first two were evaluated by varying the p T and y spectrum for each centrality interval, taking into account the correlations between the p T and y distributions.The systematic uncertainty was estimated as the largest difference between the nominal (A × ε) and the one estimated from the variations.It ranges between 0.2% and 4.1%.The systematic uncertainty on the muon tracking efficiency was estimated based on the difference between the single-muon tracking efficiency obtained in data and MC with a method that uses the redundancy of the tracking information in each station.The corresponding uncertainty for dimuons was evaluated to be 3% and constant over p T .An additional systematic uncertainty is ascribed to the loss of tracking efficiency due to occupancy effects in the most central events and was estimated to range between 0.5 and 1%, increasing towards more central events.The systematic uncertainty on the trigger efficiency has two components, one due to the intrinsic efficiency of the trigger chambers and another one due to the trigger response.The first component was estimated from the uncertainties on the single-muon trigger efficiency measured from data and used in the simulations.The second component was evaluated by comparing the p T dependence of the trigger response function of the single muon between data and MC.The two sources were added in quadrature and the obtained uncertainty ranges between 1.5 and 2%.Finally, a 1% systematic uncertainty is assigned, related to the choice of the χ 2 selection used to define the matching between the tracks reconstructed in the tracking system and the track segments reconstructed in the trigger system.The uncertainty on F norm was estimated by using two methods.First, the opposite-sign dimuon trigger condition was applied when analysing recorded MB events and, second, the counting rate of the dimuon and MB triggers were compared.The estimated uncertainty on F norm was obtained by comparing the two methods and it amounts to 0.7%.
The systematic uncertainties on ⟨T AA ⟩ were obtained as described in Ref. [55] and the values are listed in Table 1 for both analyses, at midrapidity and forward rapidity.The systematic uncertainty on the definition of the centrality interval was estimated using variations of ±1% of the V0 signal amplitude corresponding to 90% of the hadronic Pb-Pb cross section and redefining correspondingly the centrality intervals.The systematic uncertainty of the centrality limit depends on the width of the centrality classes Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration and it ranges from 1 to 6% and from 0 to 2.8%, as shown in Tables 2 and 3 for the analyses at midrapidity and forward rapidity, respectively.
The systematic uncertainties on the J/ψ reference cross section in pp collisions at √ s = 5.02 TeV as obtained in Refs.[57,69] are provided in Table 2 and Table 3.The correlations of the systematic uncertainties over centrality and p T depend on the mid and forward rapidity analysis and they are indicated in Table 2 and 3.

Inclusive J/ψ yields
The fully corrected inclusive J/ψ p T -differential yields, d 2 N/(dydp T ), were obtained according to Eq.1 in centrality intervals in Pb-Pb collisions at √ s NN = 5.02 TeV at midrapidity (|y| < 0.9) and at forward rapidity (2.5 < y < 4). Figure 3 shows the J/ψ yields obtained at midrapidity for the 0-10% and 30-50% centrality intervals, while Fig. 4 shows the yields measured at forward rapidity in the 0-20%, 20-40% and 40-90% centrality intervals.For all of the results, the statistical and systematic uncertainties Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration are indicated by the vertical error bars and the open boxes around the data points, respectively.These results are compared with calculations performed using the statistical hadronisation model (SHMc) by Andronic et al. [32], and two microscopic transport models by Rapp et al. [70] and Zhuang et al. [33].
The physics assumptions in these model calculations are discussed in more detail in the next section.The lower panels of the Figs. 3 and 4 depict the ratio between the experimental data and the different model calculations, with the width of the bands representing the model uncertainties.These uncertainties are due to uncertainties on input parameters, mainly the total charm-quark production cross section and CNM effects.The filled boxes around unity represent the uncertainties of the measured results, shown as the quadratic sum of statistical and systematic uncertainties.Both transport models describe the p T -differential yields in central and semicentral collisions, while they overestimate them at forward rapidity in peripheral collisions.The SHMc calculations coupled with a hydrodynamics inspired freezeout parameterisation are in good agreement with the data in the low-p T region, but underestimate the measurements at higher p T .Data are compared to model calculations from Refs.[32,33,70].The ratios between data and models are shown in the lower panels.The filled boxes around unity depict the quadratic sum of statistical and systematic uncertainties from the measurement, while the bands indicate model uncertainties.[32,33,70].The ratio between data and models is shown in the lower panels.The filled boxes around unity depict the quadratic sum of statistical and systematic uncertainties from the measurement, while the bands indicate model uncertainties.

The J/ψ nuclear modification factor R AA
The nuclear modification factor R AA was obtained using the measured yields, according to Eq. 2. Figure 5 shows the p T -integrated J/ψ R AA as a function of N part in Pb-Pb collisions at √ s NN = 5.02 TeV, obtained in the current analysis at midrapidity in comparison to the results at forward rapidity, previously reported by the ALICE Collaboration in Ref. [19].Global uncertainties represented the centralitycorrelated uncertainties, shown as filled boxes around unity, and are largely uncorrelated between the forward and midrapidity analyses.Both results exclude low-p T J/ψ, with a selection of p T > 0.15 GeV/c and p T > 0.3 GeV/c at midrapidity and forward rapidity, respectively, in order to reject J/ψ produced via photoproduction processes [71][72][73], which contribute significantly to the J/ψ yield in particular in peripheral collisions [74,75].The R AA is compatible with unity in the most peripheral collisions, while √ s NN = 5.02 TeV ALICE Collaboration a suppression of the J/ψ production in Pb-Pb collisions with respect to binary scaled pp collisions is observed in semicentral and central collisions, in particular at forward rapidity.At midrapidity, R AA exhibits a slightly increasing trend from approximately ⟨N part ⟩ = 100 towards the most central collisions, with slightly larger R AA values at midrapidity than at forward rapidity, which confirms previous observations reported by ALICE [21].The results at midrapidity are larger than those measured at forward rapidity, with a significance of the difference of 2.2σ when considering the data points in the centrality range 0-10%.The larger R AA in central collisions at midrapidity is expected in phenomenological models due to the larger dσ cc /dy at midrapidity which leads to larger fraction of J/ψ produced via (re)generation [32,33,76].The p T -differential R AA results are shown in Fig. 6 for midrapidity (left panel) and forward rapidity (right panel) in various centrality intervals.The main feature of these measurements is that the R AA values are relatively large at low p T (p T < 5 GeV/c), in contrast with the strong suppression in the p T > 5 GeV/c range, for central and semicentral collisions.A weaker p T dependence of the R AA values is observed from central to peripheral collisions, up to a constant R AA , within uncertainties, in the 40-90% centrality interval at forward rapidity.Such a behaviour in the data can be qualitatively understood by the dominance of hot nuclear matter effects for central and semicentral collisions, acting on top of the CNM ones, which were discussed in Refs.[39,44].In the low-p T region and in particular in central and semicentral collisions, where the density of charm quarks is larger, the coalescence of charm quark pairs has an important contribution counterbalancing the impact of quarkonium suppression in the QGP.At higher p T , the J/ψ production is dominated by effects such as dissociation and energy loss, which are expected to be stronger in the most central collisions.A comparison of the J/ψ p T -differential R AA in the most central Pb-Pb collisions at √ s NN = 5.02 TeV between midrapidity and forward rapidity is shown in Fig. 7. Neglecting the slight difference in centrality intervals, the R AA is higher at midrapidity with respect to forward rapidity at low p T (p T < 3 GeV/c) with a 2.7σ significance, highlighting the strong dependence of R AA on the local charm quark density, further supporting the picture of quarkonium production via coalescing charm quarks.The two measurements converge to similar values at higher p T suggesting a weaker dependence on rapidity for the suppression effects.
The measurements presented so far are discussed in the rest of this section in comparison to calculations which employ different approaches in modelling the collision fireball and the hot medium effects having an impact on J/ψ production in Pb-Pb collisions.
The statistical hadronisation model by Andronic et al. [32] assumes that all charm quarks are produced  during the initial hard partonic interactions and then thermalize in the QGP.The relative yields of the charmed hadrons are then determined solely by the equilibrium thermodynamical parameters at the chemical freeze-out, which are fixed from fits to the yields of light-flavoured hadrons.A recent extension of the model [30,32] implements a hydro-inspired freeze-out hypersurface which allows the calculation of p T -differential yields in addition to the integrated ones.
The microscopic transport model of Zhuang et al. [33,77] implements a real time evolution of the J/ψ, ψ(2S) and χ c production using a Boltzmann-type rate equation that includes both dissociation and coalescence terms.The dissociation term includes contributions from the temperature dependent colour screening effect and scatterings with thermal partons, i.e. gluon dissociation.The (re)generation of charmonia is implemented by exploiting the detailed balance of the gluon dissociation process.The Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration space-time evolution of the fireball is described using the equations for (2+1)D ideal hydrodynamics.This model includes also the production of non-prompt J/ψ, with the precursor beauty quarks being propagated through the QGP using the Langevin equation.
Similarly to the previously described model, the transport model proposed by Rapp et al. [76,[78][79][80][81], is also based on a kinetic rate equation to compute the time evolution of charmonium (J/ψ, ψ(2S), and χ c ) yields.The dissociation term of the rate equation employs an inelastic parton scattering cross section of charmonia in the QGP, computed using next-to-leading order perturbative QCD Feynman diagrams, and also includes the effect of in-medium reduced binding energy.The (re)generation rate depends on the charmonium dissociation temperature, which is extracted from lattice QCD calculations, and equilibrium limits computed based on the thermal model.The space-time evolution of heavy-ion collisions is simulated by a cylindrically expanding fireball model with the regenerated charmonium p T -spectra being calculated in a thermal blast-wave approximation at a temperature and flow velocity reflecting the average production time of each charmonium state [82].
All of the models described above consider the initial state of the nuclear collisions by making assumptions on the total charm quark density produced during the hard partonic collisions and modified by the CNM effects.The two transport models obtain the total charm density based on the measured total charm cross section in pp collisions [83] multiplied by the number of binary nucleon-nucleon collisions.The CNM effects are introduced via different approaches.The microscopic transport model of Rapp et al. [78,79] estimates CNM effects using fits of the measured p-A data, while the transport model of Zhuang et al. [80,81] uses the nuclear parton distribution functions and their uncertainties from EPS09 [84].The SHMc extracts the total charm cross section from the ALICE measurements of D meson production in Pb-Pb collisions [85].The large uncertainty on the estimation of CNM effects is inherited by these model calculations.
At large p T , the fragmentation of high-energy partons may become the main mechanism for J/ψ production.In that case, energy loss of partons due to multiple scattering in the QGP leads to J/ψ suppression at high p T .In the model by Arleo et al. [35] , the quenching of large-p T particles (p T > 10 GeV/c) is assumed to be mostly due to radiative parton energy loss.In this approach, the p T dependence of R AA is fully predicted from the model proposed by Baier et al. [86,87], which employed a medium-induced gluon spectrum.The R AA value is computed from the mean energy loss, which is extracted from a fit to the charged hadron R AA , measured in various collision systems, the average fragmentation function, and the colour coupling factor of the parton.At forward rapidity, the mean energy loss is further corrected for the charged-hadron multiplicity difference between midrapidity and forward rapidity.The model uncertainties arise from the uncertainties on these inputs.This model does not include the production of non-prompt J/ψ, but the R AA variation, when accounting for this contribution, is expected to lie within the theoretical uncertainties.
The p T -integrated nuclear modification factor measured in Pb-Pb collisions at √ s NN = 5.02 TeV at midrapidity is shown in Fig. 8 in comparison with results from the SHMc and the two transport-model calculations.The calculations are shown as coloured bands, illustrating the uncertainties on the initial effects, mainly CNM effects, described above.Within the model uncertainties, all three predictions agree with the data.One can note though that the data lie on the upper edge of the transport-model calculations, while they are in good agreement with the central values from the SHMc calculations for semicentral and central collisions.
Figures 9 and 10 show the p T -differential R AA measurements for various centrality intervals at midrapidity and forward rapidity, respectively, in comparison with the available model calculations.With the exception of the energy-loss calculations, available only for p T > 10 GeV/c, all of the models cover the full p T range in which these measurements were performed.The SHMc model calculations are in good agreement with the data at low p T at both midrapidity and forward rapidity.However, the R AA is under-Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration estimated for p T > 5 GeV/c in all centrality intervals in both rapidity ranges.This might be attributed to physical sources missing in this approach, such as the contributions from surviving primordial J/ψ or non-prompt J/ψ from beauty-hadron decays, but also to an underestimated amount of radial flow acquired by the charm quarks during the system evolution.A similar conclusion can also be drawn from the comparison of the p T -differential yields in Pb-Pb collisions shown in Fig. 3 where the measured spectrum is harder than the one from the SHMc calculations.The two transport models are in better quantitative agreement with data than the SHMc model.Both of the transport models provide a good description of the R AA at both low and high p T .However, the model calculations in the low-p T region, where J/ψ production is dominated by coalescence in these models, do not describe the detailed shape of the p T dependence of R AA , in particular in semicentral collisions, which points to a still not perfectly understood dynamics of charm-quark coalescence.
The energy-loss calculations by Arleo et al [35], performed in all studied centrality ranges for p T > 10 GeV/c, are in good agreement with the measurements, which, based on the model assumptions, suggests that the dominant mechanism in this kinematic regime is indeed energy loss, similar to that of the other hadrons measured at LHC energies.The data are compared with model calculations from Refs.[32,33,35,70].[32,33,35,70] 4.3 The inclusive J/ψ ⟨p T ⟩ and ⟨p 2  T ⟩ Observables that allow for a more differential study of the J/ψ p T spectrum with respect to the centrality of the collisions are the J/ψ ⟨p T ⟩ and ⟨p 2 T ⟩.The latter is typically quantified using the r AA , defined in Eq.3, which is related to the broadening or narrowing of the J/ψ p T spectrum relative to that in pp collisions.The ⟨p T ⟩ and the r AA measured at midrapidity are shown in Fig. 11 in the left and right panels, respectively, as a function of ⟨N part ⟩.Similar measurements were done at the forward rapidity as well [20].These results are compared with similar measurements in heavy-ion collisions at RHIC [11,26,88] and SPS [89] energies.The ⟨p T ⟩ results are also compared with the value reported by the ALICE Collaboration in pp collisions at midrapidity at √ s = 5.02 TeV [69], showing a good agreement with the value measured in the most peripheral Pb-Pb collisions.The data show that for a given N part , the J/ψ ⟨p T ⟩ grows with increasing collision energy.However, while the data indicate no centrality dependence of the ⟨p T ⟩ at the SPS and RHIC energies, a monotonically decreasing trend from the most peripheral to the most central collisions is observed in the ALICE measurements, which reflects the gradual increase of the low p T (re)generation component.The r AA results support the observations for the ⟨p T ⟩, showing a decrease from unity towards central collisions for the ALICE measurements.The RHIC results are compatible with unity over the whole covered centrality range, while the SPS data indicate a strong increase from peripheral to central collisions, suggesting that CNM effects such as the Cronin effect [90] have an impact on the J/ψ p T shape.
The ⟨p T ⟩ and ⟨p 2 T ⟩ results are also compared with the aforementioned transport model calculations, which show a good agreement with the trends observed in data as demonstrated in Fig. 12 4.4 The J/ψ to D 0 yield ratio A long awaited measurement which helps understanding the details of the J/ψ production in heavy-ion collisions is the ratio between the J/ψ and the D 0 yields, both measured in the same collision system.Such a measurement provides a tight constraint to models because some of the model parameters and most model uncertainties related to the cc cross section cancel in the ratio.This ratio is sensitive to the hadronisation mechanisms of the different charm hadrons.While a model independent measurement of the cc production density in heavy-ion collisions is not available, it is still useful to compare the J/ψ yield with the recently published ALICE measurements of the D 0 yield down to zero p T [85].
Figure 13 shows the measured p T -integrated J/ψ to D 0 yield ratio in central (0-10%) and semicentral (30-50%) collisions.The largest source of systematic uncertainty for both measurements comes from tracking efficiency and it is considered correlated between the D 0 and J/ψ measurements, and Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration consequently cancels in the ratio.Without this uncertainty, the numerical values of the J/ψ yields are 0.12 ± 0.005 (stat.)±0.012 (syst.)and 0.016 ± 0.0008 (stat.)±0.0004 (syst.)for (0-10%) and (30-50%) centrality intervals, respectively.The statistical (systematical) uncertainty of the ratio is the quadratic sum of the statistical (systematical) uncertainties of the two measurements.The results suggest a higher value for this ratio in central compared to semicentral collisions.This is supported by the SHMc calculations [32], which suggests both the J/ψ and D 0 are produced via the coalescence of charm quarks at the phase boundary, the ratio being determined by the charm fugacity.The SHMc model gives a good description of the data.The model uncertainty from the SHMc model is due to uncertainties on the charm fugacity parameter, which is fitted to the ALICE D 0 data [85].

Conclusions
The J/ψ p T -differential yields and nuclear modification factors R AA measured from p T = 0.15 GeV/c up to 15 GeV/c in the 0-10% and 30-50% centrality ranges at midrapidity and from p T = 0.3 GeV/c up to p T = 20 GeV/c in the 0-20%, 20-40% and 40-90% centrality ranges at forward rapidity, and the centrality dependent J/ψ ⟨p T ⟩ and r AA measured at midrapidity, are reported and discussed in comparison with model calculations.
The centrality dependent p T -integrated R AA in peripheral collisions shows similar values at both midrapidity and forward rapidity, while a hint of an increasing trend of the R AA is observed at midrapidity towards central collisions.When looking at the p T -differential R AA , relatively large values are observed at low p T , which are compatible with unity for p T < 3 GeV/c in the 0-10% centrality interval at |y| < 0.9, while a strong nuclear suppression is seen at higher p T in central and semicentral collisions.In addition, R AA is higher at midrapidity than at forward rapidity for p T < 3 GeV/c in the most central collisions.A weaker p T dependence of the R AA values is observed for more peripheral collisions.Such a behaviour, for the central and semicentral collisions, can be explained by a large contribution from (re)generation to the J/ψ yields.This is supported by the statistical hadronisation model and by two microscopic transport model calculations when compared to the data.However, the large model uncertainties, which are mainly due to the assumptions on the collision initial conditions, prevent from drawing a clear conclusion on the phenomenology of the J/ψ production in heavy-ion collisions at LHC energies.The J/ψ nuclear modification factor at high p T is well described by the transport models and also by a J/ψ energy loss model, while it is largely underestimated in the hydro-inspired freeze-out approach implemented together with Inclusive J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration the SHMc model.For most central events, the R AA converge to similar values at high p T at mid and forward rapidity suggesting a weaker dependence on rapidity for the suppression effects.The centrality dependent J/ψ ⟨p T ⟩ and r AA measurements are compared with similar results at lower energies from RHIC [11,88,91] and SPS [89], with the centrality trends showing an opposite behaviour between the LHC and the lower energy results.This behaviour is compatible with a strong contribution from the regeneration component which tends to soften the p T distributions.The two microscopic transport models describe the data within the uncertainties.
The ratio of inclusive p T -integrated J/ψ to D 0 yields measured by ALICE at midrapidity in Pb-Pb collisions is presented for the first time in the 0-10% and 30-50% centrality ranges.The data shows a larger value of this ratio in 0-10% compared to 30-50% collisions, which is in good agreement with the expectation from the SHMc model that fixes the charm fugacity parameter based on the D 0 yields measured by ALICE.
The large improvements in experimental accuracy expected for the LHC Run 3 and 4 for charmonium measurements and more in general for heavy-quark production, the improved measurements of total charm quark cross section and CNM effects are critical to constrain the phenomenological model calculations, will allow to settle the longstanding questions regarding the mechanisms behind charmonium production in heavy-ion collisions. Inclusive

Figure 1 :
Figure 1: Upper panels: invariant-mass distribution of opposite-sign lepton pairs from the same event (black points) and mixed events (red histograms) at midrapidity (left) and forward rapidity (right) in Pb-Pb collisions at √ s NN = 5.02 TeV.Lower panels: invariant-mass distribution after the background subtraction with the eventmixing technique.The fit curves, shown in red, represent the sum of the signal and background shapes and the blue curves correspond to the residual background.

Figure 2 :
Figure 2: J/ψ ⟨p T ⟩ extraction in Pb-Pb collisions at √ s NN = 5.02 TeV at midrapidity for the 0-5% (left panel) and the 70-90% (right panel) centrality interval.The data points correspond to opposite-sign e + e − pairs from the same event, the blue line to the e + e − pairs from mixed events, and the red line is the combined fit that includes the mixed events and residual background which is described by the polynomial function.

Figure 3 :
Figure 3: J/ψ p T -differential production yields in Pb-Pb collisions at √ s NN = 5.02 TeV at midrapidity in the 0-10% (left panel) and 30-50% (right panel) centrality intervals.The statistical and systematic uncertainties are indicated, respectively, by the vertical error bars and the open boxes.The horizontal bars indicate the p T intervals.Data are compared to model calculations from Refs.[32,33,70].The ratios between data and models are shown in the lower panels.The filled boxes around unity depict the quadratic sum of statistical and systematic uncertainties from the measurement, while the bands indicate model uncertainties.

Figure 4 :
Figure 4: J/ψ p T -differential production yields in Pb-Pb collisions at √ s NN = 5.02 TeV at forward rapidity in the 0-20%, 20-40%, and 40-90% centrality intervals.The statistical and systematic uncertainties are indicated, respectively, by the vertical error bars and the open boxes.The horizontal bars indicate the p T intervals.Data are compared to model calculations from Refs.[32,33,70].The ratio between data and models is shown in the lower panels.The filled boxes around unity depict the quadratic sum of statistical and systematic uncertainties from the measurement, while the bands indicate model uncertainties.

Figure 5 :
Figure 5: Inclusive J/ψ nuclear modification factor at midrapidity and forward rapidity [19], integrated over p T , as a function of the number of participants in Pb-Pb collisions at √ s NN = 5.02 TeV.The statistical and systematic uncertainties are indicated, respectively, by the vertical error bars and the open boxes around the data points.The filled boxes around unity show the global uncertainties.

40 Figure 6 :
Figure 6: Inclusive J/ψ p T -differential R AA in Pb-Pb collisions at √ s NN = 5.02 TeV in various centrality intervals.The left panel shows the comparison of R AA measured in central (0-10%) and semicentral (30-50%) collisions at midrapidity.For the data point in the p T bin 10 < p T < 15 GeV/c an open symbol is used to highlight the usage of pp reference from the extrapolation approach.The right panel shows the measured R AA in three centrality classes, 0-20%, 20-40%, and 40-90%, at forward rapidity.The statistical and systematic uncertainties are indicated, respectively, by the vertical error bars and the open boxes around the data points.The filled boxes around unity show the global uncertainties.

Figure 7 :
Figure 7: Inclusive J/ψ R AA as a function of p T in Pb-Pb collisions at √ s NN = 5.02 TeV at midrapidity and forward rapidity, in the 0-10% centrality class and 0-20% centrality class, respectively.For the data point in the p T bin 10 < p T < 15 GeV/c an open symbol is used to highlight the usage of pp reference from the extrapolation approach.The statistical and systematic uncertainties are indicated, respectively, by the vertical error bars and the open boxes around the data points.The filled boxes around unity show the global uncertainties.

Figure 8 :Figure 9 :
Figure 8: Inclusive J/ψ R AA at midrapidity, integrated over p T , as a function of N part in Pb-Pb collisions at √ s NN = 5.02 TeV and compared to model calculations from Refs.[32, 33, 70].

Figure 11 :Figure 12 :
Figure 11: Left panel: Inclusive J/ψ ⟨p T ⟩ as a function of the mean number of participants in Pb-Pb collisions at √ s NN = 5.02 TeV at midrapidity.Right panel: Inclusive J/ψ r AA as a function of centrality at √ s NN = 5.02 TeV and compared with measurements at lower energies from RHIC [11, 88, 91] and the SPS [89].The statistical and systematic uncertainties are indicated, respectively, by the vertical error bars and the open boxes around the data points.The filled box around unity on the right panel shows the global uncertainty.

Table 1 :
Average nuclear overlap function ⟨T AA ⟩ and the average number of participants ⟨N part ⟩ in Pb-Pb collisions at √ s NN = 5.02 TeV for the centrality classes used in the analyses at midrapidity (upper) and forward rapidity (lower).

Table 2 :
Systematic uncertainties on the p T -integrated (0.15 < p T < 15 GeV/c) measurement at midrapidity for different centrality intervals.The individual contributions and the total uncertainties are given in percentage.It is considered that all the uncertainties are correlated over centrality and p T to a very large extent.

Table 3 :
Systematic uncertainties on the p T -differential measurement at forward rapidity for various centrality intervals.The individual contributions are given in percentage.When a range is given, it corresponds to the minimum and maximum values obtained in the p T interval.Values marked with an asterisk correspond to the uncertainties correlated over p T .These uncertainties are considered as global ones for R AA and added in quadrature to the uncorrelated uncertainties for the yields.