# Charmonium Triangular Flow in High Energy Nuclear Collisions

Jiaxing Zhaoa Baoyi Chenb Pengfei Zhuanga aPhysics Department, Tsinghua
University, Beijing 100084, China  bDepartment of Physics, Tianjin University,
Tianjin 300350, China

###### Abstract

We calculate, for the first time, $J/\psi$ triangular flow in high energy
nuclear collisions. The charmonium motion in medium is controlled by a
transport equation with loss and gain terms, and the evolution of the hot
medium is governed by a single-shot hydrodynamic equation with a triangularly
deformed initial condition. In comparison with the elliptic flow $v_{2}$, the
triangular flow $v_{3}$ comes almost completely from the charmonium
regeneration in the quark gluon plasma and therefore is more sensitive to the
medium properties.

##  I Introduction

Quarks and gluons can be deconfined at high temperature and form a new phase
of strongly interacting matter, called quark-gluon plasma (QGP) Aoki:2006we ;
Bazavov:2011nk . The only way to generate the new state of matter in
laboratory is through high energy nuclear collisions STAR:2005gfr ;
PHENIX:2004vcz ; ALICE:2010suc . From the study on the collisions at the
Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC), one of
the characteristic observables to describe the collective properties of the
QGP medium is the anisotropic flow of the final state hadrons
Ollitrault:1992bk , which is defined as the coefficients $v_{n}$ of a Fourier
decomposition of the hadron distribution Voloshin:1994mz ,

$\displaystyle{dN\over d^{2}{\bm{p}}_{T}dy}={1\over\pi
dp_{T}^{2}dy}\left[1+2\sum_{n=1}^{\infty}v_{n}\cos\left(n\left(\phi-\Psi_{n}\right)\right)\right],$
(1)

where $y,{\bm{p}}_{T}$ and $\phi$ are respectively the hadron longitudinal
rapidity, transverse momentum and azimuth angle, and $\Psi_{n}$ is the $n$-th
harmonic symmetry plane angle. The momentum anisotropy of the produced hadrons
in the final state originates from the spatial anisotropy of the collision
overlap region in the initial state Ollitrault:1992bk ; Voloshin:1994mz ;
Qiu:2011iv ; Luzum:2008cw . For instance, the elliptic flow $v_{2}$ is a
response to the ellipsoidal shape of the overlap region in non-central
nucleus-nucleus collisions, and the triangular flow $v_{3}$ is mostly from the
initial fluctuation of the energy density in the overlap region Alver:2010dn ;
Alver:2010gr ; Teaney:2010vd . On the other hand, the collective flow develops
in the medium and is therefore sensitive to the QGP properties, like the
equation of state and shear and bulk viscosity McDonald:2016vlt ; Shen:2015qta
.

Charmonia, the bound state of $c$ and $\bar{c}$ quarks, have long been
considered as a sensitive signature of the QGP formation in nuclear collisions
Matsui:1986dk . Different from light hadrons which are all produced through
the decay of the medium at the phase transition boundary Andronic:2006ky ,
there are two sources for charmonium production: one is the initial production
through hard processes which is then largely suppressed in the medium
Gerschel:1988wn ; NA50:1996lag , and the other is the regeneration in the
medium through coalescence mechanism Thews:2000rj ; Grandchamp:2002wp ;
Andronic:2003zv ; Yan:2006ve ; Zhao:2017yan . As colorless particles, the
initially created quarkonia interact weakly with the QGP and their collective
flow is very small. On the contrary, colored charmed quarks interact strongly
with the medium and the regenerated charmonia can carry sizeable flow
Zhou:2014kka ; Chen:2018kfo . The competition between the suppression and
regeneration can explain well Zhao:2010nk ; Zhao:2011cv ; Liu:2009wza ;
Zhou:2014kka ; Du:2015wha ; Chen:2018kfo  the experimentally measured $J/\psi$
nuclear modification factor $R_{AA}$ and elliptic flow $v_{2}$ at RHIC and LHC
energies, especially in the low momentum region.

Recently, the inclusive $J/\psi\ v_{3}$ at forward rapidity and $v_{2}$ at
middle and forward rapidity are measured in Pb-Pb collisions at colliding
energy $\sqrt{s_{NN}}$= 5.02 TeV by ALICE Collaboration at LHC ALICE:2020pvw ;
ALICE:2018bdo . The first measurement of the charmonium triangular flow
provides us an opportunity to probe the initial energy density fluctuation of
the medium via heavy flavors and to understand the interaction mechanism of
heavy quarks with the medium. In this paper we study the $J/\psi$ triangular
flow in the frame of a transport approach.

##  II Medium evolution

