Centrality dependence of J/$\psi$ and $\psi$(2S) production and nuclear modification in p-Pb collisions at $\sqrt{s_{\rm NN}} =$ 8.16 TeV

The inclusive production of the J/$\psi$ and $\psi$(2S) charmonium states is studied as a function of centrality in p-Pb collisions at a centre-of-mass energy per nucleon pair $\sqrt{s_{\rm NN}} = 8.16$ TeV at the LHC. The measurement is performed in the dimuon decay channel with the ALICE apparatus in the centre-of-mass rapidity intervals $-4.46<y_{\rm cms}<-2.96$ (Pb-going direction) and $2.03<y_{\rm cms}<3.53$ (p-going direction), down to zero transverse momentum ($p_{\rm T}$). The J/$\psi$ and $\psi$(2S) production cross sections are evaluated as a function of the collision centrality, estimated through the energy deposited in the zero degree calorimeter located in the Pb-going direction. The $p_{\rm T}$-differential J/$\psi$ production cross section is measured at backward and forward rapidity for several centrality classes, together with the corresponding average $\langle p_{\rm T} \rangle$ and $\langle p^{2}_{\rm T} \rangle$ values. The nuclear effects affecting the production of both charmonium states are studied using the nuclear modification factor. In the p-going direction, a suppression of the production of both charmonium states is observed, which seems to increase from peripheral to central collisions. In the Pb-going direction, however, the centrality dependence is different for the two states: the nuclear modification factor of the J/$\psi$ increases from below unity in peripheral collisions to above unity in central collisions, while for the $\psi$(2S) it stays below or consistent with unity for all centralities with no significant centrality dependence. The results are compared with measurements in p-Pb collisions at $\sqrt{s_{\rm NN}} = 5.02$ TeV and no significant dependence on the energy of the collision is observed. Finally, the results are compared with theoretical models implementing various nuclear matter effects.


Introduction
Quarkonia, bound states of a heavy quark and its antiquark, are prominent probes of the properties of the strong interaction, which is described by quantum chromodynamics (QCD). In high-energy hadronic collisions, the production of quarkonia is usually factorised in a two-step process: the creation of a heavy-quark pair, mainly by gluon fusion at LHC energies, followed by its evolution and binding into a colour-singlet state. The former is described using perturbative QCD calculations, while the latter involves non-perturbative processes and is described using effective models [1][2][3].
In high-energy heavy-ion collisions, the creation of a deconfined state of nuclear matter made of quarks and gluons, the so-called quark-gluon plasma (QGP), modifies the production rates of the various quarkonium states. On the one hand, the production of quarkonium states is expected to be suppressed by the large density of colour charges in the QGP [4], with the suppression increasing with decreasing binding energy of the resonance [5]. Such sequential suppression has been observed, most notably in the bottomonium (bb) sector in Pb-Pb collisions at the LHC by the CMS [6][7][8][9] and ALICE [10] collaborations. On the other hand, quarkonia could also be regenerated during the QGP phase [11] or at its late boundary [12] by recombination of deconfined heavy quarks. Strong indications supporting such a regeneration mechanism, which (partially) compensates the aforementioned suppression, have been reported by the ALICE Collaboration in the charmonium (cc) sector for the J/ψ in Pb-Pb collisions at the LHC [13][14][15][16][17].
However, to fully exploit those experimental results for the understanding of the inner-workings of the QGP, other nuclear effects, not related to the presence of the QGP, must be addressed. These are typically referred to as cold nuclear matter (CNM) effects, as opposed to those related to the hot medium, and include the effects described below. A significant contribution involves the nuclear modification of the parton distribution function (PDF) of the nucleons inside the nucleus [18], i.e. the modification of the probability for a parton (quark or gluon) to carry a fraction x of the momentum of the nucleon. The gluon nuclear parton distribution function (nPDF) includes, most notably, a shadowing region at low x (x 0.01) corresponding to a suppression of gluons and an antishadowing region at intermediate x (0.01 x 0.3) corresponding to an enhancement of gluons [18]. The modification of the initial state of the nucleus with respect to an incoherent superposition of free nucleons can also be described in terms of the saturation of low-x gluons as implemented in the Colour Glass Condensate (CGC) effective field theory [19]. In addition, coherent energy-loss effects involving the initial-and final-state partons can modify the production of heavy-quark pairs and thus of quarkonium states [20]. The pre-resonant quarkonium state could also interact with the surrounding spectator nucleons. This nuclear absorption is expected to be negligible at LHC energies due to the short crossing time of the colliding nuclei [21]. The CNM effects discussed above are expected to affect similarly all states of the same quarkonium family, as they act on the production cross section of heavy-quark pairs or on the pre-resonant quarkonium state. On the contrary, final-state interactions with the co-moving medium [22] or with a medium including a short-lived QGP and a hadron resonance gas [23] could affect differently the various states of the same family. Soft-color exchanges between the hadronising cc pair and long-lived co-moving partons [24] could also affect differently the various charmonium states.
Cold nuclear matter effects are typically investigated using proton-nucleus collisions, where the formation of the QGP is not expected. At the LHC, the production of quarkonia was extensively studied in p-Pb collisions at a centre-of-mass energy per nucleon pair √ s NN = 5.02 TeV by the ALICE [25-31], ATLAS [32], CMS [33,34], and LHCb [35][36][37] collaborations. In the charmonium sector, a significant suppression of J/ψ yields is observed at forward rapidity y, i.e. in the p-going direction, at low transverse momentum p T , with the effect vanishing with increasing p T . The suppression at midrapidity is compatible with the one at forward y, while at backward y, i.e. in the Pb-going direction, no suppression of the J/ψ yields is observed [25,28,31]. Interestingly, the ψ(2S) appears to be more suppressed than the J/ψ at both forward and backward rapidity [26]. This observation cannot be explained by the first group of CNM effects discussed above and seems to indicate the need to consider additional final-state effects. The centrality dependence of the J/ψ and ψ(2S) suppression was also measured in p-Pb collisions at √ s NN = 5.02 TeV [29,30]. The difference between the ψ(2S) and J/ψ suppression increases with increasing centrality, especially at backward rapidity, indicating, once again, that shadowing or coherent parton energy-loss mechanisms are not enough to explain the ψ(2S) suppression [30]. Complementarily, the ALICE Collaboration also studied the J/ψ production at forward, mid, and backward rapidity as a function of the multiplicity of charged particles measured at midrapidity [38]. Such study does not require the interpretation of the centrality classes in terms of collision geometry and allows for the investigation of rare events with the highest charged particle multiplicities. An increase of the relative J/ψ yields with the relative charged-particle multiplicity is observed. At forward rapidity the increase saturates towards the highest multiplicities, while at backward rapidity a hint of a faster-than-linear increase with multiplicity is seen.
Also, in high-multiplicity p-Pb events, long-range angular correlations between the J/ψ at large rapidity and charged particles at midrapidity are observed [39]. These correlations are reminiscent of those observed in Pb-Pb collisions, which are often interpreted as signatures of the collective motion of the particles during the hydrodynamic evolution of the hot and dense medium.
More recently, the J/ψ and ψ(2S) production cross sections were also measured in p-Pb collisions at √ s NN = 8. 16 TeV as a function of transverse momentum and rapidity [40][41][42] confirming, with better statistical precision, the earlier findings. Namely, a significant suppression of the J/ψ is observed at forward rapidity but not at backward rapidity, and a stronger suppression of the ψ(2S) is seen, especially at backward rapidity. The J/ψ production at forward and backward rapidity as a function of the multiplicity of charged particles measured at midrapidity was also studied at √ s NN = 8.16 TeV [43] confirming the earlier observations.
In the bottomonium sector, a significant suppression of the ϒ(1S) yield is observed at mid and forward rapidity, vanishing from low to high transverse momentum, while at backward rapidity the yields are consistent with the expectations from pp collisions [27,32,36,44,45]. Interestingly, the excited ϒ(2S) state at midrapidity [32, 46] and ϒ(3S) state at backward rapidity [45] appear to be more suppressed than the fundamental ϒ(1S) state, which is similar to the comparison of the ψ(2S) and J/ψ discussed above. This paper presents the centrality dependence of the production of inclusive J/ψ and ψ(2S) in p-Pb collisions at √ s NN = 8.16 TeV. The inclusive ψ(nS) production contains contributions from direct ψ(nS), from decays of higher-mass excited states in the case of the J/ψ (mainly ψ(2S) and χ c ), as well as from non-prompt ψ(nS), from weak decays of beauty hadrons. Section 2 briefly presents the experimental setup and event selection, Section 3 describes the data analysis procedure, while the results are presented and discussed in Section 4. A summary is given in Section 5.

