#  $J/\psi$ and $\psi$(2S) polarization in proton-proton collisions at
energies available at the CERN Large Hadron Collider using PYTHIA8

Bhagyarathi Sahoo1<EMAIL_ADDRESS>Dushmanta Sahu1
<EMAIL_ADDRESS>Suman Deb2<EMAIL_ADDRESS>Captain R. Singh1
<EMAIL_ADDRESS>Raghunath Sahoo1 Corresponding Author:
<EMAIL_ADDRESS>1Department of Physics, Indian Institute of Technology
Indore, Simrol, Indore 453552, India  2Laboratoire de Physique des 2 infinis
Irène Joliot-Curie, Université Paris-Saclay, CNRS-IN2P3, F-91405 Orsay, France

###### Abstract

The production mechanisms of charmonium states in both hadronic and heavy-ion
collisions hold great significance for investigating the hot and dense QCD
matter. Studying charmonium polarization in ultra-relativistic collisions can
also provide insights into the underlying production mechanisms. With this
motivation, we explore the $J/\psi$ and $\psi$(2S) polarization in
proton+proton collisions at $\sqrt{s}$ = 7, 8, and 13 TeV using a pQCD-
inspired Monte-Carlo event generator called PYTHIA8. This work considers
reconstructed quarkonia through their dimuons decay channel in the ALICE
forward rapidity acceptance range of $2.5<y_{\mu\mu}<4$. Further, we calculate
the polarization parameters $\lambda_{\theta}$, $\lambda_{\phi}$,
$\lambda_{\theta\phi}$ from the polar and azimuthal angular distributions of
the dimuons in helicity and Collins-Soper frames. This study presents a
comprehensive measurement of the polarization parameters as a function of
transverse momentum, charged-particle multiplicity, and rapidity at the LHC
energies. Our findings of charmonium polarization are in qualitative agreement
with the corresponding experimental data.

##  I Introduction

Despite being discovered nearly five decades ago, heavy quarkonia states
remain a challenging puzzle for QCD-based models Brambilla:2010cs , due to
their non-relativistic nature and the complex multi-scale dynamics involved in
ultra-relativistic heavy-ion collisions (HICs). Several theoretical models
have been developed to comprehend the quarkonium production mechanisms; the
non-relativistic quantum chromodynamics (NRQCD) is one such model
Bodwin:1994jh . The color octet NRQCD Kramer:2001hh  explains the quarkonium
production cross section and matches with the experimental data of high energy
collider experiments such as Tevatron CDF:1997ykw ; CDF:1999fhr ; D0:1996awi ,
RHIC PHENIX:2019ihw ; STAR:2019vkt ; Trzeciak:2014cma , and the LHC
ALICE:2017leg ; ALICE:2019pid ; LHCb:2015foc ; CMS:2017dju ; ATLAS:2015zdw ;
ATLAS:2012lmu ; LHCb:2018yzj ; ALICE:2015pgg . In comparison, the color
singlet model (CSM) of NRQCD Baier:1981uk ; Berger:1980ni ; Chang:1979nn ;
Einhorn:1975ua ; Kartvelishvili:1978id  predicts $J/\psi$ and $\psi$(2S)
production cross-sections 50 times smaller than the experimental observation
at CDF collaboration in proton+proton ($pp$) collisions at $\sqrt{s}$ = 1.8
TeV CDF:1997ykw . There are other various factorization techniques to predict
the production cross-section of quarkonium, such as NRQCD factorization
Bodwin:1994jh ; Bodwin:2013nua , leading power fragmentation, next-to-leading-
power fragmentation, the color singlet model, the color evaporation model
(CEM) Gluck:1977zm ; Fritzsch:1977ay ; Halzen:1977rs , and the
$k_{T}$-factorization approach Baranov:2016clx ; Braaten:2014ata . The NRQCD
calculations use color-octet matrix elements to account for the non-
perturbative long-distance physics in heavy quarkonium systems. These matrix
elements are adjusted to explain experimental data, improving the agreement
between theory and experiment. However, this adjustment introduces model
dependence and should be interpreted cautiously. Including next-to-leading-
order (NLO) QCD corrections in CSM, the quarkonium production rates have shown
a significant increase in the large transverse momentum ($p_{\rm T}$) region,
i.e., $p_{\rm T}>$ 20 GeV Lansberg:2008gk . This increase has led to a notable
reduction in the required contributions from color octet components to match
the measured quarkonium production cross-section at the Tevatron
Lansberg:2008gk .

However, for a comprehensive understanding of the quarkonium production
mechanism, it is crucial to explore the dynamics responsible for its
polarization. Studying the quarkonium polarization provides valuable
information, e.g., quarkonium production mechanisms in $pp$ collisions, the
effect of the deconfined medium on the formation of a bound state of two heavy
quarks, and the role of spin-vorticity coupling in a thermal rotating medium
Faccioli:2022peq ; DeMoura:2023jzz , etc. Polarization refers to the alignment
of the quarkonium spin with respect to a chosen axis in a reference frame. The
details about the chosen reference frames are discussed in Sec. II.