We first discuss the initial energy density of the medium. Since the medium is
dominated by light partons, the feedback from charmonia to medium can be
safely neglected, and the medium can be considered as a background for the
charmonium motion. In this case we can take a well established hydrodynamic
model for light hadron production to describe the evolution of the background.
The initial energy profile of the medium without fluctuations can be
factorized as Schenke:2010nt

$\displaystyle\epsilon({\bm{x}},\tau_{0}|{\bm{b}})=\epsilon_{0}f(\eta)W({\bm{x}}_{T}|{\bm{b}})/W(0|0),$
(2)

where $\epsilon_{0}$ is the maximum energy density at the center of the
fireball created in central collisions, which is related to the maximum
temperature $T_{0}$ via the equation of state and can be fixed by matching to
the experimentally measured charged hadron multiplicity Schenke:2010nt . In
our numerical calculation $T_{0}$ is taken to be $510$ MeV for Pb-Pb
collisions at $\sqrt{s_{NN}}$=5.02 TeV. The longitudinal and transverse space
dependence is characterized by the distributions $f(\eta)$ and
$W({\bm{x}}_{T}|{\bm{b}})$, where $\eta=1/2\ln[(t+z)/(t-z)]$ is the space-time
rapidity and we take $f(\eta)$ as the one in Ref. Schenke:2010nt . The initial
state of matter in the transverse plane is normally described by two
components Kharzeev:2000ph ; Bozek:2010wt ; Schenke:2010nt , the number of
wounded nucleons $n_{\text{w}}$ and the number of binary nucleon-nucleon
collisions $n_{\text{b}}$, which respectively control the initial soft and
hard processes. While how the deposited energy density or entropy density
precisely scales with the two numbers is not clear from the first principle,
we know that the soft process is the main source of the initial energy
deposition at SPS energy, and the contribution from the hard process increases
at RHIC and LHC energy. In general, the spatial dependence of the initial
energy distribution in the transverse plane can be parameterized as
Schenke:2010nt ,

$\displaystyle
W({\bm{x}}_{T}|{\bm{b}})=(1-\gamma)n_{\text{w}}({\bm{x}}_{T}|{\bm{b}})+\gamma
n_{\text{b}}({\bm{x}}_{T}|{\bm{b}}),$ (3)

where ${\bm{x}}_{T}$ is the transverse coordinate of the colliding nucleon,
${\bm{b}}$ the impact parameter of the nuclear collisions, and the parameter
$\gamma$ the fraction of binary collisions. At LHC energy, one takes
$\gamma=0.16$ for Pb-Pb collisions at $\sqrt{s_{NN}}$=5.02 TeV Bozek:2010wt .
The numbers of wounded nucleons and binary collisions are controlled by the
nuclear geometry,

$\displaystyle n_{\text{w}}({\bm{x}}_{T}|{\bm{b}})$ $\displaystyle=$
$\displaystyle
T_{A}({\bm{x}}_{T}-{\bm{b}}/2)\left(1-e^{-\sigma_{\text{inel}}T_{B}({\bm{x}}_{T}+{\bm{b}}/2)}\right)$
$\displaystyle+$ $\displaystyle
T_{B}({\bm{x}}_{T}+{\bm{b}}/2)\left(1-e^{-\sigma_{\text{inel}}T_{A}({\bm{x}}_{T}-{\bm{b}}/2)}\right),$
$\displaystyle n_{\text{b}}({\bm{x}}_{T}|{\bm{b}})$ $\displaystyle=$
$\displaystyle\sigma_{\text{inel}}T_{A}({\bm{x}}_{T}+{\bm{b}}/2)T_{B}({\bm{x}}_{T}-{\bm{b}}/2),$
(4)

where $T_{A}$ and $T_{B}$ are the thickness functions Miller:2007ri , and
$\sigma_{\text{inel}}$ is the inelastic scattering cross-section between
nucleons. For Pb-Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV one takes
$\sigma_{\text{inel}}=68\ \text{mb}$ ALICE:2012fjm .

In relativistic heavy ion collisions, the initial condition of the medium
evolution fluctuates event by event, and the fluctuations can be simulated
through the Monte Carlo Glauber model (MC-Glauber) Miller:2007ri  and Monte
Carlo fKLN model (MC-KLN) Kharzeev:2001gp ; Drescher:2006ca . The QCD-based
IP-Glasma model can also describe well the initial stage of the medium
Schenke:2012wb . With these initial conditions including event-by-event
fluctuations, the anisotropic flows for charged hadrons measured at RHIC and
LHC are explained very well in the frame of relativistic viscous hydrodynamics
followed by a transport cascade in hadronic phase McDonald:2016vlt ;
Shen:2015qta .

In principle, aiming to describe the $J/\psi$ triangular flow $v_{3}$, one
needs to evolve the event-by-event hydrodynamics with fluctuating initial
conditions McDonald:2016vlt ; Qiu:2011iv . However, to simplify the
calculations, one can take, as an approximation, a single-shot hydrodynamics
where the initial condition is given by an averaged, smooth but deformed
profile Qiu:2011iv ; Alver:2010dn ; Alver:2010gr . We follow the way shown in
Ref. Alver:2010dn  to add a deformation factor in the energy density through
the transformation