Experimental apparatus and event selection
A detailed description of the ALICE apparatus and its performance can be found in Refs. [47,48]. The main detectors used in this analysis are briefly discussed below.
The ALICE muon spectrometer is used to detect muons in the pseudorapidity region −4 < η lab < −2.5. It includes five tracking stations each having two planes of cathode pad chambers, with the third station being placed inside a dipole magnet with a field integral of 3 T · m. Two trigger stations, each composed of two planes of resistive plate chambers, provide the trigger for single muon as well as dimuon events with a programmable single-muon p T threshold. The setup is completed by a set of absorbers. A front absorber made of carbon, concrete, and steel is placed between the nominal interaction point (IP) and the first tracking station, to remove hadrons coming from the interaction vertex. An iron filter is positioned between the tracking and trigger stations and absorbs the remaining hadrons escaping the front absorber and the low p T muons originating from the decay of pions and kaons. Finally, a conical absorber surrounding the beam pipe protects the muon spectrometer against secondary particles produced by the primary particles emerging at large pseudorapidities and interacting with the beam pipe.
The Silicon Pixel Detector (SPD), corresponding to the two innermost layers of the Inner Tracking System [49] and covering the pseudorapidity ranges |η lab | < 2 (first layer) and |η lab | < 1.4 (second layer), is used to reconstruct the primary vertex of the collision. The two V0 hodoscopes [50] have 32 scintillator tiles each, are placed on each side of the IP, and cover the pseudorapidity ranges 2.8 < η lab < 5.1 and −3.7 < η lab < −1.7. The coincidence of signals from the two hodoscopes defines the minimum bias (MB) trigger condition and a first luminosity signal during van der Meer scans [51]. The V0s are also used to remove beam-induced background. A second luminosity signal in van der Meer scans is defined by the coincidence of signals from the two T0 arrays, which are located on opposite sides of the IP (4.6 < η lab < 4.9 and −3.3 < η lab < −3.0). Each array consists of 12 quartz Cherenkov counters read out by photomultiplier tubes [52]. Finally, two Zero Degree Calorimeters (ZDC) [53] are placed along the beam axis at ±112.5 m from the IP. Each ZDC is composed of a neutron calorimeter (ZN), positioned between the two beam pipes downstream of the first machine dipole that separates the beams, and a proton calorimeter (ZP), installed externally to the outgoing beam pipe. The ZN are used to estimate the centrality of the collision (described in Section 3) and to remove beam-induced background events.
The data were collected in 2016 with two beam configurations obtained by reverting the direction of the proton and lead ion beams. The corresponding acceptance ranges of the muon spectrometer, in terms of dimuon centre-of-mass rapidity, are −4.46 < y cms < −2.96 and 2.03 < y cms < 3.53. The backward and forward rapidity intervals correspond to the muon spectrometer being located in the Pb-going and p-going direction, and are denoted as Pb-p and p-Pb, respectively.
The non-symmetric rapidity ranges arise from the energy-per-nucleon asymmetry of the p and Pb beams, which shifts the rapidity of the nucleon-nucleon centre-of-mass system with respect to the laboratory system by 0.465 units of rapidity in the direction of the proton beam. The events were collected using an opposite-sign dimuon trigger, which requires the coincidence of the MB trigger condition and two opposite-sign track segments in the muon trigger chambers. For the data samples used here, the programmable online p T threshold for each muon track was set to 0.5 GeV/c. This threshold is not sharp in p T and the single-muon trigger efficiency is about 50% at p µ T = 0.5 GeV/c and reaches a plateau value of about 96% at p µ T 1.5 GeV/c. Beam-induced background was removed using the timing information provided by the V0 and the ZDC. The events are classified in classes of centrality according to the energy deposited in the ZN located in the direction of the Pb beam, as will be discussed in Section 3. Events in which two or more interactions occur in the same colliding bunch (in-bunch pile-up) or during the readout time of the SPD (out-of-bunch pile-up) are removed using the information from the SPD and V0. The integrated luminosity for the two beam configurations is L int = 12.8 ± 0.3 nb −1 for Pb-p and L int = 8.4 ± 0.2 nb −1 for p-Pb collisions.

