Multiplicity dependence of inclusive J/ ψ production at midrapidity in pp collisions at

Measurements of the inclusive J/ ψ yield as a function of charged-particle pseudorapidity density d N ch / d η in pp collisions at √ s = 13 TeV with ALICE at the LHC are reported. The J/ ψ meson yield is measured at midrapidity ( | y | < 0 . 9) in the dielectron channel, for events selected based on the charged-particle multiplicity at midrapidity ( | η | < 1) and at forward rapidity ( − 3 . 7 < η < − 1 . 7 and 2 . 8 < η < 5 . 1); both observables are normalized to their corresponding averages in minimum bias events. The increase of the normalized J/ ψ yield with normalized d N ch / d η is signiﬁcantly stronger than linear and dependent on the transverse momentum. The data are compared to theoretical predictions, which describe the observed trends well, albeit not always quantitatively.


Introduction
Hadronic charmonium production at collider energies is a complex and not yet fully understood process, involving hard-scale processes, i.e. the initial heavy-quark pair production, which can be described by means of perturbative quantum chromodynamics (pQCD), as well as soft-scale processes, i.e. the subsequent binding into a color-neutral charmonium state. The latter stage is addressed via models which assume that it factorizes with respect to the perturbative early stage. The widely used non-relativistic QCD (NRQCD) effective theory [1] incorporates contributions from several hadronization mechanisms, like color-singlet or color-octet models (see Ref. [2] for a recent review on models and Ref. [3] for a comparison with data of Run 1 at the LHC). The NRQCD formalism combined with a Color Glass Condensate (CGC) description of the incoming protons [4] is a recent example of a comprehensive treatment of the transverse momentum p T and rapidity dependent production, in particular extended down to zero transverse momentum. Measurements of inclusive J/ψ production, as reported in this publication, contain a non-prompt contribution from bottomhadron decays and the production of bottom quarks can be calculated in QCD pertubatively.
The event-multiplicity dependent production of charmonium and open charm hadrons in pp and p-Pb collisions are observables having the potential to give new insights on processes at the parton level and on the interplay between the hard and soft mechanisms in particle production and is widely studied at the LHC. ALICE has studied the multiplicity dependence in pp collisions at E-mail address: alice -publications @cern .ch. √ s = 7 TeV of inclusive J/ψ production at mid-and forward rapidity [5], and prompt J/ψ (including feed down from ψ(2S) and χ c ), non-prompt J/ψ and D-meson production at midrapidity [6].
The general observation is an increase of open and hidden charm production with charged-particle multiplicity measured at midrapidity. For the J/ψ production, multiplicities of about 4 times the mean value were reached. The results are consistent with an approximately linear increase of the normalized yield as a function of the normalized multiplicity (both observables are normalized to their corresponding averages in minimum bias events). For the D-meson production, normalized event multiplicities of about 6 were reached; a stronger than linear increase of D-meson production was observed at the highest multiplicities. Observations made by the CMS Collaboration for ϒ(nS) production at midrapidity at √ s = 2.76 TeV indicate a linear increase with the event activity, when measuring it at forward rapidity, and a stronger than linear increase with the event activity measured at midrapidity [7]. At RHIC, a measurement of J/ψ production as a function of multiplicity was recently performed by the STAR Collaboration [8] for √ s = 0.2 TeV, showing similar trends as observed in the LHC data.
The J/ψ production as a function of charged-particle multiplicity was studied also in p-Pb collisions, exhibiting significant differences for different ranges of rapidity of the J/ψ meson [9,10]. A clear correlation with the event multiplicity (and event shape) was experimentally established for the inclusive charged-particle production [11] as well as for identified particles, including multistrange hyperons [12]. Several theoretical models, described briefly in Section 4, predict a correlation of the normalized J/ψ production with the normalized event multiplicity which is stronger than linear. These include a coherent particle production model [13] [14], the EPOS3 event generator [15], a CGC-complemented NRQCD model [16], the PYTHIA 8.2 event generator [17,18], and the 3-Pomeron CGC model [19]. While for instance multiparton interactions (as implemented in PYTHIA) play an important role in charm(onium) production, it is important to notice that the predicted correlation is, in all the models to first order, the result of a (N ch -dependent) reduction of the charged-particle multiplicity. Well known is the color string reconnection mechanism implemented in PYTHIA, but initial-state effects as in CGC models lead, with very different physics, similarly to a reduction in particle multiplicity.
In this Letter, the measurements of the inclusive J/ψ yield as a function of charged-particle pseudorapidity density in pp collisions at √ s = 13 TeV are presented. The measurements are performed in the dielectron channel at midrapidity with the ALICE detector at the LHC. The p T -integrated and differential results are presented for minimum bias events as well as for events triggered on high multiplicity, which extend the multiplicity range up to 7 times the average multiplicity, and on the electromagnetic calorimeter signals, which allow to access p T values up to 15-40 GeV/c. Section 2 outlines the experimental setup and the data sample; Section 3 describes the analysis, while Section 4 presents the results; a brief summary and outlook are given in Section 5.