Quarkonium polarization is predominantly investigated through the dilepton
decay channel in experimental studies. The polarization of quarkonium states
is obtained by analyzing the angular distributions of decay products. So far,
from the experimental side, quarkonium polarization is observed at Tevatron
CDF:2000pfk ; CDF:2001fdy ; D0:2008yos ; CDF:2007msx ; CDF:2011ag , RHIC
PHENIX:2009ghc ; STAR:2013iae ; STAR:2020igu , and LHC ALICE:2011gej ;
ALICE:2018crw ; ALICE:2020iev ; ALICE:2023jad ; Etzion:2009zz ; LHCb:2013izl ;
LHCb:2014brf ; CMS:2013gbz ; CMS:2016xpm ; CMS:2012bpf  at collider
experiments, as well as at fixed target experiments such as E866 (NuSea)
NuSea:2003fkm ; NuSea:2000vgl  and HERA-B HERA-B:2009iab . From a theoretical
standpoint, the polarization study has been discussed in Ref. Gong:2012ug ;
Chao:2012iv ; Ma:2018qvc ; Butenschoen:2012px ; Beneke:1998re ;
Haberzettl:2007kj ; Shao:2014yta ; Faccioli:2010kd ; Shao:2012fs ;
Faccioli:2010ji ; Faccioli:2010ej ; Faccioli:2008dx ; Beneke:1996yw ;
Braaten:1999qk ; Kniehl:2000nn ; Leibovich:1996pa ; Gong:2008sn . The color-
octet model of NRQCD successfully explains the quarkonium production cross
section but fails to account for the polarization results of $J/\psi$ obtained
by the CDF experiment at $\sqrt{s}$ = 1.96 TeV CDF:2007msx . It predicts that
at very high momenta, quarkonia are produced from the fragmentation of gluons,
preserving their natural spin alignment Braaten:1993rw ; Faccioli:2010kd .
Therefore a large transverse polarization of $J/\psi$ is estimated with
respect to their momentum direction Faccioli:2014cqa . In addition, the
leading order (LO) calculation in CSM of NRQCD predicts a strong transverse
polarization for $J/\psi$ Gong:2008sn , which also fails to explain the
polarization data of $J/\psi$ at the LHC, RHIC and Tevatron Butenschoen:2012px
; Chung:2009xr . While the inclusion of NLO calculation to the CSM predicts a
strong longitudinal polarization Gong:2008sn ; LHCb:2013izl ; LHCb:2014brf ;
Lansberg:2008gk . The effective theory NRQCD includes all the important
contributions from both color-octet and color singlet intermediate states and
gives a detailed overview of $J/\psi$ polarization Butenschoen:2012px ;
Gong:2012ug ; Chao:2012iv . Furthermore, the quarkonium polarization
parameters are estimated in $pp$, $p-A$, and $A-A$ collisions using an
improved CEM employing the $k_{T}$ factorization Cheung:2018tvq ;
Cheung:2018upe  and the collinear factorization Cheung:2021epq ;
Cheung:2022nnq  approach. It predicts either zero or slightly transverse
polarization at high-$p_{\rm{T}}$ and a small longitudinal polarization at
low-$p_{\rm{T}}$ depending on the polarization frame Cheung:2021epq ;
Cheung:2022nnq ; Vogt:2019zmr . On the other hand, with the current
statistics, the experimental data of ALICE show a little or zero polarization
for $J/\psi$ within uncertainty in both hadronic and nucleus-nucleus
collisions, although LHCb predicts a small longitudinal polarization in the
helicity frame. This discrepancy between theory and experiment is commonly
called the “J/$\psi$ polarization puzzle” ALICE:2011gej ; Ma:2018qvc .
However, a Color Glass Condensate (CGC)+NRQCD approach Ma:2018qvc  provides a
qualitative description of the experimental data towards high $p_{\rm T}$ (
$p_{\rm T}>4$ GeV)at the LHC and STAR in minimum bias $pp$ collisions at
forward rapidities STAR:2020igu .

In this work, we attempt to understand the charmonia polarization in $pp$
collisions by studying the $\lambda$-polarization parameters;
$\lambda_{\theta}$, $\lambda_{\phi}$, and $\lambda_{\theta\phi}$. These
parameters are obtained using PYTHIA8 simulation by taking the angular
distribution of dimuons produced from $J/\psi$ and $\psi$(2S). Further, these
$\lambda$-polarization parameters are studied as functions of $p_{\rm{T}}$,
charged-particle multiplicity ($N_{ch}$), and rapidity ($y_{\mu\mu}$)
corresponding to $\sqrt{s}$ = 7, 8, and 13 TeV collision energies. The
experimental measurement of quarkonium decay angular distribution is
challenging because it demands a large number of event samples and a high
level of accuracy in the subtraction of various kinematic correlations induced
by detector acceptance. Thus, it is difficult to analyze the charged particle
multiplicity and rapidity dependence of polarization parameters from the
angular distribution of decay muons. Such difficulties can be easily overcome
in the Monte-Carlo (MC) simulation studies. With the present understanding of
$pp$ collisions dynamics at the LHC energies, the charged-particle
multiplicity dependence study of quarkonium polarization would be an
interesting topic to investigate the medium effect, and it may serve as a
benchmark for heavy-ion collisions as well. The $p_{\rm{T}}$-dependence of the
polarization parameter may help us to understand the dynamics of the particle
production. Similarly, the rapidity dependence of polarization studies may
reveal the phase-space analysis of the particles produced in ultra-
relativistic collisions.

This paper is organized as follows. The brief details of dimuon angular
distribution and event generation are described in Sec. II. Section III
includes the results obtained by analyzing the angular distribution of dimuons
in $pp$ collisions at the LHC energies. Section III consists of three
subsections. The transverse momentum dependence of polarization parameters is
discussed in Sec. III.1. In Sec. III.2, the charged-particle multiplicity
dependence of polarization parameters is studied. The rapidity dependence
study of polarization parameters is discussed in Sec. III.3. Finally, the
important findings are summarized, and the future scopes are presented in Sec.
IV.

Figure 1: (Color online) Illustration of the three different definitions of
the polarization axis, z, in the helicity (HF), Collins-Soper (CS), Gottfried-
Jackson (GJ) reference frames, with respect to the direction of motion of the
colliding beams (Beam 1 and Beam 2) and of the Quarkonia.

##  II  Analysis Methodology and Event generation

###  II.1 Dimuon decay angular distribution

The polarization of the $J^{PC}=1^{--}$ quarkonium states can be measured
through the study of the angular distribution W($\theta,\phi$) of decay
daughters in the dimuons decay channel and can be parameterized as
Faccioli:2010kd ; ALICE:2020iev ; Braaten:2014ata ;

$W(\theta,\phi)\propto\frac{1}{3+\lambda_{\theta}}\left(1+\lambda_{\theta}\cos^{2}\theta+\lambda_{\phi}\sin^{2}\theta\cos
2\phi+\lambda_{\theta\phi}\sin 2\theta\cos\phi\right)$ (1)

where, $\theta$ and $\phi$ are the polar and azimuthal angles of the $\mu^{+}$
with respect to the spin-quantization axis (say, z-axis) of the chosen
polarization frames and $\lambda_{\theta}$, $\lambda_{\phi}$,
$\lambda_{\theta\phi}$ are the polarization parameters. In particular, the two
cases ($\lambda_{\theta}$ = 1, $\lambda_{\phi}$ =0 , $\lambda_{\theta\phi}$ =
0 ) and ($\lambda_{\theta}$ = -1, $\lambda_{\phi}$ =0, $\lambda_{\theta\phi}$
= 0 ) correspond to the transverse and longitudinal polarization,
respectively. The case ($\lambda_{\theta}$ = 0, $\lambda_{\phi}$ =0,
$\lambda_{\theta\phi}$ = 0) correspond to zero polarization ALICE:2020iev .

There are three different conventions to define the polarization reference
frames (definition of the z-axis), which are illustrated in Fig. 1
Faccioli:2010kd .

  * •

