J/$\psi$ production as a function of charged-particle multiplicity in p-Pb collisions at $\sqrt{\textit{s}_{\rm NN}}~=~8.16$ TeV

Inclusive J/$\psi$ yields and average transverse momenta in p-Pb collisions at a center-of-mass energy per nucleon pair $\sqrt{s_{\rm NN}}$ = 8.16 TeV are measured as a function of the charged-particle pseudorapidity density with ALICE. The J/$\psi$ mesons are reconstructed at forward $(2.03<y_{\rm cms}<3.53)$ and backward ($-4.46<y_{\rm cms}<-2.96$) center-of-mass rapidity in their dimuon decay channel while the charged-particle pseudorapidity density is measured around midrapidity. The J/$\psi$ yields at forward and backward rapidity normalized to their respective average values increase with the normalized charged-particle pseudorapidity density, the former showing a weaker increase than the latter. The normalized average transverse momenta at forward and backward rapidity manifest a steady increase from low to high charged-particle pseudorapidity density with a saturation beyond the average value.

Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8. 16 TeV ALICE Collaboration

Experimental setup and data samples
In this section, the detector subsystems relevant for this analysis are presented. A complete description of the ALICE detector and its performance can be found in [42,43].
The muon spectrometer [42,43] covers the pseudorapidity window of −4.0 < η < −2.5 and consists of: a 4 m long composite front absorber, corresponding to about 10 interaction lengths (10 λ int ), starting at 90 cm from the nominal interaction point, ten layers of muon tracking chambers (MCH), coupled to a dipole magnet with a 3 Tm field integral, and four layers of muon trigger chambers (MTR). The MCH and MTR systems are separated by an additional iron wall of about 7.2 λ int that absorbs the remaining hadronic and low-momentum particle contamination. A rear absorber positioned downstream of the MTR filters out the background from beam-gas interactions. A conical absorber surrounds the beam pipe and protects the spectrometer against secondary particles produced mainly by large-η primary particles interacting with the beam pipe.
The Silicon Pixel Detector (SPD) [44] is the innermost part of the Inner Tracking System (ITS). It consists of two cylindrical silicon pixel layers at radial distances of 3.9 and 7.6 cm from the beam line.
The respective pseudorapidity coverage of the two layers are |η| < 2 and |η| < 1.4. The SPD is used to reconstruct the primary vertex and to measure the charged-particle pseudorapidity density at midrapidity.
The V0 scintillator arrays [45] are located at each side of the interaction point, covering the pseudorapidity ranges of −3.7 < η < −1.7 and 2.8 < η < 5.1. In this analysis, the V0 provides an online trigger and helps to reject contamination from beam-gas events.
The neutron Zero Degree Calorimeter (ZDC) [42] located at about 112.5 m on either side from the interaction point are used to reject electromagnetic interactions and beam-induced background.
The results presented in this letter are obtained with data recorded during the p-Pb run at √ s NN = 8. 16 TeV in 2016. The J/ψ are reconstructed in the dimuon channel with data taken in two different beam configurations. Due to the asymmetry of the beam energy per nucleon in p-Pb collisions at the LHC, the nucleon-nucleon center-of-mass rapidity frame is shifted by ∆y = 0.465 in the direction of the proton beam. As a consequence, the J/ψ are measured in the forward rapidity range of 2.03 < y cms < 3.53 (with protons going in the direction of the muon spectrometer, p-going direction) and in the backward rapidity region −4.46 < y cms < −2.96 (Pb-going direction). Events used in this analysis were collected with a dedicated dimuon trigger which requires the coincidence of signals in both V0 arrays (minimum bias trigger, MB) with at least two opposite-sign muons registered in the MTR. The trigger has an adjustable online threshold, which for this data sample was set to only accept muons with transverse momenta p T > 0.5 GeV/c (p T for which an efficiency of 50% is reached). The p T differential single-muon trigger efficiency reaches a plateau of ∼ 96% at p T ∼ 1.5 GeV/c. In this data-taking period, the maximum pile-up probability was about 4%. A dedicated event-selection strategy-exploiting the signals in the V0 and the ZDC, the correlation of the number of clusters and track segments reconstructed in the SPD, as well as an algorithm to tag events with multiple vertices-allowed us to keep the pile-up below 0.5% for the analysed events, even at large multiplicities. The data sample analyzed corresponds to an integrated luminosity of L int = 7.2 ± 0.2 nb −1 (10.6 ± 0.3 nb −1 ) for the p-going (Pb-going) configuration [46].
3 Charged-particle multiplicity measurement The charged-particle pseudorapidity density (dN ch /dη) is measured at midrapidity exploiting the information provided by the SPD detector [47,48]. It is evaluated by counting the number of tracklets (N tracklet ), i.e. track segments joining pairs of hits in the two layers of the SPD pointing to the primary vertex. The primary vertex is also computed with the SPD information. To minimize non-uniformities in the SPD acceptance, only events with a z-vertex position determined within |z vtx | < 10 cm are considered, and tracklets are counted within |η| < 1.
Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8. 16 TeV ALICE Collaboration The raw N tracklet counts are corrected (N corr tracklet ) for the variation of the detector conditions with time (fraction of active SPD channels) and its limited acceptance as a function of z vtx using a data-driven eventby-event correction [29,30]. This correction ensures a uniform response as a function of z vtx . In this analysis, the correction is done by renormalising the N tracklet (z vtx ) distributions to the overall maximum with a Poissonian smearing to account for the fluctuations. The events are sliced in N corr tracklet intervals. Monte Carlo (MC) simulations using the DPMJET [49] event generator and the GEANT3 transport code [50] are used to estimate dN ch /dη from N corr tracklet . A second order polynomial correlation is assumed between these two quantities for the full N corr tracklet interval. Several sources of systematic uncertainty were taken into account. Possible deviations from the second order polynomial correlation were estimated by using other functions to quantify the correlation or MC averages in each interval, with values ranging from 0.1% at intermediate multiplicities to 6.9% (5.8%) at the lowest (highest) multiplicity intervals. The systematic uncertainty on the residual z vtx dependence due to differences between data and MC amounts to 3%. Finally, the event generator influence was considered and evaluated by comparing the DPMJET simulations with events generated in EPOS [51], resulting in a 2% uncertainty. The average charged-particle pseudorapidity density, dN ch /dη , in non-single diffractive (NSD) events was obtained from an independent analysis and amounts to dN ch /dη = 20.33 ± 0.83 (20.32 ± 0.83) in p-Pb (Pb-p) collisions for |η| < 1 [48], where the quoted uncertainty is systematic. Table 1 summarizes the contributions to the normalized charged-particle multiplicity uncertainty. The total uncertainty is evaluated assuming that the different sources are uncorrelated.

