1 Introduction

Quarkonium production in heavy ion collisions has a rich history. In their original article [1], Matsui and Satz proposed that Debye color screening of the heavy-quark potential in a hot medium prevents the production of \({\mathrm {J}/\psi }\) mesons (and this applies also to other heavy-quark bound states such as \(\psi \text {(2S)}\), and \(\mathrm {\varUpsilon (1S)}\) mesons [2]). Consequently, the suppression of quarkonium yields in heavy ion collisions, relative to those in \(\mathrm {p}\mathrm {p}\) collisions, has long been considered to be a sensitive probe of deconfinement and quark-gluon plasma formation. The \({\mathrm {J}/\psi }\) meson suppression observed in \(\text {PbPb}\) collisions at the CERN SPS [3] and \(\text {AuAu}\) collisions at the BNL RHIC [4] is compatible with this picture. Similarly, the disappearance of \(\mathrm {\varUpsilon }\) resonances in \(\text {PbPb}\) collisions at the CERN LHC [5, 6] is consistent with the Debye screening scenario.

When produced abundantly in a single heavy ion collision, uncorrelated heavy quarks may combine to form quarkonia states in the medium [7, 8]. This additional source of quarkonium, commonly referred to as recombination, would enhance its production in heavy ion collisions, in contradistinction with the Debye screening scenario. Signs of this effect can be seen in the recent results from the ALICE Collaboration at the LHC [9, 10], which measured a weaker \({\mathrm {J}/\psi }\) meson suppression than at RHIC [4, 11], despite the higher medium energy density. Note that recombination is only expected to affect charmonium production at low transverse momenta (\(p_\text {T}\)), typically for values smaller than the charmonium mass (\(p_{\mathrm {T}} \lesssim m_{\psi } \, c\)), where the number of charm quarks initially produced in the collision is the largest [8].

At large \(p_{\mathrm {T}}\), other mechanisms may contribute to charmonium suppression. Until recently, no quarkonium results were available at high \(p_{\mathrm {T}} \), because of kinematic constraints at the SPS and too low counting rates at RHIC. At the LHC, a strong \({\mathrm {J}/\psi }\) suppression has been measured up to \(p_{\mathrm {T}} = 30\) \({\,\text {Ge}\text {V}/}\text {c}\) by the CMS Collaboration [12] in \(\text {PbPb}\) collisions at a centre-of-mass energy per nucleon pair of \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 2.76\,\text {Te}\text {V} \). Results at \(5.02\,\text {Te}\text {V} \) have also been reported, up to \(p_{\mathrm {T}} = 10\) \({\,\text {Ge}\text {V}/}\text {c}\), by the ALICE Collaboration [10]. According to Refs. [13, 14], quarkonium suppression by Debye screening may occur even at high \(p_{\mathrm {T}} \). At the same time, when \(p_{\mathrm {T}} \gg m_{\psi } \, c\), heavy quarkonium is likely to be produced by parton fragmentation, hence it should rather be sensitive to the parton energy loss in the quark-gluon plasma. The similarity of \({\mathrm {J}/\psi }\) meson suppression with the quenching of jets, light hadrons, and D mesons supports this picture [12, 15, 16].

At the LHC, the inclusive \({\mathrm {J}/\psi }\) meson yield also contains a significant nonprompt contribution coming from \(\mathrm {b}\) hadron decays [17,18,19]. The nonprompt \({\mathrm {J}/\psi }\) component should reflect medium effects on \(\mathrm {b}\) hadron production in heavy ion collisions, such as b quark energy loss. Measuring both prompt and nonprompt \({\mathrm {J}/\psi }\) meson production in \(\text {PbPb}\) collisions thus offers the opportunity to study both hidden charm and open beauty production in the same data sample.

In this paper we report on a new measurement of the prompt and nonprompt \({\mathrm {J}/\psi }\) and \(\psi \text {(2S)}\) nuclear modification factors (\(R_{\text {AA}}\)) using \(\text {PbPb}\) data, collected at the end of 2015 with the CMS experiment at \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 5.02\,\text {Te}\text {V} \). The analysis is performed via the dimuon decay channel. The results are compared to those obtained at 2.76\(\,\text {Te}\text {V}\)  [12]. The larger integrated luminosities allow for more precise and more differential measurements of \(R_{\text {AA}}\), as functions of centrality, rapidity (y), and \(p_{\mathrm {T}}\) up to 50\({\,\text {Ge}\text {V}/}\text {c}\).

2 The CMS detector

The central feature of the CMS apparatus is a superconducting solenoid of 6\(\text {\,m}\) internal diameter, providing a magnetic field of 3.8\(\text {\,T}\). Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Forward calorimeters extend the coverage provided by the barrel and endcap detectors. Muons are measured in the pseudorapidity range \(|\eta | < 2.4\) in gas-ionisation detectors embedded in the steel flux-return yoke outside the solenoid, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers. The hadron forward (HF) calorimeters use steel as an absorber and quartz fibres as the sensitive material. The two HF calorimeters are located 11.2\(\text {\,m}\) from the interaction region, one on each side, and together they provide coverage in the range \(2.9< |\eta | < 5.2\). They also serve as luminosity monitors. Two beam pick-up timing detectors are located at 175\(\text {\,m}\) on both sides of the interaction point, and provide information about the timing structure of the LHC beam. Events of interest are selected using a two-tiered trigger system [20]. The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimised for fast processing. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [21].

For \(\mathrm {p}\mathrm {p}\) data the vertices are reconstructed with a deterministic annealing vertex fitting algorithm using all of the fully reconstructed tracks [22]. The physics objects used to determine the primary vertex are defined based on a jet finding algorithm [23, 24] applied to all charged tracks associated with the vertex, plus the corresponding associated missing transverse momentum. The reconstructed vertex with the largest value of summed physics object \(p_{\mathrm {T}} ^2\) is taken to be the primary \(\mathrm {p}\mathrm {p}\) interaction vertex. In the case of \(\text {PbPb}\) data, a single primary vertex is reconstructed using a gap clustering algorithm [22], using pixel tracks only.

3 Data selection

3.1 Event selection