Helicity frame: In the direction of $J/\psi$ (or $\psi$(2S)) momentum in the
center of the mass frame of the colliding beams.

  * •

Collins-Soper frame: The bisector of the angle between the momentum of one
beam and the successive direction of the other beam Collins:1977iv .

  * •

Gottfried-Jackson frame: The direction of the momentum of one of the colliding
beams Gottfried:1964nx .

The quantization axis of the Gottfried-Jackson frame lies in between the
helicity and Collins-Soper reference frame Faccioli:2010kd ; Sarkar:2010zza .
So, in the current work, we solely cover the helicity and Collins-Soper frames
as two extreme cases that are physically relevant. It is noteworthy to mention
that the default setting for quarkonium production in hadronic collisions in
all existing MC generators use an isotropic dilepton distribution, which is
discussed in Ref. Faccioli:2010kd . In NA38, NA50, NA51, and NA60 experiments
have measured a flat $\cos\theta_{\rm CS}$ ($\cos\theta$ distribution in
Collins-Soper frame) angular distribution in the window $|\rm\cos\theta_{\rm
CS}|<$ 0.5 covering 50 % of the phase space and assumed that $J/\psi$ is
unpolarized Sarkar:2010zza . But, the recent global analysis of $J/\psi$
polarization measurement indicates that the $J/\psi$ is significantly
polarized and its polarization changes longitudinal to transverse from low
$p_{\rm{T}}$ to high $p_{\rm{T}}$ ALICE:2011gej . We use PYTHIA8 to visualize
the angular distribution of decay muons in Fig. 2, in order to validate the
mentioned assumption. From Fig. 2, it is observed that the $\cos\theta$
angular distribution is almost isotropic in the angular range
$|\rm\cos\theta|<$ 0.5, however, the deviation from isotropic distribution
starts at $|\rm\cos\theta|>$ 0.5 for both frames of reference. This is the
consequence of physics processes involved in PYTHIA8, such as the production
and decay of higher excited resonances, and the emission of gluons in the
final state radiations, etc. manual . The distribution of dimuons with uniform
acceptance and efficiency over $\cos\theta$ and $\phi$ distribution at the
generation level allows us to determine the polarization parameters from the
observed angular distribution Faccioli:2010kd ; Faccioli:2022peq .

$<\cos^{2}\theta>=\frac{1+\frac{3}{5}\lambda_{\theta}}{3+\lambda_{\theta}}$
(2)

$<\cos 2\phi>=\frac{\lambda_{\phi}}{3+\lambda_{\theta}}$ (3)

$<\sin
2\theta\cos\phi>=\frac{4}{5}\frac{\lambda_{\theta\phi}}{3+\lambda_{\theta}}$
(4)

The polarization parameters $\lambda_{\theta}$, $\lambda_{\phi}$, and
$\lambda_{\theta\phi}$ are obtained for helicity and Collins-Soper frames by
taking the average over $\cos^{2}\theta$, $\cos 2\phi$, and $\sin
2\theta\cos\phi$. This approach provides an alternative method over the
multiparameter fit of Eq. (1) to the dimuon angular distribution
Faccioli:2010kd .

###  II.2 Event Generation in PYTHIA8

For modeling ultra-relativistic collisions between particles such as electron-
electron, electron-positron, proton-proton, and proton-antiproton, one of the
commonly used event generators is PYTHIA8. It is quite effective at explaining
the LHC results Deb:2018qsl ; Thakur:2017kpv . Numerous physical processes are
involved in PYTHIA8, including hard and soft scattering, parton distributions,
initial and final state parton showers, multi-partonic interaction (MPI),
string fragmentation, color reconnection, resonance decays, rescattering, and
beam remnants Sjostrand ; manual . In this study, we have used PYTHIA8 to
generate $pp$ collisions at $\sqrt{s}$ = 7, 8 and 13 TeV with 4C Tune (Tune:pp
= 5) Corke . One of the key benefits of PYTHIA8 is the subsequent MPI
processes, which, combined with impact parameter dependence of collisions,
enables the generation of heavy-flavor quarks through 2 $\rightarrow$ 2 hard
subprocesses. A detailed explanation of all physics processes involved in
PYTHIA8 can be found in Ref. manual .

This analysis is performed by generating 1.5 billion events for $pp$
collisions at $\sqrt{s}$ = 7, 8, and 13 TeV. For our study, we contemplate
inelastic and non-diffractive simulated events. So in the total scattering
cross section, only the non-diffractive component of all hard QCD processes
(HardQCD:all = on) will contribute. Hard processes involve the production of
heavy quarks. We have considered color reconnection on
(ColourReconnection:mode = on) along with MPI (PartonLevel:MPI = on). To avoid
the divergences of QCD processes in the limit $p_{\rm{T}}\rightarrow 0$ a
transverse momentum cut of 0.5 GeV (PhaseSpace:pTHatMinDiverge = 0.5) is
taken. For the production of $J/\psi$ and $\psi$(2S), we use Charmonium:all
flag (Charmonium:all = on) in the simulation Shao ; Caswell ; Bodwin through
NRQCD framework. The polarization study of $J/\psi$ and $\psi$(2S) has been
performed in the dimuon channel by forcing a $J/\psi$ and $\psi$(2S) to decay
into dimuons ($\mu^{+}\mu^{-}$) in the MC simulation. The $J/\psi$ and
$\psi$(2S) yields are then obtained through invariant mass reconstruction
considering the detector acceptance. This helps in comparing the observations
directly with the experimental data.

To check the compatibility of PYTHIA8 with experimental data, we have used the
same tuning as used in our previous works Deb:2018qsl ; Thakur:2017kpv , where
we have compared the production cross section obtained from PYTHIA8 as a
function of transverse momentum and rapidity with the ALICE experimental data
of $J/\psi$, and found them to be comparable within uncertainties.

##  III Results and Discussion