J/ψ measurement
The normalized J/ψ yield, i.e. the yield in each multiplicity interval i normalized to the multiplicityintegrated value, is evaluated as from the reconstructed number of J/ψ, N J/ψ , the number of minimum bias (MB) events equivalent to the analysed dimuon sample, N eq MB , the J/ψ acceptance and efficiency correction, (Aε) J/ψ , and the NSD event selection efficiency in the minimum bias sample, ε MB .
The J/ψ are reconstructed for each multiplicity interval by combining opposite-sign muons and computing the invariant mass of the pairs. The muon identification is ensured by requiring that the track candidates reconstructed in the MCH have a matching track segment in the MTR. Furthermore, the individual tracks must fulfill the following criteria to make sure they are within the acceptance of the spectrometer: their radial distance from the beam axis at the end of the front absorber is within 17.6 < R abs < 89.5 cm and their pseudorapidity in the detector reference frame is within −4 < η < −2.5.
To extract the signal, the invariant-mass distributions are first corrected for the J/ψ acceptance times efficiency (Aε), differentially in p T and y. The resulting distributions are then fitted with a superposition of J/ψ and ψ(2S) signals and a background lineshape. Various combinations of lineshapes are Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8. 16 TeV ALICE Collaboration used in order to evaluate the signal counts and their uncertainties. The two charmonium resonances are parametrized by a sum of either two Crystal Ball or two pseudo-Gaussian functions with power-law tails [52]. The tail parametrizations are fixed to the values determined from either fits of the J/ψ signal from MC simulations or to values taken from fits to the multiplicity-integrated distribution in p-Pb data at √ s NN = 8. 16 TeV [13] and in pp data at √ s = 13 TeV [53]. The tails obtained from fitting the multiplicity-integrated distributions using the Crystal Ball function are also considered, and fixed in the binned fits. The J/ψ peak mean position and width are left free in the multiplicity-integrated fit, whilst the ψ(2S) ones are bound to those of the J/ψ following the same procedure as in [54]. Note that the ψ(2S) yields obtained are not physical values, as the invariant-mass spectrum is corrected by the Aε correction for the J/ψ. In the multiplicity-differential fits, the mass and width of the J/ψ peak are fixed to the integrated values to ensure the convergence of the fits in the few cases where statistical significance is low. The background is parameterized by either a sum of two exponentials or the product of an exponential and a fourth-order polynomial. Two fit mass ranges are taken into account when computing the average number of J/ψ and its uncertainty: 1.7 < m µ µ < 4.8 GeV/c 2 and 2.0 < m µ µ < 5.0 GeV/c 2 . Examples of fits at low, intermediate, and high multiplicity for data in the rapidity range 2.03 < y cms < 3.53 are shown in Fig. 1. The signal lineshape is found to be independent of multiplicity, while the background does change with multiplicity. Therefore, in order to minimize the uncertainty on the signal extraction, the same signal lineshape is used in the fit function for both the numerator and denominator in Eq. 1.
The number of equivalent MB events N eq MB is computed from the number of dimuon triggered events, N µ µ , and the normalization factor of dimuon triggered to MB events (calculated as explained in next section) as N eq MB = F norm · N µ µ . The number needs to be corrected for by the NSD event selection efficiency, ε MB = (97 ± 1)% [48], to take into account the fraction of events without a reconstructed SPD vertex that are rejected. This factor ε MB is found to be independent of the charged-particle multiplicity in all the intervals studied, with the exception of the lowest multiplicity interval, where it decreases by 1%.
The J/ψ acceptance and efficiency correction is obtained from MC simulations as a function of p T and y cms . The J/ψ are generated using p T and y cms distributions tuned to data [13]. They are simulated to decay into a muon pair using EvtGen [55]. The final state radiation is described with PHOTOS [56]. The acceptance and efficiency correction is independent of multiplicity in the measurement intervals. Therefore, when estimating the uncertainty on the MC input, only the possible variation of the input p T and y cms distributions is taken into account by using as input a subsample of the lower/higher multiplicity events.
To extract the J/ψ mean transverse momentum p J/ψ T , the Aε-corrected transverse momentum of the dimuon pair is fitted with the following function [26]: where the ratios of signal over the sum of signal and background of the two charmonium states α J/ψ = S J/ψ /(S J/ψ + S ψ ′ + B) and α ψ ′ = S ψ ′ /(S J/ψ + S ψ ′ + B) are fixed to the value extracted from fitting the invariant-mass spectrum corrected by the J/ψ Aε. The background is described by a function p bkgd T (m µ µ ). Two functional forms are used: either a sum of two exponentials or the product of an exponential and a fourth-order polynomial. Note that the p ψ ′ T does not represent a physical mean transverse momentum of the ψ(2S) as the spectra are corrected by the Aε for J/ψ.