$\displaystyle\epsilon({\bm{x}},\tau_{0}|{\bm{b}})\to\epsilon(\tilde{\bm{x}},\tau_{0}|{\bm{b}})$
(5)

with the deformed coordinate
$\tilde{\bm{x}}=(x_{T}\sqrt{1+\varepsilon_{3}\cos[3(\phi_{s}-\Psi_{3})]},\phi_{s},\eta)$,
where $x_{T}=\sqrt{x^{2}+y^{2}}$ and $\phi_{s}=\arctan(|y|/|x|)$ are the
transverse radius and azimuth angle, $\Psi_{3}$ is the reference angle which
can be the reaction plane angle in smooth hydrodynamics (we take $\Psi_{3}=0$
in the following calculation). The magnitude of the deformation
$\varepsilon_{3}$, which is called triangularity and depends on the collision
centrality, can be given by the MC-Glauber or MC-KLN model Alver:2010dn ;
Qiu:2011iv . We take $\varepsilon_{3}$= 0.08, 0.1 and 0.2 for Pb-Pb collisions
at $\sqrt{s_{NN}}$=5.02 TeV with impact parameters $b=$ 3.2, 6.8 and 9.6 fm,
corresponding to the centrality bin 0-10%, 10-30% and 30-50% Alver:2010dn ;
Qiu:2011iv . The initial energy density profile with and without triangular
deformation is shown in Fig. 1.

Figure 1: The initial energy density $\epsilon(x,y)$ in the transverse plane
for Pb-Pb collisions at $\sqrt{s_{NN}}=$5.02 TeV with impact parameters
$b=6.8$ fm. The left and right panels are the profiles without and with
triangular deformation.

The quark and gluon medium produced in high energy nuclear collisions is a
strongly coupled system and therefore can be locally thermalized fast. With
the fluctuating energy density (5) as the initial condition at time
$\tau_{0}=0.6$ fm/c Schenke:2010nt , the evolution of the medium can be
described by hydrodynamic equations. In this paper, we employ a
(2+1)-dimensional hydrodynamic model, the MUSIC package Schenke:2010nt ;
McDonald:2016vlt , to characterize the space and time dependence of the
temperature and velocity of the hot medium. To close the hydrodynamic
equations, the equation of state for both QGP and hadron phases is needed. We
will adopt ”s95p-v1” which matches Lattice QCD data at high temperature and
the hadron resonance gas at low temperature Huovinen:2009yb . The two phases
are connected with a smooth crossover. The critical temperature $T_{c}$
between QGP and hadron gas is taken to be 170 MeV. We choose an effective
shear viscosity $\eta/s=0.08$ Bernhard:2016tnd  and a zero bulk viscosity.

##  III Charmonium transport

We now focus on the charmonium motion in hot medium. Since quarkonia are
colorless and very heavy, they are unlikely thermalized with the medium.
Therefore, their phase space distribution
$f_{\psi}({\bm{p}},{\bm{x}},\tau|{\bm{b}})$ for
$\psi=J/\psi,\chi_{c},\psi^{\prime}$ should be governed by a transport
equation including both initial production as well as regeneration. The
distribution satisfies the relativistic Boltzmann equation Zhao:2020jqu ,

$\displaystyle\left[\cosh(y-\eta)\partial_{\tau}+{\sinh(y-\eta)\over\tau}\partial_{\eta}+{\bm{v}}_{T}\cdot\nabla_{T}\right]f_{\psi}$
$\displaystyle=$ $\displaystyle-\alpha f_{\psi}+\beta,$ (6)

where $y=1/2\ln[(E+p_{z})/(E-p_{z})]$ is the charmonium momentum rapidity, and
${\bm{v}}_{T}={\bm{p}}_{T}/E_{T}$ is the charmonium transverse velocity with
transverse energy $E_{T}=\sqrt{m_{\psi}^{2}+{\bm{p}}_{T}^{2}}$. 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
direction. The anomalous suppression and regeneration mechanisms in QGP medium
are reflected in the loss term $\alpha$ and gain term $\beta$.

