Elucidating the multiplicity dependence of J/$\psi$ production in proton-proton collisions with PYTHIA8

A study of prompt and non-prompt J/$\psi$ production as a function of charged-particle multiplicity in inelastic proton--proton (pp) collisions at a centre-of-mass energy of $\sqrt{s}$ = 13 TeV based on calculations using the PYTHIA8 Monte Carlo is reported. Recent experimental data shows an intriguing stronger-than-linear increase of the self-normalized J/$\psi$ yield with multiplicity; several models, based on initial or final state effects, have been able to describe the observed behaviour. In this paper, the microscopic reasons for this behaviour, like the role of multiple parton interactions, colour reconnections and auto-correlations are investigated. It is observed that the stronger-than-linear increase and the transverse momentum ($p_{\rm T}$) dependence, contrary to what is predicted by the other available models, can be attributed to auto-correlation effects only. In absence of auto-correlation effects, the increase of the yield of J/$\psi$ with multiplicity -- and in general for all hard processes -- is weaker than linear for multiplicities exceeding about three times the mean multiplicity. The possibility of disentangling auto-correlation effects from other physical phenomena by measuring the charged-particle multiplicity in different pseudo-rapidity and azimuthal regions relative to the J/$\psi$ direction is investigated. In this regard, it is suggested to extend the experimental measurements of J/$\psi$ production as a function of the charged-particle multiplicity by determining the multiplicity in several azimuthal regions and in particular in the Transverse region with respect to the direction of the J/$\psi$ meson.


