# J/ψ polarization in relativistic heavy ion collisions

Jiaxing Zhao¹,ᵃ, Baoyi Chen²

¹ SUBATECH, Université de Nantes, IMT Atlantique, IN2P3/CNRS, 4 rue Alfred Kastler, 44307 Nantes cedex 3, France

² Department of Physics, Tianjin University, Tianjin 300354, China

Received: 16 January 2024 / Accepted: 20 August 2024 / Published online: 30 August 2024
© The Author(s) 2024

**Abstract** Understanding the polarization property of $J/ψ$ is critical to constrain its production mechanism. In addition, the polarization of $J/ψ$ can reveal the impact of strong electromagnetic and vorticity fields in relativistic heavy ion collisions. In this study, we analyzed the yield and polarization of $J/ψ$ in relativistic heavy ion collisions at different centrality and transverse momentum regions, using three different reference frames: the Collins–Soper frame, the helicity frame, and the event plane frame. The polarization of initially produced $J/ψ$ is determined by the NRQCD calculation and is similar to that of pp collisions. However, both unpolarization and transverse polarization are considered for the regenerated $J/ψ$. Our results indicate that the polarization at high $p_T$ is similar to that observed in pp collisions. However, at low $p_T$, where regenerated $J/ψ$ dominates, it is likely that the polarized charm quarks in the rotational QGP medium are responsible for this phenomenon. Our study supplies a baseline for future research on the effects of strong electromagnetic and vorticity fields on $J/ψ$ polarization.

## 1 Introduction

Quarkonium is one of the most ideal probes for quantum chromodynamics (QCD), which provides significant progress in understanding of the QCD in both perturbative and non-perturbative ways. Quarkonium production can be dealt with the factorization theory. Heavy quark pairs are produced via perturbative QCD, due to $m_Q \gg \Lambda_{QCD}$, with heavy quark mass $m_Q$ and the QCD cutoff $\Lambda_{QCD}$. The formation of colorless quarkonium from the heavy quark pairs $Q\bar{Q}$ is a non-perturbative process, which is usually encoded in some parameters, such as the long-distance matrix elements in non-relativistic QCD (NRQCD), and it can be fixed by the experimental data. Quarkonium production and momentum

spectrum in pp collisions have been well described by many theoretical models, such as the Color-Evaporation Model (CEM) [1,2], the Color-Singlet Model (CSM) [3,4] and the Color-Octet Model (COM) [5], the latter two are encompassed in the NRQCD [6–9].

However, the yield and spectrum are insufficient to rule out anything so far, and it's hard to get a further understanding of the production mechanism. The quarkonium polarization introduces a new dimension to constrain the production mechanism. For the polarization of $J/ψ$, the CSM gives a transversal polarization at leading order (LO) but a longitudinal polarization at next to leading order (NLO) [10]. The COM shows a small longitudinal polarization at small $p_T$ regions [11], while a transverse polarization at large $p_T$ regions [12]. NRQCD includes all important contributions from both color-singlet and color-octet intermediate states and gives a detailed interview of the $J/ψ$ polarization [7–9]. Further, NRQCD plus color glass condensate (CGC) makes the prediction down to very low $p_T$ regions [13]. Recently, the polarization of $J/ψ$ is also studied via the improved color evaporation model [14,15]. In experiments, $J/ψ$ polarization has been observed in different colliding systems, energies, and different kinematic ranges [16–19]. However, due to the large uncertainty, it is hard to get a conclusion so far.

In relativistic heavy ion collisions, there are two sources for $J/ψ$ production: one is the initial production through hard processes which is largely suppressed in the medium [20,21], and the other one is the regeneration in the quark–gluon plasma (QGP) via the recombination from uncorrelated charm and anticharm [22–24]. The competition between the suppression and regeneration can explain well the experimentally measured $J/ψ$ nuclear modification factor $R_{AA}$ and collective flows [25–32] at RHIC and LHC energies. Except for the yield, the polarization of $J/ψ$ concerning the event plane in $PbPb$ collisions at $\sqrt{s_{NN}} = 5.02$ TeV at the LHC is observed in the first time [33]. Initial produced $J/ψ$ may carry different polarization information compared

<footnote>ᵃ e-mail: zhao-jx15@tsinghua.org.cn (corresponding author)</footnote>
to the regenerated one, as discussed before [34]. So, this will give new constraints on the regeneration component and also reveal the influence of strong electromagnetic and vorticity fields on the polarization of $J/\psi$.

In this paper, we want to study the yield and polarization of $J/\psi$ in relativistic heavy ion collisions. The paper is organized as follows. The $J/\psi$ polarization in pp collisions is shown in Sect. 2. Next the production and transverse spectrum of $J/\psi$ in PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV is investigated in Sect. 3. In Sect. 4, we will show the polarization of $J/\psi$ in PbPb collisions. A conclusion will be given in Sect. 5.

## 2 Polarization in proton–proton collisions

The geometrical shape of the angular distribution of the two decay products reflects the polarization of the quarkonium state in dilepton decays. The distribution of the $J/\psi$ decay products can be expressed in general as

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