Charmonia in hot QGP medium suffers Debye screening Matsui:1986dk . With
increasing temperature, the interaction distance between a pair of heavy
quarks becomes shorter and shorter. When it is less than the charmonium size,
the charmonium is dissociated. Considering different sizes for the ground and
excited states of $c\bar{c}$, the dissociation temperatures are
$T_{d}\sim(2.3,1.2,1.1)T_{c}$ for $J/\psi$, $\chi_{c}$ and $\psi^{\prime}$.
This is the so-called sequential suppression Satz:2005hx . The above analyses
based on the Debye screening effect is typically based on the assumption of a
constant temperature in connection with a sharp transition of the inelastic
charmonium widths from zero below $T_{d}$ to infinity above $T_{d}$. However,
the volume of the produced fireball in relativistic heavy ion collisions is
relatively small and expands rapidly, implying rather fast temperature changes
and short fireball lifetimes. In this case, the conclusion from the static
Debye screening effect may deviate from the real system, and it becomes
essential to include the concrete interactions between partons and charmonia,
leading to sizable inelastic reaction rates comparable to the fireball
expansion or cooling rate. In particular, the charmonia can be destroyed below
$T_{d}$ and survivable above $T_{d}$. The Debye screening is still operative,
by controlling the binding energy which in turn determines the phase space and
thus the width of the dynamic dissociation reactions Satz:2005hx ;
Chen:2017jje . An important such process in the QGP is the gluon dissociation
process Peskin:1979va ; Bhanot:1979vb  $g\psi\to c\bar{c}$, an analogy to the
photon dissociation process of electromagnetic bound states. 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 Peskin:1979va ; Bhanot:1979vb . For the
excited states $\chi_{c}$ and $\psi^{\prime}$, the cross sections can be
obtained through their geometric relation to the ground state Chen:2018kfo .
Taking only the gluon dissociation as the loss term and its inverse process
$c\bar{c}\to g\psi$ as the gain term (both are $2\to 2$ processes), $\alpha$
and $\beta$ can be explicitly expressed as Zhao:2020jqu

$\displaystyle\alpha({\bm{p}},{\bm{x}},\tau|{\bf b})$ $\displaystyle=$
$\displaystyle{1\over
2E_{T}}\int{d^{3}{\bm{p}}_{g}\over(2\pi)^{3}2E_{g}}W_{g\psi}^{c\bar{c}}(s)f_{g}({\bm{p}}_{g},{\bm{x}},\tau)$
$\displaystyle\times\Theta(T({\bm{x}},\tau|{\bm{b}})-T_{c}),$
$\displaystyle\beta({\bm{p}},{\bm{x}},\tau|{\bm{b}})$ $\displaystyle=$
$\displaystyle{1\over
2E_{T}}\int{d^{3}{\bm{p}}_{g}\over(2\pi)^{3}2E_{g}}{d^{3}{\bm{p}}_{c}\over(2\pi)^{3}2E_{c}}{d^{3}{\bm{p}}_{\bar{c}}\over(2\pi)^{3}2E_{\bar{c}}}$
(7) $\displaystyle\times
W_{c\bar{c}}^{g\psi}(s)f_{c}({\bm{p}}_{c},{\bm{x}},\tau|{\bm{b}})f_{\bar{c}}({\bm{p}}_{\bar{c}},{\bm{x}},\tau|{\bm{b}})$
$\displaystyle\times(2\pi)^{4}\delta^{(4)}(p+p_{g}-p_{c}-p_{\bar{c}})$
$\displaystyle\times\Theta(T({\bm{x}},\tau|{\bm{b}})-T_{c}),$