Hadronic collisions are selected offline using information from the HF calorimeters. In order to select \(\text {PbPb}\) collisions, at least three towers with energy deposits above 3 GeV are required in each of the HF calorimeters, both at forward and backward rapidities. A primary vertex reconstructed with at least two tracks is also required. In addition, a filter on the compatibility of the silicon pixel cluster width and the vertex position is applied [25]. The combined efficiency for this event selection, including the remaining non-hadronic contamination, is \((99 \pm 2)\)%. Values higher than 100% are possible, reflecting the possible presence of ultra-peripheral (i.e. non-hadronic) collisions in the selected event sample.

The \(\text {PbPb}\) sample is divided into bins of collision centrality, which is a measure of the degree of overlap of the colliding nuclei and is related to the number of participating nucleons (\(N_{\text {part}}\)). Centrality is defined as the percentile of the inelastic hadronic cross section corresponding to a HF energy deposit above a certain threshold [26]. The most central (highest HF energy deposit) and most peripheral (lowest HF energy deposit) centrality bins used in the analysis are 0–5% and 70–100% respectively. Variables related to the centrality, such as \(N_{\text {part}}\) and the nuclear overlap function (\(T_{\text {AA}}\)) [27], are estimated using a Glauber model simulation described in Ref. [28].

The \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) data sets correspond to integrated luminosities of 28.0\(\,\text {pb}^\text {-1}\) and 464\(\,\mu \mathrm {b}^{-1}\), respectively. Both \({\mathrm {J}/\psi }\) and \(\psi \text {(2S)}\) mesons are reconstructed using their dimuon decay channel. The dimuon events were selected online by the L1 trigger system, requiring two tracks in the muon detectors with no explicit momentum threshold, in coincidence with a bunch crossing identified by beam pick-up timing detectors. No additional selection was applied by the HLT. Because of the high rate of the most central dimuon events, a prescale was applied at the HLT level during part of the PbPb data taking: as a consequence only 79% of all the dimuon events were recorded, resulting in an effective luminosity of 368\(\,\mu \mathrm {b}^{-1}\). For peripheral events we were able to sample the entire integrated luminosity of 464\(\,\mu \mathrm {b}^{-1}\). This was done by adding an additional requirement that events be in the centrality range of 30–100% to the dimuon trigger. The prescaled data sample is used for the results integrated over centrality and those in the centrality range 0–30%, while for the results in the 30–100% range the data sample with 464\(\,\mu \mathrm {b}^{-1}\) was used instead. The results reported in this paper are unaffected by the small number of extra collisions potentially present in the collected events: the mean of the Poisson distribution of the number of collisions per bunch crossing (pileup), averaged over the full data sample, is approximately 0.9 for the \(\mathrm {p}\mathrm {p}\) data and less than 0.01 for \(\text {PbPb}\) collisions.

Simulated events are used to tune the muon selection criteria and the signal fitting parameters, as well as for acceptance and efficiency studies. These samples, produced using pythia 8.212 [29], and decaying the \(\mathrm {b}\) hadrons with evtgen 1.3.0 [30], are embedded in a realistic \(\text {PbPb}\) background event generated with hydjet 1.9 [31] and propagated through the CMS detector with Geant4  [32]. The prompt \({\mathrm {J}/\psi }\) is simulated unpolarised, a scenario in good agreement with pp measurements [33,34,35]. For nonprompt \({\mathrm {J}/\psi }\), the polarisation is the one predicted by evtgen, roughly \(\lambda _{\theta } = 0.4\). The resulting events are processed through the trigger emulation and the event reconstruction sequences. The assumptions made on the quarkonium polarisation affect the computation of the acceptance. Quantitative estimates of the possible effect evaluated for several polarisation scenarios can be found in Refs. [36, 37]. While there are no measurements on quarkonium polarisations in \(\text {PbPb}\) collisions, a study in \(\mathrm {p}\mathrm {p}\) collisions as a function of the event activity [38] has not revealed significant changes. Therefore the effects of the \({\mathrm {J}/\psi }\) polarisation on the acceptance are not considered as systematic uncertainties.

3.2 Muon selection

The muon reconstruction algorithm starts by finding tracks in the muon detectors, which are then fitted together with tracks reconstructed in the silicon tracker. Kinematic selections are imposed to single muons so that their combined trigger, reconstruction and identification efficiency stays above 10%. These selections are: \(p_{\mathrm {T}} ^{\mu }>3.50{\,\text {Ge}\text {V}/}\text {c} \) for \(|\eta ^{\mu } |<1.2\) and \(p_{\mathrm {T}} ^{\mu }>1.89{\,\text {Ge}\text {V}/}\text {c} \) for \(2.1<|\eta ^{\mu } |<2.4\), linearly interpolated in the intermediate \(|\eta ^{\mu } |\) region. The muons are required to match the ones selected by the dimuon trigger, and soft muon selection criteria are applied to global muons (i.e. muons reconstructed using the combined information of the tracker and muon detectors), as defined in Ref. [39]. Matching muons to tracks measured in the silicon tracker results in a relative \(p_{\mathrm {T}}\) resolution for muons between 1 and 2% for a typical muon in this analysis [39]. In order to remove cosmic and in-flight decay muons, the transverse and longitudinal distances of approach to the measured vertex of the muons entering in the analysis are required to be less than 0.3 and 20 cm, respectively. The probability that the two muon tracks originate from a common vertex is required to be larger than 1%, lowering the background from \(\mathrm {b}\) and \(\mathrm {c}\) hadron semileptonic decays.

4 Signal extraction

