Forward rapidity J/ψ production as a function of charged-particle multiplicity in pp collisions at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 5.02 and 13 TeV

The production of J/ψ is measured as a function of charged-particle multiplicity at forward rapidity in proton-proton (pp) collisions at center-of-mass energies s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 5.02 and 13 TeV. The J/ψ mesons are reconstructed via their decay into dimuons in the rapidity interval (2.5 < y < 4.0), whereas the charged-particle multiplicity density (dNch/dη) is measured at midrapidity (|η| < 1). The production rate as a function of multiplicity is reported as the ratio of the yield in a given multiplicity interval to the multiplicity-integrated one. This observable shows a linear increase with charged-particle multiplicity normalized to the corresponding average value for inelastic events (dNch/dη/〈dNch/dη〉), at both the colliding energies. Measurements are compared with available ALICE results at midrapidity and theoretical model calculations. First measurement of the mean transverse momentum (〈pT〉) of J/ψ in pp collisions exhibits an increasing trend as a function of dNch/dη/〈dNch/dη〉 showing a saturation towards high charged-particle multiplicities.


Introduction
The study of charmonia, bound states of a charm and anti-charm quark (cc), production in hadronic collisions represents a stringent test for theory of the strong interaction, quantum chromodynamics (QCD). While the production of heavy quark-antiquark pair can be calculated within pQCD, their evolution into a bound colorless qq pair is a non-perturbative process. Different theoretical approaches exist, which mainly differ in the treatment of the (non-perturbative) bound state formation [1][2][3][4]. Inclusive differential studies of J/ψ production as a function of transverse momentum (p T ) and rapidity (y) have been performed at the LHC in pp collisions at center-of-mass energies ranging from √ s = 2.76 TeV up to 13 TeV [5-10]. The measurements are described by the sum of non-relativistic quantum chromodynamics (NRQCD) calculations for the prompt component, and fixed-order nextto-leading logarithm calculations for the contribution arising from beauty-hadron decays (referred as "non-prompt") [10]. Charmonium polarization measurements add additional constraints to models [11][12][13][14]. While the polarization of prompt J/ψ measured by LHCb in the low-p T region is described by the NRQCD models, the predicted transverse polarization at high-p T is not observed in any of the LHC experiments [14].
At the LHC energies, several parton interactions may occur in a single pp collision. Multi-parton interactions (MPI) influence the production of light quarks and gluons, affecting the total event multiplicity [15][16][17]. If the MPI also affect heavy-flavour production, this -1 -

JHEP06(2022)015
could introduce a dependence on the charged-particle multiplicity. A faster-than-linear correlation has been observed between the production of charged particles and that of prompt D mesons, as well as that of inclusive, prompt, and non-prompt J/ψ in pp collisions at √ s = 7 TeV [18], suggesting that the heavy-flavour quark production mechanism is at the origin of this trend, while hadronization does not have a dominant effect. Recent results on inclusive J/ψ (at midrapidity) in pp collisions at √ s = 13 TeV [19] confirm the observed correlation and extend the measurement multiplicity reach. In addition, measurements of the associated production of charmonia, bottomonia and/or open-charm hadrons in pp collisions at √ s = 7 and 8 TeV indicate that double-parton scatterings have a significant contribution to their production [20,21].
The structures seen in two-particle angular correlations in high-multiplicity pp and p-Pb collisions at the LHC are similar to those in Pb-Pb data [22][23][24][25][26]. These structures in Pb-Pb collisions are interpreted as signatures of the collective motion of particles in the quark-gluon plasma. In the heavy-flavour sector, the correlation of J/ψ and charged particles in p-Pb collisions also shows similar features to those observed with charged hadrons [27,28]. The comparison with Pb-Pb measurements [29][30][31] suggests that there is a similar effect responsible for azimuthal asymmetries in both collision systems. The results at √ s = 13 TeV extend the charged-particle multiplicity reach up to about seven times minimum bias multiplicity, corresponding to about 50 charged particles per unit of rapidity, values similar to those of p-Pb collisions (about 40 [32]) where collective-like effects have been observed.
In this publication, the measurement of the inclusive yield and the first moment of the J/ψ p T distribution as a function of the charged-particle multiplicity in pp collisions at √ s = 5.02 and 13 TeV is presented. J/ψ mesons are reconstructed via their dimuon decay channel at forward rapidity (2.5 < y < 4), 1 whereas the charged-particle multiplicity is measured at midrapidity (|η| < 1). These measurements complement those performed for J/ψ at forward rapidities in pp collisions at √ s = 7 TeV [33].