Systematic uncertainties
The following sources of systematic uncertainty on the J/ψ yields in multiplicity classes are considered: (i) the signal extraction, (ii) the normalisation, (iii) the effect of resolution and pile-up, (iv) the event-by-event N tracklet to N corr tracklet correction, and (v) the event selection efficiency of the NSD event class. For the measurement of the yields in each multiplicity interval normalized to the event average, the systematic uncertainties are estimated directly for this ratio. Details on the signal extraction uncertainty were addressed in the previous section. The values are estimated by varying the signal and background shapes of the fit function, as well as by varying the invariant-mass range of the fit. The systematic uncertainty is computed as the root-mean-square of the uncertainties on the ratio for each of these fits, ranging between 0.8-2.3% (0.5-1.9%) at forward (backward) rapidity, being larger at large multiplicities where the number of events is smaller. The normalisation factor of the dimuon triggered to MB events F norm is studied using three alternative methods [13]. The first method evaluates the probability of a coincidence of a dimuon-and a MB-triggered event in a MB-triggered data set. The second method exploits the higher probability of occurrence of a single-muon trigger by looking at the product of the probability of coincidence of a single-muon-and MB-triggered event and of the probability of finding a dimuon event in the single-muon triggered data. The third method is based on information from the trigger scalers. The run-by-run spread of the F norm /F i norm values, ratio of the normalisation values in the integrated and specific multiplicity intervals, computed for these three methods determines a 2.5% systematic uncertainty, independent of multiplicity. The effect of the method of choice for the event-by-event correction from N tracklet to N corr tracklet on the J/ψ yield is also studied [29,30]. Both the randomisation function (Poisson or binomial) and the reference normalisation of the correction are varied. The Poissonian smearing is applied when the maximum is selected as normalisation reference, while the binomial correction should be used when considering all other possible reference values (in our case the minimum). The influence of these modifications on the yield ranges from 0.1% to 2.6% (4.3%) at forward (backward) rapidity, as a function of multiplicity. The uncertainty coming from pile up and multiplicity axis resolution is estimated as a single contribution by repeating the analysis multiple times with a different randomisation seed for the event-by-event correction, or introducing a small shift of the N corr tracklet intervals, or varying the pile-up rejection criteria. The uncertainty amounts to 2%, independent of multiplicity. The uncertainty on the event selection efficiency for the NSD event class is estimated as in Ref. [48]. The uncertainty amounts to 1% and is correlated in all multiplicity intervals. Table 2 summarizes all contributions to the systematic uncertainty on the normalized yield.
Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8. 16 TeV ALICE Collaboration For the p T , the effects of the uncertainty on the p T extraction procedure and of the Aε are considered. Similar to the yields, the signal extraction uncertainty is estimated by varying the fit function and its range. In addition, as the S/(S + B) terms in Eq. 2 are fixed in the fit to the p T invariant-mass spectrum, the influence of the statistical uncertainty on the J/ψ signal S is introduced via a Gaussian smearing of S (with respect to its statistical uncertainty) to prevent artificially minimising the uncertainty. It ranges from 0.2% to 3.0% (1.2%) at forward (backward) rapidity, increasing with multiplicity as a consequence of the smaller number of events. The uncertainty on the absolute p T also takes into account the uncertainty on: (i) the MC input shapes as a function of p T and y cms , ranging from < 0.1 to 6% (< 0.1 to 11%) at forward (backward) rapidity, (ii) the tracking efficiency, 1% [13], (iii) the trigger efficiency, 2.6% (3.1%) [13], and (iv) the matching efficiency between the tracks in the MCH and the MTR, 1% [13]. To evaluate the uncertainty on the MC input, the data are divided into two multiplicity classes at the mean of the N corr tracklet distribution for each rapidity interval. For each of these bins, the p T is estimated using a modified Aε correction, which was re-weighted to better describe the p T -and y-dependent distributions of J/ψ in given bin. The systematic uncertainty is taken as the difference of the original p T value, computed with the initial Aε correction, and the new p T estimated with re-weighted correction. The uncertainty on all the measured multiplicity intervals is extrapolated from these two values assuming that in each class the uncertainty is proportional to the p T . The contributions of the tracking, the trigger and their matching to the uncertainty are correlated between multiplicity intervals. The normalized p T values are only affected by the uncertainty on the signal extraction procedure and the MC input, which is partly correlated in multiplicity and ranges from < 0.1% to 2% (< 0.1% to 4%) in the forward (backward) rapidity interval. Table 3 summarizes all contributions to the average, p T , and normalized average, p T p int T , p T measurements. The correlated uncertainties are added in quadrature and quoted in the plot as a text.