Because of the long lifetime of \(\mathrm {b}\) hadrons compared to that of \({\mathrm {J}/\psi }\) mesons, the separation of the prompt and nonprompt \({\mathrm {J}/\psi }\) components relies on the measurement of a secondary \(\mu ^+ \mu ^- \) vertex displaced from the primary collision vertex. The \({\mathrm {J}/\psi }\) mesons originating from the decay of \(\mathrm {b}\) hadrons can be resolved using the pseudo-proper decay length [40] \(\ell _{{\mathrm {J}/\psi }} = L_{xyz} \, m_{{\mathrm {J}/\psi }} \, c / |p_{\mu \mu }|\), where \(L_{xyz} \) is the distance between the primary and dimuon vertices, \(m_{{\mathrm {J}/\psi }}\) is the Particle Data Group [41] world average value of the \({\mathrm {J}/\psi }\) meson mass (assumed for all dimuon candidates), and \(p_{\mu \mu }\) is the dimuon momentum. Note that due to resolution effects and background dimuons the pseudo-proper decay length can take negative values. To measure the fraction of \({\mathrm {J}/\psi }\) mesons coming from \(\mathrm {b}\) hadron decays (the so-called nonprompt fraction), the invariant mass spectrum of \( \mu ^+ \mu ^- \) pairs and their \(\ell _{{\mathrm {J}/\psi }}\) distribution are fitted using a two-dimensional (2D) extended unbinned maximum-likelihood fit. In order to obtain the parameters of the different components of the 2D probability density function (PDF), the invariant mass and the \(\ell _{{\mathrm {J}/\psi }}\) distributions are fitted sequentially prior to the final 2D fits, as explained below. These fits are performed for each \(p_{\mathrm {T}}\), rapidity and centrality bin of the analysis, and separately in \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) collisions.

The sum of two Crystal Ball functions [42], with different widths but common mean and tail parameters, is used to extract the nominal yield values from the \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) invariant mass distributions. The tail parameters, as well as the ratio of widths in the \(\text {PbPb}\) case, are fixed to the values obtained from simulation. The background is described by a polynomial function of order N, where N is the lowest value that provides a good description of the data, and is determined by performing a log-likelihood ratio test between polynomials of different orders, in each analysis bin, while keeping the tail and width ratio parameters fixed. The order of the polynomial is chosen in such a way that increasing the order does not significantly improve the quality of the fit. The typical order of the polynomial is 1 for most of the analysis bins. The invariant mass signal and background parameters are obtained in an initial fit of the invariant mass distribution only and then fixed on the 2D fits of mass and \(\ell _{{\mathrm {J}/\psi }}\) distributions, while the number of extracted \({\mathrm {J}/\psi }\) mesons and background dimuons are left as free parameters.

The prompt, nonprompt, and background components of the \(\ell _{{\mathrm {J}/\psi }}\) distributions are parameterised using collision data and Monte Carlo (MC) simulated events, and the signal and background contributions unfolded with the \({}_{s}\mathcal {P}lot\) technique [43]. In the context of this analysis, this technique uses the invariant mass signal and background PDFs to discriminate signal from background in the \(\ell _{{\mathrm {J}/\psi }}\) distribution. The \(\ell _{{\mathrm {J}/\psi }}\) per-event uncertainty distributions of signal and background, provided by the reconstruction algorithm of primary and secondary vertices, are extracted from data and used as templates. The \(\ell _{{\mathrm {J}/\psi }}\) resolution is also obtained from the data by fitting the distribution of events with \(\ell _{{\mathrm {J}/\psi }} < 0\) with a combination of three Gaussian functions. The resolution varies event-by-event, so the per-event uncertainty is used as the width of the Gaussian function that describes the core. To take into account the difference on the per-event uncertainty distributions of signal and background dimuons, the resolution PDF is multiplied by the per-event uncertainty distribution of signal and background dimuons separately. All the resolution parameters are fixed in the 2D fits. The \(\mathrm {b}\) hadron decay length is allowed to float freely in the fit, and it is initialised to the value extracted by fitting the \(\ell _{{\mathrm {J}/\psi }}\) distribution of nonprompt \({\mathrm {J}/\psi }\) mesons from a MC sample with an exponential decay function, at generator level. The \(\ell _{{\mathrm {J}/\psi }}\) distribution of background dimuons is obtained from fits to the data, using an empirical combination of exponential functions. The parameters of the \(\ell _{{\mathrm {J}/\psi }}\) background distribution are also fixed in the 2D fits. Finally, the number of extracted \({\mathrm {J}/\psi }\) mesons, the number of background dimuons and the nonprompt fraction are extracted from the 2D fits. An example of a 2D fit of the invariant mass and pseudo-proper decay length for the \(\text {PbPb}\) data is shown in Fig. 1 for a representative analysis bin.

Fig. 1
figure 1

Invariant mass spectrum of \(\mu ^+ \mu ^- \) pairs (upper) and pseudo-proper decay length distribution (lower) in \(\text {PbPb}\) collisions for \(1.8<|y |<2.4\), \(4.5<p_{\mathrm {T}} <5.5{\,\text {Ge}\text {V}/}\text {c} \), for all centralities. The result of the fit described in the text is also shown

5 Acceptance and efficiency corrections

Correction factors are applied to all results to account for detector acceptance, trigger, reconstruction, and selection efficiencies of the \(\mu ^+ \mu ^- \) pairs. The corrections are derived from prompt and nonprompt \({\mathrm {J}/\psi }\) meson MC samples in \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\), and are evaluated in the same bins of \(p_{\mathrm {T}}\), centrality, and rapidity used in the \(R_{\text {AA}}\) and cross section analyses. The prompt and nonprompt \({\mathrm {J}/\psi }\) meson \(p_{\mathrm {T}}\) distributions in bins of rapidity in MC samples are compared to those in data, and the ratios of data over MC are used to weight the MC \({\mathrm {J}/\psi }\) distributions to describe the data better. This weighting accounts for possible mis-modelling of \({\mathrm {J}/\psi }\) kinematics in MC. The acceptance in a given analysis bin is defined as the fraction of generated \({\mathrm {J}/\psi }\) mesons in that bin which decay into two muons entering the kinematic limits defined above, and reflects the geometrical coverage of the CMS detector. The value of the acceptance correction ranges from 4 to 70%, depending on the dimuon \(p_{\mathrm {T}}\), both for prompt and nonprompt \({\mathrm {J}/\psi }\) mesons in \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) collisions. The efficiency in a given analysis bin is defined as the ratio of the number of reconstructed \({\mathrm {J}/\psi }\) mesons in which both muons pass the analysis selection and the number of generated \({\mathrm {J}/\psi }\) mesons in which both muons pass the analysis selection. The efficiency correction depends on the dimuon \(p_{\mathrm {T}}\), rapidity and event centrality, and ranges from 20 to 75% (15 to 75%) for prompt (nonprompt) \({\mathrm {J}/\psi }\) mesons in \(\text {PbPb}\) data, and from 40 to 85% for both prompt and nonprompt \({\mathrm {J}/\psi }\) mesons in \(\mathrm {p}\mathrm {p}\) data. The efficiency is lower at low than at high \(p_{\mathrm {T}}\), and it decreases from mid to forward rapidity; it is also lower for central than peripheral events. The individual components of the efficiency (tracking reconstruction, standalone muon reconstruction, global muon fit, muon identification and selection, and triggering) are also measured using single muons from \({\mathrm {J}/\psi }\) meson decays in both simulated and collision data, using the tag-and-probe (T&P) technique [36, 44]. The values obtained from data and simulation are seen to differ only for the muon trigger efficiency and the ratio of the data over simulated efficiencies is used as a correction factor for the efficiency. The correction factor for dimuons is at most 1.35 (1.38) for the \(\mathrm {p}\mathrm {p}\) (\(\text {PbPb}\)) efficiency in the \(3< p_{\mathrm {T}} < 4.5\) \({\,\text {Ge}\text {V}/}\text {c}\) and forward rapidity bin, but the \(p_{\mathrm {T}}\) and rapidity integrated value of the correction is about 1.03. The other T&P efficiency components are compatible, hence only used as a cross-check, as well as to estimate systematic uncertainties.