Experimental apparatus and data sample
The ALICE apparatus is described in detail in refs. [34,35]. The detectors employed in this measurement, namely V0, Silicon Pixel Detector (SPD), and Muon Spectrometer (MS), are described below. The V0 detector consists of two scintillator hodoscopes located on each side of the interaction point (2.8 < η < 5.1 and −3.7 < η < −1.7) [36]. It provides a minimum bias (MB) trigger which requires a signal in both hodoscopes. The charged-particle multiplicity at midrapidity is measured using the SPD [37]. The SPD consists of two layers of silicon pixel detectors covering pseudorapidity ranges |η| < 2 and |η| < 1.4, respectively. It is used to reconstruct the primary vertex as well as tracklets, short two-point track segments covering the pseudorapidity region |η| < 1.4. Tracklets are required to point to the primary interaction vertex within ±1 cm in the transverse plane and ±3 cm in the beam (z) direction.

JHEP06(2022)015
The muons originating from J/ψ decays are detected at forward rapidity (2.5 < y < 4.0) in the MS. The MS is composed of five tracking stations with two planes of Cathode Pad Chambers for the first two stations and two planes of Cathode Strip Chambers for the rest. The third station is placed inside a dipole magnet with a field integral of 3 T m. Two trigger stations, each with two planes of Resistive Plate Chambers, are positioned downstream of the tracking system and provide a single muon as well as a dimuon trigger. A 4.1 m long front absorber of 10 interaction lengths (λ int ) is placed between the interaction point and the first tracking station to stop the high hadron flux. Hadrons which escape this front absorber are further filtered out by a second absorber, a 1.2 m long (7.2 λ int ) thick iron wall, placed between the tracking and the triggering system, which also removes lowmomentum muons originating from pion and kaon decays. A conical absorber shields the muon spectrometer against the secondary particles produced by the interaction of primary particles in the beam pipe throughout the entire length of the MS.
The reported results are based on data collected in pp collisions at √ s = 5.02 and 13 TeV by ALICE during 2015 and 2016, respectively. The J/ψ production is measured using dimuon triggered events. A dimuon trigger is obtained as the coincidence of a MB trigger and at least one pair of track segments reconstructed in the muon trigger system with low p T threshold of 0.5 GeV/c. The analysis presented in this publication utilizes ∼1.2 million and ∼121.1 million dimuon trigger events for pp collisions at √ s = 5.02 and 13 TeV, respectively. This corresponds to an integrated luminosity of ≈ 109.1 nb −1 (≈ 4.9 pb −1 ) in pp collisions at √ s = 5.02 (13) TeV. Events containing more than one distinct vertex were tagged as pileup and discarded from the analysis. The maximum probability of the pileup was about 5×10 −3 for the data collected at √ s = 13 TeV, while at √ s = 5.02 TeV the pileup was below 2.5% [10].

Data analysis
J/ψ yield as well as dN ch /dη are measured for INEL > 0 events, which are defined as inelastic collisions for which at least one charged-particle track is recorded within |η| < 1.0. Beam-induced background events are removed by using the time information from the V0 detectors and the correlation between the number of clusters and track segments reconstructed in the SPD. Pileup events, i.e. events with multiple collision vertices are removed as follows: using the V0 time information to select events within a bunch crossing time window (out-of-bunch pilepup); using a vertex finding algorithm based on the SPD [38] to remove events with multiple vertices (in-bunch pileup). The impact of a possible residual contamination from in-bunch pile-up was estimated by dividing the data sample in two groups, according to the average pileup probability per bunch crossing below/above 0.6%. A negligible effect was observed on the final results by analysing the sub-sample below 0.6% average probability per bunch crossing.