Data analysis
In this section the various elements involved in the cross section and the nuclear modification factor measurements are discussed.
In p-Pb collisions, a centrality determination based on the charged-particle multiplicity can be biased by fluctuations related to the variation of the event topology, which are unrelated to the collision geometry. In contrast, an event selection depending on the energy deposited in the ZDC by nucleons emitted in the nuclear de-excitation process after the collision or knocked out by the nucleons participating in the collision (participant or wounded nucleons) should not be affected by this kind of bias. In this analysis the centrality estimation is based on a hybrid method, as described in detail in Refs. [54,55]. In this approach, the centrality classes are determined using the ZN detector, while the average number of binary nucleon-nucleon collisions N coll and the average nuclear overlap function T pPb for each centrality class are obtained assuming that the charged-particle multiplicity measured at midrapidity scales with the number of participant nucleons N part = N coll + 1. The centrality classes used in this analysis and the corresponding N coll and T pPb as well as their uncertainties, which reflect possible remaining biases (as discussed in Refs. [54,55]), are shown in Table 1. Monte Carlo (MC) simulations reproducing the LHC running conditions indicate that a residual pile-up may be present in the 2% most central collisions. The 0-2% centrality interval is therefore excluded and a 2% systematic uncertainty is conservatively assigned to the results in the other centrality classes. Furthermore, the 90-100% centrality interval is also excluded as the dimuon trigger may suffer from residual background contamination. It is worth noting that the previous analysis at √ s NN = 5.02 TeV was performed in the wider 80-100% centrality class where such possible contamination was not apparent. Charmonium candidates are built by forming pairs of opposite-sign charged tracks that were reconstructed by the tracking chambers of the muon spectrometer satisfying the following criteria. Each muon track candidate should be within −4 < η µ lab < −2.5 to avoid the edges of the acceptance. The tracks crossing the thicker part of the absorber are removed with the condition that the radial transverse position of the muon track at the end of the front absorber must be in the range 17.6 < R abs < 89.5 cm. The tracks must match a track segment in the muon trigger chambers above the aforementioned p T threshold of 0.5 GeV/c. The rapidity of the muon pair should be within the fiducial acceptance of the muon spectrometer, namely 2.03 < y cms < 3.53 and −4.46 < y cms < −2.96, for the p-Pb and Pb-p data samples, respectively.
The charmonium signal is estimated with a binned maximum likelihood fit to the dimuon invariant mass distribution. The J/ψ and ψ(2S) mass shapes are described with a Crystal Ball function with asymmetric tails on both sides of the peak (denoted as extended Crystal Ball) or a pseudo-Gaussian function [56]. The J/ψ mass and width are free parameters of the fit, while the other parameters, which correspond to the non-Gaussian tails of the signal shape, are fixed to those extracted from MC simulations. In addition, other sets of tails obtained from fits to the centrality-integrated invariant mass distribution in p-Pb at √ s NN = 8.16 TeV and in pp collisions at √ s = 8 TeV are used to test the stability of the fit and are included in the evaluation of the charmonium signal and its systematic uncertainty. The ψ(2S) fit parameters, apart from the amplitude, are constrained to those of the J/ψ, since its signal-to-background ratio is rather small. For the position of the mass peak, the following relation is used m ψ(2S) = m J/ψ + m PDG ψ(2S) − m PDG J/ψ , where the value obtained from the J/ψ fit is shifted by the difference between the two mass poles reported by the PDG [57]. The ψ(2S) width is fixed to the J/ψ one, applying a correction factor given by the ratio of the widths obtained in MC simulations (σ ψ(2S) = σ J/ψ × σ MC ψ(2S) σ MC J/ψ ). The background continuum is parameterised by either a Gaussian having a mass-dependent width or the product of a fourth degree polynomial function and an exponential. Finally, to test the background description, the signal is extracted using different fit ranges (2 < m µ µ < 5 GeV/c 2 and 2.2 < m µ µ < 4.5 GeV/c 2 ). The number of J/ψ and ψ(2S) and their statistical uncertainties are evaluated as the averages of the results of each test, i.e. the aforementioned signal extraction variations, and of their statistical uncertainty, respectively. The systematic uncertainty is given by the root-mean-square of the distribution of the results. For the ψ(2S), an additional uncertainty of 5% is added in quadrature. It corresponds to the uncertainty on the ψ(2S) width obtained from the large pp data sample used to validate the assumption on the relative widths for J/ψ and ψ(2S) from the MC [58]. In Fig. 1 the fits to the dimuon invariant mass distribution for the forward and the backward rapidity ranges are shown for two centrality classes.  Figure 1: Fit to the dimuon invariant mass distribution for the p-Pb (top panels) and Pb-p (bottom panels) data sets, for the 20-40% (left panels) and 60-80% (right panels) ZN centrality classes. The extended Crystal Ball function is used to describe the J/ψ and ψ(2S) signals, while a Variable Width Gaussian function is used for the background. The red line represents the total fit. realistically describe the J/ψ and ψ(2S) spectra, the MC input p T and y shapes are tuned directly on data performing an iterative procedure [40]. The decay products of the generated charmonia are then propagated inside a realistic description of the ALICE detector, based on GEANT 3.21 [64]. The p Tand y-integrated (A × ε) values are 0.264 ± 0.001 (0.235 ± 0.001) in the p-Pb (Pb-p) data sample for the J/ψ and 0.280 ± 0.008 (0.250 ± 0.004) for the ψ(2S). The larger (A × ε) in the p-Pb than Pb-p data taking period is due to different running conditions. The quoted uncertainties are the systematic uncertainties on the input p T and y shapes used for the MC generation, which are evaluated comparing the (A × ε) values obtained using different input MC distributions. For the J/ψ, these were obtained by adjusting the input MC distributions to the data in various p T and y intervals. For the ψ(2S), due to the larger statistical uncertainties of the data, the input p T and y shapes used for the J/ψ were considered in addition to the ones tuned directly on the ψ(2S) data. For the J/ψ the uncertainty on the p T -integrated (A × ε) is 0.5% for both p-Pb and Pb-p, varying with p T from 1% to 3%, while for the ψ(2S) the uncertainty amounts to 3% and 1.5% in p-Pb and Pb-p collisions, respectively. The same values of (A × ε) are used for all centrality classes since no dependence on the detector occupancy is observed within the multiplicities reached in p-Pb collisions. Possible changes of (A × ε) due to shape variations of the p T -and y-differential cross sections with centrality are accounted for in the systematic uncertainties by using different p T shapes extracted from different centrality intervals as inputs to the MC simulations. The corresponding systematic uncertainty on the p T -integrated J/ψ (A × ε) varies from 1.6% (2.5%) to 1.7% (2.7%) as a function of centrality, while as a function of p T in different centrality classes it varies from 1.2% (1.4%) to 4.4% (2.2%) in Pb-p (p-Pb) collisions.
The normalisation of the J/ψ and ψ(2S) yields is obtained following the prescription described in Ref. [54]. It is based on the evaluation of the number of minimum bias events N i MB for each centrality class i as N i is the inverse of the probability of having a dimuontriggered event in a MB-triggered one, and N i 2µ is the number of analyzed dimuon-triggered events. The value of F i 2µ/MB depends on the centrality class and increases from central to peripheral events, passing from 384 ± 3 to 1855 ± 18 in p-Pb and from 161 ± 1 to 2036 ± 16 in Pb-p collisions, for the 2-10% and 80-90% centrality classes, respectively. The quoted systematic uncertainties, which vary between 1% and 1.4%, contain two contributions. The first one, which is correlated in p T and centrality and amounts to 1%, is estimated by comparing the centrality-integrated F 2µ/MB obtained with the method described above with the one obtained using the information of the online trigger counters, as described in Ref.
[25]. The second one, which is not correlated in centrality, is obtained comparing F i 2µ/MB evaluated with two different methods, as detailed in Ref. [29]. Namely, F i 2µ/MB can be evaluated directly in each centrality class, or derived from the centrality integrated F 2µ/MB factor normalised by the ratio of N i MB /N MB to N i 2µ /N 2µ . The resulting systematic uncertainty varies from 0.1% (0.1%) to 0.8% (1%) in Pb-p (p-Pb) collisions.
The inclusive cross section for J/ψ and ψ(2S) for centrality class i is calculated using the expression where N i ψ(nS)→µ + µ − is the raw yield for the given resonance, (A × ε) ψ(nS)→µ + µ − is the corresponding product of the detector acceptance and reconstruction efficiency, and B.R. ψ(nS)→µ + µ − is the branching ratio of the corresponding dimuon decay channel as reported in Ref.
[57]. The integrated luminosity L int of the analyzed data sample is given by the ratio of the equivalent number of minimum bias events N MB to the cross section for events satisfying the minimum bias trigger condition σ MB . The latter is evaluated through a van der Meer scan and results in a value of 2.09 ± 0.04 b for p-Pb collisions and 2.10 ± 0.04 b for Pb-p [51], where the quoted uncertainties are the systematic uncertainties. The integrated luminosity can be independently calculated using the luminosity signal provided by the T0 detector. The difference between the integrated luminosity obtained with the V0 and T0 detectors amounts to 1.1% (0.6%) [51] in the p-Pb (Pb-p) data sample and is assigned as a further systematic uncertainty of σ MB . The correlated uncertainty on σ MB for the p-Pb and Pb-p data samples are 0.5% and 0.7%, respectively [51].
The relative modification between the two charmonium states in proton-nucleus collisions can be firstly observed through the evaluation of the ratio B.R. ψ(2S)→µ + µ − σ ψ(2S) /B.R. J/ψ→µ + µ − σ J/ψ , where the systematic uncertainties on trigger, tracking, and matching efficiencies, as well as on the luminosity, which are common for the J/ψ and ψ(2S), cancel out. The only remaining systematic uncertainties are those related to the signal extraction and to the shape of the input p T and y distribution used for the MC simulations. In turn this ratio can be normalised to the same quantity evaluated in pp collisions, providing a direct access to the relative ψ(2S) production modification with respect to J/ψ moving from a pp to a p-Pb collision system. Since there is no measurement available at √ s = 8.16 TeV in pp collisions, the ratio ψ(2S)/J/ψ is evaluated through an interpolation procedure using ALICE data at √ s = 5, 7, 8, and 13 TeV in the interval 2.5 < y < 4 [58, 65, 66]. The uncertainty associated to the interpolated value contains a contribution of 6% due to the energy-interpolation procedure and a further 1% contribution due to the rapidity-extrapolation procedure [42]. In addition an extra 1% is included due to the assumption of non-flat dependence of the ratio as a function of The nuclear modification factor as a function of centrality is calculated using the following expression where T i pPb is the nuclear overlap function for the centrality class i, while σ pp ψ(nS) is the ψ(nS) production cross section in proton-proton collisions. The notation Q pPb is used instead of the usual R pPb in order to point out the possible bias in the centrality determination, which depends on the loose correlation between the centrality estimator and the collision geometry [54]. The J/ψ cross section in pp collisions at √ s = 8.16 TeV is obtained from the available results in the interval 2.5 < y < 4 for inclusive J/ψ production at √ s = 8 TeV from ALICE [66] and LHCb [70] using the energy and rapidity extrapolation procedure described in Ref. [40]. A resulting first contribution of 7.1% to the systematic uncertainty of the extrapolation procedure is correlated in p T , y, and centrality. A second contribution of 1.8% (1.5%) for the p T -integrated cross section and ranging from 3.0% to 4.6% (2.9% to 4.7%) for the p Tdifferential cross section at backward (forward) rapidity, correlated with centrality, arises from the energy and rapidity interpolation procedures (see Ref. [40] for details). The ψ(2S) cross section in pp collisions at √ s = 8.16 TeV is obtained from the extrapolated J/ψ cross section and the interpolated ψ(2S)/J/ψ ratio. The related total systematic uncertainty is 9.4% and is correlated in p T , y, and centrality. The resulting extrapolated cross sections are reported in Ref. [69].
In addition to the various contributions to the systematic uncertainty discussed above, the following sources, which are common for the J/ψ and ψ(2S) states, are also taken into account. The systematic uncertainty of the trigger efficiency includes two contributions, one related to the intrinsic efficiency of each trigger chamber and one related to the muon trigger response function. The former is calculated from the uncertainties on the trigger chamber efficiencies measured from data and applied to simulations and it amounts to 1%. The latter is obtained from the difference between the (A × ε) obtained using the response function in data or in MC simulations and for the p T -integrated case this uncertainty is 2.9% for Pb-p and 2.4% for p-Pb, and it varies between 1% and 4% as a function of p T . The total systematic uncertainty of the trigger efficiency, obtained by adding in quadrature the aforementioned contributions, is 3.1% for Pb-p and 2.6% for p-Pb, varying as a function of p T from 1.4% up to 4.1%. The evaluation of the systematic uncertainty on the tracking efficiency follows a similar approach as reported in Ref. [26]. The discrepancy between the efficiencies in data and MC corresponds to 2% in Pb-p and 1% in p-Pb, without any appreciable dependence on the dimuon kinematics and event centrality. Finally, the choice of the χ 2 selection applied for the definition of the matching between tracks in the trigger and tracking chambers leads to a 1% systematic uncertainty.
In Table 2, a summary of all the sources of systematic uncertainty which contribute to the cross section and nuclear modification factor measurements is reported. Table 2: Summary of the systematic uncertainties (in percentage) of the quantities associated to the measurements of the differential J/ψ cross section and Q pPb of J/ψ and ψ(2S). The uncertainties for the p T -differential case are indicated in parentheses if the values are different from the p T -integrated case. When appropriate, a range of variation (for centrality, rapidity, or p T intervals) of the uncertainty is given. Type I, II, and III stands for uncertainties correlated over centrality, rapidity, or p T , respectively. 4 Results 4.1 p T -differential cross section of inclusive J/ψ for various centrality classes Figure 2 shows the p T -differential cross section of inclusive J/ψ at backward (left) and forward (right) rapidity measured in six centrality classes: 2-10%, 10-20%, 20-40%, 40-60%, 60-80%, and 80-90%. The vertical error bars represent the statistical uncertainties and the open boxes the uncorrelated systematic uncertainties. A global systematic uncertainty, which is correlated over centrality, rapidity, and p T and is obtained as the quadratic sum of the systematic uncertainty of the branching ratio and the correlated systematic uncertainty of σ MB amounts to 0.9% (0.7%) at backward (forward) rapidity.