Figure 2: (Color online) The cosine of the polar angle (upper panel) and the
azimuthal angle (lower panel) distribution in $pp$ collisions for $J/\psi$ and
$\psi$(2S) at $\sqrt{s}$ = 13 TeV in helicity (left panel) and Collins-Soper
(right panel) reference frame. Figure 3: (Color online) The $J/\psi$ and
$\psi$(2S) polarization parameters as a function of transverse momentum for
$pp$ collisions at $\sqrt{s}$ = 7, 8, and 13 TeV using PYTHIA8. The obtained
results are compared with the $J/\psi$ polarization measurement in pp
collisions from ALICE at $\sqrt{s}$ = 7, 8 TeV, LHCb at $\sqrt{s}$ = 7 TeV and
the NLO-NRQCD model predictions in color singlet (CS) and color singlet +
color octet (CS+CO) states at $\sqrt{s}$ = 7 TeV in both helicity and Collins-
Soper reference frames. Figure 4: (Color online) Polarization of $J/\psi$
(left) and $\psi$(2S) (right) as a function of transverse momentum obtained
from PYTHIA8 compared with the LHCb experimental data in $pp$ collisions at
$\sqrt{s}$ = 7 TeV in the rapidity interval $2.5<y<4.0$ in helicity frame. The
MC simulation results are compared with NLO CSM from Butenschoen:2012px  and
NLO NRQCD calculation from (1) Butenschoen:2012px , (2) Gong:2012ug , and (3)
Chao:2012iv ; Shao:2012fs .

In this section, we discuss the polar and azimuthal angular distribution of
dimuons (corresponding to $J/\psi$ and $\psi$(2S)) obtained from the PYTHIA8
simulation of $pp$ collisions at $\sqrt{s}$ = 7, 8, 13 TeV. The parameters
$\lambda_{\theta}$, $\lambda_{\phi}$ and $\lambda_{\theta\phi}$ are obtained
from the averaged angular distribution of dimuons using Eq. (2), Eq. (3), and
Eq. (4), respectively. The left and right columns of Fig. 2 correspond to the
helicity and Collins-Soper frames, respectively. In both the frames, the top
panel shows the cosine of the polar angle distribution and the bottom panel
represents the azimuthal angular distribution for $J/\psi$ and $\psi$(2S) in
$pp$ collisions at $\sqrt{s}$ = 13 TeV. To obtain these angular distributions,
we have considered a $p_{\rm T}$ bin of $2.0<p_{\rm T}<3.0$ GeV/c and a
rapidity window $2.5<y_{\mu\mu}<4.0$. Further, the present section is divided
into three subsections. The $p_{\rm T}$, charged particle multiplicity, and
rapidity dependence of $\lambda$-polarization parameters (i.e.,
$\lambda_{\theta}$, $\lambda_{\phi}$, and $\lambda_{\theta\phi}$ ) are
discussed in consecutive Secs. III.1, III.2, and III.3 respectively.

Figure 5: (Color online) The $J/\psi$ and $\psi$(2S) polarization parameters
as a function of charged-particle multiplicity for $pp$ collisions at
$\sqrt{s}$ = 7, 8, and 13 TeV using PYTHIA8 in both helicity and Collins-Soper
reference frames.

###  III.1 Transverse momentum dependence of $\lambda_{\theta}$,
$\lambda_{\phi}$, $\lambda_{\theta\phi}$

We explore the polarization parameters of $J/\psi$ and $\psi$(2S) as a
function of $p_{\rm{T}}$ in $pp$ collisions using PYTHIA8. The variation of
these parameters at $\sqrt{s}$ = 7, 8, and 13 TeV in both helicity and
Collins-Soper reference frame are shown in Fig. 3. The $J/\psi$ polarization
parameters obtained using PYTHIA8 are compared with the corresponding
experimental data for $\sqrt{s}$ = 7 TeV at LHCb LHCb:2013izl  and $\sqrt{s}$
= 7 and 8 TeV at ALICE ALICE:2011gej ; ALICE:2018crw . In addition, the
obtained results are compared with the color singlet (CS) and color singlet +
color octet (CS+CO) mechanism-based NRQCD, which include the NLO corrections
Butenschoen:2012px . The $p_{\rm T}$-interval in PYTHIA8 is chosen in
accordance with the ALICE measurement of $J/\psi$ at $\sqrt{s}$ = 8 TeV in
both reference frames for all energies. The rapidity cut is set to
$2.5<y_{\mu\mu}<4.0$ in accordance with ALICE detector acceptance. From Fig.
3, it is observed that the $\lambda_{\theta}$ parameter indicates a
longitudinal polarization at low-$p_{\rm{T}}$ regime, and a transverse
polarization at high-$p_{\rm{T}}$, in both the frames of references for
$J/\psi$ and $\psi$(2S). This trend qualitatively agrees with the $J/\psi$
polarization measured by ALICE for $\sqrt{s}$ = 7 TeV ALICE:2018crw  in the
helicity frame. In the obtained result, longitudinal polarization is observed
at low-$p_{\rm{T}}$, which decreases towards high-$p_{\rm{T}}$. At lower
energy, the HERA-B experiment HERA-B:2009iab  predicts a longitudinal
polarization at low-$p_{\rm{T}}$ in the Collins-Soper frame. However, the LO
calculation of the NRQCD approach predicts a sizable transverse polarization
at high $p_{\rm T}$ Beneke:1996yw ; Braaten:1999qk ; Kniehl:2000nn ;
Leibovich:1996pa .

Further, it is observed that at low $p_{\rm{T}}$, $\psi$(2S) has a
comparatively lower longitudinal polarization than $J/\psi$ while the
polarization of $\psi$(2S) increases at high $p_{\rm{T}}$. This result seems
to be apparent because, at low-$p_{\rm{T}}$, the formation of $J/\psi$ through
$c-\bar{c}$ is more favourable than $\psi$(2S). As a consequence of this,
$\psi$(2S) yield reduces at low $p_{\rm{T}}$, which affects its polarization.
This study shows no clear dependence on the center of mass energy on quarkonia
polarization in PYTHIA8 for both reference frames. Next, the $\lambda_{\phi}$
parameter indicates a transverse polarization in the helicity and a
longitudinal polarization in the Collins-Soper reference frame. This
dissimilarity arises due to differences in the azimuthal angle distribution of
dimuons around the chosen reference axis in both frames. The finite value of
$\lambda_{\theta}$ and $\lambda_{\phi}$ parameters indicates that the
probability of finding the $J/\psi$ (or $\psi$(2S)) vector mesons in the three
spin states are not equal and hence the emission of their daughter particles
is not intrinsically isotropic. We found that the PYTHIA8 predicts a
relatively higher value of $\lambda_{\phi}$ compared to the experimental data.
Although, the $\lambda_{\theta\phi}$ parameter is almost zero for $J/\psi$ and
$\psi$(2S) in PYTHIA8, which is consistent with the LHC results as displayed
in Fig. 3.