Charged-particle multiplicity density
The charged-particle multiplicity per unit of pseudorapidity (dN ch /dη) is estimated using information from the SPD. A pair of hits in the two SPD layers that are pointing towards -3 -JHEP06(2022)015 the interaction vertex are chosen to form tracklets (track segments). The dN ch /dη is proportional to the measured number of tracklets (N trk ). Primary charged particles, defined as all particles produced in the collision (including products of strong and electromagnetic decays) except those from the weak decay of strange hadrons [39], are considered for the analysis.
In order to minimize non-uniformities in the SPD acceptance, only events with a zvertex position within |z SPD vtx | < 10 cm are considered, and the tracklet multiplicity is measured within |η| < 1.
In addition, the variation of the SPD efficiency and its number of inactive channels throughout the data taking period, causes the measured N trk to vary with z SPD vtx . About 4 to 15% dead channels were registered for the analyzed data sets. In order to correct for these detector effects, a data-driven correction is applied to equalize N trk variation as a function of z SPD vtx and time on an event-by-event basis [18,33]. The N trk for each event reconstructed for a given z-vertex position is corrected by the average fraction of missing tracklets at the interaction z-vertex with respect to a reference multiplicity. The maximum of the average SPD tracklet multiplicity as a function of z SPD vtx position was chosen as the reference value or N corr trk (z SPD vtx ) profile, similarly as done in previous analyses [18,33]. However, minimum reference multiplicity has also been tested and the result was found to be consistent with the maximum one. The correction factors are smeared randomly using a Poisson (or Binomial depending on the reference value) distribution [33,40].
Monte Carlo (MC) simulations based on PYTHIA6 [41,42], PYTHIA8 [43] and EPOS-LHC [44] event generators and GEANT3 [45] transport code were used for evaluating the correlation between N ch and N corr trk . Assuming a second order polynomial function (f ) to describe the correlation between N corr trk and N ch , the average charged-particle multiplicity density in a given multiplicity interval i, is estimated as [46] where ∆η = 2 is the pseudorapidity window in which charged particles are measured. Possible deviations from this assumption are evaluated by either assuming a linear correlation between N corr trk and N ch or using a Bayesian unfolding procedure [47,48]. The values of N ch obtained using these methods are consistent within systematic uncertainties. The dN ch /dη is measured by ALICE in the same acceptance and is found to be 5.48 ± 0.05 (uncorr. syst.) ± 0.05 (corr. syst.) and 6.93 ± 0.07 (uncorr. syst.) ± 0.06 (corr. syst.) for pp collisions at √ s = 5.02 TeV and 13 TeV, respectively [49]. A summary of dN ch /dη/ dN ch /dη values for the INEL > 0 class is shown in table 1. Statistical uncertainties are negligible.

J/ψ reconstruction and signal extraction
In order to provide muon identification, each candidate track reconstructed in the muon tracking chamber (MCH) of the MS is required to match with a corresponding track segment in the muon trigger (MTR). A selection on the pseudorapidity (−4 < η < −2.5) is applied on both J/ψ daughter tracks to reject muons at the edges of the spectrometer acceptance. In addition, a selection on the distance from the z axis of the track at the end of the absorber within 17.6 to 89.5 cm is applied. The selection removes tracks which cross the high density part of the absorber, since they are significantly affected by multiple Coulomb scattering, which results in poor mass resolution of the corresponding dimuon pair. A condition on the rapidity of the dimuon pair, 2.5 < y < 4, is applied.
The J/ψ raw signal yield is obtained by fitting the opposite sign dimuon invariant mass distribution with a superposition of J/ψ and ψ(2S) signals and a function to account for the background. The estimation of the J/ψ signal and background is carried out using several fit functions. The J/ψ signal is extracted using two functions, an extended Crystal Ball function, consisting of Gaussian core and asymmetric power-law tails, and a pseudo-Gaussian with mass-dependent asymmetric tails [50]. The J/ψ peak position and width are left free in the fit, while the tail parameters are fixed to the values obtained from MC simulations anchored to the data with injected J/ψ signals. Similarly, as a part of the checks related to the systematic uncertainty of the signal extraction, tail parameters from a previous analysis at √ s = 13 TeV are also used [10]. The mass position and width of the ψ(2S) are bound to the mass and width of J/ψ using the method described in ref.
[10]. Similar approaches are considered for the signal extraction at √ s = 5.02 and 13 TeV, except the different mass ranges and background functions. The invariant mass fit is performed within [2,5] GeV/c 2 at both energies, and this range was varied in order to compute the systematic uncertainty.