where $E_{g},E_{c}$ and $E_{\bar{c}}$ are the gluon, charm quark and anti-
charm quark energies, ${\bm{p}}_{g},{\bm{p}}_{c}$ and ${\bm{p}}_{\bar{c}}$ are
their momenta, and $s$ is the $g\psi$ interaction energy in their center-of-
mass frame. $W_{g\psi}^{c\bar{c}}$ is the dissociation probability, and
$W_{c\bar{c}}^{g\psi}$ the regeneration probability which can be obtained from
$W_{g\psi}^{c\bar{c}}$ through considering detailed balance between the two
processes Yan:2006ve ; Liu:2009wza . The step function $\Theta$ is used to
guarantee the calculation in QGP phase above the critical temperature $T_{c}$,
and the local temperature of the medium comes from the hydrodynamic
calculation discussed in Section II. Since gluons are constituents of the
medium, their distribution $f_{g}$ is the Bose-Einstein distribution. Again
the local temperature and velocity of the medium in the distribution are from
the hydrodynamics. Charm quarks are not constituents of the medium, their
motion in the medium should in principle be controlled by a transport
equation, and the distribution is in between two limits: the pQCD limit
without interaction with the medium and the equilibrium limit with strong
interaction with the medium. The experimentally measured large charmed meson
flow indicates that charm quarks with small transverse momentum might be
thermalized with the medium STAR:2017kkh ; ALICE:2020pvw . As a first
approximation, we take a thermally equilibrated distribution for charm (anti-
charm) quarks $f_{c}=\rho_{c}N_{c}/(e^{p_{c}^{\mu}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 is governed by the
nuclear geometry of the colliding system,

$\rho_{c}({\bm{x}},\tau_{0}|{\bm{b}})={T_{A}(\tilde{\bm{x}}_{T}+{\bf
b}/2)T_{B}(\tilde{\bm{x}}_{T}-{\bm{b}}/2)\cosh\eta\over\tau_{0}}{d\sigma_{pp}^{c\bar{c}}\over
d\eta},$ (8)

where $d\sigma^{c\bar{c}}_{pp}/d\eta$ is the rapidity distribution of charm
quark production cross section in p-p collisions. We take 0.56 mb
$<d\sigma^{c\bar{c}}_{pp}/d\eta<$ 0.78 mb in the forward rapidity of Pb-Pb
collisions at $\sqrt{s_{NN}}=5.02$ TeV LHCb:2016ikn ; ALICE:2021dhb . Since
charm quarks are produced through initial hard scattering processes controlled
by the number of binary collisions $n_{\text{b}}$, their initial density
should also be triangularly deformed. This is introduced via the change from
${\bm{x}}$ to $\tilde{\bm{x}}$ in $\rho_{c}$.

Besides the hot nuclear matter effect which affects charmonium motion through
the above discussed suppression and regeneration terms ($\alpha$ and $\beta$),
there is also the cold nuclear matter effect which changes the initial
condition of the transport equation (III). If the cold effect is completely
omitted, the initial condition is just a superposition of the charmonium
distribution in p-p collisions which can be parameterized as Zhou:2014kka ;
Chen:2018kfo ; Zhao:2017yan ,

$\displaystyle{d^{2}\sigma^{J/\psi}_{pp}\over 2\pi p_{T}dp_{T}dy}$
$\displaystyle=$ $\displaystyle{n-1\over n-2}{1\over\pi\langle
p_{T}^{2}\rangle_{pp}}\left(1+{1\over n-2}{p_{T}^{2}\over\langle
p_{T}^{2}\rangle_{pp}}\right)^{-n}{d\sigma^{J/\psi}_{pp}\over dy},$ (9)

where the mean transverse momentum square $\langle p_{T}^{2}\rangle_{pp}=$
12.5 (GeV/c)2, the parameter $n=3.2$, and the rapidity distribution
$d\sigma^{J/\psi}_{pp}/dy=$ 3.25 $\mu$b in forward rapidity $2.5<|y|<4$ of Pb-
Pb collisions at $\sqrt{s_{NN}}$=5.02 TeV are extracted from the experimental
data ALICE:2012vup ; ALICE:2012vpz . The cold nuclear matter effect on
charmonium production includes mainly nuclear shadowing Mueller:1985wy ,
Cronin effect Cronin:1974zm ; Hufner:1988wz  and nuclear absorption
Gerschel:1988wn . For heavy ion collisions at LHC energy, the collision time
is very short, the nuclear absorption can be safely neglected. Before two
gluons fuse into a charmonium, they 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.
Therefore, when we do the superposition of the p-p distribution (III), the
averaged transverse momentum square $\langle p_{T}^{2}\rangle_{pp}$ should be
replaced by

$ \langle p_{T}^{2}\rangle_{pp}+a_{gN}l, $ (10)

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 Pb-Pb collisions at
$\sqrt{s_{NN}}$=5.02 TeV. In practice, we take a Gaussian smearing
Huefner:2002tt ; Zhao:2010nk ; Zhou:2014kka  for the modified transverse
momentum distribution. The shadowing effect comes from the parton collective
motion in a nucleus Norton:2003cb  and can be parameterized as a modification
of the parton distribution function $f_{i}(x,\mu_{F})$ in a free nucleon by a
factor $R_{i}=\bar{f}_{i}/(Af_{i})$, where $\bar{f}_{i}(x,\mu_{F})$ is the
distribution for parton $i$ in a nucleus, and the factorization scale is taken
as $\mu_{F}=\sqrt{m_{\psi}^{2}+p_{T}^{2}}$. The modification factor can be
taken from the EPS09 package Helenius:2012wd . Including both the shadowing
and Cronin effects, the initial charmonium distribution can be written as
Zhou:2014kka

$\displaystyle $ $\displaystyle f_{\psi}({\bm{p}},{\bm{x}},\tau_{0}|{\bm{b}})$
$\displaystyle=$ $\displaystyle{(2\pi)^{3}\over E_{T}\tau_{0}}\int
dz_{A}dz_{B}\rho_{A}(\tilde{\bm{x}}_{T}+{{\bm{b}}/2},z_{A})\rho_{B}(\tilde{\bm{x}}_{T}-{{\bm{b}}/2},z_{B})$
(11)
$\displaystyle\times\mathcal{R}_{g}(x_{1},\mu_{F},\tilde{\bm{x}}_{T}+{\bm{b}}/2)\mathcal{R}_{g}(x_{2},\mu_{F},\tilde{\bm{x}}_{T}-{\bm{b}}/2)$
$\displaystyle\times
f_{pp}^{\psi}(\tilde{\bm{x}}_{T},{\bm{p}},z_{A},z_{B}|{\bm{b}}), $

where $\rho_{A}$ and $\rho_{B}$ are the nucleon distributions in the two
colliding nuclei, the charmonium momentum distribution $f_{pp}^{\psi}$ is with
the Cronin-effect modified averaged transverse momentum square (10). The local
shadowing effect is embedded in the factor $\mathcal{R}_{g}$. Again, the
triangular deformation is considered in the initial distributions via the
replacement of ${\bm{x}}_{T}$ by the deformed coordinate $\tilde{\bm{x}}_{T}$.

##  IV Numerical result

The transport equation (III) with the initial distribution (III) can be solved
analytically Yan:2006ve ; Liu:2009wza . To compare the result with the
experimentally measured inclusive $J/\psi$ data, one should consider not only
the above discussed direct production but also the feed-down from the excited
charmonium states and the B-decay contribution. The feed-down branch ratio
from $\chi_{c}$ and $\psi^{\prime}$ to $J/\psi$ are taken as 22% and 61%
ParticleDataGroup:2020ssz , and the momentum dependence of the B-decay
fraction in p-p collisions can be fitted as $f_{B}(p_{T})=N_{pp}^{B\to
J/\psi}/(N_{pp}^{\text{prompt}}+N_{pp}^{B\to
J/\psi})=0.04+0.023p_{T}/(\text{GeV/c})$ which depends on rapidity and
colliding energy weakly Chen:2013wmr . We take the same fraction in heavy ion
collisions. For bottom quarks, their initial momentum distribution can be
given by FONNL FONLL  and the evolution in QGP can be simulated by the
Langevin equation with a spatial diffusion coefficient $2\pi TD_{s}=2$
Zhao:2020jqu . After the evolution bottom quarks will carry an anisotropic
flow $v_{n}^{\text{B}}(n=2,3)$. With both prompt (direct + feed-down) and non-
prompt (B-decay) contributions, the $J/\psi$ anisotropic flow $v_{n}$ in the
final state can be expressed as,

$\displaystyle
v_{n}(p_{T})=v_{n}^{\text{prompt}}(p_{T})\left(1-f_{B}(p_{T})\right)+v_{n}^{\text{B}}(p_{T})f_{B}(p_{T}).
$ (12)

Figure 2: The $J/\psi$ elliptic flow $v_{2}$ as a function of transverse
momentum $p_{T}$ in different centrality bins in Pb-Pb collisions at
$\sqrt{s_{NN}}$= 5.02 TeV. Dotted and dashed lines are calculations with only
initial production and regeneration, slash and dark bands are total results
without and with B-decay contribution, and the lower and higher limits of the
bands correspond to charm quark cross section $d\sigma^{c\bar{c}}_{pp}/d\eta=$
0.56 and 0.78 mb. The experimental inclusive data are from ALICE Collaboration
ALICE:2020pvw . Figure 3: The $J/\psi$ triangular flow $v_{3}$ as a function
of transverse momentum $p_{T}$ in different centrality bins in Pb-Pb
collisions at $\sqrt{s_{NN}}$= 5.02 TeV. Dotted and dashed lines are
calculations with only initial production and regeneration, slash and dark
bands are total results without and with B-decay contribution, and the lower
and higher limits of the bands correspond to charm quark cross section
$d\sigma^{c\bar{c}}_{pp}/d\eta=$ 0.56 and 0.78 mb. The experimental inclusive
data are from ALICE Collaboration ALICE:2020pvw .

The numerical calculations for $J/\psi\ v_{2}$ and $v_{3}$ and the comparison
with inclusive $J/\psi$ data in Pb-Pb collisions at $\sqrt{s_{NN}}$= 5.02 TeV
are shown in Figs. 2 and 3. Dotted and dashed lines are the calculations with
only initial production and regeneration, slash and dark bands are the total
results for prompt and inclusive $J/\psi$s, and the higher and lower limits of
the bands correspond to the charm quark cross section
$d\sigma_{pp}^{c\bar{c}}/d\eta=0.78$ and $0.56$ mb. Considering the fact that,
the initially produced charmonia via hard processes carry high momentum, and
the regeneration via coalescence of thermalized charm quarks happens at low
momentum, their contributions to charmonium production are mainly in the
higher and lower momentum regions. Since bottom quark is too heavy, its decay
contribution to $J/\psi$ yield is small, and its interaction with the medium
is much weaker than a charm quarks which leads to a rather small
$v_{n}^{\text{B}}$. That is the reason why the difference between the slash
and dark bands is very small. Our calculated $v_{2}$ agrees reasonably well
with the experimental data at $p_{T}\lesssim 3$ GeV/c, the sizeable $v_{2}$ at
high momentum may come from charm quark interaction with the magnetic field
created in non-central collisions Fukushima:2015wck  and the fragmentation
mechanism which dominates the high-momentum charmonium production Bain:2016rrv
; Kang:2017yde .

Different from the elliptic flow $v_{2}$ which is initiated from the
anisotropic nuclear geometry and thus increases from central to peripheral
collisions, the source of the triangular flow $v_{3}$ is the initial
fluctuation of the medium which is not directly related to the collision
centrality. This difference leads to the following three characteristics of
the $J/\psi$ triangular flow: 1) $v_{3}$ is sensitive to the hot medium
formation, but its dependence on centrality is weak provided the fireball is
formed; 2) The initial production is before the medium formation and the
following interaction between colorless charmonia and the medium is weak, the
contribution from initial production to $v_{3}$ is extremely small and $v_{3}$
is almost fully from the regeneration; 3) High $p_{T}$ charmonia are not from
the regeneration in medium and therefore not sensitive to the initial medium
fluctuation, $v_{3}$ approaches to zero fast at high momentum.