where $\theta$ and $\phi$ are the polar and azimuthal angles in a given reference frame. The polarization parameters are related to the spin states of decay products, so the spin state of $J/\psi$. The polarization is completely longitudinal if the set of polarization parameters $(\lambda_{\theta}, \lambda_{\phi}, \lambda_{\theta\phi})$ takes the values $(-1, 0, 0)$ and it is completely transverse if it takes the values $(1, 0, 0)$. In the unpolarization scenario, the parameters are $(0, 0, 0)$.

The polarization parameters are frame-dependent, which is defined with respect to a particular polarization axis, normally taken as z-axis. There are three main reference frames. There are the Collins–Soper (CS) frame, the helicity frame (HX), and the newly used event plane frame (EP). In the CS frame, the z-axis is defined as the bisector of the angles between the direction of one beam and the opposite of another beam in the rest frame of the $J/\psi$. In the HX frame, z-axis is defined as the direction of the $J/\psi$ in the centre of mass frame of two colliding nucleus, which is same as the Laboratory frame for the colliding system. In the EP frame, z-axis is defined as the direction orthogonal to the event plane in the centre of mass frame of two colliding nucleus. This is also the direction of the strongest net-magnetic field and global vorticity field in heavy ion collisions.

Although the polarization parameters $(\lambda_{\theta}, \lambda_{\phi}, \lambda_{\theta\phi})$ depend on the frame, the values in different frame are connected via the geometrical transformation. The polarization parameters between two reference frames are given by [35],

$$\lambda'_{\theta} = \frac{\lambda_{\theta} - 3\Lambda}{1 + \Lambda},$$

$$\begin{aligned} \lambda'_{\phi} &= \frac{\lambda_{\phi} + \Lambda}{1 + \Lambda}, \\ \lambda'_{\theta\phi} &= \frac{\lambda_{\theta\phi} \cos 2\delta - (\lambda_{\theta} - \lambda_{\phi})/2 \sin 2\delta}{1 + \Lambda}, \end{aligned} \quad (2)$$

where $\Lambda = (\lambda_{\theta} - \lambda_{\phi})/2 \sin^2 \delta - \lambda_{\theta\phi}/2 \sin 2\delta$ with the angle $\delta$ between two polarization axis. Usually, this transformation is with a strong kinematic dependence, which is encoded in the angles $\delta$. The angle between the CS and HX frame can be given,

$$\delta_{\text{CS} \to \text{HX}} = \arccos \left( \frac{m_{J/\psi} \sinh y}{\sqrt{p_T^2 + E_T^2 \sinh^2 y}} \right), \quad (3)$$

where $y$ is the rapidity and $p_T$ is transverse momentum. The transverse energy $E_T = \sqrt{m_{J/\psi}^2 + p_T^2}$. Based on the definition, the angle between the HX and EP frame can be given by

$$\delta_{\text{HX} \to \text{EP}} = \arccos \left( \frac{p_T \sin \phi}{\sqrt{p_T^2 + E_T^2 \sinh^2 y}} \right). \quad (4)$$

So it depends not only on the $J/\psi$ momentum $p_T$ but also the angle $\phi$ in the transverse plane. The angle $\phi$ dependent distribution can be expressed as the anisotropic flow. Such as the elliptic flow, its definition is $v_2 \equiv \langle \cos 2\phi \rangle$. So, $\langle \phi \rangle = 1/2 \arccos v_2$. The angles $\delta$ are related to the value of the $v_2$. If the $v_2 = 0$ for any $p_T$ regions, we obtain $\langle \phi \rangle = \pi/4$. So, the $\delta_{\text{HX} \to \text{EP}} \approx \pi/2$ in the forward rapidity.

The experimental data of $J/\psi \lambda_{\theta}, \lambda_{\phi}$, and $\lambda_{\theta\phi}$ are shown in Fig. 1. The solid lines in the CS and HX frames are from the NRQCD calculation [13], which are related to each other via the transformation shown in Eq. (2). We can see that at small $p_T$, the NRQCD result shows a small transverse polarization, $\lambda_{\theta} > 0$, while the data on the other hand seem to favor a weak longitudinal polarization $\lambda_{\theta} < 0$ (especially the data from pp 7 TeV) but with larger uncertainty. For higher transverse momentum, e.g. $p_T > 6$ GeV, the NRQCD result agrees with the data. For the $\lambda_{\phi}$ and $\lambda_{\theta\phi}$, the experimental data can be described very well by the NRQCD. The values in the EP frame, which is shown with blue lines in the right panel, are obtained from the HX frame quantities via the above transformation and the angle $\delta_{\text{HX} \to \text{EP}}$.

## 3 J/ψ transport in heavy ion collisions