Inclusive J/ψ average transverse momentum and p T broadening
A first insight into the modification of J/ψ production in p-Pb collisions can be obtained by studying the average transverse momentum p T and the average squared transverse momentum p 2 T as a function of the collision centrality. The p T and p 2 T are extracted for each centrality class by performing a fit of the p T -differential cross section with a widely used function proposed in Ref.
[71] and defined as where C, p 0 , and n are free parameters of the fit. The central values of p T and p 2 T are obtained from the fit using the quadratic sum of statistical and uncorrelated systematic uncertainties of the data points. The uncertainties on the free parameters obtained from the fit are propagated to the values of p T and ))   Figure 2: Inclusive J/ψ p T -differential cross section for different centrality classes at backward (left) and forward (right) rapidity in p-Pb collisions at √ s NN = 8.16 TeV. The vertical error bars, representing the statistical uncertainties, and the boxes around the points, representing the uncorrelated systematic uncertainties, are smaller than the marker. The global systematic uncertainty, which is correlated over centrality, rapidity, and p T and is obtained as the quadratic sum of the systematic uncertainty of the branching ratio and the correlated systematic uncertainty of σ MB , amounts to 0.9% (0.7%) at backward (forward) rapidity and is shown as text.
p 2 T . The statistical and systematic uncertainties on p T and p 2 T are obtained by performing the fit using, respectively, only the statistical or the uncorrelated systematic uncertainties on the data points. The range of integration on p T for this calculation is limited to the p T interval 0 < p T < 16 GeV/c. Extending the integration range to infinity has a negligible effect with respect to the quoted uncertainties. Table 3 shows the values of p T and p 2 T of inclusive J/ψ for each centrality class. Both p T and p 2 T increase with increasing centrality, which indicates a hardening of the J/ψ p T distribution from peripheral to central collisions in both rapidity intervals.
The p T broadening defined as the difference between the average squared transverse momentum in p-Pb and pp collisions (∆ p 2 T = p 2 T pPb − p 2 T pp ) can be used to quantify the nuclear effects on the J/ψ production [72-74]. The value of p 2 T pp is evaluated from the p T -differential cross section in pp collisions at √ s = 8.16 TeV obtained with the interpolation procedure described in Ref.
[40], and using the same p T integration range as for p-Pb collisions. The resulting values are reported in Ref. [69]. Figure 3 shows ∆ p 2 T as a function of the number of binary collisions at backward and forward rapidity. In all cases, ∆ p 2 T is larger than zero, indicating a broadening of the p T distribution in p-Pb collisions compared to  pp collisions. For the most peripheral collisions, corresponding to N coll ∼ 2.5, the ∆ p 2 T measured at backward y is compatible, within uncertainties, with that at forward y. In both backward and forward rapidity ranges the p T broadening increases with increasing centrality. However, the increase of ∆ p 2 T is stronger in the p-going direction than in the Pb-going direction. Thus, nuclear effects appear to increase with the centrality of the collision and to be stronger in the p-going than in the Pb-going direction. Here, it is worth noting that under the naive assumption of a 2 → 1 production process (gg → ψ(nS)), the sampled x ranges of the lead nuclei correspond to the shadowing and anti-shadowing regions for the p-going and lead-going direction measurements, respectively. Also shown in Fig. 3 are the results at √ s NN = 5.02 TeV [29]. The same trend of ∆ p 2 T as a function of N coll is seen at both collision energies in the two rapidity ranges. Overall, ∆ p 2 T slightly increases with the collision energy. The ∆ p 2 T as function of N coll is also compared in Fig. 3 to the results of an energy loss model, which is based on a parameterisation of the prompt J/ψ pp cross section and includes coherent energy loss effects from the incoming and outgoing partons [20]. The band in this model represents the uncertainty on the parton transport coefficient and the parameterisation used for the pp reference cross section. The model describes the centrality dependence of ∆ p 2 T at forward rapidity reasonably well, but it underestimates the data at backward rapidity. Figure 4 shows the p T -integrated Q pPb of J/ψ as a function of N coll in p-Pb collisions at √ s NN = 8.16 TeV at backward and forward rapidity. At forward y, the production of inclusive J/ψ in p-Pb collisions is suppressed with respect to expectations from pp collisions for all centrality classes. Furthermore, Q pPb decreases with increasing collision centrality from a value of 0.85 ± 0.02(stat.) ± 0.03(syst.) for the 80-90% centrality class to 0.69 ± 0.01(stat.) ± 0.04(syst.) for the 2-10% centrality class. At backward y, on the contrary, a significant suppression is seen for the most peripheral collisions (Q 80−90% pPb = 0.80 ± 0.02(stat.) ± 0.03(syst.)) with Q pPb increasing with increasing centrality and reaching values above unity for the most central collisions (Q 2−10% pPb = 1.16 ± 0.01(stat.) ± 0.07(syst.)). The Q pPb as a function of N coll is compared with the results at √ s NN = 5.02 TeV [29]. No strong dependence with the energy of the collision is observed in the two rapidity intervals.