Fig. 4 shows the rapidity-integrated ($2.5<y<4.0$) $p_{\rm T}$ dependence
polarization measurement of the angular observable $\lambda_{\theta}$ values
for $J/\psi$ and $\psi$(2S) in helicity frame, compared with NLO CSM
Butenschoen:2012px  and different NLO NRQCD calculation Butenschoen:2012px ;
Gong:2012ug ; Chao:2012iv ; Shao:2012fs (the NRQCD calculations differ in the
experimental data samples used to evalute the long-distance matrix elements).
The $p_{\rm T}$ dependence calculations based on the CSM disagree with the
experimental data and PYTHIA8 simulation result for $J/\psi$ and $\psi$(2S).
For $J/\psi$ resonance, a qualitative agreement between the experimental data
and theory is achieved by the NRQCD calculation Chao:2012iv ; Shao:2012fs .
However, the increasing trend of $\lambda_{\theta}$ parameter with $p_{\rm T}$
obtained from PYTHIA8 is in agreement up to some extent with NRQCD(1)
Butenschoen:2012px  and NRQCD(2) Gong:2012ug  calculations. Furthermore, for
$\psi$(2S), the LHCb data and PYTHIA8 favours the calculation from Ref.
Chao:2012iv ; Shao:2012fs  and from Ref.Gong:2012ug  at low $p_{\rm T}$. The
main difference between the polarization measurement of $J/\psi$ and
$\psi$(2S) in PYTHIA8 could be due to their different masses. In experimental
data, the $J/\psi$ polarization can come from directly produced $J/\psi$
mesons and from those produced in the decay of heavier charmonium states.
However, the polarization of $\psi$(2S) will remain unaffected by feed-down
decays of heavier charmonia.

Figure 6: (Color online) The $J/\psi$ and $\psi$(2S) polarization parameters
as a function of rapidity for $pp$ collisions at $\sqrt{s}$ = 7, 8, and 13 TeV
using PYTHIA8 in both helicity and Collins-Soper reference frames.

###  III.2 Charged-particle multiplicity dependence of $\lambda_{\theta}$,
$\lambda_{\phi}$, $\lambda_{\theta\phi}$

The charged particle multiplicity-dependent study of charmonia polarization
may reveal the underlying dynamics associated with the particle density of the
system produced in $pp$ collisions. Figure 5 shows the charged particle
multiplicity dependence of polarization parameters for $J/\psi$ and $\psi$(2S)
mesons in $pp$ collisions at $\sqrt{s}$ = 7, 8, and 13 TeV. The charged-
particle multiplicity classes used in the present analysis are taken from Ref.
Deb:2018qsl . The experimental study in this regard is reported by CMS
collaboration for three $\Upsilon(nS)$ states in $pp$ collisions at $\sqrt{s}$
= 7 TeV CMS:2016xpm . From the theoretical front, the relative multiplicity
dependence
$\left(\frac{dN_{ch}}{d\eta}/\left<\frac{dN_{ch}}{d\eta}\right>\right)$ study
of polarization parameters for $J/\psi$ in helicity and Collins-Soper frame is
studied in CGC+NRQCD approach Stebel:2021bbn  in $pp$ and $p$-Pb collisions at
$\sqrt{s}$ = 13 and 8.16 TeV, respectively. From Fig. 5, we observe that the
polarization parameter $\lambda_{\theta}$ indicates that the degree of
longitudinal polarization increases towards high multiplicity for $\psi$(2S)
meson, while for $J/\psi$ the longitudinal polarization remains almost
constant from low- to high multiplicity in both reference frames. The
$\lambda_{\phi}$ shows a transverse polarization in the helicity frame and a
longitudinal polarization in the Collins-Soper reference frame for $J/\psi$
and $\psi$(2S) for all multiplicity classes. The polarization parameter,
$\lambda_{\theta\phi}$, is negligible for $J/\psi$ and $\psi$(2S) at low
multiplicities, while at higher multiplicities, it has non-zero values. The
present study indicates that the charmonia polarization weakly depends on the
center of mass collision energy ($\sqrt{s}$).

###  III.3 Rapidity dependence of $\lambda_{\theta}$, $\lambda_{\phi}$,
$\lambda_{\theta\phi}$

In this section, we explore the rapidity dependence study of
$\lambda$-parameters for $J/\psi$ and $\psi$(2S) mesons in $pp$ collisions at
$\sqrt{s}$ = 7, 8, and 13 TeV, as shown in Fig. 6. The $J/\psi$ polarization
for various rapidity bins is reported in the LHCb experiment for $pp$
collisions at $\sqrt{s}$ = 7 TeV LHCb:2013izl  and data show a small
polarization within the uncertainties. In Fig. 6, the $p_{\rm{T}}$-integrated
polarization is obtained for the rapidity range 2.5 to 4.0, with a step of
0.3. The top panel of Fig. 6 shows that the degree of longitudinal
polarization for $J/\psi$ is larger than $\psi$(2S) in the helicity frame. It
observed that in the Collins-Soper frame, $\psi$(2S) is transversely polarized
while $J/\psi$ is longitudinally polarized. However, as observed in the
present study, there is no clear dependence of $\lambda_{\theta}$ parameter
with rapidity in both reference frames. Similar to $p_{\rm{T}}$ and charged-
particle multiplicity dependence, the $\lambda_{\phi}$ as a function of
rapidity shows a positive value for polarization in the helicity frame and a
negative value of polarization in the Collins-Soper frame for $J/\psi$ and
$\psi$(2S). The $\lambda_{\theta\phi}$ parameter shows almost zero
polarisation with rapidity irrespective of the particles under consideration,
the chosen center of mass energy, and the reference frame.

The finite polarization of $J/\psi$ and $\psi$(2S) in PYTHIA8 may have several
possible sources, a few of which are discussed below. In HICs (Au+Au, Pb+Pb,
etc.), the potential sources of hadrons polarization are the vorticity fields,
electromagnetic fields, and a strong vector meson force field STAR:2017ckg ;
Han:2017hdi ; Gao:2012ix ; Sheng:2020ghv ; Sheng:2022wsy . However, the
formation of vorticity due to the initial orbital angular momentum and the
electromagnetic fields is almost negligible in pp collisions. Therefore, the
charmonium polarization in ultra-relativistic pp collisions is an expected
consequence of the involved production mechanisms Faccioli:2022peq ;
Faccioli:2010kd . As PYTHIA8 considers LO NRQCD mechanisms for quarkonia
production, it can be anticipated that the finite polarization of charmonium
states is due to the QCD processes involved in PYTHIA8. In this work, we use
the flag HardQCD:all = on, which favors the production of heavy flavors
through all hard-QCD $2\rightarrow 2$ scattering processes. It involves pair
creation [q$\bar{q}$ (gg) $\rightarrow$ c$\bar{c}$], flavour excitation, and
gluon splitting processes manual . Moreover, a few basic principles, such as:
i) helicity conservation in the dilepton decays; ii) rotational covariance of
angular momentum eigenstates; and iii) conservation of parity and its
violation, govern the polarization of charmonium states Faccioli:2022peq ;
Faccioli:2010kd . In addition, we use the flag Charmonium:all = on, which
includes the production of charmonium and its higher excited states via the
color-singlet and color-octet mechanism of NRQCD, requiring information on
long-distance NRQCD matrix elements for the various wavefunctions involved in
PYTHIA8 manual . Further, the role of MPI on the decay angular distribution of
dimuon is investigated in Appendix A.