We now focus on charmonium transport and production in heavy ion collisions, which can be described by a classical transport equation [25–31]. The charmonium phase space distribution $f_{\psi}(p, x, \tau|b)$, where $\psi = J/\psi, \chi_c, \psi'$, is governed by the Boltzmann equation [36],
<figure><img src="image_1.png" /><figcaption>**Fig. 1** $p_T$ dependent $\lambda_\theta$, $\lambda_\phi$, and $\lambda_{\theta\phi}$ in forward rapidity $2.5 < y < 4.0$ of pp collisions in the Collins-Soper frame (left), the helicity frame (middle), and the event plane frame (right). The experiment data is from LHCb $\sqrt{s} = 7$ TeV [17], ALICE $\sqrt{s} = 7$ TeV [18], and ALICE $\sqrt{s} = 8$</figcaption></figure>

TeV [19]. Solid lines in the left and middle panels are from the NRQCD calculation and taken from Ref. [13] directly, while the lines in the right panel are from the transformation as shown in Eq. (2)

$$ \left[ \cosh(y - \eta)\partial_{\tau} + \frac{\sinh(y - \eta)}{\tau} \partial_{\eta} + \mathbf{v}_T \cdot \nabla_T \right] f_{\psi} = -\alpha f_{\psi} + \beta, \qquad (5) $$

where $y = 1/2 \ln[(E + p_z)/(E - p_z)]$ is the charmonium momentum rapidity, $\eta = 1/2 \ln[(t+z)/(t-z)]$ is the pseudorapidity of charmonium, and $v_T = p_T/E_T$ is the charmonium transverse velocity with transverse energy $E_T$. The second and third terms on the left-hand side arise from the free streaming of $\psi$, which leads to the leakage effect in the longitudinal and transverse directions. The charmonium suppression and regeneration processes in the QGP medium are reflected in the loss term $\alpha$ and gain term $\beta$ in the right-hand side, respectively.

In the hot QGP medium, the interaction potential between charm and anticharm suffers Debye screening [37]. This may leads to the suppression of the charmonium states. The dissociation temperatures $T_d$ of charmonium states can be given by solving the Schrödinger equation with the temperature-dependent heavy quark potential. In this study, we take $T_d \approx 2.3T_c$ for $J/\psi$, which can be obtained by using the potential from lattice QCD [38] as well as the inter-

nal energy [39]. Besides the Debye screening, charmonium may also suffer the dynamical dissociation in the QGP via scattering with thermal partons, such as the gluon dissociation [40,41] $g\psi \to c\bar{c}$. For the ground state $J/\psi$, the gluon dissociation cross-section $\sigma_{g\psi}^{c\bar{c}}$ in vacuum is derived using the operator-production-expansion (OPE) method [40,41]. For the excited states $\chi_c$ and $\psi'$, the cross sections can be obtained through their geometric relation to the ground state [29]. In this study, we take only the gluon dissociation as the loss term and its inverse process $c\bar{c} \to g\psi$ as the gain term, $\alpha$, and $\beta$ can be explicitly expressed as [36]

$$ \begin{aligned} \alpha(\mathbf{p}, \mathbf{x}, \tau|\mathbf{b}) &= \frac{1}{2E_T} \int \frac{d^3\mathbf{p}_g}{(2\pi)^3 2E_g} W_{g\psi}^{c\bar{c}}(s) f_g(\mathbf{p}_g, \mathbf{x}, \tau) \\ & \times \Theta(T(\mathbf{x}, \tau|\mathbf{b}) - T_c), \\ \beta(\mathbf{p}, \mathbf{x}, \tau|\mathbf{b}) &= \frac{1}{2E_T} \int \frac{d^3\mathbf{p}_g}{(2\pi)^3 2E_g} \frac{d^3\mathbf{p}_c}{(2\pi)^3 2E_c} \frac{d^3\mathbf{p}_{\bar{c}}}{(2\pi)^3 2E_{\bar{c}}} \\ & \times W_{c\bar{c}}^{g\psi}(s) f_c(\mathbf{p}_c, \mathbf{x}, \tau|\mathbf{b}) f_{\bar{c}}(\mathbf{p}_{\bar{c}}, \mathbf{x}, \tau|\mathbf{b}) \\ & \times (2\pi)^4 \delta^{(4)}(p + p_g - p_c - p_{\bar{c}}) \\ & \times \Theta(T(\mathbf{x}, \tau|\mathbf{b}) - T_c), \end{aligned} \qquad (6) $$
where $E_g$, $E_c$, and $E_{\bar{c}}$ are the gluon, charm quark, and anti-charm quark energies, $p_g$, $p_c$ and $p_{\bar{c}}$ are their momenta, and $s$ is the invariant mass of $g\psi$ system. $b$ is impact parameter of two nucleus. $W_{g\psi}^{cc}(W_{c\bar{c}}^{g\psi})$ is the dissociation(regeneration) probability. The regeneration probability and dissociation probability are connected via the detailed balance [24,30]. The step function $\Theta$ is used to guarantee the calculation in the QGP phase above the critical temperature $T_c$, and the local temperature of the medium $T(x, \tau)$ is given by the hydrodynamic simulation. In this paper, we employ a (2+1)-dimensional hydrodynamic model, the MUSIC package [42,43], to characterize the space and time dependence of the temperature and velocity of the hot medium.