Experiment and data sample
The reconstruction of J/ψ in the e + e − decay channel at midrapidity is performed using the ALICE central barrel detectors, described in detail in Refs. [20,21]. The setup is located in a solenoidal magnet providing a field of 0.5 T oriented along the beam direction.
For this analysis, a minimum bias (MB) trigger, a high multiplicity (HM) trigger, and two triggers based on the deposited energy in the combined Electromagnetic Calorimeter (EMCal) and the Di-jet Calorimeter (DCal) [22][23][24] are employed. Both the MB and HM triggers are provided by the V0 detector, that consists of two forward scintillator arrays [25] covering the pseudorapidity ranges −3.7 < η < −1.7 and 2.8 < η < 5.1. The MB trigger signal consists of a coincident signal in both arrays, while the HM trigger requires a signal amplitude in the V0 arrays above a threshold which corresponds to the 0.1% highest multiplicity events. The EMCal and DCal are located back-to-back in azimuth and form a two-arm electromagnetic calorimeter. While the EMCal detector covers |η| < 0.7 over an azimuthal angle of 80 • < ϕ < 187 • , the DCal covers 0.22 < |η| < 0.7 for 260 • < ϕ < 320 • and |η| < 0.7 for 320 • < ϕ < 327 • . As a consequence of identical construction, both have identical granularity and intrinsic energy resolution. In this paper, EMCal and DCal will be referred to together as EMCal. The EMCal trigger consists of the sum of energy in a sliding window of 4 × 4 towers above a given threshold (a tower is the smallest segmentation of the EMCal). In this data set, the trigger requires the presence of a cluster with a minimum energy of 9 GeV (EG1) or 4 GeV (EG2) in coincidence with the MB trigger condition.
Tracks are reconstructed in the pseudorapidity range |η| < 0.9 using the Inner Tracking System (ITS) [26], which consists of six layers of silicon detectors around the beam pipe, and the Time Projection Chamber (TPC) [27], a large cylindrical gas detector providing tracking and particle identification via specific ionization energy loss dE/dx. The first two layers of the ITS (covering |η| < 2.0 and |η| < 1.4), the Silicon Pixel Detector (SPD), are used for the charged-particle multiplicity measurement at midrapidity by counting tracklets, reconstructed from pairs of hits in the two SPD layers pointing to the collision vertex.
The results presented in this Letter are obtained using data recorded by ALICE during the LHC Run 2 data taking period for pp collisions at √ s = 13 TeV. The number of selected events and the corresponding integrated luminosities [28] are listed in Table 1 for the different triggers used in this analysis. For the analyzed data set, the maximum interaction rate was 260 kHz, and the maximum pileup probability was about 5 × 10 −3 .

Analysis
In this work the inclusive production of J/ψ mesons is studied as a function of the pseudorapidity density of charged particles at midrapidity, dN ch /dη. The J/ψ yield in a given multiplicity interval and in a given rapidity ( y) range dN J/ψ /dy is normalized to the J/ψ yield in the INEL>0 event class, dN J/ψ /dy . The INEL>0 event class contains all events with at least 1 charged particle in |η| < 1. In this ratio, most of the systematic uncertainties related to tracking and particle identification cancel.