6 Systematic uncertainties

The systematic uncertainties in these measurements arise from the invariant mass signal and background fitting model assumptions, the parameterisation of the \(\ell _{{\mathrm {J}/\psi }}\) distribution, the acceptance and efficiency computation, and sample normalisation (integrated luminosity in \(\mathrm {p}\mathrm {p}\) data, counting of the equivalent number of minimum bias events in \(\text {PbPb}\), and nuclear overlap function). These systematic uncertainties are derived separately for \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) results, and the total systematic uncertainty is computed as the quadratic sum of the partial terms.

The systematic uncertainty due to each component of the 2D fits is estimated from the difference between the nominal value and the result obtained with the variations of the different components mentioned below, in the extracted number of prompt and nonprompt \({\mathrm {J}/\psi }\) mesons, or nonprompt fraction separately. In the following, the typical uncertainty is given for the observable on which each source has the biggest impact.

In order to determine the uncertainty associated with the invariant mass fitting procedure, the signal and background PDFs are independently varied, in each analysis bin. For the uncertainty in the signal, the parameters that were fixed in the nominal fits are left free with a certain constraint. The constraint for each parameter is determined from fits to the data, by leaving only one of the parameters free, and it is chosen as the root mean square of the variations over the different analysis bins. A different signal shape is also used: a Crystal Ball function plus a Gaussian function, with the CB tail parameters, as well as the ratio of widths in the PbPb case, again fixed from MC. The dominant uncertainty comes from the variation of the signal shape, yielding values for the number of extracted nonprompt \({\mathrm {J}/\psi }\) mesons ranging from 0.1 to 2.9% (0.3–5.5%) in \(\mathrm {p}\mathrm {p}\) (\(\text {PbPb}\)) data. For the background model, the following changes are considered, while keeping the nominal signal shape. First, the log-likelihood ratio tests are done again with two variations of the threshold used to choose the order of the polynomial function in each analysis bin. Also the fitted mass range is varied. Finally, an exponential of a polynomial function is also used. The dominant uncertainty in the background model arises from the assumed shape (invariant mass range) in \(\mathrm {p}\mathrm {p}\) (\(\text {PbPb}\)) data. The corresponding uncertainty ranges from 0.1 to 2.1% (0.1–2.8%). The maximum difference of each of these variations, in each analysis bin and separately for the signal and the background, is taken as an independent systematic uncertainty.

Fig. 2
figure 2

Fraction of \({\mathrm {J}/\psi }\) mesons coming from the decay of b hadrons, i.e. nonprompt \({\mathrm {J}/\psi }\) meson fraction, as a function of dimuon \(p_{\mathrm {T}}\) (upper) and rapidity (lower) for \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) collisions, for all centralities. The bars (boxes) represent statistical (systematic) point-by-point uncertainties

Fig. 3
figure 3

Differential cross section of prompt \({\mathrm {J}/\psi }\) mesons (left) and \({\mathrm {J}/\psi }\) mesons from b hadrons (nonprompt \({\mathrm {J}/\psi }\)) (right) decaying into two muons as a function of dimuon \(p_{\mathrm {T}}\) (upper) and rapidity (lower) in \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) collisions. The \(\text {PbPb}\) cross sections are normalised by \(T_{\text {AA}}\) for direct comparison. The bars (boxes) represent statistical (systematic) point-by-point uncertainties, while global uncertainties are written on the plots

Fig. 4
figure 4

Nuclear modification factor of prompt \({\mathrm {J}/\psi }\) mesons as a function of dimuon rapidity (upper left), \(N_{\text {part}}\) (upper right) and dimuon \(p_{\mathrm {T}}\) (lower) at \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 5.02\) \(\,\text {Te}\text {V}\). For the results as a function of \(N_{\text {part}}\) the most central bin corresponds to 0–5%, and the most peripheral one to 70–100%. Results obtained at 2.76\(\,\text {Te}\text {V}\) are overlaid for comparison [12]. The bars (boxes) represent statistical (systematic) point-by-point uncertainties. The boxes plotted at \(R_{\text {AA}} =1\) indicate the size of the global relative uncertainties

