J/ψ, a charmonium bound state made of a charm and an anti-charm quark, was discovered in the 1970s and confirmed the quark model. Because the mass of charm quarks is significantly above the quantum chromodynamics (QCD) scale ΛQCD, charmonia are considered excellent probes to test perturbative quantum chromodynamics (pQCD) calculations. In recent decades, they have been studied extensively at different high-energy colliders. However, their production mechanisms, which involve multiple scales, are still not very well understood. Recently, in high-multiplicity p+p collisions at RHIC and at the LHC, a significant enhancement of J/ψ production yield has been observed, which suggests a strong contribution of multi-parton interaction (MPI). This is different from the traditional pQCD picture, where charm quark pairs are produced from a single hard scattering between partons in p+p collisions. In this work, we will report the J/ψ normalized production yield as a function of normalized charged particle multiplicity over a board range of rapidity and event multiplicity in the J/ψ→μ+μ− channel with PHENIX Run 15 p+p data at s=200 GeV. The results are compared with PYTHIA 8 simulations with the MPI option turned on and off. Finally, the outlooks of J/ψ in p+Au and Au+p collisions, along with color glass condensate (CGC) predictions and the multiplicity-dependent ψ(2S)/J/ψ ratio in p+p data, will be briefly discussed.