At
√ s = 5.02 TeV, the background is parameterized with either a Gaussian with a width that varies linearly with the mass, called Variable Width Gaussian function, or a ratio of a first-order to a second-order polynomial function. In the case of √ s = 13 TeV, the background functions are the Variable Width Gaussian function, and the product of a fourth order polynomial and an exponential function. Examples of J/ψ signal extraction in the lowest and highest multiplicity intervals at both the collision energies are shown in figure 1. Several tests are performed by taking different combinations of signal and background functions, invariant mass ranges as well as tail parameters. The raw yields of J/ψ and the corresponding statistical uncertainties are calculated by taking the mean of the values obtained in these tests. The square root of variance (RMS) of different tests is assigned as systematic uncertainty. The dimuon invariant-mass distribution is sliced in various N corr trk intervals as discussed in section 3.1. The signal extraction in each multiplicity interval is performed in the same way as in the multiplicity-integrated case.
The J/ψ yield is extracted from the uncorrected invariant mass distribution and corrected for J/ψ acceptance and efficiency (A × ) of the MS. The A × factor is estimated based on MC simulations. A J/ψ (p T ,y) distribution corresponding to the data is used in the simulations. The generated (p T ,y) distributions are modified via an iterative process to reproduce the reconstructed ones. The MC simulation is performed enforcing J/ψ decays into two opposite-sign muons along with radiative photon emissions using EVTGEN [51] and PHOTOS [52]. The simulated decay muons are transported through the detector using GEANT3 [45].

Relative J/ψ yield measurement
The relative yield of J/ψ in a given multiplicity interval i is computed by normalizing to the corresponding multiplicity integrated value as where N J/ψ and N MB eq are the number of J/ψ reconstructed candidates and the equivalent number of MB events. N MB eq is obtained from the number of dimuon triggered events, N µµ , as N MB eq = F Norm × N µµ . The normalization factor F Norm (see section 3.5) was estimated as in previous analyses [40,46]. The A × correction for N J/ψ is taken into account for the multiplicity integrated case, (A × ) J/ψ , as well as in multiplicity bins, (A × ) J/ψ,i . The A × is observed to increase by around 4% from low to high multiplicity at both center-ofmass energies. This is because A × (p T ,y) increases with p T and that p T increase with multiplicity.
The 1/ MB and 1/ J/ψ represent correction factors applied on the number of MB selected events and number of J/ψ candidates, respectively. These factors account for possible event and signal losses due to the event selections applied at the analysis level. Each of these corrections include contributions from the minimum bias trigger efficiency for    MC by evaluating the fraction of INEL > 0 events that satisfy the MB trigger condition. The correction factors for measuring J/ψ (and MB events) that assure good vertex quality within the satisfied trigger condition were estimated from data by taking the ratio between the number of J/ψ (and MB triggered events) with and without vertex selection criteria. The effect of the vertex range (|z SPD vtx | < 10 cm) selection is the same for MB and J/ψ, hence this efficiency correction will cancel out in the ratio. The correction for the events that are rejected due to the pileup is evaluated from data by comparing N MB i / N MB between the high and low pileup rate runs in each multiplicity interval with and without pileup selection criteria. A similar technique was applied for J/ψ as well, using dimuon events. In addition, the integrated number of MB events is used to normalize the J/ψ yield, which includes events with zero number of charged particles (INEL = 0 events). This contamination is taken into account by a correction factor (1/ INEL=0 ) estimated through MC simulations. The values of all efficiency correction factors for the integrated case over multiplicity, as well as for the lowest multiplicity interval, are summarized in ) as a function of invariant mass of the dimuon (m µ + µ − ). This method allows us to extract the p J/ψ T in multiplicity intervals with a low number of reconstructed J/ψ mesons. A two dimensional map of acceptance and efficiency, A× (p T ,y) is produced using simulated events, as described in is studied as a function of the charged-particle multiplicity. A phenomenological function is used to extract the signal. It is based on weighting the p T of J/ψ entering in the spectrum and the p T of the background, by the ratio of the signal over signal plus background of each particle (J/ψ and ψ (2S)). Similar technique was previously used in refs. [40,46]. A second peak around the signal region, already observed in previous analyses [40], can be seen for the multiplicity integrated case, as shown in figure 3, which is due to the fact that the p J/ψ T is not constant with mass, rather it is a function of m µ + µ − . It is worth noticing that the dimuon invariant mass resolution depends on the precision in the measurement of the momentum of each muon and the angle between them. Several effects like scattering and energy-loss fluctuations in the muon absorber, and the precision of the tracking chambers may affect the precision of the measurements. These effects may induce a variation of the p µ + µ − T as a function of the m µ + µ − . The variation of p J/ψ T with dimuon invariant mass was quantified through MC simulations using pure signal, in particular a closure test was performed and it is found to be successful when the "piece-wise" function  is used to parametrize this variation. This function is given by, where,