For the \(\ell _{{\mathrm {J}/\psi }}\) distribution fitting procedure, four independent variations of the different components entering in the 2D fits are considered. For the \(\ell _{{\mathrm {J}/\psi }}\) uncertainty distribution, instead of using the distributions corresponding to signal and background, the total distribution is assumed. The contribution to the systematic uncertainty in the number of extracted nonprompt \({\mathrm {J}/\psi }\) mesons ranges from 0.3 to 2% (0.3–9.5%) in \(\mathrm {p}\mathrm {p}\) (\(\text {PbPb}\)) data. The \(\ell _{{\mathrm {J}/\psi }}\) resolution obtained from prompt \({\mathrm {J}/\psi }\) meson MC is used instead of that evaluated from data. The corresponding uncertainty in the nonprompt fraction ranges from 1 to 5% (1–11%) in \(\mathrm {p}\mathrm {p}\) (\(\text {PbPb}\)) data. A nonprompt \({\mathrm {J}/\psi }\) meson MC template replaces the exponential decay function for the \(\mathrm {b}\) hadron decay length. In this case, the contribution of this source to the systematic uncertainty in the nonprompt \({\mathrm {J}/\psi }\) yield ranges from 0.2 to 8% (0.2–20%) in \(\mathrm {p}\mathrm {p}\) (\(\text {PbPb}\)) data. A template of the \(\ell _{{\mathrm {J}/\psi }}\) distribution of background dimuons obtained from the data is used to describe the background, instead of the empirical combination of exponential functions. This variation has an impact on the nonprompt \({\mathrm {J}/\psi }\) yield ranging from 0.1 to 1.3% (0.2–22%) in \(\mathrm {p}\mathrm {p}\) (\(\text {PbPb}\)) data. Therefore the dominant sources of uncertainty in the \(\ell _{{\mathrm {J}/\psi }}\) fitting are the background parameterisation and the MC template for the nonprompt signal. They have an important impact on the nonprompt \({\mathrm {J}/\psi }\) meson yield, especially at the lowest \(p_{\mathrm {T}}\) reached in this analysis for the most central events in \(\text {PbPb}\) collisions. The reason for this is that the background dimuons largely dominate over the nonprompt \({\mathrm {J}/\psi }\) signal.

The uncertainties in the acceptance and efficiency determination are evaluated with MC studies considering a broad range of \(p_{\mathrm {T}}\) and angular spectra compatible with the pp and \(\text {PbPb}\) data within their uncertainties. These variations yield an uncertainty about 0.2% (<1.7%) in \(\mathrm {p}\mathrm {p}\) (\(\text {PbPb}\)) collisions, both for prompt and nonprompt \({\mathrm {J}/\psi }\) acceptance and efficiency. The statistical uncertainty of the weighting of the MC distributions, reflecting the impact of the limited knowledge on the kinematic distribution of \({\mathrm {J}/\psi }\) mesons on the acceptance and efficiency corrections, is used as systematic uncertainty. This uncertainty is at most 6% (11%) in \(\mathrm {p}\mathrm {p}\) (\(\text {PbPb}\)) collisions at the largest \(p_{\mathrm {T}}\) but it usually ranges from 1 to 3% in both collision systems. In addition, the systematic uncertainties in the T&P correction factors, arising from the limited data sample available and from the procedure itself, are taken into account, covering all parts of the muon efficiency: inner tracking and muon reconstruction, identification, and triggering. The dominant uncertainty in the T&P correction factors arises from muon reconstruction and ranges from 2 to 10% for both collision systems.

The global uncertainty in the pp luminosity measurement is 2.3% [45]. The number of minimum bias events corresponding to our dimuon sample in PbPb (\(N_{\mathrm {MB}}\)) comes from a simple event counting in the events selected by the Minimum Bias triggers, taking into account the trigger prescale. The corresponding uncertainty arises from the inefficiency of trigger and event selection, and is estimated to be 2%. Finally, the uncertainty in the \(T_{\text {AA}}\) is estimated by varying the Glauber model parameters within their uncertainty and taking into account the uncertainty on the trigger and event selection efficiency, and ranges from 3 to 16% from the most central to the most peripheral events used in this analysis.

7 Results

In this section, the results obtained for nonprompt \({\mathrm {J}/\psi }\) fractions, prompt and nonprompt \({\mathrm {J}/\psi }\) cross sections for each collision system, and nuclear modification factors \(R_{\text {AA}}\) are presented and discussed. In addition, a derivation of the \(\psi \text {(2S)}\) \(R_{\text {AA}}\) is also presented and discussed. For all results plotted versus \(p_{\mathrm {T}}\) or \(|y |\), the abscissae of the points correspond to the centre of the respective bin, and the horizontal error bars reflect the width of the bin. The lower \(p_{\mathrm {T}}\) thresholds in the different rapidity intervals reflect the detector acceptance. In the range \(1.8< |y | < 2.4\) \({\mathrm {J}/\psi }\) are measured down to 3\({\,\text {Ge}\text {V}/}\text {c}\), while for the bins with \(|y | < 1.8\) they are measured down to 6.5\({\,\text {Ge}\text {V}/}\text {c}\). When plotted as a function of centrality, the abscissae are the average \(N_{\text {part}}\) values for minimum bias events within each centrality bin. The weighted average \(N_{\text {part}}\) values (weighted for the number of binary nucleon-nucleon collisions) correspond in most cases to the average \(N_{\text {part}}\) values for minimum bias events, with the exception of the most peripheral bin (50–100%) where \(N_{\text {part}}\) changes from 22 to 43. The centrality binning used is 0–5–10–15–20–25–30–35–40–45–50–60–70–100% for the results in \(|y |<2.4\), and 0–10–20–30–40–50–100% for the results differential in rapidity.

7.1 Nonprompt \({\mathrm {J}/\psi }\) meson fractions

The nonprompt \({\mathrm {J}/\psi }\) meson fraction is defined as the proportion of measured \({\mathrm {J}/\psi }\) mesons coming from b hadron decays, corrected for acceptance and efficiency. It is presented in Fig. 2 for \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) collisions, as a function of \(p_{\mathrm {T}}\) and rapidity, in the full \(|y |<2.4\) and \(6.5<p_{\mathrm {T}} <50{\,\text {Ge}\text {V}/}\text {c} \) range. No significant rapidity dependence is observed, while there is a strong \(p_{\mathrm {T}}\) dependence, from about 20% at low \(p_{\mathrm {T}}\) to 60% at high \(p_{\mathrm {T}}\), reflecting the different \(p_{\mathrm {T}}\) distributions of prompt and nonprompt \({\mathrm {J}/\psi }\) mesons, which highlights the necessity of separating the two contributions.

Fig. 5
figure 5