Centrality dependence of the inclusive J/ψ nuclear modification factor
Three model calculations are also shown in Fig. 4 for comparison. First, a next-to-leading order (NLO) Colour Evaporation Model (CEM) [75] using the EPS09 parameterisation of the nuclear modification of the gluon PDF at NLO is shown and denoted as "EPS09s NLO + CEM". The band represents the systematic uncertainty of the calculation, which is dominated by the uncertainty of the EPS09 parameterisation.
The second one is the energy loss model that was described in the Section 4.2. Finally, the third one is a transport model [23] based on a thermal-rate equation framework, which implements the dissociation of charmonia in a hadron resonance gas. The fireball evolution implemented in this model includes the transition from a short QGP phase into the hadron resonance gas, through a mixed phase. The model uses a cc production cross section dσ cc /dy = 0.57 mb and a prompt J/ψ production cross section in pp collisions of dσ pp J/ψ /dy = 3.35 µb. Shadowing effects are included through the EPS09 parameterisation. In this case, the upper (lower) limit of this calculation corresponds to a 10% (25%) contribution of nuclear shadowing. The three models provide a satisfactory description of the centrality dependence of the inclusive J/ψ Q pPb at forward rapidity. However, at backward rapidity, all three calculations show a slightly decreasing trend of Q pPb with increasing centrality that appears opposite to the one indicated by the data.
It is worth noting that the model calculations discussed above are for prompt J/ψ while the inclusive measurements contain a contribution from non-prompt J/ψ too. The Q and Q incl pPb are found to be below 9% and 5% at backward and forward rapidity, respectively. Thus, the conclusions outlined above, and also in the following, are expected to remain valid also for prompt J/ψ.     4.4 Centrality-differential inclusive J/ψ Q pPb as a function of p T Figure 5 shows the inclusive J/ψ Q pPb as a function of p T at backward and forward rapidity for all centrality classes considered in this analysis. At backward rapidity, a slight suppression is seen at low p T for all centralities. However, while almost no p T dependence is observed for the most peripheral collisions, for all other centralities Q pPb increases with p T reaching a plateau for p T 5 GeV/c, with the value of the plateau being largest for more central collisions. For the three most central classes, Q pPb is above unity for p T 2 GeV/c. A similar behavior is also observed for prompt D mesons at midrapidity (−0.96 < y cms < 0.04) measured in p-Pb collisions at √ s NN = 5.02 TeV [76]. In contrast, at forward rapidity, Q pPb is below or consistent with unity for all p T in all centrality classes. At low p T , a centrality dependent hierarchy of Q pPb is observed, showing a stronger suppression in central collisions compared to peripheral ones. For all centralities, Q pPb smoothly increases towards unity at high p T .
The different shapes of the evolution of Q pPb with p T for the various centralities can be better appreciated by forming the ratio Q PC of the Q pPb in peripheral to that in central collisions. Figure 6 shows the inclusive J/ψ Q PC as a function of p T at backward and forward rapidity. The centrality-correlated systematic uncertainties cancel when calculating the ratio. The Q PC could, therefore, provide stronger constraints to the theoretical calculations. Transport model calculations by Du et al. [23] are also shown in Fig. 6 for comparison. At backward rapidity, the model calculations tend to overestimate the measured Q PC for all centrality classes. The centrality dependent hierarchy of the measured Q PC is also not reproduced by the model calculation. At forward rapidity, the transport model calculations qualitatively describe the p T and centrality dependence of the inclusive J/ψ Q PC , but do systematically overestimate the measurements.
The J/ψ Q pPb as a function of p T is shown separately for the six centrality classes in Figs. 7 and 8 for the backward and forward rapidity regions and is compared with the results at √ s NN = 5.02 TeV [29] and the same model calculations discussed previously. The results are similar at both collision energies in the two rapidity ranges, indicating that the mechanisms behind the modification of the J/ψ production in p-Pb collisions do not depend strongly on the collision energy. It is worth noting that the p T range is extended up to 16 GeV/c at √ s NN = 8.16 TeV and that the most peripheral centrality is 80-90% at the highest energy while it was 80-100% at the lowest one.
At backward rapidity, the EPS09s NLO + CEM [75] calculations show a mild increase of Q pPb with p T for all centralities, but more pronounced towards more central collisions. The EPS09s NLO + CEM Q pPb is above unity for all centralities but the strength of the anti-shadowing effect is stronger the more central the collisions are. The description of the data by the EPS09s NLO + CEM calculations is rather poor, except for the 40-60% centrality class. For more central collisions the calculations underestimate the data, but overestimate them for more peripheral collisions. Similar observations can be drawn from the energy loss [20] calculations, which in the common p T region are compatible with the EPS09s NLO + CEM calculations. Only for the more central collisions the p T dependence appears steeper for the energy loss model and closer to the data, but the overall magnitude is lower than the measured Q pPb . The transport model [23] calculations, which are in general terms lower than the EPS09 + CEM and quite similar to the energy loss ones, only describe the inclusive J/ψ Q pPb in the 40-60% centrality class, while underestimating it for more central collisions and overestimating it for more peripheral ones.
At forward rapidity, the differences between the EPS09 + CEM and the energy loss calculations are more pronounced. On the contrary, the transport model calculations are rather similar to the EPS09 + CEM ones, though on the lower edge. The uncertainties of the model calculations are also larger at forward than at backward rapidity, especially for the most central collisions. The description of the data by the EPS09 + CEM calculations is fair for all centralities, especially for p T 4 GeV/c. Below 4 GeV/c, the model tends to overestimate the measured Q pPb . The p T dependence of the energy loss calculation appears steeper than that in data, except for the most peripheral class. The model tends to underestimate the measured Q pPb at low p T and to overestimate it at high p T in all the other centrality classes. The transport model describes the data fairly well in all centrality classes for p T 8 GeV/c but tends to overestimate the Q pPb at higher p T .