Evaluation of systematic uncertainties
This section presents the summary of the assessment and implementation of the systematic uncertainties for the observables discussed in this analysis. The N ch is calculated using different MC event generators as mentioned in section 3.1. The variations of N ch from the various MC event generators are taken as a contribution to the systematic uncertainty on N ch . The dependence of N ch on the reconstructed z SPD vtx is considered to estimate another contribution to the systematic uncertainty. In particular, the correlation between N ch and N corr trk is studied in various z SPD vtx intervals and the differences are taken as systematic uncertainty for dN ch /dη. The variation of N corr trk as a function of z SPD vtx for reconstructed MC events is not identical to that in data. Hence, the N corr trk (z SPD vtx ) profile from both data and MC are used for the correction of tracklets using a data driven method. The deviation across data and MC profiles acts as another source for estimating the systematic uncertainty on N ch . All the mentioned sources are taken into account to evaluate N ch . The average over all the combinations is considered as central value, and the root mean -10 -

JHEP06(2022)015
square deviation from this value as the systematic uncertainty. The total uncertainty on N ch ranges within 0.5%−2.4% and 0.3%−2.3% for √ s = 5.02 and 13 TeV, respectively. The different contributions of systematic uncertainties on N ch are presented in table 3. The MB trigger efficiency associated to the condition INEL > 0 enters in the estimate of N ch in the lowest multiplicity interval and the systematic uncertainty assigned due to the correction is listed in table 3 for both collision energies. The systematic uncertainty of the integrated dN ch /dη is taken from ref. [49].
The systematic uncertainty due to signal extraction is determined by the RMS value of the N J/ψ i /N J/ψ tot distribution, and represents the dominant source of uncertainty on the relative J/ψ yields, as shown in table 4. The relative ratios of multiplicity dependent J/ψ yields to the multiplicity integrated one are estimated by considering different fit mass ranges, tail parameters and background functions. F Norm is computed using different methods, for the integrated multiplicity as well as in the N corr trk intervals [46]. In the first method, F Norm is obtained as the ratio of the number of MB-triggered events to the number of MB-triggered events that also satisfy the dimuon trigger condition. The second method benefits from the probability of occurrence of the single-muon triggered events by looking at the product of the probability of finding dimuon events in the single-muon triggered events and the probability of coincidence of single-muon and MB-triggered events (two-step offline method). The third method is based on the information regarding the dimuon and online MB trigger counters. Alternatively, F Norm is also obtained by re-scaling the two-step offline method with the ratio of N MB i /N MB to the N i µµ /N µµ . Run-by-run variations of F i Norm /F Norm values computed for the mentioned methods produce a maximum variation of 1.5% assigned as systematic uncertainty at √ s = 5.02 TeV, which is found to be independent of the multiplicity.
The influence of various efficiency correction factors on the relative yield of J/ψ are described in section 3.3. Different MC generators are used to evaluate the efficiency factors and the difference across MC generators is assumed as the systematic uncertainty. The corresponding systematic uncertainties for MB INEL>0 and MB vtx,QA are 1% and 2%, respectively, for pp at √ s = 5.02 TeV. However, the systematic uncertainty for the same correction factor in pp at √ s = 13 TeV can be safely ignored as the difference in obtained values between PYTHIA8 and EPOS generators is negligible. The systematic uncertainties for INEL=0 and MB INEL>0 in the lowest interval are 1% (1.5%) and 1.4 % (0.6%) at √ s = 5.02 (13) TeV, respectively. Table 4 shows the systematic uncertainty of each correction factor. The systematic uncertainty associated with pu is estimated from MC and found to be negligible at both the collision energies. The discussed systematic uncertainties are considered uncorrelated.
In the p J/ψ T analysis, the systematic uncertainty related to the signal extraction is estimated by taking the influence of various background fit functions, tail parameters, and invariant mass ranges. The uncertainty values which are determined by these tests are listed in table 5. Another source of systematic uncertainty due to A × is obtained by varying the input p T and y shapes to the MC which is used to determine the A × (p T ,y) map. Weight factors have been estimated using normalized p T and y distributions from low, high and integrated multiplicity data samples. Data were divided in sub-samples according to their multiplicity. The systematic uncertainty is evaluated by considering the -11 -   -12 -