Another source of polarization is the inhomogeneous expansion of the medium,
resulting in an anisotropic flow in the transverse plane. The second
anisotropic flow-coefficient (i.e. elliptic flow) induced polarization is
observed for $\Lambda$-hyperons in HICs Becattini:2017gcx ; STAR:2019erd . In
$pp$ collisions, a similar flow-like phenomenon is obtained in PYTHIA8 due to
the considered multi-parton interactions (MPI) and color reconnection
mechanism OrtizVelasquez:2013ofg . In the present work, we use the same tune
MPI = on and CR = on; as an outcome, the flow-like effect indicates another
possible source of $J/\psi$ and $\psi(2S)$ polarization.

It is observed that PYTHIA8 serves as an effective model for jet studies in
$pp$ collisions with MPI and CR mode on ALICE:2023oww . The jet-like
fluctuation in PYTHIA8 is also a possible source of polarization for hadrons.
Because, if a jet is produced in the medium it deposits energy and creates a
smoke-loop-type vortex in the jet plane. As a consequence, it induces
polarization to the particles associated with the jets Betz:2007kg . Likewise,
charmonia can also get polarized due to such phenomena involved in the
production processes.

Along with the discussed polarization mechanisms, there are other possible
sources of polarization for charmonium states, but it is still a matter of
investigation whether the spin orientation or polarization occurs at the
partonic level and then is transferred to the bound state or directly to the
bound state. Based on the present study, we may infer that NRCQD-based
production mechanism give rise to the possible source for the partonic level
polarization and then is transferred to the bound state. However, the flow-
like effect is attributed to the direct polarization of the bound state.

##  IV Summary

In this work, we have studied the polarization parameters for $J/\psi$ and
$\psi$(2S) from the angular distribution of dimuons in $pp$ collisions at LHC
energies using PYTHIA8. The important observations of this paper are
summarized below:

  1. 1.

The $\lambda_{\theta}$, $\lambda_{\phi}$, $\lambda_{\theta\phi}$ are obtained
in the helicity and Collins-Soper reference frames in the rapidity interval
$2.5<y_{\mu\mu}<4.0$.

  2. 2.

It is observed from the $\lambda_{\theta}$ parameter that $J/\psi$ and
$\psi$(2S) are longitudinally polarized at low-$p_{\rm{T}}$ and transversely
polarized at high-$p_{\rm{T}}$ in both the reference frames.

  3. 3.

The $\lambda_{\phi}$ parameter indicates the longitudinal polarization in
helicity and transverse polarization in the Collins-Soper frame for $J/\psi$
and $\psi$(2S) across all energies. The $\lambda_{\theta\phi}$ parameter
values are close to zero as a function of $p_{\rm{T}}$.

  4. 4.

The multiplicity dependence study of $\lambda_{\theta}$ parameter shows the
degree of longitudinal polarization increases with charged particle
multiplicity for $\psi$(2S), while the behavior of longitudinal polarization
of $J/\psi$ stays constant with charged particle multiplicity, which needs
attention from an experimental point of view.

  5. 5.

In this study, we observe no clear dependence of $\lambda_{\theta}$ parameter
with rapidity. However the $\lambda_{\phi}$, $\lambda_{\theta}$ show almost
constant polarization with rapidity both for $J/\psi$ and $\psi$(2S) in
helicity and Collins-Soper reference frames. In the future, ALICE 3 set up
with a wider kinematics acceptance range of muon spectrometer, the rapidity
dependence study of polarization parameters would be an interesting topic.

  6. 6.

It is essential to mention that the polarization results obtained in this
analysis as a function $p_{\rm{T}}$, $N_{ch}$, and $y_{\mu\mu}$ consider only
the production of $J/\psi$ and $\psi$(2S), without taking the feed-down from
higher excited states. The investigation of charmonia polarization by taking
the feed-down from higher resonances is a future scope. In other words, the
newly upgraded Muon Forward Tracker (MFT) detector in the Run 3 and Run 4
setup of ALICE 2 to the muon Spectrometer, having better vertexing and precise
measurement capabilities, helps us to separate for prompt and non-prompt
charmonium states. Therefore, the investigation of polarization parameters for
prompt and non-prompt charmonium states in both the experiment and MC
simulation would be an intriguing subject to test.

  7. 7.

Since we used PYTHIA8, which incorporates pQCD and NRQCD-based processes, our
results overestimate $J/\psi$ polarization in some $p_{\rm{T}}$ bins as
compared with experimental data. It suggests that there might be an interplay
between these fundamental processes in the realistic scenario. On the other
hand, there might be some other processes responsible for net charmonia
polarization. Therefore, charmonia polarization in ultra-relativistic $pp$
collisions requires a thorough study using theoretical models confronted with
the experimental results.

## Acknowledgement

B. S. acknowledges the financial aid from CSIR, Government of India. SD
acknowledges the financial support under the Post-Doctoral Fellowship of CNRS
at IJCLAB, Orsay (France). The authors gratefully acknowledge the DAE-DST,
Government of India funding under the mega-science project ”Indian
Participation in the ALICE experiment at CERN” bearing Project No.
SR/MF/PS-02/2021-IITI (E-37123).

##  Appendix A ROLE OF MPI ON POLARIZATION PARAMETERS

Figure 7: (Color online) The $J/\psi$ polarization parameters as a function of
transverse momentum with “MPI: ON, CR: ON” (Red Marker) and “MPI: OFF, CR: ON”
(Green Marker) for $pp$ collisions at $\sqrt{s}$ = 7 using PYTHIA8 in helicity
and Collins-Soper reference frames.

