J/ψ production as a function of charged-particle multiplicity in p-Pb collisions at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 8.16 TeV

Inclusive J/ψ yields and average transverse momenta in p-Pb collisions at a center-of-mass energy per nucleon pair sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 8.16 TeV are measured as a function of the charged-particle pseudorapidity density with ALICE. The J/ψ mesons are reconstructed at forward (2.03 < ycms< 3.53) and backward (−4.46 < ycms< −2.96) center-of-mass rapidity in their dimuon decay channel while the charged-particle pseudorapidity density is measured around midrapidity. The J/ψ 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.


Introduction
Quarkonium states have long been considered as probes of the Quark-Gluon Plasma (QGP) produced in ultra-relativistic heavy-ion collisions [1]. The large color-charge density in the plasma prevents the formation of bound states, in an analogous process to the Debye screening for electromagnetic processes [2]. The suppression of J/ψ production in nucleus-nucleus (AA) with respect to proton-proton (pp) collisions was observed by several experiments [3][4][5][6][7][8][9][10][11]. To determine whether the origin of this suppression is the influence of the QGP or of Cold Nuclear Matter (CNM), data on proton(deuteron)-nucleus collisions are also scrutinized.
The measurements in p-Pb collisions at the LHC show a suppression of J/ψ production [12][13][14], with respect to pp collisions, at low transverse momentum (p T ) and forward center-of-mass rapidity (p-going direction, positive y cms ), consistent with various combinations of CNM effects: modification of the parton distribution functions (PDFs) in nuclei, i.e. shadowing [15,16], the Color-Glass Condensate (CGC) [17,18], or coherent parton energy loss [19]. The measurement of ψ(2S) production in p-Pb collisions [20] exhibits a larger suppression, with respect to pp collisions, than the one measured for J/ψ, both at forward and backward rapidity, which was not expected from CNM predictions. This effect is reproduced by models which consider the break-up of the bound quark-anti-quark pair via interactions with the final-state comoving particles [21,22].
The p-Pb data at the center-of-mass energy per nucleon-nucleon collision of √ s NN = 5.02 TeV [23,24] showed that these effects depend on the centrality of the collision, as -1 -

JHEP09(2020)162
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 lowmomentum 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 beaminduced 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].

JHEP09(2020)162
Source |η| < 1 N corr tracklet to dN ch /dη correlation 0.1-6.9(5.8)% z-vertex dependence 3% Monte Carlo event generator 2% dN ch /dη 4% Table 1. Sources of systematic uncertainties on the normalized charged-particle multiplicity. For the N corr tracklet to dN ch /dη correlation an interval is quoted, varying with multiplicity, with a different maximum uncertainty for the Pb(p)-going configuration. 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.
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 event-by-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 multiplicity-integrated value, is evaluated as

JHEP09(2020)162
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 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 figure 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. 4.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    Figure 1. Opposite-sign muon pair invariant mass distributions for selected multiplicity intervals, corrected for the J/ψ acceptance and efficiency, at forward rapidity. The distributions are shown together with a typical fit function (solid line, see text for details). The J/ψ signal contribution is also depicted by a dot-dashed red line, and the background by a dotted line.
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/ψ. Figure 2 illustrates typical p µµ T distributions for selected multiplicity intervals.

Systematic uncertainties
The following sources of systematic uncertainty on the J/ψ yields in multiplicity classes are considered: 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-meansquare 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 -7 -

JHEP09(2020)162
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.
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. 4.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.  Table 3. Systematic uncertainty sources on the average and normalized average p T . The values in parentheses correspond to the multiplicity-integrated uncertainties related to the signal extraction. The contributions marked with an asterisk are correlated in multiplicity. The uncertainty on MC input, marked with a diamond, is partially correlated in multiplicity.
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 figure 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 grows faster than the one at forward rapidity, reaching values above those expected from a linear (with slope unity) increase at large -9 -JHEP09(2020)162  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 T -and 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 centralitydifferential 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].  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 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 figure 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 figure 5).
-11 -JHEP09(2020)162  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 on the collision energy (typically a factor of 2 in the simplified 2-body picture). Figure 8 presents a comparison of the normalized J/ψ p-Pb yields with results from pp collisions at √ s NN = 7 TeV [29] (2.5 < y cms < 4.0) and Pb-Pb collisions at √ s NN = 5.02 TeV [62] (2.5 < y cms < 4.0). 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 chargedparticle multiplicity is also shown in figure 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--12 -JHEP09(2020)162  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 -13 -JHEP09(2020)162 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 -14 -JHEP09(2020)162   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. 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 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 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.

Acknowledgments
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. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

JHEP09(2020)162
[26] ALICE collaboration, J/ψ production as a function of charged-particle pseudorapidity density in p-Pb collisions at √ s NN = 5.02 TeV, Phys. Lett [29] ALICE collaboration, J/ψ production as a function of charged particle multiplicity in pp collisions at √ s = 7 TeV, Phys. Lett [30] ALICE collaboration, Measurement of charm and beauty production at central rapidity versus charged-particle multiplicity in proton-proton collisions at