Event selection
All events selected in this analysis are required to have a reconstructed collision vertex within the longitudinal interval |z vtx | < 10 cm in order to ensure uniform detector performance and one SPD tracklet in |η| < 1. Beam-gas events are rejected using timing cuts with the V0 detector. Pileup events are rejected using a vertex finding algorithm based on SPD tracklets [21], allowing the removal of events with 2 vertices. Because of the relatively small in-bunch pileup probability and the further event selection performed in the analysis, the fraction of remaining pileup is negligible in the minimum bias events sample and at most 2% in the high multiplicity triggered sample.
Events are binned in multiplicity classes based on either the SPD or the V0 detector signals, as shown in Fig. 1. Events corresponding to the onset of the V0 HM trigger are excluded; that onset is rather sharp. The smearing seen in the distribution in the right panel of Fig. 1 is due to the different thresholds used during operation. To illustrate this, the V0-amplitude distribution for a single data taking period is included in Fig. 1 (right panel, open  squares).
For the measurement of the charged-particle pseudorapidity density dN ch /dη at midrapidity, |η| < 1, the SPD tracklets are used [29]. Given the close proximity of the SPD detector to the interaction point (the two layers are at radial distances of 3.9 and 7.6 cm), its geometrical acceptance changes by up to 50% in the z vtx interval selected for analysis. In addition, the mean number of SPD tracklets also varied during the 3-year Run 2 data taking period due to changes in the number of active SPD detector elements. In order to compensate for these detector effects, a z vtx and time-dependent correction factor is applied such that the measured average multiplicity is equalized to a reference value. This  Table 2; the first bin spans from 0 to the position of the first line). For the HM-triggered events, the V0 amplitude distribution for a single data taking period is included for illustration (open squares).
reference was chosen to be the largest mean SPD tracklet multiplicity observed over time and z vtx . This procedure is similar to what was done previously in Ref.
[5]. The correction factor for each event is randomly smeared using a Poisson distribution to take into account event-by-event fluctuations. In the case of the event selection based on the forward multiplicity measurement with the V0 detector, the signal amplitudes are equalized to compensate for detector aging and for the small acceptance variation with the event vertex position.
The overall inefficiency, the production of secondary particles due to interactions with the detector material and particle decays lead to a difference between the number of reconstructed tracklets and the true primary charged-particle multiplicity N ch (see details in Ref. [29]). Using events simulated with the PYTHIA 8.2 event generator [30] (Monash 2013 tune, Ref. [31]), the correlation between the tracklet multiplicity (after the z vtx -correction), N corr trk , and the generated primary charged particles N ch is determined. The propagation of the simulated particles is done by GEANT 3 [32] with a full simulation of the detector response, followed by the same reconstruction procedure as for real data.
The correction factor β(N corr trk ) = N ch /N corr trk to obtain the average dN ch /dη value corresponding to a given N corr trk bin is computed from the N corr trk -N ch correlation, shown in Fig. 2 for events simulated with PYTHIA 8.2 and particle transport through GEANT 3. As the generated charged-particle multiplicity in Monte Carlo differs from data, a corrected N ch distribution is constructed from the measured N corr trk distribution using Bayesian unfolding. From it, the corrected β factors are obtained. A Monte Carlo closure test in PYTHIA 8.2 with unfolding based on EPOS-LHC events is used to validate the procedure. The normalized charged-particle pseudorapidity density in each event class is calculated as: where N corr trk is the averaged value of N corr trk in each multiplicity class, corrected for the trigger and vertex finding efficiencies. The former is estimated from Monte Carlo simulations and the latter with a data driven approach. They are below unity only for the low-multiplicity events. The value corresponding to INEL > 0 events, dN ch /dη INEL>0 , was cross-checked with the published ALICE measurement [29], and is found to be in very good agreement. A similar procedure is also used for the event selection based on the V0 amplitude, measured as a sum of signals from charged particles in the intervals −3.7 < η < −1.7 and  Table 2 Average normalized charged-particle pseudorapidity density in |η| < 1 for each event class selected in N corr trk measured in SPD (|η| < 1; left part) and in V0 amplitude (−3.7 < η < −1.7 and 2.8 < η < 5.1; right part). The values correspond to the data sample used for the p T -integrated analysis. Only systematic uncertainties are shown since the statistical ones are negligible. The corresponding fraction of the INEL>0 cross section for each event class is also indicated.

