Charmonium suppression in ultra-relativistic proton-proton collisions at LHC energies: A hint for QGP in small systems

Proton-proton ($pp$) collision has been considered as a baseline to study the system produced in relativistic heavy-ion (AA) collisions with the basic assumption that no thermal medium is formed in $pp$ collisions. This warrants a cautious analysis of the system produced in $pp$ collisions at relativistic energies.In this work we investigate the charmonium suppression in $pp$ collisions at $\sqrt{s} = 5.02, 7$ and $13$ TeV energies. Further, charmonium suppression has been studied for various event multiplicities and transverse momenta by including the mechanisms of color screening, gluonic dissociation, collisional damping along with regeneration due to correlated $c\bar c$ pairs. Here we obtain a net suppression of charmonia at high-multiplicity events indicating the possibility of the formation of quark-gluon plasma in $pp$ collisions.


I. INTRODUCTION
The heavy-ion programs at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) are dedicated to understanding the rich and complex nature of the phase transition at high temperature and/or density from hadrons to a deconfined state of thermal quarks and gluons. This new state of deconfined colored quarks and gluons is called quark-gluon plasma (QGP) [1][2][3]. The suppression of quarkonia [4] has been proposed as a signal of the formation of the transient phase of QGP. The production of heavy quark pairs (cc or bb) takes place at the initial stage of the collisions which can be treated perturbatively but the evolution of the quark pairs and the formation of their bound states (quarkonia) occur through non-perturbative processes [5]. A quarkonia state traversing through QGP gets dissociated to unbound quark anti-quark pair, which results in the reduction of quarkonia yields (which will otherwise remain as a bound state) referred as quarkonia suppression. The quarkonia production gets altered by the cold nuclear matter (CNM) effects (nuclear absorption and shadowing) [6]. The asymmetric pA collision can be used to disentangle the CNM effects from the QGP medium effects [7]. The pA collision has been used to serve as an important baseline for the understanding and interpretation of the heavy-ion (AA) collision data. As the CNM effects will be present both in pA and AA collisions, and also it is expected that QGP will be formed in AA but not in pA collisions, therefore, the data from pA may be used to distinguish the effects of QGP. In fact, the pA experimental data corresponding to quarkonia suppression have been explained by considering CNM effects only at various * captainriturajsingh@gmail.com † Suman.Deb@cern.ch ‡ Raghunath.Sahoo@cern.ch (Corresponding˜Author) § jane@vecc.gov.in rapidities and collision centralities. For instance, the suppression pattern of J/ψ in d−Au collisions at RHIC is well explained by the CNM effects [8].
It is pertinent to mention here that the formation of hot QGP medium in AA collisions has been predicted by considering the ratio of hadronic spectra originating from AA to pp collisions. While calculating this ratio, the spectra from pp collisions required to be appropriately scaled up by the number of nucleon-nucleon collisions in AA interaction. The deviation of the ratio from unity has been treated as a signal of medium formation in AA interactions. The data from pp collision have been used in these studies as a benchmark. The formation of QGP and its collective nature has also been studied by evaluating various flow coefficients (like directed, elliptic, triangular, etc) [9,10]. For example, the elliptic flow coefficient has been used to extract the ratio of shear viscosity (η) to entropy density (s) and the triangular flow coefficients have been used to understand the initial-state fluctuations.
It was expected that the elliptic flow for small systems formed in pp and pA collisions will be zero due to lack of thermalization in such systems. However, data from p−Pb collisions at √ s N N = 5.02, & 8. 16 TeV and for pp collisions at √ s = 7 & 13 TeV at the LHC have provided hints towards the formation of a QGP medium [11][12][13][14][15][16][17]. The search for QGP in p−Pb collisions at √ s N N = 5.02 TeV has been carried out by the theoretical model which will be called, "Unified Model of Quarkonia Suppression (UMQS)" henceforth. This model is based on the unification of the various suppression and regeneration mechanisms of quarkonia produced in ultra-relativistic collisions [18]. In the present work, UMQS model approach is employed to investigate the formation of QGP in pp collisions through charmonium suppression. The dissociation temperature (T D ) of charmonia states is much less than in bottomonia states. Therefore, in this work, we choose arXiv:2109.07967v2 [hep-ph] 20 Jun 2022 charmonia over bottomonia because the temperature of the system formed in pp collisions will not be high enough to dissociate bottomonia (T D > 650 MeV). In previous works, it is found that bottomonia hardly gets suppressed even in p−Pb collisions at √ s N N = 5.02 TeV [18]. Therefore, charmonia is a better choice to investigate QGP formation in small systems.
The dissociation of the quarkonia (J/ψ, Υ, etc.) occurs through various mechanisms [4,[19][20][21]. As mentioned, charmonium suppression in heavy-ion collisions can be caused by the cold nuclear matter (CNM) effects as well as by the Debye screening of color interaction, gluonic dissociation and collisional damping in hot QGP [4,19,[22][23][24][25][26]. As the pp collisions lack the nuclear environment, i.e. the CNM effects, therefore, the hot QGP effects are considered to study the charmonia suppression in pp collision. Along with the suppression, there is a non-zero probability that a charmonia state can be regenerated in the medium at temperature (T ), T < T D .
In the current study, we have used the UMQS model to predict the charmonia suppression in ultra-relativistic pp collisions. This model has been used to explain the bottomonia suppression in Pb−Pb and p−Pb collisions at LHC energies [18,20,27]. UMQS includes color screening, gluonic dissociation, collisional damping, and regeneration as the hot matter effects. Due to the absence of a nuclear environment in pp collisions, the cold nuclear matter effects like shadowing need not be considered here. We have used the second-order dissipative fluid dynamics [28][29][30][31] in (1+1) dimension with the assumption of boost invariance to study the propagation of the charmonia in the expanding QGP background.
The paper is organized as follows. The space-time evolution of QGP is discussed in section II. Section III is devoted for studying the charmonium kinematics and medium effects along with brief discussion of UMQS. The model prediction and comparison with experimental data are shown in section IV. Finally, the summary and the outlook are presented in section V.

II. SPACE-TIME EVOLUTION OF QGP
The QGP formed in pp and AA collisions with high internal pressure will expand in space with time. The evolution of QGP is governed by the hydrodynamics. Here we use the second order relativistic viscous hydrodynamics. For solving the hydrodynamical equations we need to specify the initial conditions and the equation of state (EoS). In absence of any method based on the first principle to estimate the initial temperature (T 0 ), the following relation [32] has been used in the present work to con-strain T 0 by data: Although this equation is derived by assuming isentropic expansion, we have used it also for the present case. We feel that the uncertainties involved in this assumption at least for small viscosity, may not be large as compared to the uncertainties involved in the other quantities like transverse area of the system (A T ), statistical degeneracy (g k ) of the QGP phase and the thermalization time (τ 0 ). In Eq. 1, C = 2π 4 45ζ(3) ≈ 3.6. We have also taken, dN ch dy ∼ = dN ch dη which is exact in the massless limit.
Here the transverse overlap area A T = πR 2 T , is computed by using IP-Glasma model, where R T is the transverse radius of the fireball [33]. This description of the transverse area is based on the impact parameter defined for pp collisions and combined with a description of the particle production based on the color glass condensate model. It is worthwhile to mention that initial thermalization time (τ 0 ) plays a key role in estimating the T 0 . While, there is no way to estimate the τ 0 from first principle but it may be reasonable to assume that thermalization time will decrease with increasing center-of-mass collision energy, i.e. τ 0 ∝ 1/ √ s [34,35]. Here we consider τ 0 = 0.3, 0.2 & 0.1 fm respectively for √ s = 5.02, 7.0, & 13 TeV for pp collisions. Corresponding to these values of τ 0 , we obtain the values of T 0 as a function of multiplicity. The results have been depicted in Fig. 1.  At high multiplicity pp collision at √ s = 13 TeV, initial temperature reaches a value ∼ 370 MeV which is comparable to the initial temperature achieved in heavy-ion collisions. However, T 0 obtained from Eq. 1, is a slowly varying function of the multiplicity (T 0 ∼ (dN/dy) 1/3 ). Therefore, we get T 0 200 MeV at low multiplicity in pp collision at √ s = 5.02 TeV. However, we find that T 0 > T c , suggesting that QGP-like medium may be produced even at low multiplicities in pp collisions at LHC energies, T c ≈ 155 MeV is predicted in the Lattice QCD calculation [36]. It is expected that in such a hot medium formed at √ s = 13 TeV, say, the J/ψ will be suppressed.
A. Variation of temperature (T ) with proper time (τ ) The temperature of the system produced in ultrarelativistic heavy-ion collisions decrease with time. The rate of the decrease of temperature can be obtained by solving the second-order hydrodynamical equation [28]. These second order equations are derived from kinetic theory using Grad's 14-moment approximation method within the framework of the Müller-Israel-Stewart equation [29][30][31]37] which in (1+1) dimension with boost invariance reads: and The φ in Eq. 3 measures the change in the shear viscosity (η) with time τ . In short, φ indicates the characteristics of the medium formed under extreme conditions in ultra-relativistic collisions. For the first order solution of the Eq.2, the φ is defined as, φ = 4η/3τ . The constants a and b are given by, and where N f = 3, is the number of flavors and α s is the strong coupling constant. We have taken α s = 0.5 here as in Ref. [28]. We have used an EoS with conformal symmetry (P = /3). The bulk viscosity is zero for such EoS. For an ideal system, the proper-time dependence of temperature is given by, However, for a viscous system, the Eqs. 2 & 3 have been solved numerically. We have obtained the initial condition for φ with the help of quasi-particle model equation of state (QPM EoS) where QGP is considered as viscous medium. In QPM EoS, entropy density, s = c + dT 3 here c = 0.00482 GeV 3 and d = 16.46 are the fit parameters [18]. By taking the lower bound of η/s = 1/(4π) from Ads/CFT, the value of φ at τ 0 is obtained as: φ 0 = 1 3π s0 τ0 , here s 0 = c + dT 3 0 . The variation of T with τ is crucial to determine the lifetime of QGP and hence for the signal of QGP. This variation is sensitive to φ 0 as demonstrated in Ref. [38]. We have taken the low value of η/s to set the intial condition for φ(τ 0 ) = φ 0 . The lower η/s will give faster cooling and hence a smaller lifetime of the thermal QGP system. Therefore, if the suppression is observed for such a conservative scenario (QGP with a small lifetime), then the realization of thermal system in pp collisions will be better assured.
In order to predict the transverse expansion effect on the medium evolution, transverse expansion is incorporated as a correction to the (1+1) dimension cooling rate. It is assumed that transverse expansion starts at a time (τ tr ) as: τ tr = τ + 0.2928(r/c s ) [18,39,40]. Here r is the transverse distance and c s is the speed of sound in QCD medium. Expansion rate in transverse direction is estimated by replacing τ tr in coupled cooling rate Eqs. 2 & 3. In Fig. 2, the variation of temperature with propertime is displayed for viscous QGP (solid blue line), viscous QGP with the inclusion of transverse expansion (green dashed line) and ideal QGP (violet dashed line) for system produced in pp collisions at √ s = 13 TeV. The generation of heat due to dissipation makes the cooling of viscous QGP slower. However, considering the transverse expansion correction to current formulation, marginally increases cooling rate compared to (1+1)-dimension ex-pansion as the lifetime of the QGP in smaller system is shorter than bigger system produced in nuclear collisions. Transverse temperature cooling rate, T SO:T r , mentioned in Fig. 2, uses initial transverse position r = 0.2 fm (say) at τ = τ 0 . Change in the viscous cooling rate after including transverse correction is not much prominent. While cooling for ideal QGP medium is large throughout medium evolution as compare to both (with and without transverse correction) viscous cases. Solution of the second-order hydrodynamics excluding transverse expansion shows that the QGP lifetime will be enhanced (3.8 fm/c) compared to the ideal fluid (1.5 fm/c). Such an increase in lifetime will be crucially important for the signal of QGP.

B. Effective Temperature
Charmonia being a massive object, does not become a part of the thermalized medium which induces a temperature gradient between the medium and charmonium. This can be measured by using the relativistic Doppler shift effect, which comes into the existence due to the velocity difference between medium particles and the charmonium. Therefore, relative velocity (v r ) between medium and charmonia is calculated to obtain an effective temperature felt by the charmonium. The velocities of the medium and charmonium are denoted by v m and v J/ψ(nl) , respectively. This relativistic Doppler shift causes an angle dependent effective temperature (T ef f ) expressed as (see Refs. [41,42] for details): where θ is the angle between v r and incoming light partons. To calculate the v r , we have taken medium velocity, v m ∼ 0.7c, and charmonium velocity v J/ψ(nl) = p T /E T , where p T is transverse momentum of charmonia and E T = p 2 T + M 2 nl is its transverse energy, M nl is the mass of the charmonium state. We have taken the average of Eq.(7) over the solid angle and obtained the average effective temperature as given by: The formulation of the current work is based on a model, UMQS recently proposed in [18]. The UMQS is described below briefly for the sake of completeness with particular emphasize on the modifications wherever required.

A. Charmonium Kinematics in Evolving QGP
The charmonia during their propagation through evolving medium get influenced by the medium in several ways. The kinematics of the charmonium (J/ψ(nl)) and cc pairs formation and recombination with the evolution of the QGP medium, can be written by using the transport equation as: The first term in the right hand side of Eq. (9), is a formation term and second one corresponds to the dissociation. Γ F,nl and Γ D,nl are the recombination and dissociation rates corresponding to the regeneration and dissociation of J/ψ(nl), respectively. V (τ ) is the volume of the medium. We assume that initially, the number of charm (N c ) and anti-charm quarks (Nc) are produced in equal numbers, N c = Nc = N cc . The Eq.(9) can be solved analytically under the assumption of N J/ψ (nl) < N cc at τ 0 [43,44]: Here N J/ψ(nl) (τ QGP , b, p T ) is the net number of charmonium formed during the evolution of QGP, τ QGP . The N J/ψ(nl) (τ 0 , b) is the number of initially produced charmonium at various centralities [69] at time τ 0 . We have calculated the N J/ψ(nl) (τ 0 , b) corresponding to the centrality bins: where, T pp (b) is the overlap function for pp obtained by using the Glauber model as demonstrated in [45]. We have obtained the number of charm and anti-charm quarks given by, . The values of σ N N J/ψ(nl) and σ N N cc , used in the calculation, are given in Table I: cross-sections at mid-rapidity [46,47]. In Eq.(10), 1 (τ QGP ) and 2 (τ ) are the decay factors for the meson due to gluonic dissociation and collisional damping in QGP (of lifetime τ QGP ) and τ is the evolution time, respectively. These factors are obtained by using the following expressions: and Here, Γ D,nl (τ, b, p T ) is the sum of collisional damping and gluonic dissociation decay rates, discussed in the following sections. The lower limit of the above integral (τ nl = γτ nl , here γ is Lorentz factor) is taken as the charmonium dilated formation time where the dissociation due to color screening becomes negligible. In the equilibrated scenario of the QGP, these dissociation factors strongly depends on the evolution of the medium.

B. Suppression Mechanisms
The charmonium floating in medium gets dissociated into its constitutions due to QGP medium effects and as a result the initial production of charmonia gets suppressed. The input parameters used in the model for calculating the charmonium suppression in the QGP medium are given in the Table II.  [19,21].
In the following sections we describe the suppression mechanisms in brief along with the regeneration process which causes reduction in the suppression.

C. Color Screening
The suppression mechanism due to color screening was first proposed by Matsui and Satz [4]. They proposed that like electric charge screening in QED plasma, the partons in the QGP medium screens the color charges. The screening prevents the formation of the cc bound states [48]. The color screening of the real part of the quark-antiquark potential is an independent suppression mechanism that dominates in the initial phase of QGP where the medium temperature is very high. Original color screening mechanisms have gone through many refinements. In this work, the pressure is parameterized in the transverse plane instead of energy density [49][50][51]. In this model it is assumed that pressure vanishes at the phase boundary, i.e. at r = R T , where R T is the transverse radius of cylindrical QGP. where, The factor p(τ 0 , 0) is obtained in the refs. [50,51]. h(r) is the radial distribution function in transverse direction and θ is the unit step function. The exponent β in the above equation depends on the energy deposition mechanism whose dependence on pressure is shown in Ref. [18].
The variation of pressure as the function of τ is given by: 3(c 2 s −1) and D = c 3 , here c 1 , c 2 , c 3 are constants and have been calculated using different boundary conditions on energy density and pressure. Other parameters appeared above are: c s is speed of sound in QGP medium, η is shear viscosity of the medium and q = c 2 s + 1. Determining the pressure profile at initial time τ = τ 0 and at screening time τ = τ s , we get: where p QGP is the pressure of the QGP phase. Putting variation of T and P with τ we have determined the radius of the screening region, r s . The screening radius defines a region where effective medium temperature of a particle is larger than the dissociation temperature (i.e. T ef f ≥ T D ). If T ef f < T D , then r s → 0 which suggests that melting of the quarkonia due to color screening would be negligible in such a situation.
The cc pairs formed inside the screening region at a point r J/ψ , may escape the region, if | r J/ψ + v J/ψ t f | > r s . Here v J/ψ = p T /E T , is charmonium velocity. The condition for escape of cc pair is expressed as: where, φ is the angle between the velocity ( v J/ψ ) and position vector ( r J/ψ ), and m is mass of a particular charmonium state.
Based on Eq. ( 19), the allowed values of the azimuthal angle, φ max (r) for survival of charmonium is expressed as: Here r is the the radial separation betweenc and c in the QGP medium.
The integration over φ max along with radial separation r gives the escape probability of cc pair from the screening region, which is defined as the survival probability. The survival probability, S J/ψ c , for a charmonium state is expressed as: T αs (20) where α s = 0.5 [48,49]. The R T depends on the impact parameter (b). We have calculated it by using the transverse overlap area A T as; R T (b) = A T /π. The value of α s is chosen in such a way that beyond the chosen value, color screening mechanism becomes almost independent with respect to change in its values.

D. Collisional Damping
Collisional damping arises due to the inherent nature of the complex potential between (cc) located inside the QCD medium. The imaginary part of the potential in the limit of t → ∞, represents the thermal decay width induced due to the low frequency gauge fields that mediate interaction between two heavy quarks [52].
The charmonium dissociation due to collisional damping is obtained by using the effective potential models. In the model under consideration, the singlet potential for cc bound state in the QGP medium is given by [19,52,53]: In Eq. (21) the first and second term in the right hand side are the string and the Coulombic terms, respectively. The third term is the imaginary part of the heavy-quark potential responsible for the collisional damping. Here σ is the string constant for cc bound state, given as σ The collisional damping, Γ damp,nl defines the charmonium decay induced due to the imaginary part of the complex potential. It is calculated using the first-order perturbation, by folding of imaginary part of the potential with the radial wave function: (22) where g nl (r) is the charmonia singlet wave function. Corresponding to different values of n and l (here n and l are the principal and the orbital quantum numbers), the wave functions can be obtained by solving the Schrödinger equation for J/ψ, χ c , ψ(2S).

E. Gluonic Dissociation
The dissociation of J/ψ by gluons contributes to the reduced yields of J/ψ in heavy-ion collisions. The the cross-section for such dissociation process is given by [19]: where J ql nl is the probability density obtained by using the singlet and octet wave functions as follows: where m c = 1.5 GeV, is the mass of the charm quark and α u s 0.59, is the coupling constant, scaled as α u s = α s (α s m 2 c /2). The E nl is the energy eigen values corresponding to the charmonium wave function, g nl (r). The octet wave function h ql (r) has been obtained by solving the Schrödinger equation with the octet potential V 8 = α ef f /8r. The value of q is determined by using the conservation of energy, q = m c (E g + E nl ).
The gluonic dissociation rate, Γ gd,nl is obtained, by taking the thermal average of the dissociation crosssection [18], where p T is the transverse momentum of the charmonium and g d = 16 is the statistical degeneracy of the gluons. Now summing the decay rates corresponding to the collisional damping and the gluonic dissociation, we have calculated the combined effect in terms of total decay width denoted by, Γ D,nl (τ, p T , b):

F. Regeneration Factor
There are two processes by which charmonia can be reproduced in the QGP medium. The first one is through the uncorrelated cc pairs present in the medium. They can recombine within the QGP medium at a later stage [43,44,[54][55][56][57][58] to form a charmonium state. This regeneration process is found to be significant for charmonia in heavy-ion collisions at LHC energies because cc are produced abundantly in the hot QGP medium. As pp is a small system, the production of cc will be small compared to the heavy-ion collisions, therefore, the probability of regeneration due to uncorrelated cc is not considered in the present work. The regeneration due to correlated cc pairs is the regeneration mechanism which is just the reverse of gluonic dissociation. In this process, cc pairs may undergo a transition from color octet state to color singlet state in the due course of time in QGP medium. To account for the regeneration via correlated cc pairs in our current UMQS model, we have considered the de-excitation of octet state to singlet state via a gluon emission. We have calculated this de-excitation in terms of recombination cross-section σ f,nl for charmonium by using the detailed balance to the gluonic dissociation cross-section σ d,nl [18]: where s is the Mandelstam variable, related with the center-of-mass energy of cc pair, given as; s = (p c + pc) 2 , where p c and pc are four momenta of c andc, respectively.
Finally we have obtained the recombination factor, Γ F,nl =< σ f,nl v rel > pc by taking the thermal average of the product of recombination cross-section and relative velocity v rel between c andc : Γ F,nl = pc,max pc,min pc,max pc,min dp c dpc p 2 c p 2 c f c fc σ f,nl v rel pc,max pc,min pc,max pc,min dp c dpc p 2 c p 2 c f c fc , where, p c and pc are the 3-momentum of charm and anti-charm quark, respectively. The f c,c is the modified Fermi-Dirac distribution function of charm, anti-charm quark and expressed as; f c,c = λ c,c /(e Ec,c/T ef f + 1). Here E c,c = p 2 c,c + m 2 c,c is the energy of the quarks and λ c,c is their respective fugacity factors [59]. We have calculated the relative velocity of cc pair in medium as follows: Since the gluonic dissociation increases with the increase in temperature, it leads to the production of significant number of cc octet states in central collision where temperature is found to be around 300 MeV. Such that the de-excitation of cc octet states to J/ψ enhance the regeneration of J/ψ in central collisions as compared with the peripheral collisions.

G. Net Yield
The survival probability of charmonium (S J/ψ gc ) in pp collisions due to gluonic dissociation along with collisional damping is given by: It is assumed here that at τ = τ 0 the color screening is the most dominating mechanism and would not allow for the charmonium to be formed. However, as QGP cools down, its effect on quarkonia suppression becomes insignificant at the later stage of the evolution. The color screening is incorporated in the model as an independent mechanism with the other suppression mechanisms in QGP. The net yield in terms of survival probability, S P (p T , b) is expressed as: The feed-down formulation of the excited states of the charmonia, χ c and ψ(2S) into J/ψ has been adopted from the Refs. [18][19][20]. To obtain the feed-down, here we have taken the ratio between the net numbers of the initial and final J/ψ. The net initial number is obtained considering the feed-down of the higher resonances into J/ψ in the absence of the QGP medium. This is given as N in J/ψ = J≥I C IJ N (J), where C IJ is the branching ratio of state J to decay into the state I. The net final number of J/ψ includes medium effects in terms of survival probability (S P (p T , b)) along with feed-down: N f i J/ψ = J≥I C IJ N (J)S P (J). The net generalized survival probability including feed-down correction can be written as; It has been observed that only 60% of J/ψ originates from direct production whereas 30% and 10% from the decays of χ c and ψ(2S) respectively. The net survival probability of the J/ψ after considering feed-down correction is given by,

IV. RESULTS AND DISCUSSIONS
The results obtained from the calculation mentioned above is contrasted with the data available in the form of the J/ψ yield corresponding to the normalized multiplicities, defined as; where N J/ψ m and N ch m are normalized by the corresponding mean values in minimum bias pp collisions.
The normalized p T -integrated J/ψ production as a function of normalized charged-particle multiplicity is contrasted with the data from ALICE at mid-rapidity (|y| < 0.9) in pp collisions to check the feasibility of the UMQS model. Left panel of Fig 3 shows [61]. It is worth mentioning here that in Ref. [61], the case of prompt J/ψ production (predicted by PYTHIA 8.2) is also included to illustrate the effect of non-prompt J/ψ in the inclusive production. Clearly, both the results of PYTHIA 8.2 reproduce the experimental data reasonably well at low N ch but fails at higher N ch as shown in the right panel of Fig 3. It is interesting to observe that UMQS prediction explains the data within the uncertainties. From Fig 3, it is further observed that at high multiplicity the production of J/ψ is larger than low multiplicity regime. In UMQS charmonium can be produced through hard scattering and recombination. The combined effect of both of these production mechanisms dominates at high multiplicities. We have found that the regeneration plays an important role at high multiplicities while its contribution at low multiplicity bins is found to be negligible. The agreement of results displayed in Fig 3 gives confidence to proceed with further analysis of quarkonia production by using UMQS in pp collisions at the LHC energies.
Further, we have predicted the charmonium suppression in ultra-relativistic pp collisions at the LHC energies. The UMQS model is used to calculate the p T and centrality dependent survival probability of charmonium states formed in pp collisions at LHC energies at the midrapidity. The relative suppression of ψ(2S) over J/ψ is measured in the form of double ratio or yield ratio, S ψ(2S) P /S J/ψ P which is plotted as a function of multiplicity (dN ch /dη) and transverse momentum, (p T ). It is to be noted that the "S", "R" and "F" stand for "suppression", "regeneration" and "feed-down", respectively. The results are shown with and without combining all the medium effects to demonstrate the significance of an individual mechanism.

A. Centrality Dependent Suppression
We have obtained the centrality dependent survival probability for J/ψ and ψ(2S) by averaging over p T with the p T -distribution function 1/E 4 T [62]. The p Tintegrated centrality dependent survival probability is given by, The initial temperature obtained for low multiplicity exceeds the value of the QCD transition temperature (T c ) obtained from lattice QCD calculation (see Fig. 1). A marginal suppression of J/ψ is obtained at low multiplicity too and it increases with the final state multiplicities.
However, ψ(2S) is effectively more suppressed at low multiplicity as compared with J/ψ. The effects of the individual mechanisms including feed-down are also shown in this section. The higher resonances of charmonium are more suppressed than J/ψ and therefore the feed-down of these states into J/ψ enhances its suppression. However, regeneration mechanism reduces the suppression depending upon the medium size and temperature. For a relatively small QGP lifetime τ QGP for 5.02 TeV, the regeneration significantly reduces the suppression.
In Fig. 4, charmonium suppression is shown as a function of the charged-particle multiplicity and the impacts of individual mechanisms are also presented. In this figure, the effective regeneration is visible at high multiplicity while at low multiplicity, the effect of regeneration is found to be very small. Here, it is also shown that feed-correction enhances the J/ψ suppression a bit. It is an indirect observation of the fact that higher resonances of charmonium are more suppressed than J/ψ. However, at √ s = 5.02 TeV, we have observed an effective regeneration for ψ(2S) which reduces its suppression significantly at high multiplicities. Inclusion of the regeneration along with feed-down correction almost cancels out the feed-down effect and suggests the possibility of regeneration for the higher resonances at high multiplicities for √ s = 5.02 TeV. Fig. 5 depicts the relative suppression of ψ(2S) over J/ψ including all the proposed medium effects. It shows, that at mid multiplicity bins ψ(2S) is more suppressed as compared to J/ψ for √ s = 5.02. The coinciding point corresponding to "S+F" and "S+R" at the lowest multiplicity bin, suggests that the difference between ψ(2S) and J/ψ suppression is almost the same At √ s = 7 TeV, J/ψ and ψ(2S) are almost equally suppressed at the high multiplicity, as shown in Fig. 6. However, ψ(2S) is found to be more suppressed at low and mid multiplicities. Regeneration effect is also found to be negligible for J/ψ and ψ(2S) at √ s = 7 TeV. An infinitesimal regeneration at very high multiplicity can be observed for J/ψ and ψ(2S) in Fig. 6. As regeneration being insignificant at √ s = 7 TeV, feed- down correction is found to take over the regeneration. Therefore, suppression along with regeneration and feed-down correction ("S+R+F") follows the foot-step of the combined effects of the suppression and feed-down ("S+F"). Non-vanishing effect of regeneration at highest multiplicity bin, reduces the suppression a bit for J/ψ and ψ(2S) and almost takes away the effect of feed-down for J/ψ. Fig. 7, where at high multiplicity ψ(2S) to J/ψ yield ratio is approximately one, while at low multiplicities the ratio is less than one. This indicates that J/ψ is less suppressed than ψ(2S) at peripheral pp collisions, which correspond to lower final state event multiplicities. Results in Fig. 7 also suggest a higher regeneration of J/ψ than ψ(2S) at high multiplicity. The results for suppression and suppression along with feed-down, and regeneration are found to coincide at the high multiplicity domain.

Similar result is displayed in
It is expected that at √ s = 13 TeV the energy deposition is large which provides substantial suppression of J/ψ and ψ(2S), comparable to the heavy-ion collisions (see Fig. 8). Results displayed in Fig. 8 show that the effect of regeneration is almost negligible at low multiplicity (as in Fig. 6) while at high multiplicities it slightly reduces the suppression for both the charmonium states. ψ(2S) is suppressed slightly more than J/ψ at high multiplicity contrary to the results shown in Fig. 6.
As a consequence, "S+F" increases the suppression at high multiplicity, however, an effective regeneration for high resonances neutralizes feed-down effect in "S+R+F" (Fig. 8). The suppression ratio plotted in Fig. 9 shows that ψ(2S) is less suppressed at low and high multiplicities as compared to the intermediate multiplicity bins. Because for higher multiplicities, the temperature is higher and the regeneration effect creates the difference, which reduces the suppression of ψ(2S) relative to J/ψ. Results depicted in Fig. 9 also show that individual effect (S) predicts a slightly more suppression of ψ(2S) than the combined effects (S+R) at the highest multiplicity bins.

B. pT -dependent Suppression
Now we discuss the transverse momentum (p T ) dependence of charmonia production by considering the p Tdependent survival probability, (S P ) for minimum bias (0 − 100% centrality) case by taking the weighted average over all centrality bins. The weighted average for S P is given by, where i = 1, 2, 3, ..., indicate various centrality bins corresponding to different multiplicity classes. The weight function W i is given as W i = bi max bi min N coll (b)π b db. The number of binary collision, N coll is calculated by using a Glauber model for pp collisions. In the Glauber model for pp collisions, a azimuthally asymmetric and inhomogeneous density distribution of a proton is considered, which is motivated by the shape of the structure function obtained in deep inelastic scattering. Therefore, same density profile is used to obtained N coll [45]. The p T -dependent suppression of charmonium is shown in Fig. 10. It shows an increasing trend for regeneration of J/ψ for p T < 20 GeV and then it slightly decreases with increasing p T . It seems that the thermodynamic conditions created at √ s = 5.02 TeV suits the regeneration of J/ψ, as indicated by the results marked by S+R for 12 < p T < 30 GeV. While including the feed-down along with regeneration (S+R+F) predicts the suppression instead of enhancement for the mentioned p T -range. However, for p T < 13 GeV both the data sets predict a reduced suppression for J/ψ. For p T < 3 GeV, J/ψ and ψ(2S) are found to be less interactive with the medium, consequently we get a bit less suppression and regeneration at this p T -range. All the medium effects predict a reduction in suppression with increasing p T . As shown in Fig. 10, a regeneration effect is observed for ψ(2S) for p T > 6 GeV. A relative suppression of ψ(2S) over J/ψ is plotted in Fig. 11, the double ratio shows that inclusion of regeneration (S+R) for ψ(2S) reduces its suppression for p T < 3 GeV. While inclusion of feed-down and exclusion of regeneration effect (S+F) increases the J/ψ suppression relative to suppression (S) with increasing p T . A dip in double ratio at p T = 6 GeV, corresponds to the suppression of ψ(2S) when regeneration is almost insignificant while there is a regeneration for J/ψ. The charmonia suppression for √ s = 7 TeV is shown in Fig. 12 as a function of p T . We find less suppression at very low p T range: 1 ≤ p T ≤ 3 GeV because the multiple scattering of the charmonium with medium constituents is less effective for very low p T . While for the intermediate p T range: 3 ≤ p T ≤ 6 GeV, the scattering with medium will be at the peak, therefore J/ψ receives maximum suppression at the intermediate p T -range. Further, suppression reduces with increase in p T , because the charmonia scattering with medium constituents will be comparably less for high-p T particles as they quickly traverse through the medium and their abundances are low. The regeneration of charmonia increases up to a certain range of the p T and at very high-p T it gets saturated as the probability of the recombination of fast moving cc reduces. The ψ(2S) feels maximum suppression for intermediate p T : 2 ≤ p T ≤ 10 GeV and gets a large suppression than J/ψ due to the low dissociation temperature (T D ) as compared to J/ψ (see Table II). The p T -dependent suppression of ψ(2S) over J/ψ corresponding to minimum bias is plotted in Fig. 13. It shows that ψ(2S) is suppressed more over the entire p T -range considered here. It depicts varying ψ(2S) suppression from intermediate to high-p T as the difference in the suppression between J/ψ and ψ(2S) reduces at high-p T regime. The results indicate that the higher resonances are more suppressed than the J/ψ and that is why the inclusion of feed-down correction (S+F) always predicts a higher suppression for J/ψ. Fig. 14 for √ s = 13 TeV show that the charmonium states are largely suppressed at low-p T and follows a trend similar to the results obtained for √ s = 7 TeV. But it shows a larger regeneration for J/ψ for p T > 5. The regeneration effect gets saturated for p T > 20 GeV. The regeneration of ψ(2S) is found to be negligible over the chosen p T -range, however, its suppression reduces at higher transverse momenta: p T > 10 GeV. Similar behavior for double ratio is also obtained in Fig. 15, where a large suppression for J/ψ reduces the relative suppression for ψ(2S), i.e. "S+F". However, the relative suppression of ψ(2S) over J/ψ comes out to be approximately constant, a slightly reduced suppression may be observed at high-p T regime (Fig. 15).  It is worthwhile to mention that in this work an interplay between regeneration and feed-down mechanisms for √ s = 5.02 to 13 TeV has been demonstrated. In UMQS, not a single parameter is freely varied in order to explain the experimental data or predict the charmonium suppression.

V. CONCLUSIONS
The present analysis shows a QGP-like medium formation for all the multiplicity bins in ultra-relativistic pp collisions at LHC energies. We have observed that ψ(2S) is suppressed more and the effect of regeneration is also marginal for ψ(2S) as compared to J/ψ. This makes ψ(2S) a relatively cleaner probe to investigate QGP-like medium in the small collision systems. As the cold nuclear environment is absent in the pp collision systems any modification in the charmonium yield can be considered as a pure effect of hot QCD medium. In contrast to this, the system formed in pA collision will show the combined effects of the hot and cold nuclear matter, therefore, pp collisions can be used as a scale to distinguish CNM effects. The pp collision results which are routinely used as a benchmark to understand the AA collision system is now become questionable because of the possibility of the creation of QGP-like medium in pp interactions. Thus, the formation of QGP in pp collisions will create a new puzzle and inevitably a new baseline will be required to analyze AA collisions.