Results and discussion
The normalized J/ψ yield, at forward and backward rapidities, is presented in Fig. 3 as a function of the normalized charged-particle pseudorapidity density, measured at midrapidity (|η| < 1). The normalized yield increases with increasing multiplicity in both rapidity intervals. The yield at backward rapidity Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8.16 TeV ALICE Collaboration grows faster than the one at forward rapidity, reaching values above those expected from a linear (with slope unity) increase at large multiplicities. On the other hand, at forward rapidity the values show a slower-than-linear increase at large multiplicities. The forward and backward rapidity yields cross a linear increase estimate (and each other) at around 1.5 times the average multiplicity. The underlying mechanism remains unclear. The forward (p-going) rapidity region probes the Pb-nucleus low Bjorken-x regime (x Pb ∼ 10 −5 in a naive 2-body calculation for p T = 0), whereas the backward (Pb-going) rapidity is sensitive to the intermediate-to-large values (x Pb ∼ 10 −2 ). The observed suppression of the p Tand multiplicity-integrated J/ψ yield at forward rapidity, with respect to pp collisions, is described by different cold nuclear matter models considering the probed shadowing/saturation domain [13]. The centrality-differential measurements at √ s NN = 5.02 TeV [23] of the nuclear modification factor, p T and p 2 T , corresponding to relative multiplicities of at most 2.5 times the average one, can also be described by these models. The contribution from beauty-quark decays to the inclusive J/ψ yield amounts to ∼ 10% [57]. It is not expected to affect significantly these results, since a similar trend was observed for prompt and non-prompt J/ψ as a function of the charged-particle pseudorapidity density in pp collisions [30]. Moreover, the autocorrelations influence is negligible in this analysis due to the large rapidity gap between the measurement of the charged-particle multiplicity and the J/ψ yield [58]. Figure 4 presents p T as a function of the relative charged-particle pseudorapidity density, in p-Pb collisions at √ s NN = 8. 16 TeV. The measured p T is systematically smaller at backward than at forward rapidity. This is also true for the multiplicity-integrated value, which is consistent with the observed decrease of p T with increasing |y cms | in pp collisions [59]. The p T increases steadily for multiplicities below the average, and saturates above the average multiplicity. Two naive scenarios are typically considered to explain high-multiplicity events: the incoherent superposition of multiple parton-parton collisions, or single parton interactions with higher energy transfer. One would expect the latter to be Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8.16 TeV ALICE Collaboration characterized by a higher p T of the produced J/ψ. Reality is probably somewhere in between these two simplified scenarios. The simultaneous increase of the p-Pb yield together with the saturation of p T at large multiplicities may point to J/ψ production from an incoherent superposition of parton-parton collisions. The measured yield in p-Pb collisions can be described with the EPOS 3 event generator [60, 61] (see Fig. 5) based on a combination of Gribov-Regge theory and pQCD: where the individual scatterings are identified with parton ladders emerging as flux tubes, the existence of multiple nucleon-nucleon collisions in p-Pb collisions is accounted for, the initial conditions of the collision are modified due to CNM effects including parton saturation, and slow string segments (far from the surface) can be further mapped to fluid dynamic fields using a core-corona description. The J/ψ bound-state formation in EPOS 3 assumes a color-evaporation approach, i.e. it is associated to a charm quark-anti-quark pair in a given mass range. The influence of the 3D+1 viscous hydrodynamic evolution of the bulk (starting from flux tube initial conditions) in the EPOS 3 calculation is small (see Fig. 5). However the number of simulated events at large multiplicities is limited and does not allow us to elucidate possible hydrodynamic effects. EPOS 3 description of the measurement suggests J/ψ production from an incoherent superposition of parton-parton collisions.
The normalized J/ψ yield and p T are compared with the results in p-Pb collisions at √ s NN = 5.02 TeV [26] in Fig. 6 and 7, respectively. The measurements are in remarkable agreement, within the uncertainties, at both energies and rapidities. These results extend the probed charged-particle pseudorapidity density interval, both at low and high multiplicity, examining events of up to almost six times the average value. The more precise √ s NN = 8. 16 TeV data evidence a continuous increase of the normalized yield with multiplicity up to the largest multiplicities attained. The similarities at √ s NN = 8. 16 TeV and √ s NN = 5.02 TeV suggest a common origin of the multiplicity trend, with a mechanism whose effect varies with rapidity, but might have a small dependence on the collision energy, in the explored interval. This is consistent with the large variation of the probed x Pb with rapidity and its relative slow evolution Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8.16 TeV ALICE Collaboration   The corresponding dN ch /dη in |η| < 1 for those measurements is 6.01 ± 0.01 (stat.) +0.20 −0.12 (syst.) [29] and 544.7 ± 0.2 (stat.) ± 7.3 (syst.) for 0-90% centrality [63], respectively. The ratio of the yields over the corresponding charged-particle multiplicity is also shown in Fig. 8. The trend exhibited by the pp data is similar to the one observed in the backward (Pb-going) direction. It should be noted that the pp results are normalized to the inelastic 'INEL' event class, whereas the p-Pb measurements are normalized to the non-single-diffractive 'NSD' one. In p-Pb collisions these two event classes mostly overlap when comparing MB results [48]. The Pb-Pb data also show a faster-than-linear increase with the normalized charged-particle pseudorapidity density. They are compatible within uncertainties with the p-Pb backward rapidity result in the restricted multiplicity interval of the measurement. Whereas the pp and p-Pb data include J/ψ with p T > 0, the Pb-Pb data points include J/ψ with 0.3 < p T < 12 GeV/c to reduce the low-p T contribution from photoproduction, which is significant only in more peripheral collisions [62,64,65]. The overall increase of multiple parton-parton collisions with the colliding system (from pp up to Pb-Pb collisions) described in [66] is expected to cancel in the relative quantities reported in this publication, sensitive to the relative evolution with charged-particle density in a given colliding system. Model calculations are needed to interpret the similarities of pp, p-Pb, and Pb-Pb normalized J/ψ yields at large rapidity as a function of the normalized charged-particle pseudorapidity density at midrapidity.