J/ψ signal extraction
The J/ψ meson is measured in the dielectron decay channel at midrapidity. Electrons and positrons are reconstructed in the central barrel detectors by requiring a minimum of 70 out of maximally 159 track points in the TPC and a value of the track fit χ 2 over the number of track points smaller than 4 [27]. Only tracks with at least two associated hits in the ITS, and one of them in the two innermost layers, are accepted. This requirement ensures both a good tracking resolution and the rejection of electrons and positrons produced from photons converting in the detector material. In the MB and HM trigger analysis, a further veto on the tracks belonging to identified photon conversion topologies is applied. The electron identification is achieved by the measurement of the specific energy loss of the track in the TPC, which is required to be compatible with that expected for electrons within 3 standard deviations. Tracks with a specific energy loss being consistent with that of the pion or proton hypothesis within 3.5 standard deviations are rejected. For the analysis of the EMCal-triggered events, the energy deposition of the track in the TPC is required to be in a range between −2.25 to +3 standard deviations around the mean expected value for the electrons. In addition, at least one of the J/ψ decay electrons is required to be matched to a cluster in the EMCal, with a cluster energy above the trigger threshold and an energy-to-momentum ratio in the range 0.8 < E/p < 1.3.
Electrons and positrons are selected in the pseudorapidity range |η| < 0.9 and in the transverse momentum range p T > 1 GeV/c.
The number of reconstructed J/ψ is obtained from the invariant mass distribution of all the opposite-sign (OS) pairs, which contains e + e − pairs from J/ψ decays as well as combinatorics and other sources. In the MB and HM trigger analysis, the combinatorial background is estimated using a track rotation procedure in which one of the tracks is rotated by a random azimuthal angle multiple times to obtain a high statistics invariant mass distribution. This distribution is then normalized such that its integral over a range of the invariant mass well above the J/ψ mass peak matches the one of real OS pairs, and is subtracted from the latter distribution. The remaining residual background, which can be attributed to physical sources, e.g. correlated semileptonic decays of heavy-quark pairs, is estimated using a second-order polynomial function. For the analysis of the EMCal-triggered events, a fit to the OS invariant mass distribution is performed using a MC shape for the signal added to a polynomial to describe the background. A second-or third-order polynomial function is used, depending on the p T range. The number of J/ψ is extracted by summing the dielectron yield in the background-subtracted invariant mass distribution in the mass interval 2.92 < m ee < 3.16 GeV/c 2 , which contains approximately 2/3 of the total reconstructed yield. The yield falling outside of the counting window at low invariant mass is due to the electron bremsstrahlung in the detector material and to the radiative J/ψ decay, and is corrected for using Monte Carlo simulations. Also, a correction for the yield loss due to the limited trigger and vertex finding efficiencies at low multiplicities is applied.
Due to the trigger enhancement, the yields obtained using the EMCal-triggered events were corrected by the trigger scaling factor, which is observed to be identical for all event classes. This correction is necessary to convert the yield per EMCal-triggered events into a yield per MB-triggered event and is accomplished by a datadriven method using the ratio of the cluster energy distribution in triggered data divided by the cluster energy distribution in minimum bias data. The ratio flattens above the trigger threshold and the scaling factor is then obtained by fitting a constant to the flat interval.
In the top panels of Fig. 3 are shown the OS invariant mass distribution for MB events (left), a high multiplicity interval from the HM-(middle) and EMCal-triggered events (right), together with the estimated background distribution. The combinatorial background distribution from the track rotation method is shown in the left and middle panels with the blue lines, while the total background is shown as black squares in all the panels. The signal obtained after background subtraction is described well by the signal shape obtained from Monte Carlo simulations (discussed below); these MC templates have been scaled and overlaid on the data points in the bottom panels of Fig. 3.
The J/ψ measurement is performed integrated in transverse momentum and in the p T intervals 0 < p T < 4 GeV/c and 4 < p T < 8 GeV/c, using the MB and HM triggers. At higher p T , the J/ψ mesons are reconstructed using the EMCal triggered events in the transverse momentum intervals 8 < p T < 15 GeV/c and 15 < p T < 40 GeV/c. It was checked that the acceptance and efficiency for J/ψ reconstruction are not dependent on the event multiplicity. This was performed using pp collisions simulated with the PYTHIA 8.2 event generator with an injected J/ψ signal. The dielectron decay is simulated with the EvtGen package [33] using PHOTOS [34] to describe the final-state radiation. The J/ψ mesons are assumed to be unpolarized consistent with measurements in pp collisions at the LHC [35].
To account for the multiplicity dependence of the p T spectrum of the J/ψ mesons, a correction for the relative efficiency, namely the efficiency for a given p T value relative to the p T -integrated value, is applied to each dielectron pair. This is contained in the invariant mass distributions shown in Fig. 3.