Introduction
Hadronic charmonium production at collider energies is a complex and not yet fully understood process. It involves partonic interactions with large momentum transfer (hard 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 charmonium state. A comprehensive description of the transverse momentum and rapidity dependent production down to transverse momentum p T =0 was recently achieved within the non-relativistic QCD (NRQCD) formalism combined with a Colour Glass Condensate (CGC) description of the incoming protons [1]; polarization is also well described [2]. The correlation of charmonium (and also of heavy quarks in general) production with the charged particle multiplicity is of high interest, as it could give new insight into the Send offprint requests to: interplay between hard and soft mechanisms in particle production, both at parton level and at hadronization.
In pp collisions at √ s = 7 TeV, ALICE has performed multiplicity dependent measurements of inclusive J/ψ production at mid and forward rapidity [3], and prompt J/ψ, non-prompt J/ψ and D-mesons production at mid-rapidity [4]. The general observation is an increase of open and hidden charm production with multiplicity. For J/ψ production, multiplicities of about 4 times the mean value observed in minimum bias events are reached. The results are consistent with a linear, or stronger than linear increase. For D-meson production, multiplicities of about 6 times the mean value are reached, with a stronger than linear increase at the highest multiplicities. Similar observations have been made by CMS for Υ (nS) mesons at mid-rapidity: a linear increase with the event activity measured at forward rapidity and a stronger-than-linear increase with the event activity measured at mid-rapidity was observed [5].
Different theoretical models attribute the observed behaviour to different underlying processes, such as the per-colation mechanism [6], higher Fock states in the proton [7], effects from the colour glass condensate EFT [8], and multi-pomeron interaction combined with high-density effects (EPOS3 event generator [9]).
In this article, we study prompt and non-prompt J/ψ production as a function of the charged-particle multiplicity in proton-proton (pp) collisions at a centre-of-mass energy of √ s = 13 TeV. The study is based on Monte Carlo (MC) event samples generated using PYTHIA8. It is worth noting that the presented observations for nonprompt J/ψ production (J/ψ from heavy flavor hadron decays) are equally valid for open heavy-flavour hadrons in general, and the findings for prompt J/ψ apply also to bottomonium production. 2 The PYTHIA event generator PYTHIA [10,11] is an event generator for collisions of protons, leptons and nuclei. It has a complex physics model with a multitude of different processes implemented at different stages of the collision. Proton-proton collisions contain one or more perturbative scattering processes between the incoming partons implemented within the MPI framework [12]. A typical event at LHC energies contains roughly between four and ten partonic interactions (PI) [13]. The number of PI per collision depends on the matter overlap in the collisions and, hence, on the pp impact parameter. The perturbative scattering processes are accompanied by Initial-State Radiation (ISR) and Final-State Radiation (FSR).
Hadronization is implemented according to the Lund string fragmentation model [14]. The created partons and the beam-remnants are connected via colour fluxtubes, or strings storing potential energy. As the partons move apart the string breaks, producing light quark-antiquark pair. The process repeats itself, until small enough pieces of strings remain, which are then identified with on-shell hadrons. In the Colour Reconnection (CR) scenario [15], strings can be rearranged between partons, so as to reduce the total string length. Partons from different PI can become connected to each other. The reduction of the total string length leads to a reduction of the total multiplicity, since the bulk of particles are produced from the string breaking mechanism.
The results reported in this paper are obtained from simulated non-diffractive events using PYTHIA version 8.230 [16] with the default Monash 2013 tune [17].

Charged-particle multiplicity
The probability distributions of the charged particle at mid-rapidity (|η| < 1) for the default settings, for CR switched-off, and for MPI switched-off are shown in Fig. 1. Multiplicity is defined as the number of primary chargedparticles, according to the ALICE definition [18]. Note that the dependence of J/ψ production on the multiplicity at mid-rapidity was investigated, as well as at forward rapidity, i.e. at the pseudo-rapidity of the V0 [19] detectors (2.8 < η < 5.1 and −3.7 < η < −1.7) of the ALICE apparatus.   The inclusion of MPI increases the average multiplicity by about a factor of 3 and the distribution becomes much wider. On the other hand the colour reconnection mechanism reduces the average multiplicity by about 30% and also makes the distribution narrower. From earlier analysis it is known that the full simulation including both MPI and CR reproduces the charged-particle multiplicity distribution measured at LHC reasonably well [20]. Table  1  Without the CR mechanism the individual PIs are independent of each other, thus the charged-particle multiplicity N ch increases roughly linearly with the number of MPI per collision, N MPI , as seen in Fig. 2 (left). The increase is observed to be somewhat weaker than linear, since the total available momentum transfer in one pp collision is fixed, such that for a higher number of MPI the individual PI most likely are softer and produce less particles. With CR switched on, the charged-particle multiplicity distribution for a given N MPI gets wider and the mean grows slower as a function of N MPI . Fig. 2 (right) shows the average N MPI as a function of N ch at mid-rapidity. Since the N MPI probability distribution is steeply falling for large N MPI upwards fluctuation of the multiplicity for a given N MPI strongly contribute to high multiplicity events. For this reason the mean N MPI grows weaker as linearly   Figure 2. Left: Mean self-normalized charged-particle multiplicity in |η| < 1 as a function of the self-normalized number of MPI for activated and deactivated colour reconnections (CR) in PYTHIA8. Right: Mean self-normalized number of MPI as a function of self-normalized charged-particle multiplicity in |η| < 1 in PYTHIA8.      with N ch . It is interesting to note that CR also reduce the increase of N MPI with multiplicity. Without the influence of fluctuations one expects the opposite behaviour and for this reason CR have been put forward as an explanation for a stronger than linear increase. In PYTHIA8, the baseline for the dependence of the yield of hard probes on multiplicity is a function (MPI-CR-baseline) which is approximately linear in the range N ch /N ch < 3 and increases weaker than linearly above this value.
Further insight into this deviation from linearity can be obtained by investigating the impact parameter dependence of MPI. As mentioned earlier, in PYTHIA the number of MPI per event is related to the matter overlap in the pp collisions and, hence, to the impact parameter b [21].

Heavy quark and quarkonium production in PYTHIA8
Heavy quark pair production is implemented in PYTHIA as perturbative scattering processes in the form of gluon fusion gg → QQ (Q denotes either charm or bottom) or light quark-antiquark annihilation qq → QQ. At the Q 2 of the hard scattering, heavy quarks can be also present in the parton distribution function leading to heavy quark production via Qg → Qg (flavour excitation). Moreover, in the parton shower heavy quarks can be produced by gluon splitting g → QQ.  For quarkonia, different production mechanisms are implemented. First, in the perturbative scattering processes, leading order NRQCD channels via colour-singlet and colour-octet pre-resonant states are included [22]. For colour-octet states, one additional gluon is emitted in the transition to the physical colour-singlet state. Secondly, quarkonia can be produced from the cluster collapse mechanism [23]. This occurs when at the hadronization stage, a heavy quark is connected to a corresponding heavy antiquark. If they are close in phase space the potential energy in the string is insufficient to create a light quarkantiquark pair, and instead the two heavy quarks bind into a quarkonium bound state. Finally, in so called nonprompt charmonium production, charmonia can be produced from the weak decay of a hadron containing a beauty quark.
In the left panel of Fig. 4 the relative contributions of the different sources of J/ψ production to the total J/ψ yield at mid-rapidity are shown as a function of p T . The largest fraction of J/ψ is produced in NRQCD processes. Non-prompt J/ψ production has a harder p T spectrum than prompt J/ψ and its relative contribution rises from about ∼10% at p T = 0 to above ∼40% at p T ∼16 GeV/c. These observations are in agreement with experiment, see e.g. [24]. The cluster collapse contribution amounts to between ∼8% and ∼4% of the total J/ψ yield, depending on p T .
For non-prompt J/ψ the origin of the initial beauty quark is investigated. This is shown in the right panel of Fig. 4, which depicts the relative contributions of the different sources of beauty quark production to the total yield of non-prompt J/ψ at mid-rapidity as a function p T . About ∼65% of all beauty quarks are produced in a primary perturbative scattering process, about ∼20% in the splitting of a gluon, which was in turn produced in a hard scattering, and the remaining ∼15% in the splitting of a gluon from initial or final-state radiation. In the latter case, the p T distribution is slightly softer than for leading order processes. The contribution of beauty quarks from beam remnants is negligible at mid-rapidity.

Results
In our MC simulations, the J/ψ was forced to decay in the dielectron channel in order to have the same conditions as in experiments that typically reconstruct J/ψ either from this or the dimuon decay channel. The role of MPI, CR and auto-correlations, are investigated in the following in order to understand the origin of the experimentally observed stronger than the MPI-CR-baseline increase and the p T dependence [25].

Multi-Parton Interactions and Colour Reconnections
In the PYTHIA model each PI has a certain probability to produce a J/ψ. Consequently the self-normalized J/ψ yield rises approximately linearly with the self-normalized number of MPI, as shown in the left panel of Fig. 5.
Here CR was deactivated, as its influence will be discussed later. In this scenario, the J/ψ yields from NRQCD vs N MPI is closest to linear. For non-prompt J/ψ and J/ψ yields from cluster collapse, the increase tends to saturate at very high multiplicities. This can be explained by the fact that the total available energy in one collision is limited. Hence, for a large number of PIs in one collision, the momentum transfer per single PI is on average smaller, so the cross section for hard processes such as heavy quark production is reduced.
The right panel of Fig. 5 shows the same dependence with CR activated. J/ψ production from NRQCD and non-prompt J/ψ are not affected by CR. This is as expected, since CR acts on the colour strings that are responsible for light particle production but does not enter in these J/ψ production processes. On the other hand J/ψ from cluster collapse shows a completely different behaviour, as the yield grows quadratically with N MPI . In this process, the charm and the anti-charm quark typically originate from independent pairs, since heavy quarks produced as a pair usually have a large opening angle. With     Figure 6. Mid-rapidity J/ψ production as a function of charged-particle multiplicity at mid (left) and forward (right) rapidity split into the different production processes implemented in PYTHIA8, i.e. non-prompt J/ψ, prompt J/ψ from NRQCD and prompt J/ψ from cluster collapse.

|<1)
CR activated, also a charm quark from one PI can bind with an anticharm quark from a different PI. The probability to produce one heavy quark pair in a collision increases linearly with N MPI : P (cc) ∝ N MPI . Since the processes are independent, the probability to produce a second one under the condition that a first one was produced also increases linearly with N MPI : P (2cc|cc) ∝ N MPI , so the total probability to produce two charm-anticharm pairs increases quadratically P (2cc) = P (2cc|cc) · P (cc) ∝ N 2 MPI . In Fig. 6 the self-normalized yields of J/ψ, for the different production mechanism, are shown as a function of the charged-particle multiplicity estimated at mid-rapidity (left panel) and forward rapidity (right panel). The J/ψ yield from NRQCD grows linearly with multiplicity, independent of whether the latter was measured at mid or forward rapidity. The J/ψ yield from cluster collapse grows stronger than linearly independently from the rapidity where the multiplicity was measured. Also the nonprompt J/ψ yield grows stronger than linearly, however, in this case the increase is stronger as a function of the mid-rapidity than the forward rapidity multiplicity. This observation hints to the possibility that for non-prompt J/ψ auto-correlation effects are important. Note that also the linear behaviour of the J/ψ yield from NRQCD lies above the MPI-CR-baseline and, hence, also in this case auto-correlation effects might have an influence. The p T dependence of the increase is shown in Fig. 7. Both as a function of mid-rapidity (left) and forward-rapidity (right) multiplicity, the increase is steeper with rising p T , in agreement with what has been experimentally measured [25]. The effect is slightly stronger for the mid-rapidity multiplicity and it is mostly due to the non-prompt J/ψ contribution, where the dependence on p T is most pronounced.