Conclusions
The production of inclusive J/ψ at large rapidities in p-Pb collisions at √ s NN = 8. 16 TeV is reported as a function of the charged-particle pseudorapidity density at midrapidity. The normalized J/ψ yield shows an increase with increasing normalised charged-particle pseudorapidity density. The yield at backward rapidity grows faster than the forward rapidity one, reaching values above those of the linear (with slope unity) increase estimate at large normalised multiplicity, whereas the values at forward rapidity show a slower-than-linear increase. The trends of the normalised yield are reproduced by the EPOS 3 [60,61] event generator. The p T is smaller at backward than at forward rapidity, consistent with the expected Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8.16 TeV ALICE Collaboration  Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8.16 TeV ALICE Collaboration    Figure 8: Top: Normalized yield of inclusive J/ψ as a function of the normalized charged-particle pseudorapidity density, measured at midrapidity, in various collision systems. Bottom: Ratio of the normalized yields to the corresponding normalized charged-particle pseudorapidity density. The pp results are normalized to INEL collisions [29], whereas p-Pb ones refer to the NSD event class; all for p T > 0. The Pb-Pb data points include J/ψ with 0.3 < p T < 12 GeV/c to reduce the low-p T contribution from photoproduction, which is significant only in more peripheral collisions [62,64,65]. The vertical bars represent the statistical uncertainties, the boxes the systematic ones. The dashed line indicates the one-to-one correlation, to guide the eye. √ s NN = 8. 16 TeV ALICE Collaboration softening of the spectra with increasing |y cms |. The p T increases steadily for multiplicities below the average, and saturates above the average multiplicity. The simultaneous increase of the yield together with the saturation of p T may point to J/ψ production from an incoherent superposition of parton-parton collisions. These measurements show trends compatible with those observed at √ s NN = 5.02 TeV [26] in p-Pb collisions, but in this work an improved precision and extended multiplicity coverage were reached. The similarities suggest a common origin, with a mechanism whose effect varies with rapidity, but with only a small dependence (if any) on the collision energy.
Multiplicity dependent J/ψ production in p-Pb at √ s NN = 8. 16 TeV ALICE Collaboration

Acknowledgements
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The    [13] ALICE Collaboration, S. Acharya et al., "Inclusive J/ψ production at forward and backward rapidity in p-Pb collisions at √ s NN = 8. 16          [46] ALICE Collaboration, "ALICE luminosity determination for p-Pb collisions at √ s NN = 8.16 TeV", ALICE-PUBLIC-2018-002 (Apr, 2018) . http://cds.cern.ch/record/2314660.