Systematic uncertainties
Normalized multiplicity: The systematic uncertainty on the normalized multiplicity contains contributions from the trigger, vertex finding, and SPD efficiencies. The first two contributions are estimated using alternative approaches: the trigger efficiency is calculated in a data-driven way, and for the vertex finding efficiency Monte Carlo simulations are used. The differences to the values obtained with the default methods are taken as the systematic uncertainties. The contribution from the vertex finding efficiency is below 1% (relative uncertainty) in all event classes, the one from the trigger efficiency reaches a maximum value of 1.3% for the lowest multiplicity class.
In order to estimate uncertainties due to the SPD tracklet reconstruction efficiency, the number of corrected tracklets is scaled up and down by 3%, which is the maximum observed discrepancy of the average number of SPD tracklets between data and Monte Carlo simulations. This uncertainty amounts to 3.6% in the lowest multiplicity class, and to less than 1.5% in all other classes. The uncertainty from the unfolding of the charged-particle multiplicity distribution is estimated by varying the number of iterations used in the Bayesian unfolding, as well as by using an alternative unfolding method [36]. The uncertainty is found to be negligible. All the aforementioned uncertainty sources are added in quadrature, leading to a total uncertainty on the normalized multiplicity of 3.7% for the lowest multiplicity class, and to less than 2% for all other classes.

Normalized J/ψ yield:
The systematic uncertainties of the normalized J/ψ yield are due to the signal extraction, bin-flow caused by the Poissonian smearing applied for the z vtx -dependent correction of the SPD acceptance and vertex finding, trigger and SPD efficiencies. For the analysis of the EMCal-triggered events, there is an additional component due to the matching of tracks to EMCal clusters and the electron identification via the E/p measurement, which has a non-negligible multiplicity dependence. The E/p interval and the value of E used to select only electrons above the EMCal trigger threshold are varied to determine the systematic uncertainty of the electron identification with the EMCal, leading to values from 1% to 12%, depending on the multiplicity bin.
The uncertainty of the J/ψ signal extraction is determined by varying the functions used to fit the background (first-or seconddegree polynomials or exponential) and the fitting ranges, with the RMS of the distribution of normalized yields obtained from these variations being taken as a systematic uncertainty (the normalized yield corresponds to the default selection). The bin-flow effect is estimated from the RMS of the results obtained by repeating the analysis several times with different seeds for the random number generator. The uncertainties from the signal extraction and the bin-flow effect are the dominant ones. They are of comparable size, with values between 1% and 8% depending on the multiplicity and p T interval. The uncertainties of the vertex finding, trigger and The uncertainties of the vertex finding and SPD efficiencies are below 1% in most classes, while the uncertainty due to the trigger efficiency reaches up to 4%, depending on the multiplicity class. All the mentioned sources are added in quadrature to obtain the total systematic uncertainty, which, for the p T -integrated results, varies between 3% and 7% with the multiplicity class. For the selected p T intervals, the uncertainties are larger, varying between 3% and 10% with multiplicity and p T interval, mainly due to the signal extraction, which is affected by statistical fluctuations of the background. The results at high p T , employing the EMCal, have uncertainties up to 13%, which are larger because of the additional selection requirements on the track-cluster matching and the EMCal electron identification selections.