Auto-correlation effects
The charged-particle multiplicity and the J/ψ yield are not independent quantities, the latter influences the former. This influence comes from the following mechanisms:   Figure 7. Mid-rapidity J/ψ production as a function of charged-particle multiplicity at mid (left) and forward (right) rapidity in different pT intervals.
-The J/ψ decay daughters enter the charged-particle multiplicity. The simulations, in order to have the same feature of the experimental measurement are done in the dielectron decay channel. In this case two additional charged particles are produced in events containing a J/ψ, if the multiplicity is measured in the same rapidity as the J/ψ.
-In NRQCD processes the J/ψ is typically produced together with a gluon, e.g. via g g → [QQ] g which will in turn hadronize and increase the multiplicity. Additionally, if the pre-resonant state is a colour-octet, an additional gluon is emitted in the transition to the physical J/ψ state. Since the mass difference between the colour-octet and the colour-singlet state is small, the gluon will typically be emitted under a small opening angle, so the multiplicity in the flight direction of the J/ψ will be affected most.
-In the case of J/ψ from cluster collapse, the charm quark and antiquark are both produced together with another charm antiquark and charm quark, which in turn will produce additional particles.
-In the case of non-prompt J/ψ, the mother particle of the J/ψ decays into several particles and the decay daughters can decay further. Furthermore, the initially produced beauty quark can be accompanied by final state radiation, enhancing the multiplicity in the region around it. Finally, beauty quarks are always produced as pairs, mostly back-to-back in hard interactions. Thus, the non-prompt J/ψ will typically be accompanied by a high p T parton going in the opposite direction, fragmenting into a jet of particles. This recoil jet is at an azimuthal angle of 180 • with respect to the initial b quark, but can be at a different rapidity.
These auto-correlation effects can be best studied in events with only one hard interaction, which is with the MPI mechanism switched off. Then, J/ψ and charged particle production both originate from the same process and the described effects should be clearly visible.  Figure 9. Definition of the Toward, Transverse and Away region in ϕ with respect to the J/ψ direction.  Figure 10. Multiplicity dependence of J/ψ produced at mid-rapidity in PYTHIA8.230 with multiparton interactions switched off, split into the different production processes. The different panels are for the multiplicity evaluated in different kinematic regions. Left: multiplicity at mid-rapidity, right: multiplicity at forward rapidity; top row: multiplicity in Toward region in ϕ w.r.t. J/ψ, middle row: multiplicity in Transverse region, bottom row: multiplicity in Away region. Fig. 8 shows the self-normalized J/ψ yield at midrapidity as a function of the self-normalized charged-particle multiplicity at mid-rapidity for events without MPI. For non-prompt J/ψ, a strong increase of the yield with the charged-particle multiplicity can be observed, likewise but weaker for J/ψ from cluster collapse. For J/ψ from NRQCD a different picture emerges: at low multiplicity, the yield increases with multiplicity, up to around 1.5 times the mean value, afterwards it decreases again slightly with multiplicity. This behavior can be understood in the following way: at low multiplicity the additional particle production from the gluons produced alongside the J/ψ leads to an increase of J/ψ production with multiplicity. At higher multiplicity, the competition between J/ψ and charged particle production for the limited total phase space in the collision becomes relevant and leads to the observed decrease of the self-normalized J/ψ yield.
To further investigate the different auto-correlation effects, the J/ψ yield can be studied as a function of the charged-particle multiplicity in different angular regions with respect to the direction of the J/ψ, i.e. in different regions of the azimuthal angle ϕ and at different rapidities. The ϕ direction is split into three regions as indicated in Fig. 9: -Toward region: ∆ϕ ≡ |ϕ − ϕ J/ψ | < π/3 -Transverse region: π/3 < ∆ϕ < 2π/3 -Away region: 2π/3 < ∆ϕ.
Furthermore, the charged-particle multiplicity was determined either at mid-rapidity, or at forward rapidity. The dependence of mid-rapidity J/ψ production on the multiplicity in the different regions is shown in Fig. 10. In the left panels, the multiplicity was determined at midrapidity, in the right ones at forward rapidity; the top, middle, and bottom panels are for the multiplicity determined in the Toward, Transverse, and Away regions, respectively.  Transverse region Away region  18  Figure 11. Multiplicity dependence of non-prompt J/ψ in different pT intervals produced at mid-rapidity in PYTHIA8.230 with multiparton interactions switched off, for the beauty quark produced in a hard pQCD process. The panels are as in Fig. 10.
The yield of J/ψ from NRQCD does not depend strongly on the charged-particle multiplicity. It is independent of the multiplicity in the Away region, and decreases slightly with increasing multiplicity in the regions separated from the flight direction of the J/ψ. As a function of the multiplicity in the flight direction of the J/ψ, first an increase of the yield with multiplicity can be observed, which then saturates and changes into a decrease. The yield of nonprompt J/ψ increases strongly as a function of the multiplicity in the Away region, and of the multiplicity in the Toward region at mid-rapidity. It is weakly dependent on the multiplicity in the other regions. For J/ψ from cluster collapse the increase with multiplicity is qualitatively similar to the one of non-prompt J/ψ, but slightly weaker. These observations are in line with the expectations from the described effects: an auto-correlation of the multiplicity in the region around the J/ψ from the additional decay daughters, and an auto-correlation from the recoil jet, especially for non-prompt J/ψ, which is spread out in rapidity.
Since the auto-correlation effects are strongest for nonprompt J/ψ, the multiplicity dependence is further inves-tigated in transverse momentum intervals (Fig. 11). The dependence on the multiplicity in regions where no autocorrelation effects are expected, the Transverse region, and the Toward region at forward rapidity, is largely p T independent. On the other hand, in the region affected by auto-correlations, the increase with multiplicity is strongly p T dependent, i.e. the increase is much stronger for high p T . This is again in line with expectations, since higher-p T beauty quarks should fragment into more particles.
Based on the auto-correlation arguments, in principle non-prompt J/ψ production should be independent of the multiplicity in the Transverse region, and of the multiplicity in the Toward region at different rapidity. However, also in these cases a slight increase with multiplicity is observed. This fact can be further investigated by splitting the non-prompt J/ψ yield depending on the production process of the initial beauty quark. As discussed earlier, this can be a hard interaction, a gluon splitting -either of a gluon produced in a hard interaction or in ISR or FSR -or it can come from the beam remnant. The multiplicity dependence for these different cases is shown in Fig. 12. If the beauty quark is produced in a hard inter-  Figure 12. Multiplicity dependence of non-prompt J/ψ produced at mid-rapidity in PYTHIA8.230 with multiparton interactions switched off, for different production mechanisms of the beauty quark. The panels are as in Fig. 10.
action, the multiplicity dependence is as expected from auto-correlation effects only. The increase with multiplicity is present in the Toward region at the same rapidity of the J/ψ and in the Away region. In the latter case, the increase is also observed as a function of the multiplicity at forward rapidity. In these hard processes, the beauty quark and antiquark are produced back-to-back so the influence of the recoil jet is very clear, leading to the discussed auto-correlation pattern. For beauty quark-antiquark pairs from gluon splitting, where the gluon is produced in a hard process, the multiplicity dependence is also spread out into the other regions. In this case, the J/ψ yield also increases as a function of the multiplicity in the Transverse region. This is not unexpected, since the topology of the produced partons is different from a back-to-back leading order topology. The beauty quark is produced together with a beauty antiquark in the gluon splitting process with a small opening angle. The non-zero opening angle results in the slight dependence on the multiplicity in the Transverse region. Additionally, the gluon is produced back-to-back with another parton in the hard scattering process, producing the recoil jet signature. In the case of the splitting of a gluon from ISR or FSR, the increase with multiplicity is spread over all regions. This is expected, since the ISR/FSR gluon can have a large opening angle to the partons produced in the hard interaction which should produce the bulk of the multiplicity in the event. For non-prompt J/ψ from a beauty quark of the beam remnant, the multiplicity dependence is much weaker. The reason is the lower transverse momentum of this contribution.
From these observations in events without MPI it is concluded that the multiplicity dependence of J/ψ production is affected by auto-correlation effects. In the full simulation, i.e. with inclusion of the MPI mechanism, the observed auto-correlation patterns are preserved, as shown in Fig. 13. However, in the Transverse region the prompt J/ψ yield follows approximately the MPI-CR-baseline and it grows slightly stronger than linear for non-prompt J/ψ production.
A stronger-than-linear increase for non-prompt J/ψ is observed when measured as a function of the multiplicity in the Toward region (top panel), or the multiplicity in   Figure 13. Self-normalized J/ψ yield at mid-rapidity as a function of self-normalized charged-particle multiplicity at midrapidity in PYTHIA8.230. The panels are as in Fig. 10.
the Away region (bottom panel). In these cases also the stronger increase for higher p T inclusive J/ψ is observed. An interesting recent experimental observation is the production of prompt and non-prompt J/ψ inside jets at forward rapidity [26] and mid-rapidity [27] in pp collisions by the LHCb and CMS collaborations, respectively. It was found that the momentum fraction carried by prompt J/ψ inside jets is significantly lower than what is predicted from PYTHIA8. In other words prompt J/ψ are observed to be much less isolated than predicted. As a consequence, auto-correlation effects for prompt J/ψ are likely underestimated in PYTHIA8, which might explain the observed disagreement with the self-normalized yield of inclusive J/ψ as a function of multiplicity [25].
We also studied the self-normalized yield of high-p T charged hadrons as a function of the charged-particle multiplicity at mid-rapidity. These results are shown in the left panel of Fig. 14 for different p T intervals. They reflect the behaviour observed in experimental measurements and the one previously discussed for inclusive J/ψ production. A stronger than linear increase is observed and the nonlinearity rises with the p T of the hadron. Moreover, the self-normalized hadron yields have been computed as a function of the charged-particle multiplicity measured in the Transverse region as shown in the right panel of Fig. 14. As observed for the J/ψ, also for the hadrons the increase with multiplicity becomes closer to linear and the p T dependence is fully removed. However, the deviations from the MPI-CR-baseline show that auto-correlations have still an influence. These studies confirm our interpretation of the results for the J/ψ, and hence our conclusion that the stronger-than-linear increase (consequently also stronger than the MPI-CR-baseline), and its p T dependence, is fully driven by auto-correlation effects.  Figure 14. Mid-rapidity hadron production as a function of charged-particle multiplicity at mid rapidity in different pT intervals (left) and as a function of the charge-particle multiplicity in the Transverse region (right).