In this study, our investigation centered on charmonia state polarization
within proton-proton collisions simulated through PYTHIA8. Multiple parton
interactions (MPI) in PYTHIA8 primarily affect parton distribution and
activity. However, their indirect influence potentially leads to alterations
in the angular momentum distribution among final-state particles due to
modified partonic interactions and subsequent hadronization processes. These
changes can manifest as variations in polarization parameters, particularly
$\lambda_{\theta}$, which are sensitive to the angular momentum states of
generated particles. Fig.7 highlights our investigation of the impact of MPI
(with on/off settings) on polarization parameters. Notably, our focus rests on
the observed behavior of $\lambda_{\theta}$, revealing distinctive trends:
when MPI is enabled (on), $\lambda_{\theta}$ displays an initial negative
trend, transitioning to positive values with increasing transverse momentum.
Conversely, in the absence of MPI, $\lambda_{\theta}$ consistently maintains
negativity with increasing transverse momentum. These trends persist
consistently in both the helicity frame and the Collins-Soper frame. The
intriguing shift from negative to positive $\lambda_{\theta}$ values with MPI
compared to the sustained negativity without MPI suggests a significant
influence of multiple parton interactions on polarization phenomena. Possible
explanations for these observations might involve variations in partonic
interactions and their consequential indirect effects on angular momentum
dynamics. This underscores the imperative need for further theoretical and
experimental investigations to unravel the underlying mechanisms driving
polarization in high-energy collisions.