Inclusive ψ(2S) to J/ψ ratio and double ratio
The relative production of the excited ψ(2S) state compared to that of the J/ψ state can be quantified by the ψ(2S)/J/ψ ratio, which is defined here as B.R. ψ(2S)→µ + µ − σ ψ(2S) /B.R. J/ψ→µ + µ − σ J/ψ . The relative modification of the production of the two states in p-Pb collisions with respect to pp collisions is then obtained by comparing the ψ(2S)/J/ψ ratio in the two collision systems. Several systematic uncertainties cancel in the ψ(2S)/J/ψ ratio, and the remaining ones are due to the signal extraction and the MC input shapes. The centrality dependence of the ψ(2S)/J/ψ ratio at backward and forward rapidity in p-Pb collisions at √ s NN = 8.16 TeV are shown in Fig. 9. The results are compared with the same ratio in p-Pb collisions at √ s NN = 5.02 TeV [30] as well as in pp collisions at √ s = 7 TeV [65]. Firstly, the ψ(2S)/J/ψ ratio does not exhibit any significant dependence on the collision energy. Secondly, the ratio appears to be smaller in p-Pb than in pp collisions, in both explored rapidity regions and for all centralities, except the most peripheral, where the uncertainty is considerably large, and the most central ones. Here, it is important to note that also no significant energy dependence is observed in the ψ(2S)/J/ψ ratio in pp collisions [58]. Thus, the production of the ψ(2S) in p-Pb collisions appears to be suppressed compared to that of the J/ψ with respect to the expectation from pp collisions. Thirdly, given the current experimental uncertainties, no clear trend of the ratio as a function of centrality can be drawn. Finally, the suppression of the ψ(2S) relative to the J/ψ in p-Pb compared to pp collisions appears to be stronger in the Pb-going (backward rapidity) than in the p-going direction (forward rapidity).
The same conclusions can be also drawn from the so-called double ratio, i.e. the ratio of the ψ(2S) to the J/ψ cross section in p-Pb collisions divided by the same ratio in pp collisions, Figure 10 shows the double ratio [σ ψ(2S) /σ J/ψ ] pPb /[σ ψ(2S) /σ J/ψ ] pp as a function of centrality in the backward and forward rapidity regions for p-Pb collisions at √ s NN = 8.16 TeV. For the cross section ratio in pp collisions, the energy and rapidity interpolated value discussed in Section 3 is used. The double ratio is also compared with the one measured in p-Pb collisions at √ s NN = 5.02 TeV [30]. Calculations from the Comovers + EPS09LO model [22] are also shown in Fig. 10 for comparison. In the Comovers + EPS09LO model, resonances may be dissociated via interactions with "comoving particles" (their nature, partonic or hadronic, not being defined in the model) produced in the same rapidity region. The dissociation is governed by the comover interaction cross sections, σ co−J/ψ = 0.65 mb and σ co−ψ(2S) = 6 mb, which are fixed from fits to low-energy experimental data. The main source of uncertainty in this model is the nPDF parameterisation, which is strongly correlated between the J/ψ and the ψ(2S) and thus cancels out when calculating the cross section ratio. Overall, the agreement between the model calculations and the measurements is good at both collision energies. The decrease of the double ratio with increasing collision energy in the model is due to the increase of the comover density. The measurement uncertainties do not allow for the experimental confirmation of such decrease of the double ratio.