Results and discussion
The top panel of Fig. 4 shows the normalized J/ψ yield as a function of the normalized charged-particle pseudorapidity density at midrapidity, dN ch /dη/ dN ch /dη . The dashed line also shown in the figure is a linear function with the slope of unity.
These results include both the MB and HM triggered events, which allow for a coverage of up to 7 times the average chargedparticle multiplicity, when events are selected based on the measured midrapidity multiplicity. This is a significant extension with respect to our previous results in pp collisions at where only the range up to 4 was covered and with larger uncertainties. Using the event selection based on the V0 forward multiplicity (green squares), which should largely remove a possible auto-correlation bias, the measurement extends only up to 5 times the dN ch /dη . The results for the two event selection methods are in very good agreement. In both cases, the normalized J/ψ yield grows significantly faster than linear with the normalized multiplicity. Included in Fig. 4 is also the double ratio of the normalized J/ψ yield to the normalized multiplicity (bottom panel). Two regimes could be identified, with a stronger increase of the double ratio for events with small multiplicity and a weaker increase for high-multiplicity events. It is noted that the "energy cost" for the production of a J/ψ meson, characterized by a transverse mass m T = m 2 J/ψ + p 2 T /c 2 5 GeV/c 2 , is similar to the one for particle production per unit rapidity of the underlying MB event, estimated as dN ch /dη · p T . A linear (diagonal) correlation with multiplicity is then expected to first order and observed in PYTHIA 8.2 simulations [18]. As seen in Fig. 4, the data exhibit richer features than this baseline expectation. The data in intervals of p T of the J/ψ meson are shown in Fig. 5. The data exhibit a significant increase of the normalized J/ψ yield with the normalized multiplicity between the J/ψ p T intervals 0-4 and 4-8 GeV/c. This effect could be attributed to various contributions [18], like associated J/ψ production with other hadrons in jet fragmentation or from beauty-quark fragmentation, as the fraction of J/ψ from b-hadron decays increases with p T [37].
Measurements of the correlation with the event multiplicity for inclusive charged-particle production have identified similar trends [11] as for the J/ψ p T dependence. The strength of this correlation is similar for J/ψ and for inclusive charged particles (dominated by pions) for p T values giving a comparable m T value. The production of strange hyperons at midrapidity was also observed to exhibit a correlation with event multiplicity in proportion to their mass [38]; a strong correlation was also measured for the ϒ The theoretical models currently available attribute the observed behavior to different underlying processes. In the PYTHIA 8.2 event generator [17], multiparton interactions (MPI) are an important factor in charm-quark production. Indeed, from MPIs alone a stronger than linear scaling is expected for open-charm production, as was demonstrated in Ref.
[6] with PYTHIA 8.157. Taking into account all sources of heavy-quark production, however, a close to linear increase is predicted [18]. PYTHIA 8.2 reproduces well the observation in data with a very similar correlation with multiplicity for the two different rapidity intervals used for multiplicity measurement, as seen in the left panel of Fig. 6, although the overall slope of the trend is underestimated. To illustrate the effect of non-prompt J/ψ in the inclusive production, in Fig. 6 the case of prompt J/ψ meson production as predicted by PYTHIA 8.2 is shown in addition. A significant reduction of the correlation is observed, which appears to be stronger for the SPD event selection case.
In the EPOS3 event generator [15,39], initial conditions are generated according to the parton-based Gribov-Regge formalism [40]. Sources of particle production in this framework are parton ladders, each composed of a pQCD hard process with initial-and final-state radiation. This model already predicted the stronger than linear increase with multiplicity observed for open-charm mesons [6], originating from a collective (hydrodynamical) evolution of the system. The predictions from EPOS3, here without the hydrodynamic component, are similar in magnitude to those from PYTHIA 8. In the percolation model [14], spatially extended color strings are the sources of particle production in high-energy hadronic collisions. In a high-density environment they overlap; such a decrease in the effective number of strings leads to a reduction in particle production. Since the transverse size of a string is determined by its transverse mass, lighter particles are affected in a stronger way than heavier ones. This results in a linear increase of heavy-particle production at low multiplicities, gradually changing to a quadratic one at high multiplicities. The coherent particle production (CPP) model [13,41] employs phenomenological parametrizations for light hadrons and J/ψ derived from p-Pb collisions, and predicts a stronger than linear relative increase of J/ψ production with the normalized event multiplicity. In the Color Glass Condensate (CGC) approach, the NRQCD framework is employed for J/ψ production. This effective field theory predicts, both for J/ψ and D mesons, a relative increase with the normalized multiplicity that is faster than linear, both for pp and p-Pb collisions [16]. In a CGC saturation model, a 2, the case of prompt J/ψ meson production is included for illustration. Right: comparison of data (with SPD event selection) with model predictions from the coherent particle production model [13], the percolation model [14], the EPOS3 event generator [15], the CGC model [16], the 3-Pomeron CGC model [19], and PYTHIA 8.2 predictions. Except for the latter, none of the models include the non-prompt component. Fig. 7. Normalized inclusive J/ψ yield at midrapidity as a function of normalized charged-particle pseudorapidity density at midrapidity for different p T intervals; the data are compared to theoretical model predictions from PYTHIA 8.2. faster than linear trend generically arises from the Bjorken-x dependent saturation scale which would suppress more the softparticle multiplicity, produced at low-x, compared to J/ψ production which is sensitive to larger values of x. In the 3-Pomeron fusion model [19], the correlation arises as J/ψ production via 3gluon fusion processes from various Pomeron configurations are considered. The larger configuration space for the particular case of the overlapping rapidity interval for J/ψ and charged particles leads to a significantly stronger correlation. Gluon saturation is implemented in this model; its effect, interestingly a reduced correlation, becomes significant for normalized multiplicities above 7.
All models predict an increase which is faster than linear, as shown in the right panel of Fig. 6. In all models this is effectively the result of a (N ch -dependent) reduction of the charged-particle multiplicity, realized through rather different physics mechanisms in the various approaches (color string reconnection or percolation, gluon saturation, coherent particle production, 3-gluon fusion in gluon ladders/Pomerons). The PYTHIA 8.2 and EPOS3 models underpredict the data, while the percolation model slightly overpredicts them at high multiplicity; good agreement is seen for the CGC, the coherent particle production, and the 3-Pomeron models.
These observations need to be considered having in mind that in all models except PYTHIA 8.2 only the prompt J/ψ production is included. As illustrated in Fig. 6 for PYTHIA 8.2, the prompt J/ψ meson production exhibits a weaker relative increase with multiplicity compared to the inclusive production. The agreement with data will improve in case of EPOS3 and will degrade for all the other models, in a consistent comparison. That could be realized either once the data for the prompt component will become available or as soon as the non-prompt component will be added to the current model predictions.
The contribution from decays of beauty hadrons increases significantly with p T [37] and might also have a different dependency on multiplicity; the existing measurement of charm and beauty production [6] is not precise enough to be conclusive, but a study in PYTHIA 8.2 [18] showed that the feed-down from beauty hadrons influences the result. The trend of stronger increase in the p T intervals above 4 GeV/c seen in the data is qualitatively reproduced by PYTHIA 8.2, which, however, underestimates the data for p T < 8 GeV/c, as shown in Fig. 7.