The gluon distribution $f_g$ is chosen as the Bose-Einstein distribution in the loss term. The experimentally measured large charmed meson flow indicates that charm quarks with small transverse momentum might be thermalized in the medium [44]. So, as a first approximation, we take a thermal distribution for charm (anti-charm) quarks in the momentum space $f_c(p_c, x) = \rho_c(x)N_c/(e^{p_c^2 u_\mu/T} + 1)$, where $N_c$ is the normalization factor, and the density in coordinate space $\rho_c$ is controlled by the charm conservation equation $\partial_\mu(\rho_c u^\mu) = 0$. The initial density at time $\tau_0$ is governed by the nuclear geometry of the colliding system,

$$ \rho_c(\mathbf{x}, \tau_0 | \mathbf{b}) = \frac{T_A(\mathbf{x}_T + \mathbf{b}/2)T_B(\mathbf{x}_T - \mathbf{b}/2) \cosh \eta \, d\sigma_{pp}^{cc}}{\tau_0 \, d\eta}, \quad (7) $$

where $T_A$ and $T_B$ are the thickness functions [45]. $x_T$ is transverse radius. $d\sigma_{pp}^{cc}/d\eta$ is the rapidity distribution of charm quark production cross section in $pp$ collisions. We take $d\sigma_{pp}^{cc}/d\eta = 0.61$ mb as the upper limit in the forward rapidity of $PbPb$ collisions at $\sqrt{s_{NN}} = 5.02$ TeV, which is almost the upper limit of the experimental data [46]. Considering the shadowing effect in the high energy $PbPb$ collisions, the charm cross section can be suppressed. In this study, we assume a constant shadowing factor of 0.7 for the charm production, which is close to the value given by the EPS09 [52]. This leads to a lower limit $d\sigma_{pp}^{cc}/d\eta = 0.43$ mb.

To solve the transport equation, the initial charmonium distribution is needed, which can be written as

$$ \begin{aligned} f_{\psi}(\mathbf{p}, \mathbf{x}, \tau_0 | \mathbf{b}) &= \frac{(2\pi)^3}{E_T \tau_0} \int dz_A dz_B \rho_A(\mathbf{x}_T + \mathbf{b}/2, z_A) \rho_B(\mathbf{x}_T - \mathbf{b}/2, z_B) \\ & \quad \times \mathcal{R}_g(x_1, \mu_F, \mathbf{x}_T + \mathbf{b}/2) \mathcal{R}_g(x_2, \mu_F, \mathbf{x}_T - \mathbf{b}/2) \\ & \quad \times f_{pp}^{\psi}(\mathbf{x}_T, \mathbf{p}, z_A, z_B | \mathbf{b}). \end{aligned} \quad (8) $$

It is just a superposition of the charmonium distribution in $pp$ collisions with the cold nuclear effects, such as nuclear shadowing [49] and Cronin effect [50,51]. For heavy ion

collisions at LHC energy, the collision time is very short, the nuclear absorption can be safely neglected [20]. $\rho_A$ and $\rho_B$ are the nucleon distributions in the two colliding nuclei. The shadowing effect is embedded in the factor $R_g$, which is taken from the EPS09 package [52]. Before two gluons fuse into a charmonium, they may acquire additional transverse momentum via multi-scattering with the surrounding nucleons, and this extra momentum would be inherited by the produced charmonium, this is the so-called Cronin effect. In this study, we replace the averaged transverse momentum square by $\langle p_T^2 \rangle_{pp} + a_{gN}l$, where the Cronin parameter $a_{gN}$ is the averaged charmonium transverse momentum square obtained from the gluon scattering with a unit of length of nucleons, and $l$ is the mean trajectory length of the two gluons in the projectile and target nuclei before $c\bar{c}$ formation. We take $a_{gN} = 0.15 \text{ GeV}^2/\text{fm}$ for $PbPb$ collisions at $\sqrt{s_{NN}} = 5.02 \text{ TeV}$. $\langle p_T^2 \rangle_{pp}$ is the mean transverse momentum square, which can be obtained together with the charmonium momentum distribution $f_{pp}^{\psi}$ by fitting the $pp$ collisions, see detail [27–29].

The transport equation with the initial distribution (8) can be solved analytically, as shown in detail [27–29]. To obtain the prompt $J/\psi$, which is normally measured in the experiment, the feed-down from the excited charmonium states contribution. The excited states decay into the ground state with the branching ratio $B_{\chi_c \to J/\psi} = 22\%$ and $B_{\psi(2S) \to J/\psi} = 61\%$ respectively after they moves out of the hot medium [53]. The nuclear modification factor $R_{AA}$ as a function of the number of participants $N_{part}$ and transverse momentum $p_T$ are studied and shown in Fig. 2. We can see the $J/\psi$ suffers a large suppression in central collisions, especially the low $p_T$ $J/\psi$. While the regeneration play an important role in low $p_T$ and central collisions, which compensates the suppression and together give a good description of the experimental data.