Centrality dependence of the inclusive ψ(2S) nuclear modification factor
The nuclear modification factor of the ψ(2S) is calculated using Eq. 2. Figure 11 shows the inclusive ψ(2S) Q pPb as a function of N coll , for the backward and forward rapidity intervals, compared with the inclusive J/ψ Q pPb . At forward rapidity, the suppression and its centrality dependence are similar for the ψ(2S) and the J/ψ. At backward rapidity, on the contrary, a systematically stronger suppression of the ψ(2S) relative to the J/ψ is observed, except for the most peripheral and most central collisions, where the large uncertainties prevent a firm conclusion. The ψ(2S) Q pPb at √ s NN = 8.16 TeV shows the same dependence with the centrality of the collision than at √ s NN = 5.02 TeV [29].
Also shown in Fig. 11 are the results of model calculations. The EPS09s NLO + CEM calculations [75] of Q pPb are very similar for both ψ(2S) and J/ψ. The model fails to describe ψ(2S) results at forward rapidity, while the J/ψ results lie in the lower edge of the model calculation. At backward rapidity, the model calculation is close to the J/ψ data, although exhibiting different centrality trends, but fails at explaining the stronger ψ(2S) suppression. The transport model [23] calculations yield significantly smaller Q pPb for the ψ(2S) than for the J/ψ, with the difference being more pronounced in the Pb-going direction, where this difference increases with increasing centrality. The description of the forward rapidity results is fair for both charmonium states. At backward rapidity, the model tends to overestimate the ψ(2S) measurement in the most peripheral centrality classes. In this model, the lower Q pPb for the ψ(2S) than for the J/ψ is caused by a larger suppression of the ψ(2S) in the short QGP and the hadron resonance gas phases. Finally, the Comovers + EPS09LO model [22] predicts a significantly lower Q pPb for the ψ(2S) than for the J/ψ in the backward rapidity region. In the forward rapidity region the model uncertainties are too large to draw any firm conclusion. It is worth noting that the model uncertainties are largely correlated between the J/ψ and ψ(2S), as they are dominantly due to the nPDF parameterisation, and thus mostly cancel when calculating the double ratio as shown in Fig. 10. Nuclear shadowing is included using the EPS09 LO parameterisation [18] and the uncertainties of this parameterisation dominate the uncertainties of the model. The effect of the comovers, responsible for the stronger suppression of the ψ(2S) compared to the J/ψ, is stronger at backward rapidity due to the larger density of comovers in the Pb-going direction [22]. This model provides a fair description of ψ(2S) Q pPb at backward rapidity. However, the trend with centrality exhibited for the J/ψ does not reproduce the one observed in the data. Although not shown in the figure, the energy loss model [20] predicts sensibly the same Q pPb for the two reported charmonium states. Only models including final-state interactions are able to describe, at least qualitatively, a stronger suppression of the less bound ψ(2S) state than of the more tightly bound J/ψ state.
As for the J/ψ, it is also possible to estimate the Q prompt pPb of ψ(2S). In this case, the value of f B is about 0.18 and it is calculated using the LHCb measurements in pp collisions at √ s = 7 TeV for p T < 16 GeV/c and 2 < y < 4.5 [79]. Since the non-prompt ψ(2S) Q pPb has not been measured yet as a function of centrality, it is conservatively assumed to vary between 0.4 and 1 in all centrality classes for both forward and backward rapidity. That variation range for non-prompt ψ(2S) Q pPb englobes the centrality-integrated non-prompt ψ(2S) R pPb measured by LHCb at backward and forward rapidity in p-Pb collisions at √ s NN = 5.02 TeV [37] as well as all the inclusive ψ(2S) Q pPb reported here. The Q prompt pPb calculated under these assumptions is compatible within uncertainties with the inclusive one, showing a maximum difference of 25% with respect to the latter.