JHEP06(2022)015
difference observed on the A × corrected p J/ψ T when re-weighted A × (p T ,y) maps are used instead of the default ones. The systematic uncertainty due to A × is found to be negligible in low multiplicity intervals, whereas 2% and 3% are taken at the highest multiplicity interval for √ s = 5.02 and 13 TeV with a conservative approach. The results are obtained with p J/ψ T as a constant parameter, and with the piece-wise function of p J/ψ T (m µ + µ − ) (see section 3.4). The difference in values obtained from both the assumptions are taken as systematic uncertainty. The uncertainty on p J/ψ T as a function of invariant mass is 0.4% (0.7%) at √ s = 5.02 (13) TeV. The systematic uncertainties due to the efficiencies of tracking, trigger, and matching of track to trigger are correlated among the multiplicity classes [46] and the multiplicity integrated values are obtained from ref.
[10]. The systematic uncertainty due to signal extraction of p J/ψ T and A× variation with event multiplicity are considered as uncorrelated, hence these affect the relative p J/ψ T measurement.

Results and discussion
The measured relative J/ψ yield at forward rapidity (2.5 < y < 4.0) as a function of dN ch /dη/ dN ch /dη at √ s = 5.02 and 13 TeV is presented in the upper panel of figure 4. At both energies, the relative J/ψ yield shows approximately linear increase with midrapidity relative charged-particle multiplicity. Also in figure 4, these results are compared with those from a previous measurement at √ s = 7 TeV [33], in which a close-to-linear trend was observed as well, albeit with significantly higher uncertainties. It is to be noted here that J/ψ yield measurements at √ s = 7 TeV are for INEL events and the effect of INEL to INEL > 0 is found to be below 1%. The similarity of the results at various collision energies suggests that in a same final state multiplicity domain, the J/ψ production at forward rapidity depends to a lesser extent of √ s. In order to assess better possible deviations from a linear behaviour, the ratio of the relative yield of J/ψ to the relative charged-particle multiplicity density is shown as a function of relative charged-particle multiplicity in the bottom panel of figure 4. The data points at the three collision energies are compatible with each other and the data points at √ s = 5.02 and 7 TeV are also compatible with unity within the uncertainties. The points at √ s = 13 TeV are above unity for dN ch /dη/ dN ch /dη ≥ 1.6. This hints at a 4.9σ departure from the linear behaviour for dN ch /dη/ dN ch /dη ≥ 3 in the case of √ s = 13 TeV, although the larger uncertainties at √ s = 5.02 and 7 TeV do not allow to exclude a similar behaviour at the lower collision energies.
In figure 5, the results are also compared with midrapidity measurements at √ s = 13 TeV. The midrapidity relative J/ψ yield exhibits faster than linear increase as a function of midrapidity relative charged-particle multiplicity. The results using midrapidity multiplicity selection based on the SPD detector (|η| < 1) and forward-rapidity multiplicity selection based on the V0 detector (−3.7 < η < −1.7 and 2.8 < η < 5.1) are found to be compatible within the uncertainties. Therefore, the different trends in the multiplicity dependence of the J/ψ production observed at midrapidity and forward rapidity are not due to a possible auto-correlation bias, arising from the multiplicity selection.
The p J/ψ T as a function of the relative charged-particle multiplicity measured in pp collisions at √ s = 5.02 and 13 TeV is shown in figure 6. At both the collision energies,  an increase of p J/ψ T with dN ch /dη/ dN ch /dη is observed with a possible saturation towards high multiplicity. A similar increase is observed in case of charged-particle p T in pp collisions [53]. Within the PYTHIA event generator, this increase is due to the Color Reconnection (CR) mechanism, which represents a fusion of hadronizing final-state strings produced in MPI [54,55]. It is worth noting that in the absence of CR, the incoherent superposition of MPI leads to a flat behaviour of p T at high multiplicities. The EPOS event generator [56,57] is also able to qualitatively reproduce the charged particle p T as a function of charged-particle multiplicity [44]. The generator is based on a combination of Gribov-Regge theory and perturbative QCD, and takes into account multiple parton interactions, non-linear effects (via saturation scales), as well as a 3D+1 viscous hydrodynamical evolution starting from flux tube initial conditions. Within EPOS, the charged-particle p T increases as a function of the multiplicity as a consequence of the col--14 -JHEP06(2022)015 lective hadronization of the high-density core produced in the collisions and the increasing relative contribution of the core to the particle production at high multiplicities. Figure 6 also shows that the absolute p J/ψ T at √ s = 13 TeV is higher than that at √ s = 5.02 TeV, as expected from the observed hardening of the corresponding transverse momentum distributions with increasing √ s [10]. However, the relative p J/ψ T , defined as the ratio of p J/ψ T in a multiplicity interval to that of minimum bias, is found to be consistent between both collision energies (figure 6). Figure 7 shows the comparison between the measured relative J/ψ yield as a function of multiplicity and theoretical calculations for √ s = 5.02 and 13 TeV at forward rapidity. The Coherent Particle Production (CPP) [58], CGC with ICEM (improved color evaporation model) [59], 3-Pomeron CGC [60], Percolation [61], EPOS3 event generator [57], and PYTHIA 8.2 (Monash 2013) [43] models are represented by the blue shaded region, red dotted line, dashed magenta line, grey shaded region, green shaded region, and red shaded region, respectively. The EPOS3 [57] and PYTHIA 8.2 [43] event generators, already introduced in section 4, predict slightly lower than linear increase of the relative J/ψ yield at forward rapidity as a function of relative multiplicity. Thus, they are able to describe the data at low multiplicity, but significantly underestimate the data at high multiplicity. It is worth noting that similar underestimation of the data is present at midrapidity [19].