## 4 Polarization in heavy ion collisions

The polarization of $J/\psi$ produced in relativistic heavy ion collisions can be expressed as

$$ \lambda_{\theta}(p_T) = \frac{N_{\text{ini}}(p_T)\lambda_{\theta}^{\text{ini}}(p_T) + N_{\text{reg}}(p_T)\lambda_{\theta}^{\text{reg}}(p_T)}{N_{\text{ini}}(p_T) + N_{\text{reg}}(p_T)}, \quad (9) $$

where $\lambda_{\theta}^{\text{ini}}$ and $\lambda_{\theta}^{\text{reg}}$ are for initial and regenerated $J/\psi$, respectively. The same relationship holds for $\lambda_{\phi}$ and $\lambda_{\theta\phi}$.

First, we assume the regenerated $J/\psi$ without polarization information, $\lambda_{\theta}^{\text{reg}} = 0$. The polarization of final $J/\psi$ comes completely from the initially produced ones, which is the same as the $pp$. Due to the small interaction cross section and the heavy flavor symmetry, the spin flips of the heavy quark in the hot medium are largely suppressed. The
<figure><img src="image_2.png" /><figcaption>**Fig. 2** $J/\psi$ $R_{AA}$ in forward rapidity $2.5 < y < 4$ of PbPb collisions with $\sqrt{s_{NN}} = 5.02$ TeV as a function of number of participant $N_{part}$ (upper) and transverse momentum $p_T$ (lower) in 0–20% centrality regions. The dotted and dashed lines represent initial and regeneration contributions, respectively. The thick and thin dashed lines are the regeneration fraction with the charm cross section $d\sigma_{pp}^{cc}/d\eta = 0.61$ mb and 0.43 mb. The bands are the full result with different charm cross sections. The experimental data are from ALICE [47,48]</figcaption></figure>

hot medium may melt the initial produced $J/\psi$, but doesn't wash out the polarization of the survived $J/\psi$. With the initial production fraction, we calculate the $p_T$ and number of collisions $N_{part}$ dependent polarization parameters $\lambda_\theta$, $\lambda_\phi$, and $\lambda_{\theta\phi}$ in PbPb collisions within three different frames as shown in Fig. 3. Because high $p_T$ $J/\psi$, e.g. $p_T > 5$ GeV, almost comes from the initial production, as shown in Fig 2, their polarization should be the same as the initial production. It's indicated by the experimental data (blue points) and shown in Fig. 3 that the $\lambda_\theta$ at $p_T > 5$ GeV is close to the theoretical calculation (blue line). For the centrality dependence, we can see the polarization disappear at central collisions in this case due to the regeneration dominates. In the very peripheral collisions, the polarization is consistent with the $pp$ collisions.

For the low $p_T$, we can see the larger difference between the theoretical result and experimental data. This indicates the contribution of the regenerated $J/\psi$ cannot be ignored. The polarization of the regenerated $J/\psi$ is from the polarized and uncorrelated charm and anticharm quarks in the

<figure><img src="image_3.png" /><figcaption>**Fig. 3** $p_T$ (left) and $N_{part}$ (right) dependent $J/\psi$ $\lambda_\theta$, $\lambda_\phi$, and $\lambda_{\theta\phi}$ in forward rapidity $2.5 < y < 4$ of PbPb collisions with $\sqrt{s_{NN}} = 5.02$ TeV. The results are shown in the Collins–Soper frame (orange), the helicity frame (green), and the event plane frame (blue). The band is caused by the uncertainty of the nuclear shadowing effect. The experimental data is from ALICE [33] and obtained concerning the event plane</figcaption></figure>

QGP. The charm quark can be polarized via coupling to the electromagnetic and vorticity field as light quark [54–56]. The spin of charm and anticharm quarks will tend towards the axis of rotation due to the coupling $\omega \cdot s$. In a strong magnetic field, the charm move towards the axis of the magnetic field while the anticharm move towards the opposite direction due to the coupling $\mu \cdot B$ with the magnetic moment $\mu$. In the limit, the vorticity field will induce a transverse polarization of $J/\psi$ with $\lambda_\theta = 1$, which is combined by the polarized charm and anticharm quarks, while the magnetic field will induce a longitudinal polarization with $\lambda_\theta = -1$. Because the strength of the magnetic field decays fast and becomes very small when $J/\psi$ regenerated [57–61], we only consider the polarization caused by the vorticity field. We can see if the regenerated $J/\psi$ is transversely polarized, the $\lambda_\theta$ at low $p_T$ has a large and positive value in the event plane frame, as shown in Fig. 4. So, the large $\lambda_\theta$ at 3 GeV observed in the experimental data indicates the largely polarized charm quarks in the QGP. The high $p_T$ behavior is taken over by the initial production and the same as Fig. 3.

## 5 Summary

In summary, we investigated the effects of the hot medium on the production and polarization of $J/\psi$ in relativistic heavy ion collisions across various transverse momentum
<figure><img src="image_4.png" /><figcaption>**Fig. 4** Same as the Fig. 3, but considering the polarization of the regenerated $J/\psi$</figcaption></figure>