##  V Summary

Charmonium collective flow is often used as a probe of the hot medium
properties in high energy nuclear collisions. We calculated $J/\psi$ elliptic
flow $v_{2}$ and triangular flow $v_{3}$ in the frame of a transport approach
for charmonium motion and an initially triangularly deformed hydrodynamics for
the medium evolution. Different from $v_{2}$ which is originated from the
anisotropic nuclear geometry in the beginning of the collisions, $v_{3}$ comes
from the fluctuations of the initial condition of the hot medium and therefore
is more sensitive to the medium properties. While there are different
mechanisms for charmonium production in different momentum regions, like
regeneration at low momentum and initial production and fragmentation at high
momentum, the medium-fluctuation controlled $v_{3}$ exists only in low and
intermediate momentum region and can be well described by the regeneration
mechanism.

Acknowledgement: The work is supported by the NSFC grant Nos. 11890712,
12047535, 12075129 and 12175165.

## References

  * (1) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Nature 443, 675-678 (2006).
  * (2) A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar, H. T. Ding, S. Gottlieb, R. Gupta, P. Hegde, U. M. Heller and F. Karsch, et al. Phys. Rev. D 85, 054503 (2012).
  * (3) J. Adams et al. [STAR], Nucl. Phys. A 757, 102-183 (2005).
  * (4) K. Adcox et al. [PHENIX], Nucl. Phys. A 757, 184-283 (2005).
  * (5) K. Aamodt et al. [ALICE], Phys. Rev. Lett. 105, 252302 (2010).
  * (6) J. Y. Ollitrault, Phys. Rev. D 46, 229-245 (1992).
  * (7) S. Voloshin and Y. Zhang, Z. Phys. C 70, 665-672 (1996).
  * (8) Z. Qiu and U. W. Heinz, Phys. Rev. C 84, 024911 (2011).
  * (9) M. Luzum and P. Romatschke, Phys. Rev. C 78, 034915 (2008) [erratum: Phys. Rev. C 79, 039903 (2009)]
  * (10) B. H. Alver, C. Gombeaud, M. Luzum and J. Y. Ollitrault, Phys. Rev. C 82, 034913 (2010).
  * (11) B. Alver and G. Roland, Phys. Rev. C 81, 054905 (2010) [erratum: Phys. Rev. C 82, 039903 (2010)].
  * (12) D. Teaney and L. Yan, Phys. Rev. C 83, 064904 (2011).
  * (13) S. McDonald, C. Shen, F. Fillion-Gourdeau, S. Jeon and C. Gale, Phys. Rev. C 95, no.6, 064913 (2017).
  * (14) C. Shen, Z. Qiu and U. Heinz, Phys. Rev. C 92, no.1, 014901 (2015).
  * (15) T. Matsui and H. Satz, Phys. Lett. B 178, 416-422 (1986).
  * (16) A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, Nucl. Phys. A 789, 334-356 (2007).
  * (17) C. Gerschel and J. Hufner, Phys. Lett. B 207, 253-256 (1988).
  * (18) M. Gonin et al. [NA50], Nucl. Phys. A 610, 404C-417C (1996).
  * (19) R. L. Thews, M. Schroedter and J. Rafelski, Phys. Rev. C 63, 054905 (2001).
  * (20) L. Grandchamp and R. Rapp, Nucl. Phys. A 709, 415-439 (2002).
  * (21) A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, Phys. Lett. B 571, 36-44 (2003).
  * (22) L. Yan, P. Zhuang and N. Xu, Phys. Rev. Lett. 97, 232301 (2006).
  * (23) J. Zhao and B. Chen, Phys. Lett. B 776, 17-21 (2018).
  * (24) K. Zhou, N. Xu, Z. Xu and P. Zhuang, Phys. Rev. C 89, no.5, 054911 (2014).
  * (25) B. Chen, Chin. Phys. C 43, no.12, 124101 (2019).
  * (26) X. Zhao and R. Rapp, Phys. Rev. C 82, 064905 (2010).
  * (27) X. Zhao and R. Rapp, Nucl. Phys. A 859, 114-125 (2011).
  * (28) Y. Liu, Z. Qu, N. Xu and P. Zhuang, J. Phys. G 37, 075110 (2010).
  * (29) X. Du and R. Rapp, Nucl. Phys. A 943, 147-158 (2015).
  * (30) S. Acharya et al. [ALICE], JHEP 02, 012 (2019).
  * (31) A. Polleri, T. Renk, R. Schneider and W. Weise, Phys. Rev. C 70, 044906 (2004).
  * (32) S. Acharya et al. [ALICE], JHEP 10, 141 (2020).
  * (33) B. Schenke, S. Jeon and C. Gale, Phys. Rev. C 82, 014903 (2010).
  * (34) D. Kharzeev and M. Nardi, Phys. Lett. B 507, 121-128 (2001).
  * (35) P. Bozek, M. Chojnacki, W. Florkowski and B. Tomasik, Phys. Lett. B 694, 238-241 (2011).
  * (36) M. L. Miller, K. Reygers, S. J. Sanders and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 57, 205-243 (2007).
  * (37) B. Abelev et al. [ALICE], Eur. Phys. J. C 73, no.6, 2456 (2013).
  * (38) D. Kharzeev and E. Levin, Phys. Lett. B 523, 79-87 (2001).
  * (39) H. J. Drescher and Y. Nara, Phys. Rev. C 75, 034905 (2007).
  * (40) B. Schenke, P. Tribedy and R. Venugopalan, Phys. Rev. Lett. 108, 252301 (2012).
  * (41) P. Huovinen and P. Petreczky, Nucl. Phys. A 837, 26-53 (2010).
  * (42) J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu and U. Heinz, Phys. Rev. C 94, no.2, 024907 (2016).
  * (43) J. Zhao, K. Zhou, S. Chen and P. Zhuang, Prog. Part. Nucl. Phys. 114, 103801 (2020).
  * (44) H. Satz, J. Phys. G 32, R25 (2006).
  * (45) S. Chen and M. He, Phys. Rev. C 96, no.3, 034901 (2017).
  * (46) M. E. Peskin, Nucl. Phys. B 156, 365-390 (1979).
  * (47) G. Bhanot and M. E. Peskin, Nucl. Phys. B 156, 391-416 (1979).
  * (48) L. Adamczyk et al. [STAR], Phys. Rev. Lett. 118, no.21, 212301 (2017).
  * (49) R. Aaij et al. [LHCb], JHEP 06, 147 (2017).
  * (50) S. Acharya et al. [ALICE], [arXiv:2105.06335 [nucl-ex]].
  * (51) B. Abelev et al. [ALICE], Phys. Lett. B 718, 295-306 (2012) [erratum: Phys. Lett. B 748, 472-473 (2015)].
  * (52) B. Abelev et al. [ALICE], JHEP 11, 065 (2012).
  * (53) A. H. Mueller and J. w. Qiu, Nucl. Phys. B 268, 427-452 (1986).
  * (54) J. Cronin, H. J. Frisch, M. Shochet, J. Boymond, R. Mermod, P. Piroue and R. L. Sumner, Phys. Rev. D 11, 3105-3123 (1975).
  * (55) J. Hufner, Y. Kurihara and H. Pirner, Acta Phys. Slov. 39, 281-287 (1989).
  * (56) J. Huefner and P. f. Zhuang, Phys. Lett. B 559, 193-200 (2003).
  * (57) P. Norton, Rept. Prog. Phys. 66, 1253-1297 (2003).
  * (58) I. Helenius, K. J. Eskola, H. Honkanen and C. A. Salgado, JHEP 07, 073 (2012).
  * (59) P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020).
  * (60) B. Chen, Y. Liu, K. Zhou and P. Zhuang, Phys. Lett. B 726, 725-728 (2013).
  * (61) M. Cacciari, M. Greco and P. Nason, JHEP 9805 (1998) 007. M. Cacciari, S. Frixione and P. Nason, JHEP 0103 (2001) 006.
  * (62) K. Fukushima, K. Hattori, H. U. Yee and Y. Yin, Phys. Rev. D 93, no.7, 074028 (2016).
  * (63) R. Bain, Y. Makris and T. Mehen, JHEP 11, 144 (2016).
  * (64) Z. B. Kang, J. W. Qiu, F. Ringer, H. Xing and H. Zhang, Phys. Rev. Lett. 119, no.3, 032001 (2017).