Comparison with models
The CPP model [58,62] relies on the phenomenological parametrization for mean multiplicities of light hadrons and J/ψ. The parametrization is derived from p-Pb collisions assuming a linear dependence on binary nucleon-nucleon collisions (N coll ). Within the model, the nuclear suppression is found to depend weakly on energy [63] and therefore the used parametrization is independent of center-of-mass energy. An excellent agreement is found between the model and the measurements in high-multiplicity events at both collision energies.
-16 -JHEP06(2022)015   Figure 7. Comparison of the relative J/ψ yield as a function of the relative charged-particle density with model predictions: CPP [58], CGC with ICEM [59], 3-Pomeron CGC [60], Percolation [61], EPOS3 event generator [57], and PYTHIA 8.2 [43] at forward rapidity in pp collisions at The Color Glass Condensate (CGC) approach with ICEM employs the NRQCD framework in order to describe the J/ψ hadronisation process [59]. According to this model, the heavy-flavour production in high-multiplicity events in pp and pA collisions is sensitive to strongly correlated gluons in the colliding protons and nuclei. The dynamics of such configurations is controlled by semihard saturation scale Q s (x) in each of the colliding hadrons, where x is the longitudinal momentum fraction carried by a parton in the hadron. Q s increases with decreasing x and increasing nuclear size. The model predicts a faster than linear increase of J/ψ yield with multiplicity in both pp and p-Pb collisions. Although the model describes well the midrapidity J/ψ results [19], it overestimates the forward rapidity measurements, which can be clearly seen from figure 7. The 3-gluon fusion model [60] assumes that the dominant contribution comes from gluon-gluon fusion and the heavy quarks formed in the process might emit soft gluons in order to form quarkonium states. It is found that two-gluon fusion mechanism significantly underestimates the multiplicity dependence of J/ψ production [60], whereas the 3-gluon fusion model reproduces well the measured J/ψ yields at both midrapidity [19] and forward rapidity (figure 7) in pp collisions at √ s = 13 TeV.
The percolation model [61] considers high-energy hadronic collisions as driven by the exchange of color sources (strings) between the projectile and the target. The strings have finite spatial extent and thus can interact reducing their effective number. One can distinguish between soft and hard strings, depending on their transverse masses, i.e. their quark compositions and their transverse momenta. The string transverse size (r T ) is determined -17 -JHEP06(2022)015 by its transverse mass (m T ), since r T = 1/m T . The softness of the source maximizes its possibility of interaction, as its transverse size is larger. At high densities, the coherence among the sources (partons or strings) leads to a reduction of their effective number, initially proportional to the number of parton-parton interactions. This reduction mainly affects the soft observables, such as the total multiplicity, while hard production remains unaltered. At low multiplicities, with smaller number of strings the dependence of relative J/ψ yield with relative charged-particle multiplicity is linear, whereas the linear dependence changes to a quadratic dependence for high multiplicity. The number of strings, is smaller at forward rapidity compared to midrapidity and therefore the predicted rise of relative J/ψ yield with relative charged-particle multiplicity at forward rapidity is slower compared to the midrapidity one, which is in agreement with the corresponding measurements. Figure 7 shows that the percolation model explains well the data within theoretical uncertainties at both √ s = 5.02 and 13 TeV. Figure 8 shows that within uncertainties, PYTHIA8 (Monash 2013) with color reconnection gives a reasonable description of the multiplicity dependence of p J/ψ T . p T as estimated within the framework of PYTHIA8 is closer to the data in the discussed multiplicity regime, except the lowest multiplicity for √ s = 13 TeV data. Both for data and PYTHIA8, the average p T increases at a constant rate with multiplicity for dN ch /dη/ dN ch /dη > 2.5 and dN ch /dη/ dN ch /dη > 4.5 for √ s = 5.02 and 13 TeV, respectively.

Summary
The multiplicity dependence of J/ψ production at forward rapidity has been measured as a function of charged-particle multiplicity at midrapidity in INEL > 0 pp collisions at √ s = 5.02 and 13 TeV. The present measurement extends the multiplicity reach with respect to the previous measurements at √ s = 7 TeV. The results exhibit an approximately linear dependence of the J/ψ yield as a function of the event multiplicity, independent of the collision energy. All the discussed models qualitatively describe the observed trend, however the 3-Pomeron CGC model, percolation model and the CPP model quantitatively reproduce the results.
The first moment of the transverse momentum of the J/ψ at forward rapidity as a function of the charged-particle multiplicity is explored for the first time in pp collisions at √ s = 5.02 and 13 TeV. The p J/ψ T is found to increase with the event multiplicity with a hint of a saturation towards high multiplicity. The absolute p -22 -