## References

  * (1) N. Brambilla, S. Eidelman, B. K. Heltsley, R. Vogt, G. T. Bodwin, E. Eichten, A. D. Frawley, A. B. Meyer, R. E. Mitchell and V. Papadimitriou, et al. Eur. Phys. J. C 71, 1534 (2011).
  * (2) G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51, 1125 (1995); [erratum: Phys. Rev. D 55, 5853 (1997)].
  * (3) M. Krämer, Prog. Part. Nucl. Phys. 47, 141 (2001).
  * (4) F. Abe et al., [CDF Collaboration], Phys. Rev. Lett. 79, 572 (1997).
  * (5) T. Affolder et al., [CDF Collaboration], Phys. Rev. Lett. 84, 2094 (2000).
  * (6) S. Abachi et al., [D0 Collaboration], Phys. Lett. B 370, 239 (1996).
  * (7) U. A. Acharya et al., [PHENIX Collaboration], Phys. Rev. D 101, 052006 (2020).
  * (8) J. Adam et al., [STAR Collaboration], Phys. Rev. D 100, 052009 (2019).
  * (9) B. Trzeciak, [STAR Collaboration], J. Phys. Conf. Ser. 612, 012038 (2015).
  * (10) S. Acharya et al., [ALICE Collaboration], Eur. Phys. J. C 77, 392 (2017).
  * (11) S. Acharya et al., [ALICE Collaboration], JHEP 10, 084 (2019).
  * (12) J. Adam et al., [ALICE Collaboration], Eur. Phys. J. C 76, 184 (2016).
  * (13) R. Aaij et al., [LHCb Collaboration], JHEP 10, 172 (2015); [erratum: JHEP 05, 063 (2017)].
  * (14) R. Aaij et al., [LHCb Collaboration], JHEP 07, 134 (2018); [erratum: JHEP 05, 076 (2019).]
  * (15) A. M. Sirunyan et al., [CMS Collaboration], Phys. Lett. B 780, 251 (2018).
  * (16) G. Aad et al., [ATLAS Collaboration], Eur. Phys. J. C 76, 283 (2016).
  * (17) G. Aad et al., [ATLAS Collaboration], Phys. Rev. D 87, 052004 (2013).
  * (18) R. Baier and R. Ruckl, Phys. Lett. B 102, 364 (1981).
  * (19) E. L. Berger and D. L. Jones, Phys. Rev. D 23, 1521 (1981).
  * (20) C. H. Chang, Nucl. Phys. B 172, 425 (1980).
  * (21) V. G. Kartvelishvili, A. K. Likhoded and S. R. Slabospitsky, Yad. Fiz. 28, 1315 (1978) [Sov. J. Nucl. Phys. 28 , 678 (1978)].
  * (22) M. B. Einhorn and S. D. Ellis, Phys. Rev. D 12, 2007 (1975).
  * (23) G. T. Bodwin, E. Braaten, E. Eichten, S. L. Olsen, T. K. Pedlar and J. Russ, [arXiv:1307.7425 [hep-ph]].
  * (24) H. Fritzsch, Phys. Lett. B 67, 217 (1977).
  * (25) F. Halzen, Phys. Lett. B 69, 105 (1977).
  * (26) M. Gluck, J. F. Owens and E. Reya, Phys. Rev. D 17, 2324 (1978).
  * (27) E. Braaten and J. Russ, Ann. Rev. Nucl. Part. Sci. 64, 221 (2014).
  * (28) S. P. Baranov and A. V. Lipatov, Phys. Rev. D 96, 034019 (2017).
  * (29) J. P. Lansberg, Eur. Phys. J. C 61, 693(2009).
  * (30) P. Facial and C. Lourenço, Particle Polarization in High Energy Physics: An Introduction and Case Studies on Vector Particle Production at the LHC,’ Lect.Notes Phys. 1002 (2022) (Springer, 2022).
  * (31) P. H. De Moura, K. J. Goncalves and G. Torrieri, Phys. Rev. D 108, 034032 (2023).
  * (32) A. Abulencia et al., [CDF Collaboration], Phys. Rev. Lett. 99, 132001 (2007).
  * (33) T. Affolder et al., [CDF Collaboration], Phys. Rev. Lett. 85, 2886 (2000).
  * (34) D. Acosta et al., [CDF Collaboration], Phys. Rev. Lett. 88, 161802 (2002).
  * (35) T. Aaltonen et al., [CDF Collaboration], Phys. Rev. Lett. 108, 151802 (2012).
  * (36) V. M. Abazov et al., [D0 Collaboration], Phys. Rev. Lett. 101, 182004 (2008).
  * (37) A. Adare et al., [PHENIX Collaboration], Phys. Rev. D 82, 012001 (2010).
  * (38) L. Adamczyk et al., [STAR Collaboration], Phys. Lett. B 739, 180 (2014).
  * (39) J. Adam et al., [STAR Collaboration], Phys. Rev. D 102, 092009 (2020).
  * (40) B. Abelev et al., [ALICE Collaboration], Phys. Rev. Lett. 108, 082001 (2012).
  * (41) S. Acharya et al., [ALICE Collaboration], Eur. Phys. J. C 78, 562 (2018).
  * (42) S. Acharya et al., [ALICE Collaboration], Phys. Lett. B 815, 136146 (2021).
  * (43) S. Acharya et al., [ALICE Collaboration], Phys. Rev. Lett. 131, 042303 (2023).
  * (44) E. Etzion et al., [ATLAS Collaboration], Nucl. Phys. B Proc. Suppl. 187, 136 (2009).
  * (45) R. Aaij et al., [LHCb Collaboration], Eur. Phys. J. C 73, 2631 (2013).
  * (46) R. Aaij et al. [LHCb Collaboration], Eur. Phys. J. C 74, 2872 (2014).
  * (47) S. Chatrchyan et al., [CMS Collaboration], Phys. Lett. B 727, 381 (2013).
  * (48) S. Chatrchyan et al., [CMS Collaboration], Phys. Rev. Lett. 110, 081802 (2013).
  * (49) V. Khachatryan et al., [CMS Collaboration], Phys. Lett. B 761, 31 (2016).
  * (50) T. H. Chang et al., [NuSea Collaboration], Phys. Rev. Lett. 91, 211801 (2003).
  * (51) C. N. Brown et al., [NuSea Collaboration], Phys. Rev. Lett. 86, 2529 (2001).
  * (52) I. Abt et al., [HERA-B Collaboration], Eur. Phys. J. C 60, 517 (2009).
  * (53) M. Butenschoen and B. A. Kniehl, Phys. Rev. Lett. 108, 172002 (2012).
  * (54) B. Gong, L. P. Wan, J. X. Wang and H. F. Zhang, Phys. Rev. Lett. 110, 042002 (2013).
  * (55) H. S. Shao and K. T. Chao, Phys. Rev. D 90, 014002 (2014).
  * (56) K. T. Chao, Y. Q. Ma, H. S. Shao, K. Wang and Y. J. Zhang, Phys. Rev. Lett. 108, 242004 (2012).
  * (57) Y. Q. Ma, T. Stebel and R. Venugopalan, JHEP 12, 057 (2018).
  * (58) P. Faccioli, C. Lourenco, J. Seixas and H. K. Wohri, Phys. Rev. Lett. 102, 151802 (2009).
  * (59) H. Haberzettl and J. P. Lansberg, Phys. Rev. Lett. 100, 032006 (2008).
  * (60) H. S. Shao, H. Han, Y. Q. Ma, C. Meng, Y. J. Zhang and K. T. Chao, JHEP 05, 103 (2015).
  * (61) M. Beneke, M. Kramer and M. Vanttinen, Phys. Rev. D 57, 4258 (1998).
  * (62) M. Beneke and M. Krämer, Phys. Rev. D 55, 5269 (1997).
  * (63) E. Braaten, B. A. Kniehl and J. Lee, Phys. Rev. D 62, 094005 (2000).
  * (64) B. A. Kniehl and J. Lee, Phys. Rev. D 62, 114027 (2000).
  * (65) A. K. Leibovich, Phys. Rev. D 56, 4412 (1997).
  * (66) P. Faccioli, C. Lourenco, J. Seixas and H. K. Wohri, Eur. Phys. J. C 69, 657 (2010).
  * (67) P. Faccioli, C. Lourenco and J. Seixas, Phys. Rev. D 81, 111502 (2010).
  * (68) P. Faccioli, C. Lourenco and J. Seixas, Phys. Rev. Lett. 105, 061601 (2010).
  * (69) B. Gong and J. X. Wang, Phys. Rev. Lett. 100, 232001 (2008).
  * (70) E. Braaten and T. C. Yuan, Phys. Rev. Lett. 71, 1673 (1993).
  * (71) P. Faccioli, V. Knünz, C. Lourenco, J. Seixas and H. K. Wöhri, Phys. Lett. B 736, 98 (2014).
  * (72) H. S. Chung, C. Yu, S. Kim and J. Lee, Phys. Rev. D 81, 014020 (2010).
  * (73) V. Cheung and R. Vogt, Phys. Rev. D 98, 114029 (2018).
  * (74) V. Cheung and R. Vogt, Phys. Rev. D 99, 034007 (2019).
  * (75) V. Cheung and R. Vogt, Phys. Rev. D 104, 094026 (2021).
  * (76) V. Cheung and R. Vogt, Phys. Rev. C 105, 055202 (2022).
  * (77) R. Vogt, Nucl. Phys. A 982, 751 (2019).
  * (78) J. C. Collins and D. E. Soper, Phys. Rev. D 16, 2219 (1977).
  * (79) K. Gottfried and J. D. Jackson, Nuovo Cim. 33, 309 (1964).
  * (80) S. Sarkar, H. Satz and B. Sinha, Lect. Notes Phys. 785 (2010).
  * (81) Pythia8 online manual:(https://pythia.org/latest-manual/Welcome.html).
  * (82) S. Deb, D. Thakur, S. De and R. Sahoo, Eur. Phys. J. A 56, 134 (2020).
  * (83) D. Thakur, S. De, R. Sahoo and S. Dansana, Phys. Rev. D 97, 094002 (2018).
  * (84) T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP 0605, 026 (2006).
  * (85) R. Corke and T. Sjostrand, JHEP 03, 032 (2011).
  * (86) H. S. Shao, Comput. Phys. Commun. 184, 2562 (2013).
  * (87) W. E. Caswell and G. P. Lepage, Phys. Lett. B 167, 437 (1986).
  * (88) T. Stebel and K. Watanabe, Phys. Rev. D 104, 034004 (2021).
  * (89) L. Adamczyk et al. [STAR Collaboration], Nature 548 62 (2017).
  * (90) Z. Z. Han and J. Xu, Phys. Lett. B 786, 255 (2018).
  * (91) J. H. Gao, Z. T. Liang, S. Pu, Q. Wang and X. N. Wang, Phys. Rev. Lett. 109, 232301 (2012).
  * (92) X. L. Sheng, Q. Wang and X. N. Wang, Phys. Rev. D 102, 056013 (2020).
  * (93) X. L. Sheng, L. Oliva, Z. T. Liang, Q. Wang and X. N. Wang, Phys. Rev. Lett. 131, 042304 (2023).
  * (94) F. Becattini and I. Karpenko, Phys. Rev. Lett. 120, 012302 (2018).
  * (95) J. Adam et al. [STAR Collaboration], Phys. Rev. Lett. 123, 132301 (2019).
  * (96) A. Ortiz Velasquez, P. Christiansen, E. Cuautle Flores, I. Maldon- ado Cervantes, G. Paić, Phys. Rev. Lett. 111, 042001 (2013).
  * (97) B. Betz, M. Gyulassy and G. Torrieri, Phys. Rev. C 76, 044901 (2007).
  * (98) S. Acharya et al. [ALICE Collaboration], [arXiv:2311.13322 [hep-ex]].