Conclusion
From the studies reported in this paper, it is observed that PYTHIA8 qualitatively reproduces the experimentally observed stronger-than-linear increase of J/ψ production as a function of charged-particle multiplicity. From the scenario of independent MPI one expects that multiplicity increases linearly with the number of PI. In absence of auto-correlation effects, the increase of the yield of J/ψ -and in general for all hard processes -is weaker than linear for multiplicities exceeding about three times the mean multiplicity. The trend in this region is caused by the interplay between the multiplicity fluctuations of individual PI and the steeply falling MPI probability distribution. Moreover, we show that colour reconnection, as implemented in PYTHIA8, cannot be responsible for a stronger than linear increase. A small fraction of J/ψ are formed from charm quarks produced independently. This contribution scales roughly quadratically as a function of the charged-particle multiplicity. Due to associated soft particle production in events containing non-prompt J/ψ, the multiplicity dependence is affected by auto-correlation effects leading to a strongerthan-linear increase. The auto-correlations are still present if the multiplicity is measured at a rapidity separated from the signal. However, it can be removed effectively by measuring the multiplicity in the azimuthal region transverse to the J/ψ direction. The increase with multiplicity is then approximately linear, albeit still above the MPI-CRbaseline, and independent of transverse momentum. Any additional effects such as string overlapping or CGC effects [6,7,8] are not implemented in PYTHIA, hence theoretical predictions have to take this into account by assuming that the baseline behaviour, without these effects, is not a linear increase, but already stronger than linear.
It has to be noted that the findings presented in this paper for non-prompt J/ψ production are equally valid for open heavy-flavour mesons in general, and the finding for prompt J/ψ can be applied to bottomonium production.
An experimental measurement of J/ψ production as a function of multiplicity in the Transverse region is thus highly advisable in order to disentangle between autocorrelation effects, which should vanish in this case, and true correlation effects between the hard probe and the underlying event as predicted several theoretical model calculations.