bins. The polarization observables are presented in three different reference frames: the Collins-Soper frame, the helicity frame, and the event plane frame considering two cases. One is that the polarization of $J/\psi$ in PbPb collisions is assumed to arise from the initial process, which is given by the NRQCD framework, while the charm quarks and regenerated $J/\psi$ are unpolarized. The suppression of the initial $J/\psi$ is attributed to the dissociation of the primordial polarized $J/\psi$ due to color screening effects and parton inelastic scatterings. However, the polarization information of the survived $J/\psi$ doesn't be washed out in the QGP. In the second case, the initial produced is the same as before, while the regenerated $J/\psi$ is transversely polarized, which comes from the polarized charm and anticharm quarks in the vorticity field. Compared to the experimental data, we find the polarization of high $p_T$ $J/\psi$ comes from the initial production, and a large $\lambda_\theta$ at low $p_T$ indicates the polarization of the charm quarks in a rotational QGP medium. The precise measurement at low $p_T$ will give the constraint of the contribution of the vorticity field and regeneration. Additionally, the polarization of charmonium excited states may exhibit differences compared to the ground state. Complete treatment of these processes will be performed in the future.

**Acknowledgements** This work is funded by the European Union's Horizon 2020 research and innovation program under grant agreement No. 824093 (STRONG-2020). Baoyi Chen is funded by the National Natural Science Foundation of China (NSFC) under Grants no. 12175165.

**Data Availability Statement** Data will be made available on reasonable request. [Authors' comment: The datasets generated during and/or

analysed during the current study are available from the corresponding author on reasonable request.]