Nuclear modification factor of prompt \({\mathrm {J}/\psi }\) meson as a function of dimuon \(p_{\mathrm {T}}\) (upper) and \(N_{\text {part}}\) (lower), in the mid- and most forward rapidity intervals. For the results as a function of \(N_{\text {part}}\) the most central bin corresponds to 0–10%, and the most peripheral one to 50–100%. The bars (boxes) represent statistical (systematic) point-by-point uncertainties. The boxes plotted at \(R_{\text {AA}} =1\) indicate the size of the global relative uncertainties

Fig. 6
figure 6

Nuclear modification factor of prompt \({\mathrm {J}/\psi }\) mesons. Upper: as a function of dimuon \(p_{\mathrm {T}}\) in three centrality bins. Lower: as a function of \(N_{\text {part}}\) at moderate and high \(p_{\mathrm {T}}\), in the forward \(1.8<|y |<2.4\) range. For the results as a function of \(N_{\text {part}}\) the most central bin corresponds to 0–10%, and the most peripheral one to 50–100%. The bars (boxes) represent statistical (systematic) point-by-point uncertainties. The boxes plotted at \(R_{\text {AA}} =1\) indicate the size of the global relative uncertainties

Fig. 7
figure 7

Nuclear modification factor of prompt \({\mathrm {J}/\psi }\) and \(\psi \text {(2S)}\) mesons as a function of \(N_{\text {part}}\) (left) and dimuon \(p_{\mathrm {T}}\) (right), at central (upper, starting at \(p_{\mathrm {T}} =6.5\) \({\,\text {Ge}\text {V}/}\text {c}\)) and forward (lower, starting at \(p_{\mathrm {T}} =3.0\) \({\,\text {Ge}\text {V}/}\text {c}\)) rapidity. The vertical arrows represent 95% confidence intervals in the bins where the double ratio measurement is consistent with 0 (see text). For the results as a function of \(N_{\text {part}}\) the most central bin corresponds to 0–10% (0–20%), and the most peripheral one to 50–100% (40–100%), for \(|y |<1.6\) (\(1.6<|y |<2.4\)). The bars (boxes) represent statistical (systematic) point-by-point uncertainties. The boxes plotted at \(R_{\text {AA}} = 1\) indicate the size of the global relative uncertainties

Fig. 8
figure 8

Nuclear modification factor of \({\mathrm {J}/\psi }\) mesons from b hadrons (nonprompt \({\mathrm {J}/\psi }\)) as a function of dimuon rapidity (upper left), \(N_{\text {part}}\) (upper right) and dimuon \(p_{\mathrm {T}}\) (lower) at \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 5.02\) \(\,\text {Te}\text {V}\). For the results as a function of \(N_{\text {part}}\) the most central bin corresponds to 0–5%, and the most peripheral one to 70–100%. Results obtained at 2.76\(\,\text {Te}\text {V}\) are overlaid for comparison [12]. The bars (boxes) represent statistical (systematic) point-by-point uncertainties. The boxes plotted at \(R_{\text {AA}} =1\) indicate the size of the global relative uncertainties

7.2 Prompt and nonprompt \({\mathrm {J}/\psi }\) meson cross sections in \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) collisions

The measurements of the prompt and nonprompt \({\mathrm {J}/\psi }\) cross sections can help to test the existing theoretical models of both quarkonium production and b hadron production. The cross sections are computed from the corrected yields,

$$\begin{aligned} \frac{\mathrm {d}^2 N}{\mathrm {d}p_{\mathrm {T}} \,\mathrm {d}{}y} = \frac{1}{\varDelta p_{\mathrm {T}} \, \varDelta y} \, \frac{N_{{\mathrm {J}/\psi }}}{\mathcal {A} \, \epsilon }, \end{aligned}$$
(1)

where \(N_{{\mathrm {J}/\psi }}\) is the number of prompt or nonprompt \({\mathrm {J}/\psi }\) mesons, \(\mathcal {A}\) is the acceptance, \(\epsilon \) is the efficiency, and \(\varDelta p_{\mathrm {T}} \) and \(\varDelta y\) are the \(p_{\mathrm {T}}\) and rapidity bin widths, respectively. To put the \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) data on a comparable scale, the corrected yields are normalised by the measured integrated luminosity for pp collisions (\(\sigma = N / \mathcal {L}\)), and by the product of the number of corresponding minimum bias events and the centrality-integrated nuclear overlap value for \(\text {PbPb}\) collisions (\(N/(N_{\mathrm {MB}}T_{\text {AA}})\)). Global uncertainties (common to all measurements) arise from these normalisation factors and account for the integrated luminosity uncertainty in \(\mathrm {p}\mathrm {p}\) collisions (±2.3%) and the \(N_{\text {MB}}\) and \(T_{\text {AA}}\) uncertainty for \(\text {PbPb}\) collisions \(\left( ^{+3.4\%}_{-3.9\%}\right) \), respectively.

The cross sections for the production of prompt and nonprompt \({\mathrm {J}/\psi }\) mesons that decay into two muons (\(\mathcal {B} \sigma \), where \(\mathcal {B}\) is the branching ratio of \({\mathrm {J}/\psi }\) to dimuons) are reported as a function of \(p_{\mathrm {T}}\) and rapidity in Fig. 3.

7.3 Prompt \({\mathrm {J}/\psi }\) meson nuclear modification factor

In order to compute the nuclear modification factor \(R_{\text {AA}}\) in a given bin of centrality (cent.), the above-mentioned \(\text {PbPb}\) and \(\mathrm {p}\mathrm {p}\) normalised cross sections are divided in the following way:

$$\begin{aligned} R_{\text {AA}}= & {} \frac{N^{\text {PbPb}}_{{\mathrm {J}/\psi }} (\text {cent.}) }{N^{{\mathrm {p}\mathrm {p}}}_{{\mathrm {J}/\psi }}} \times \frac{\mathcal {A}^{{\mathrm {p}\mathrm {p}}} \times \epsilon ^{{\mathrm {p}\mathrm {p}}}}{\mathcal {A}^{\text {PbPb}} \, \epsilon ^{\text {PbPb}} (\text {cent.}) }\\&\times \frac{\mathcal {L}^{{\mathrm {p}\mathrm {p}}}}{N_{\text {MB}} \, \langle T_{\text {AA}} \rangle \, (\text {cent. fraction})}, \end{aligned}$$