Summary
The study of the centrality dependence of the J/ψ and ψ(2S) production in p-Pb collisions at √ s NN = 8.16 TeV using the energy deposited in the neutron ZDC located in the Pb-going direction as the centrality estimator is presented. The J/ψ p T and p 2 T are reported for different centrality classes in the forward and backward rapidity regions covered by the ALICE muon spectrometer. The ∆ p 2 T measurement shows a p T broadening, relative to pp collisions, that increases from peripheral to central collisions, with larger values at forward than at backward y, except for the most peripheral events where similar values are seen in both rapidity intervals.
At forward rapidity, a clear suppression of J/ψ in p-Pb collisions compared to pp collisions is observed, which increases from peripheral to central collisions. At backward rapidity, the trend is opposite: the production of J/ψ relative to expectations from pp collisions is suppressed in peripheral collisions but enhanced in central collisions. The p T -and centrality-differential measurements of the J/ψ Q pPb indicate a stronger suppression in central than in peripheral collisions at low p T and forward rapidity, but with Q pPb approaching unity at high p T for all centrality classes. At backward rapidity, an enhancement is observed in central compared to peripheral collisions for p T > 3 GeV/c. The ratio B.R. ψ(2S)→µ + µ − σ ψ(2S) /B.R. J/ψ→µ + µ − σ J/ψ is compatible with the pp measurement in the most central and most peripheral collisions (within large uncertainties), whereas a decrease is observed in the semi-central and semi-peripheral events. Thus, in those centrality classes, the ψ(2S) production relative to the J/ψ is suppressed in p-Pb collisions compared to pp collisions. The nuclear modification factor of the ψ(2S) is compatible, within large uncertainties, with the one of the J/ψ in the most central and most peripheral events, but a stronger suppression of the ψ(2S) is observed in semi-central and semi-peripheral events, especially at backward rapidity.
The results presented here at √ s NN = 8.16 TeV confirm with improved statistical precision the earlier observations at √ s NN = 5.02 TeV and extend the p T reach up to 16 GeV/c for the J/ψ analysis. No significant dependence with collision energy is observed.
Theoretical models employing nPDF or energy loss mechanisms describe the centrality dependence of the J/ψ nuclear modification factor at forward rapidity but do not reproduce the shape at backward rapidity. The p T dependence of the J/ψ Q pPb in central collisions is not well described by the nPDF or energy loss based models, while the agreement is fair in peripheral collisions.
Among the three models considered, the one based only on nPDF cannot reproduce the ψ(2S) suppression. The model including final-state comover interactions describes the stronger ψ(2S) suppression at backward and forward rapidity, although the large model uncertainty prevents a firm conclusion at forward rapidity. The transport model is in good agreement at forward rapidity, but overestimates the ψ(2S) results at backward rapidity, especially in peripheral collisions.
The results presented here stress the need for a sound theoretical understanding of the production of quarkonia, including the excited states, in proton-nucleus collisions. Further experimental results expected from the future Run 3 and Run 4 of the LHC will push further our understanding of nuclear effects.