**Code Availability Statement** Code/software will be made available on reasonable request. [Authors' comment: The code/software generated during and/or analysed during the current study is available from the corresponding author on reasonable request].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Funded by SCOAP³.

## References

1. H. Fritzsch, Phys. Lett. B 67, 217 (1977). https://doi.org/10.1016/0370-2693(77)90108-3

2. J.F. Amundson, O.J.P. Eboli, E.M. Gregores, F. Halzen, Phys. Lett. B 372, 127 (1996). https://doi.org/10.1016/0370-2693(96)00035-4. arXiv:hep-ph/9512248

3. C.-H. Chang, Nucl. Phys. B 172, 425 (1980). https://doi.org/10.1016/0550-3213(80)90175-3

4. R. Baier, R. Ruckl, Phys. Lett. B 102, 364 (1981). https://doi.org/10.1016/0370-2693(81)90636-5

5. G.T. Bodwin, E. Braaten, T.C. Yuan, G.P. Lepage, Phys. Rev. D 46, R3703 (1992). https://doi.org/10.1103/PhysRevD.46.R3703. arXiv:hep-ph/9208254

6. G.T. Bodwin, E. Braaten, G.P. Lepage, Phys. Rev. D 51, 1125 (1995). https://doi.org/10.1103/PhysRevD.55.5853. arXiv:hep-ph/9407339. [Erratum: Phys. Rev. D 55, 5853 (1997)]

7. M. Butenschoen, B.A. Kniehl, Phys. Rev. Lett. 108, 172002 (2012). https://doi.org/10.1103/PhysRevLett.108.172002. arXiv:1201.1872 [hep-ph]

8. B. Gong, L.-P. Wan, J.-X. Wang, H.-F. Zhang, Phys. Rev. Lett. 110, 042002 (2013). https://doi.org/10.1103/PhysRevLett.110.042002. arXiv:1205.6682 [hep-ph]

9. K.-T. Chao, Y.-Q. Ma, H.-S. Shao, K. Wang, Y.-J. Zhang, Phys. Rev. Lett. 108, 242004 (2012). https://doi.org/10.1103/PhysRevLett.108.242004. arXiv:1201.2675 [hep-ph]

10. B. Gong, J.-X. Wang, Phys. Rev. Lett. 100, 232001 (2008). https://doi.org/10.1103/PhysRevLett.100.232001. arXiv:0802.3727 [hep-ph]

11. M. Beneke, I.Z. Rothstein, Phys. Rev. D 54, 2005 (1996). https://doi.org/10.1103/PhysRevD.54.2005. arXiv:hep-ph/9603400. [Erratum: Phys. Rev. D 54, 7082 (1996)]

12. B. Gong, X.Q. Li, J.-X. Wang, Phys. Lett. B 673, 197 (2009). https://doi.org/10.1016/j.physletb.2009.02.026. arXiv:0805.4751 [hep-ph]. [Erratum: Phys. Lett. B 693, 612–613 (2010)]

13. Y.-Q. Ma, T. Stebel, R. Venugopalan, JHEP 12, 057 (2018). https://doi.org/10.1007/JHEP12(2018)057. arXiv:1809.03573 [hep-ph]

14. V. Cheung, R. Vogt, Phys. Rev. D 98, 114029 (2018). https://doi.org/10.1103/PhysRevD.98.114029. arXiv:1808.02909 [hep-ph]

15. V. Cheung, R. Vogt, Phys. Rev. D 104, 094026 (2021). https://doi.org/10.1103/PhysRevD.104.094026. arXiv:2102.09118 [hep-ph]
16. T. Affolder et al. (CDF), Phys. Rev. Lett. 85, 2886 (2000). https://doi.org/10.1103/PhysRevLett.85.2886. arXiv:hep-ex/0004027

17. R. Aaij et al. (LHCb), Eur. Phys. J. C 73, 2631 (2013). https://doi.org/10.1140/epjc/s10052-013-2631-3. arXiv:1307.6379 [hep-ex]

18. B. Abelev et al. (ALICE), Phys. Rev. Lett. 108, 082001 (2012). https://doi.org/10.1103/PhysRevLett.108.082001. arXiv:1111.1630 [hep-ex]

19. S. Acharya et al. (ALICE), Eur. Phys. J. C 78, 562 (2018). https://doi.org/10.1140/epjc/s10052-018-6027-2. arXiv:1805.04374 [hep-ex]

20. C. Gerschel, J. Hufner, Phys. Lett. B 207, 253 (1988). https://doi.org/10.1016/0370-2693(88)90570-9

21. M. Gonin et al. (NA50), Nucl. Phys. A 610, 404C (1996). https://doi.org/10.1016/S0375-9474(96)00373-9

22. R.L. Thews, M. Schroedter, J. Rafelski, Phys. Rev. C 63, 054905 (2001). https://doi.org/10.1103/PhysRevC.63.054905. arXiv:hep-ph/0007323

23. L. Grandchamp, R. Rapp, Nucl. Phys. A 709, 415 (2002). https://doi.org/10.1016/S0375-9474(02)01027-8. arXiv:hep-ph/0205305

24. L. Yan, P. Zhuang, N. Xu, Phys. Rev. Lett. 97, 232301 (2006). https://doi.org/10.1103/PhysRevLett.97.232301. arXiv:nucl-th/0608010

25. X. Du, R. Rapp, Nucl. Phys. A 943, 147 (2015). https://doi.org/10.1016/j.nuclphysa.2015.09.006. arXiv:1504.00670 [hep-ph]

26. X. Zhao, R. Rapp, Phys. Rev. C 82, 064905 (2010). https://doi.org/10.1103/PhysRevC.82.064905. arXiv:1008.5328 [hep-ph]

27. K. Zhou, N. Xu, Z. Xu, P. Zhuang, Phys. Rev. C 89, 054911 (2014). https://doi.org/10.1103/PhysRevC.89.054911. arXiv:1401.5845 [nucl-th]

28. J. Zhao, B. Chen, P. Zhuang, Phys. Rev. C 105, 034902 (2022). https://doi.org/10.1103/PhysRevC.105.034902. arXiv:2112.00293 [hep-ph]

29. B. Chen, Chin. Phys. C 43, 124101 (2019). https://doi.org/10.1088/1674-1137/43/12/124101. arXiv:1811.11393 [nucl-th]

30. Y. Liu, Z. Qu, N. Xu, P. Zhuang, J. Phys. G 37, 075110 (2010). https://doi.org/10.1088/0954-3899/37/7/075110. arXiv:0907.2723 [nucl-th]

31. M. He, B. Wu, R. Rapp, Phys. Rev. Lett. 128, 162301 (2022). https://doi.org/10.1103/PhysRevLett.128.162301. arXiv:2111.13528 [nucl-th]

32. D.Y.A. Villar, J. Zhao, J. Aichelin, P.B. Gossiaux, Phys. Rev. C 107, 054913 (2023). https://doi.org/10.1103/PhysRevC.107.054913. arXiv:2206.01308 [nucl-th]

33. S. Acharya et al. (ALICE), Phys. Rev. Lett. 131, 042303 (2023). https://doi.org/10.1103/PhysRevLett.131.042303. arXiv:2204.10171 [nucl-ex]

34. B.L. Ioffe, D.E. Kharzeev, Phys. Rev. C 68, 061902 (2003). https://doi.org/10.1103/PhysRevC.68.061902. arXiv:hep-ph/0306176

35. P. Faccioli, C. Lourenco, J. Seixas, H.K. Wohri, Eur. Phys. J. C 69, 657 (2010). https://doi.org/10.1140/epjc/s10052-010-1420-5. arXiv:1006.2738 [hep-ph]

36. J. Zhao, K. Zhou, S. Chen, P. Zhuang, Prog. Part. Nucl. Phys. 114, 103801 (2020). https://doi.org/10.1016/j.ppnp.2020.103801. arXiv:2005.08277 [nucl-th]

37. T. Matsui, H. Satz, Phys. Lett. B 178, 416 (1986). https://doi.org/10.1016/0370-2693(86)91404-8

38. D. Lafferty, A. Rothkopf, Phys. Rev. D 101, 056010 (2020). https://doi.org/10.1103/PhysRevD.101.056010. arXiv:1906.00035 [hep-ph]

39. H. Satz, J. Phys. G 32, R25 (2006). https://doi.org/10.1088/0954-3899/32/3/R01. arXiv:hep-ph/0512217

40. M.E. Peskin, Nucl. Phys. B 156, 365 (1979). https://doi.org/10.1016/0550-3213(79)90199-8

41. G. Bhanot, M.E. Peskin, Nucl. Phys. B 156, 391 (1979). https://doi.org/10.1016/0550-3213(79)90200-1

42. B. Schenke, S. Jeon, C. Gale, Phys. Rev. C 82, 014903 (2010). https://doi.org/10.1103/PhysRevC.82.014903. arXiv:1004.1408 [hep-ph]

43. S. McDonald, C. Shen, F. Fillion-Gourdeau, S. Jeon, C. Gale, Phys. Rev. C 95, 064913 (2017). https://doi.org/10.1103/PhysRevC.95.064913. arXiv:1609.02958 [hep-ph]

44. L. Adamczyk et al. (STAR), Phys. Rev. Lett. 118, 212301 (2017). https://doi.org/10.1103/PhysRevLett.118.212301. arXiv:1701.06060 [nucl-ex]

45. M.L. Miller, K. Reygers, S.J. Sanders, P. Steinberg, Annu. Rev. Nucl. Part. Sci. 57, 205 (2007). https://doi.org/10.1146/annurev.nucl.57.090506.123020. arXiv:nucl-ex/0701025

46. R. Aaij et al. (LHCb), JHEP 06, 147 (2017). https://doi.org/10.1007/JHEP06(2017)147. arXiv:1610.02230 [hep-ex]

47. S. Acharya et al. (ALICE). (2023). arXiv:2303.13361 [nuclex]

48. S. Acharya et al. (ALICE), JHEP 02, 041 (2020). https://doi.org/10.1007/JHEP02(2020)041. arXiv:1909.03158 [nucl-ex]

49. A.H. Mueller, J.-W. Qiu, Nucl. Phys. B 268, 427 (1986). https://doi.org/10.1016/0550-3213(86)90164-1

50. J.W. Cronin, H.J. Frisch, M.J. Shochet, J.P. Boymond, R. Mermod, P.A. Piroue, R.L. Sumner, Phys. Rev. D 11, 3105 (1975). https://doi.org/10.1103/PhysRevD.11.3105

51. J. Hufner, Y. Kurihara, H.J. Pirner, Phys. Lett. B 215, 218 (1988). https://doi.org/10.1016/0370-2693(88)91423-2

52. I. Helenius, K.J. Eskola, H. Honkanen, C.A. Salgado, JHEP 07, 073 (2012). https://doi.org/10.1007/JHEP07(2012)073. arXiv:1205.5359 [hep-ph]

53. R.L. Workman et al. (Particle Data Group), PTEP 2022, 083C01 (2022). https://doi.org/10.1093/ptep/ptac097

54. Z.-T. Liang, X.-N. Wang, Phys. Lett. B 629, 20 (2005). https://doi.org/10.1016/j.physletb.2005.09.060. arXiv:nucl-th/0411101

55. F. Becattini, L. Csernai, D.J. Wang, Phys. Rev. C 88, 034905 (2013). https://doi.org/10.1103/PhysRevC.88.034905. arXiv:1304.4427 [nucl-th]. [Erratum: Phys. Rev. C 93, 069901 (2016)]

56. Y.-G. Yang, R.-H. Fang, Q. Wang, X.-N. Wang, Phys. Rev. C 97, 034917 (2018). https://doi.org/10.1103/PhysRevC.97.034917. arXiv:1711.06008 [nucl-th]

57. V. Voronyuk, V.D. Toneev, W. Cassing, E.L. Bratkovskaya, V.P. Konchakovski, S.A. Voloshin, Phys. Rev. C 83, 054911 (2011). https://doi.org/10.1103/PhysRevC.83.054911. arXiv:1103.4239 [nucl-th]

58. G. Inghirami, L. Del Zanna, A. Beraudo, M.H. Moghaddam, F. Becattini, M. Bleicher, Eur. Phys. J. C 76, 659 (2016). https://doi.org/10.1140/epjc/s10052-016-4516-8. arXiv:1609.03042 [hep-ph]

59. L. Yan, X.-G. Huang, Phys. Rev. D 107, 094028 (2023). https://doi.org/10.1103/PhysRevD.107.094028. arXiv:2104.00831 [nucl-th]

60. Z. Wang, J. Zhao, C. Greiner, Z. Xu, P. Zhuang, Phys. Rev. C 105, L041901 (2022). https://doi.org/10.1103/PhysRevC.105. L041901. arXiv:2110.14302 [hep-ph]

61. X.-L. Sheng, L. Oliva, Q. Wang, Phys. Rev. D 101, 096005 (2020). https://doi.org/10.1103/PhysRevD.101.096005. arXiv:1910.13684 [nucl-th]. [Erratum: Phys. Rev. D 105, 099903 (2022)]