where the centrality fraction is the fraction of the inclusive inelastic cross section probed in the analysis bin. Global uncertainties (indicated as boxes in the plots at \(R_{\text {AA}} =1\)) arise from the full \(\mathrm {p}\mathrm {p}\) statistical and systematic uncertainties and the \(\text {PbPb}\) \(N_{\text {MB}}\) uncertainty when binning as a function of the centrality; and from the integrated luminosity of the \(\mathrm {p}\mathrm {p}\) data, and the \(N_{\text {MB}}\) and \(T_{\text {AA}}\) uncertainties of the \(\text {PbPb}\) data, when binning as a function of rapidity or \(p_{\mathrm {T}}\).

In Fig. 4, the \(R_{\text {AA}}\) of prompt \({\mathrm {J}/\psi }\) mesons as a function of rapidity, \(N_{\text {part}}\) and \(p_{\mathrm {T}}\) are shown, integrating in each case over the other two non-plotted variables. The results are compared to those obtained at \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 2.76\) \(\,\text {Te}\text {V}\)  [12], and they are found to be in good overall agreement. No strong rapidity dependence of the suppression is observed. As a function of centrality, the \(R_{\text {AA}}\) is suppressed even for the most peripheral bin (70–100%), with the suppression slowly increasing with \(N_{\text {part}}\). The \(R_{\text {AA}}\) value for the most central events (0–5%) is measured for \(6.5<p_{\mathrm {T}} <50\) \({\,\text {Ge}\text {V}/}\text {c}\) and \(|y |<2.4\) to be \(0.219 \pm 0.005\,\text {(stat)} \pm 0.013\,\text {(syst)} \). As a function of \(p_{\mathrm {T}}\) the \(R_{\text {AA}}\) is approximately constant in the range of 5–20\({\,\text {Ge}\text {V}/}\text {c}\), but an indication of less suppression at higher \(p_{\mathrm {T}}\) is seen for the first time in quarkonia. Charged hadrons, for which the suppression is usually attributed to parton energy loss [16, 46], show a similar increase in \(R_{\text {AA}}\) at high \(p_{\mathrm {T}}\) for \(\text {PbPb}\) collisions at \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 5.02\) \(\,\text {Te}\text {V}\)  [27].

Double-differential studies are also performed. Figure 5 shows the \(p_{\mathrm {T}}\) (upper) and centrality (lower) dependence of prompt \({\mathrm {J}/\psi }\) \(R_{\text {AA}}\) measured in the mid- and most forward rapidity intervals. A similar suppression pattern is observed for both rapidities. Figure 6 (upper) shows the dependence of \(R_{\text {AA}}\) as a function of \(p_{\mathrm {T}}\), for three centrality intervals. Although the mean level of suppression strongly depends on the sampled centrality range, the general trend of the \(p_{\mathrm {T}}\) dependence appears similar in all three centrality ranges, including the increase of \(R_{\text {AA}}\) at high \(p_{\mathrm {T}}\). Finally, Fig. 6 (lower) considers the rapidity interval \(1.8<|y |<2.4\), where the acceptance goes down at lower \(p_{\mathrm {T}}\). The suppression is found to be similar in peripheral events at moderate (\(3<p_{\mathrm {T}} <6.5\) \({\,\text {Ge}\text {V}/}\text {c}\)) and high (\(6.5<p_{\mathrm {T}} <50\) \({\,\text {Ge}\text {V}/}\text {c}\)) transverse momentum ranges, but it is weaker for lower \(p_{\mathrm {T}}\) in the most central region. This is also reflected in the first bin of the most forward measurement in Fig. 5 (upper). A similarly reduced suppression at low \(p_{\mathrm {T}}\) is observed by the ALICE Collaboration, which is attributed to a regeneration contribution [9, 10].

7.4 Prompt \(\psi \text {(2S)}\) meson nuclear modification factor

Having measured the prompt \({\mathrm {J}/\psi }\) \(R_{\text {AA}}\), one can derive that of the \(\psi \text {(2S)}\) meson by multiplying it by the double ratio \((N_{\psi \text {(2S)}}/N_{{\mathrm {J}/\psi }})_{\text {PbPb}}/ (N_{\psi \text {(2S)}}/N_{{\mathrm {J}/\psi }})_{{\mathrm {p}\mathrm {p}}}\) of the relative modification of the prompt \(\psi \text {(2S)}\) and \({\mathrm {J}/\psi }\) meson yields from \(\mathrm {p}\mathrm {p}\) to \(\text {PbPb}\) collisions published in Ref. [47]. Since the \(\psi \text {(2S)}\) yield suffers from lower statistics, the current \({\mathrm {J}/\psi }\) analysis is repeated using the wider bins of Ref. [47]. The centrality binning used is 0–10–20–30–40–50–100% for the results in \(|y |<1.6\), and 0–20–40–100% for the results in \(1.6<|y |<2.4\). Since the statistical uncertainty in the \(\psi \text {(2S)}\) largely dominates, the \({\mathrm {J}/\psi }\) uncertainties are propagated by considering them to be uncorrelated to the double ratio uncertainties.

The results are presented in Fig. 7 as a function of dimuon \(p_{\mathrm {T}}\) and \(N_{\text {part}}\), in two rapidity ranges of different \(p_{\mathrm {T}}\) reach. In the bins where the double ratio is consistent with 0, 95% CL intervals on the prompt \(\psi \text {(2S)}\) \(R_{\text {AA}}\) are derived using the Feldman–Cousins procedure [48]. The procedure to obtain the CL intervals is the same as in the double ratio measurement, incorporating the \({\mathrm {J}/\psi }\) \(R_{\text {AA}}\) statistical and systematic uncertainties as a nuisance parameter. It can be observed that the \(\psi \text {(2S)}\) meson production is more suppressed than that of \({\mathrm {J}/\psi }\) mesons, in the entire measured range. The \(\psi \text {(2S)}\) meson \(R_{\text {AA}}\) shows no clear dependence of the suppression with \(p_{\mathrm {T}}\), and hints of an increasing suppression with collision centrality. These results show that the \(\psi \text {(2S)}\) mesons are more strongly affected by the medium created in \(\text {PbPb}\) collisions than the \({\mathrm {J}/\psi }\) mesons.