Summary and conclusions
We have presented a comprehensive measurement of inclusive production of J/ψ mesons as a function of the event multiplicity in pp collisions at √ s = 13 TeV performed with the ALICE apparatus. The J/ψ production at midrapidity is studied using a data sample including minimum bias, high event activity, and EMCal triggered events. The event selection is performed based on the charged-particle measurement at midrapidity and in the forward region. The J/ψ yield in a given multiplicity interval normalized to the J/ψ yield in INEL > 0 collisions is presented as a function of the charged-particle multiplicity similarly normalized. The advantage of such a representation is that most of the experimental systematic uncertainties cancel; also, some of the theoretical model uncertainties are mitigated for such normalized yields. A stronger than linear increase of the relative production of J/ψ as a function of multiplicity is observed for p T -integrated yields; this increase is stronger for high-p T J/ψ mesons. The trends are qualitatively, and for some of the models quantitatively, reproduced by theoretical models, but a critical appraisal of the similarity or difference between the physics mechanisms at play in various models is yet to be performed. More stringent tests of the models are needed too. Disentangling the feed-down from beauty hadrons, not included in most of the current theoretical predictions, will be an important step towards shedding light on the mechanism of hadronization of charm (and beauty) quarks, in particular in the environment of a high density of color strings created in high-multiplicity pp collisions. Data which will be collected in Run 3 at the LHC will be a significant addition for such studies.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements
We are grateful to E. Ferreiro