Fig. 9
figure 9

Nuclear modification factor of \({\mathrm {J}/\psi }\) mesons from b hadrons (nonprompt \({\mathrm {J}/\psi }\)) as a function of dimuon \(p_{\mathrm {T}}\) (upper) and \(N_{\text {part}}\) (lower) and in the mid- and most forward rapidity intervals. For the results as a function of \(N_{\text {part}}\) the most central bin corresponds to 0–10%, and the most peripheral one to 50–100%. The bars (boxes) represent statistical (systematic) point-by-point uncertainties. The boxes plotted at \(R_{\text {AA}} =1\) indicate the size of the global relative uncertainties

Fig. 10
figure 10

Nuclear modification factor of \({\mathrm {J}/\psi }\) mesons from b hadrons (nonprompt \({\mathrm {J}/\psi }\)). Upper: as a function of dimuon \(p_{\mathrm {T}}\) in three centrality bins. Lower: as a function of \(N_{\text {part}}\) at moderate and high \(p_{\mathrm {T}}\), in the forward \(1.8<|y |<2.4\) range. For the results as a function of \(N_{\text {part}}\) the most central bin corresponds to 0–10%, and the most peripheral one to 50–100%. The bars (boxes) represent statistical (systematic) point-by-point uncertainties. The boxes plotted at \(R_{\text {AA}} =1\) indicate the size of the global relative uncertainties

7.5 Nonprompt \({\mathrm {J}/\psi }\) meson nuclear modification factor

The procedure applied to derive the prompt \({\mathrm {J}/\psi }\) meson \(R_{\text {AA}}\) is applied to the nonprompt component. In Fig. 8, the \(R_{\text {AA}}\) of nonprompt \({\mathrm {J}/\psi }\) as a function of rapidity, centrality and \(p_{\mathrm {T}}\) are shown, integrating in each case over the other two non-plotted variables. The results are compared to those obtained at \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 2.76\) \(\,\text {Te}\text {V}\)  [12]. A good overall agreement is found, although no rapidity dependence is observed in the present analysis, while the suppression was slowly increasing towards forward rapidities in the lower-energy measurement. A steady increase of the suppression is observed with increasing centrality of the collision. The \(R_{\text {AA}}\) for the most central events (0–5%) measured for \(6.5<p_{\mathrm {T}} <50\) \({\,\text {Ge}\text {V}/}\text {c}\) and \(|y |<2.4\) is \(0.365\pm 0.009\,\text {(stat)} \pm 0.022\,\text {(syst)} \).

As for the prompt production case, double-differential studies are also performed. Figure 9 shows the \(p_{\mathrm {T}}\) (upper) and centrality (lower) dependence of nonprompt \({\mathrm {J}/\psi }\) meson \(R_{\text {AA}}\) measured in the mid- and most forward rapidity intervals. No strong rapidity dependence is observed, and a hint of a smaller suppression at low \(p_{\mathrm {T}}\) is seen in the \(1.8<|y |<2.4\) range. Figure 10 (upper) shows the dependence of \(R_{\text {AA}}\) as a function of \(p_{\mathrm {T}}\), for three centrality ranges. While the nonprompt \({\mathrm {J}/\psi }\) meson \(R_{\text {AA}}\) does not seem to depend on rapidity, the data indicates a larger \(p_{\mathrm {T}}\) dependence in peripheral events. Finally, Fig. 10 (lower) shows, for \(1.8<|y |<2.4\), \(R_{\text {AA}}\) as a function of \(N_{\text {part}}\), for two \(p_{\mathrm {T}}\) intervals. Hints of a stronger suppression are seen for \(p_{\mathrm {T}} >6.5{\,\text {Ge}\text {V}/}\text {c} \) at all centralities.

8 Conclusions

Prompt and nonprompt \({\mathrm {J}/\psi }\) meson production has been studied in \(\mathrm {p}\mathrm {p}\) and \(\text {PbPb}\) collisions at \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 5.02\,\text {Te}\text {V} \), as a function of rapidity, transverse momentum (\(p_{\mathrm {T}}\)), and collision centrality, in different kinematic and centrality ranges. Three observables were measured: nonprompt \({\mathrm {J}/\psi }\) fractions, prompt and nonprompt \({\mathrm {J}/\psi }\) cross sections for each collision system, and nuclear modification factors \(R_{\text {AA}}\). The \(R_{\text {AA}}\) results show a strong centrality dependence, with an increasing suppression for increasing centrality. For both prompt and nonprompt \({\mathrm {J}/\psi }\) mesons no significant dependence on rapidity is observed. An indication of less suppression in the lowest \(p_{\mathrm {T}}\) range at forward rapidity is seen for both \({\mathrm {J}/\psi }\) components. Double-differential measurements show the same trend, and also suggest a stronger \(p_{\mathrm {T}}\) dependence in peripheral events. An indication of less suppression of the prompt \({\mathrm {J}/\psi }\) meson at high \(p_{\mathrm {T}}\) is seen with respect to that observed at intermediate \(p_{\mathrm {T}}\). The measurements are consistent with previous results at \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 2.76\,\text {Te}\text {V} \).

Combined with previous results for the double ratio \((N_{\psi \text {(2S)}}/N_{{\mathrm {J}/\psi }})_{\text {PbPb}}/ (N_{\psi \text {(2S)}}/N_{{\mathrm {J}/\psi }})_{{\mathrm {p}\mathrm {p}}}\), the current \(R_{\text {AA}}\) values for \({\mathrm {J}/\psi }\) mesons are used to derive the prompt \(\psi \text {(2S)}\) meson \(R_{\text {AA}}\) in \(\text {PbPb}\) collisions at \(\sqrt{\smash [b]{s_{_{\text {NN}}}}} = 5.02\,\text {Te}\text {V} \), as a function of \(p_{\mathrm {T}}\) and collision centrality, in two different rapidity ranges. The results show that the \(\psi \text {(2S)}\) is more suppressed than the \({\mathrm {J}/\psi }\) meson for all the kinematical ranges studied. No \(p_{\mathrm {T}}\) dependence is observed within the current uncertainties. Hints of an increase in suppression with increasing collision centrality are also observed.