J/ψ elliptic and triangular flow in Pb-Pb collisions at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV

The inclusive J/ψ elliptic (v2) and triangular (v3) flow coefficients measured at forward rapidity (2.5 < y < 4) and the v2 measured at midrapidity (|y| < 0.9) in Pb-Pb collisions at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV using the ALICE detector at the LHC are reported. The entire Pb-Pb data sample collected during Run 2 is employed, amounting to an integrated luminosity of 750 μb−1 at forward rapidity and 93 μb−1 at midrapidity. The results are obtained using the scalar product method and are reported as a function of transverse momentum pT and collision centrality. At midrapidity, the J/ψ v2 is in agreement with the forward rapidity measurement. The centrality averaged results indicate a positive J/ψ v3 with a significance of more than 5σ at forward rapidity in the pT range 2 < pT< 5 GeV/c. The forward rapidity v2, v3, and v3/v2 results at low and intermediate pT (pT ≲ 8 GeV/c) exhibit a mass hierarchy when compared to pions and D mesons, while converging into a species-independent curve at higher pT. At low and intermediate pT, the results could be interpreted in terms of a later thermalization of charm quarks compared to light quarks, while at high pT, path-length dependent effects seem to dominate. The J/ψ v2 measurements are further compared to a microscopic transport model calculation. Using a simplified extension of the quark scaling approach involving both light and charm quark flow components, it is shown that the D-meson vn measurements can be described based on those for charged pions and J/ψ flow.


Introduction
Ultra-relativistic heavy-ion collisions are the means to create under laboratory conditions the deconfined state of strongly-interacting matter called quark-gluon plasma (QGP). This state behaves like an ideal fluid with a shear viscosity to entropy ratio approaching the conjectured lowest possible value of /(4πk B ) [1][2][3]. One of the most important observables for studying the properties of the QGP is the azimuthal dependence of particle production, also called anisotropic flow, quantified in terms of a Fourier expansion with respect to the azimuthal angle of the initial state symmetry plane for the n-th harmonic Ψ n as where v n is the n-th order harmonic coefficient and ϕ is the azimuthal angle of the particles. The initial state spatial anisotropy of the collision overlap region is transformed into a momentum anisotropy of the produced final state particles [4][5][6][7]. The medium response to the initial state anisotropy (ε n ), which is transformed into the v n coefficients, strongly depends on the macroscopic properties of the fireball, like the temperature dependent equation of state and the shear and bulk viscosity.
The dominant source of anisotropy is the ellipsoidal shape of the overlap region in non-central collisions that have a non-zero finite impact parameter (transverse distance separating the centers of the two nuclei), which gives rise to a large second order harmonic coefficient, v 2 , also known as elliptic flow. Fluctuations in the initial energy-density profile within the overlap region are thought to be the origin of the triangular flow, v 3 [8][9][10]. Higher order harmonics are strongly damped, do not depend linearly on the initial JHEP10(2020)141 anisotropy, and have significant contributions from the interplay of lower order harmonics [11][12][13][14][15]. The ALICE Collaboration published extensive studies of anisotropic flow measurements for identified light and strange particles [16,17]. Flow coefficients for all particles show, in the low p T range, an increasing trend with p T mainly attributed to the radial hydrodynamic expansion of the QGP, reach a maximum in the p T range 3-5 GeV/c depending on the particle mass and species, and finally drop towards higher p T . The behavior in the high p T region is commonly attributed to path-length dependent effects like energy loss [18][19][20]. At both RHIC and LHC energies, an approximate scaling of the flow coefficients with the number of valence quarks is observed for light and strange particles [16,17,[21][22][23]. In the low to moderate p T range (approximately 3 < p T < 8 GeV/c), this scaling is hypothesized to be the consequence of the hadronization process via quark coalescence and of a common underlying partonic flow during the hydrodynamic stage of the collision [24][25][26][27][28].
The production of charmonia, and especially of J/ψ, is one of the first proposed probes of the QGP properties, in particular the deconfinement [29]. Since charm quarks are produced during the early hard partonic collisions, they experience the entire evolution of the fireball. At the same time, their initial production cross section can be calculated in perturbative quantum chromodynamics (QCD). The suppression of the production of bound charmonium states by the free color charges of the dense deconfined medium is sensitive to both the medium bulk characteristics [30,31] and to the microscopic ones, like the charm-quark diffusion coefficient [32,33]. Measurements of the J/ψ nuclear modification factor R AA at RHIC in Au-Au collisions at √ s NN = 200 GeV [34] indicated a strong nuclear suppression especially for the most central collisions. At the LHC, in Pb-Pb collisions at √ s NN = 2.76 and 5.02 TeV, the ALICE Collaboration reported a much larger R AA compared to the one observed at RHIC [35][36][37], despite the higher energy density present in the system. This effect is concentrated in the low-p T region, which is consistent with charmonium regeneration by recombination of charm quarks, either at the QGP phase boundary via statistical hadronization [38] or continuously throughout the fireball evolution [39][40][41].
Within the statistical hadronization scenario, charm quarks thermalize in the QGP and all of the charmed bound hadrons are created at the phase boundary assuming chemical equilibration [38,42], except a small fraction created in the fireball corona that escape the medium. In transport model approaches, where charm quarks reach only a partial thermalization, roughly 50% of the produced J/ψ originate from the recombination process, while the rest comes from primordial production [39][40][41]. In both phenomenological approaches, it is expected that charm quarks will inherit some of the medium radial and anisotropic flow. Indeed, a significant D-meson [43][44][45] and J/ψ elliptic flow [46][47][48][49] was already observed at the LHC, indicating a hierarchy between the flow of charged particles, D and J/ψ mesons, with the J/ψ flow being the smallest. A positive J/ψ v 2 observed also at high p T , typically underestimated by transport model calculations, might suggest the presence of important path length dependent effects like energy loss and the survival probability in the medium [50,51]. In addition to v 2 , the ALICE Collaboration also published in ref. [48] an evidence of a positive J/ψ v 3 with a statistical significance of 3.7σ.
In this paper, the measurements of inclusive J/ψ v 2 and v 3 at forward rapidity (2.5 < y < 4) and v 2 at midrapidity (|y| < 0.9) in Pb-Pb collisions at √ s NN = 5.02 TeV are discussed. Inclusive J/ψ mesons include both a prompt component from direct J/ψ production and decays of excited charmonium states and a non-prompt component from weak decays of beauty hadrons. The results are presented as a function of p T in several collision centrality classes, expressed in percentages of the total hadronic cross section, and are compared with calculations from a microscopic transport model. The analyzed data include the full LHC Run 2 Pb-Pb data set, which improves the statistical precision with respect to the previous results by approximately a factor of two at forward rapidity [48], and a factor of nine (four) in central (semi-central) collisions at midrapidity [47] allowing the experimental evidence of a statistically significant non-zero J/ψ v 2 at midrapidity.

Experimental setup, data samples and event selection
A detailed description of the ALICE apparatus and its performance can be found in refs. [52,53]. At forward rapidity, J/ψ are reconstructed in the µ + µ − decay channel with the muon spectrometer which covers the pseudorapidity range −4 < η < −2.5. 1 The spectrometer includes five tracking stations, each composed of two planes of cathode pad chambers. The third station is placed inside a dipole magnet with a 3 Tm field integral. Two trigger stations, containing two planes of resistive plate chambers each, provide single and dimuon triggers with a programmable single-muon p T threshold. A front absorber, made of carbon, concrete, and steel, is placed in between the primary interaction point (IP) and the first tracking station to remove primary hadrons from the collision. A second absorber, made of iron, is placed in front of the trigger chambers to further reject secondary hadrons escaping the front absorber and low-p T muons, mainly from pion and kaon decays. An additional conical absorber surrounds the beam pipe to protect the muon spectrometer against secondary particles produced by the interaction of large-η particles with the beam pipe. At midrapidity, J/ψ mesons are reconstructed in the e + e − decay channel using the Inner Tracking System (ITS) [54] and the Time Projection Chamber (TPC) [55] in the rapidity range |y| < 0.9. The ITS is a cylindrical-shaped detector, consisting of 6 layers of silicon detectors used for precision tracking, reconstruction of the primary vertex of the event and event selection. The innermost two layers consists of pixels (SPD), the middle two are drift (SDD), while the two outermost layers are equipped with strip detectors (SSD). The tracklets, track segments reconstructed as pairs of hits in the SPD layers pointing to the primary vertex, are used for the determination of the event flow vector. The TPC is the main detector used for tracking and particle identification and consists of a cylindrical-shaped gas-filled active volume placed around the ITS. Radially, it extends between an inner radius of 0.85 m and an outer radius of 2.5 m, with a total length of 5 m along the beam axis. Particle identification in the TPC is performed via the measurement of the specific energy loss, dE/dx.
Besides the muon spectrometer and the central barrel detectors, a set of detectors for global event characterization are also used. Two arrays of 32 scintillator counters each, covering 2.8 < η < 5.1 (V0A) and −3.7 < η < −1.7 (V0C) [56], are used for triggering, 1 In the ALICE reference frame, the muon spectrometer covers a negative η range and consequently a negative y range. We have chosen to present our results with a positive y notation, due to the symmetry of the collision system. beam induced background rejection, and for the determination of the collision centrality. The 32 channels are arranged in four concentric rings with full azimuthal coverage allowing for the calculation of the event flow vector. The centrality of the events, expressed in fractions of the total inelastic hadronic cross section, is determined via a Glauber fit to the V0 amplitude as described in refs. [57,58]. In addition, two neutron Zero Degree Calorimeters [59], installed at ±112.5 m from the nominal IP along the beam axis, are used to remove beam induced background events and electromagnetic interactions.
The analyzed data samples were collected by ALICE during the 2015 and 2018 LHC Pb-Pb runs at √ s NN = 5.02 TeV using different trigger strategies for the forward muon spectrometer and the midrapidity detectors. At forward rapidity, data were collected requiring the coincidence of the minimum bias (MB) and unlike-sign dimuon triggers. The former is defined by the coincidence of signals in the V0A and V0C arrays while the latter requires at least a pair of opposite-sign track segments in the muon trigger stations. The programmable threshold of the muon trigger algorithm was set so that the trigger efficiency for muon tracks with p T = 1 GeV/c is 50% and reaches a plateau value of about 98% at p T ≈ 2.5 GeV/c. In order to study the background, additional samples of single muon and like-sign dimuon events were also collected by requiring, in addition to the MB condition and the low-p T threshold, at least one or a pair of same-sign track segments in the trigger system, respectively.
At midrapidity, data were collected using the MB trigger during the 2015 data taking period, and the MB, central, and semi-central triggers in the 2018 period. The central and semi-central triggers require the MB trigger to be fired but, in addition, a condition on the total signal amplitude in the V0 detectors, corresponding to collision centralities of 0-10% and 30-50%, respectively, was applied.
Both forward and midrapidity analyses require to have a primary vertex position within ±10 cm from the nominal IP along the beam axis. Events containing more than one collision (pile-up) are removed by exploiting the correlations between the number of clusters in the SPD, the number of reconstructed SPD tracklets, and the total signal in the V0A and V0C detectors. At midrapidity, events with pile-up occurring during the drift time of the TPC are rejected in the offline analysis based on the correlation between the number of SDD and SSD clusters and the total number of clusters in the TPC. The beam-induced background is filtered out offline by applying a selection based on the V0 and the ZDC timing information [60].
The integrated luminosity of the analyzed data samples is about 750 µb −1 for the dimuon analysis. For the measurements at midrapidity, the total luminosity recorded depends on the centrality range due to the centrality triggers, and amounts to 93 µb −1 , 41 µb −1 , and 20 µb −1 for the central, semi-central, and MB triggers, respectively.

Data analysis
The v n coefficients are obtained using the scalar product (SP) method [2,61]. This is a two-particle correlation technique based on the scalar product between the unit flow vector for a given harmonic n, u n = e inϕ , of the particle of interest (here a dilepton) and the -4 -

Dilepton analysis
Three sub-event technique, detectors used Corresponding gap between e + e − |y ee | < 0.9 V0C TPC V0A |∆η| >0.8 Table 1. Summary of the details concerning the dimuon and dielectron analyses, corresponding to the forward and midrapidity region, respectively. The detectors cited in this table are described in section 2, and the details concerning the three sub-event technique is presented in section 3.
complex conjugate of the event flow vector in a subdetector A, Q A * n . The flow coefficients are thus defined as where Q B n and Q C n are the n-th harmonic event flow vectors measured in two additional subdetectors, B and C, respectively, which are used to correct the event flow vector via the three sub-event technique [62]. The star ( * ) represents the complex conjugate and the bracket . . .
indicates the average over dileptons from all events in a given p T range, dilepton invariant mass (m ), and centrality interval. The brackets . . . in the denominator denote the average over all events in a narrow centrality interval containing the event under consideration. The V0A and V0C detectors are used in the analysis at both rapidities, while the analysis at forward rapidity uses the SPD as the third subdetector, and the analysis at midrapidity uses the TPC. As detector A, the SPD is chosen for the forward analysis and the V0C for the midrapidity one. The V0A and V0C event flow vectors are calculated using the energy deposition measured in the individual channels. For the SPD and TPC event flow vectors, the reconstructed tracklets and the tracks are used, respectively.
The effects of non-uniform acceptance of the detectors used for the flow vector determination are corrected through the procedure described in ref. [63]. As was discussed in section 2, the three detectors used for the event flow determination cover distinct pseudorapidity ranges, allowing for pseudorapidity gaps ∆η between the sub-events used for flow vector determination and dilepton reconstruction. The pseudorapidity gap between u n and Q A n , corresponding to |∆η| > 1.1 and |∆η| > 0.8 for the dimuon and dielectron analysis, respectively, suppresses the short-range correlations originating from resonance decays or jets (non-flow effects), not related to the global azimuthal anisotropy.
In the dimuon analysis, J/ψ candidates are formed by combining pairs of oppositesign tracks reconstructed in the geometrical acceptance of the muon spectrometer using the tracking algorithm described in ref. [64]. The same single-muon and dimuon selection criteria used in previous analyses [48,65] are applied. Namely, each muon track candidate should have −4 < η µ < −2.5, a radial transverse position at the end of the front absorber in the range 17.6 < R abs < 89.5 cm, and must match a track segment in the muon trigger chambers above the 1 GeV/c p T threshold. The rapidity of the muon pair should be within the acceptance of the muon spectrometer (2.5 < y < 4.0).

JHEP10(2020)141
At midrapidity, J/ψ mesons are reconstructed in the dielectron decay channel. Electron candidates are required to be good quality tracks matched in both the ITS and the TPC, and to have a p T > 1 GeV/c and |η| < 0.9. Tracks are selected to have at least 70 space points in the TPC, out of a maximum of 159, and a χ 2 /N dof < 2 for the track fit quality. At least one hit in either of the two SPD layers is required to reject secondary electrons from photons converted in the detector material and to improve the tracking resolution. Secondary electrons are further rejected by requiring the distance-of-closestapproach (DCA) to the collision vertex to be smaller than 1 cm and 3 cm in the transverse and longitudinal directions, respectively. Electrons are identified via their specific energy loss in the TPC gas, dE/dx, by selecting a band of ±3σ around the expectation value, with σ being the dE/dx measurement resolution. To reduce further the hadronic contamination, candidate tracks compatible within ±3.5σ with the pion or proton hypothesis are rejected.
The flow coefficients are extracted from sequential fits to the dilepton invariant mass distribution, m , and the v n as a function of m , which include the superposition of a J/ψ signal and a background contribution, using the function Here, v J/ψ n denotes the J/ψ v 2 or v 3 and α(m ) is the signal fraction defined as S/(S + B). The latter is extracted from fits to the dilepton invariant mass distribution as described below. The v bkg n (m ) corresponds to the dilepton background v 2 or v 3 . In the dimuon analysis, the J/ψ signal is parameterized using an extended Crystal Ball (CB2) function and the background with a Variable Width Gaussian (VWG) function [66]. In the fit, the J/ψ peak position and width are left free, while the CB2 tail parameters are fixed to the values reported in ref. [67]. The signal of the ψ(2S) is not included in the fit of the v n coefficients due to its marginal significance. At midrapidity, the signal fraction is obtained from the dielectron invariant mass distribution in two steps. First, the combinatorial background is estimated using an event mixing technique, where pairs are built from different events with similar collision centrality, flow-vector orientation, and longitudinal position of the event vertex, and then subtracted from the same-event dielectron invariant mass distribution. The combinatorial background normalization is obtained from the ratio of the number of same-event to mixed-event like-sign pairs. Second, the remaining distribution is fitted using a component for the signal and one for the residual background. For the J/ψ signal shape, the dielectron invariant mass distribution obtained from Monte Carlo simulations is used. The residual background, originating mainly from semileptonic decays of cc and bb pairs (correlated background) and imperfect matching between the same-event and mixed-event distributions, is parameterized using either a third order polynomial function at low p T or an exponential function at high p T .
The v n extraction method employed in this work is described in detail in ref. [48], where the v bkg n (m ) distribution is obtained using an event mixing technique. There, it was first demonstrated that the flow coefficients of the background can be obtained from where v (1) n (ϕ (1) ) and v (2) n (ϕ (2) ) are the flow coefficients (azimuthal angles) of the two leptons, respectively, and ϕ is the dilepton azimuthal angle. The brackets · · · m denote an average over all dileptons belonging to the given m interval. Here, it is worth to note that the denominator in eq. Here, u n and u (2) n are the unit vector of the two leptons, Q Examples of fits to the invariant mass distribution (top panels corresponding to a1, b1, c1) and to v n (m ) (bottom panels related to a2 for v 2 (m µµ ), b2 for v 2 (m ee ), and c2 for v 3 (m µµ )) are shown in figure 1 for the dimuon and dielectron analyses. The background, which is mostly combinatorial, especially in central events, is well reproduced with the event mixing technique. In the absence of correlated background, the background flow v bkg n is directly given by the mixed-event flow. At forward rapidity, the effect of the unknown flow contribution of the correlated background and residual mismatches between the same-event and mixed-event background flow, is considered as a systematic uncertainty and is discussed in section 4. In the default approach, the flow of the correlated background is assumed to be negligible, and thus the denominator of eq. (3.3) is given by the ratio N bkg between the number of background unlike-sign dileptons N bkg +− and the number of unlikesign dileptons from mixed events N mix +− , which is obtained after a proper normalization involving like-sign dileptons as described in ref. [48]. At midrapidity, due to the smaller signal-to-background ratio, the difference between mixed and same event background flow is taken into account by considering in the fit function an additional term which accounts for the flow of the correlated background and imperfections of the mixed event procedure. This term is parameterized using a second order polynomial and acts as a correction to the background flow obtained from the mixed event procedure.

Systematic uncertainties
The systematic uncertainties related to the v n extraction procedure, the track and event selection criteria, residual detector effects, and non-flow contributions are evaluated as described below and summarized in table 2. A quadratic sum of the systematic uncertainties from the independent sources is used as final systematic uncertainty on the measurements.
In the dimuon analysis, the signal fraction α(m µµ ) is estimated by fitting the invariant mass distribution with standard signal and background functions. The systematic uncertainty on the determination of α(m µµ ) is estimated by varying the signal and background functions, as well as the mass fit range. For the signal, in addition to a CB2, a pseudo-Gaussian with a mass-dependent width [66] is also used. The tail parameters were fixed to the values obtained in Monte Carlo simulations or in other analyses with better signal significance [37,67]. For the background, the VWG function was changed to a fourth order Chebyshev polynomial. The invariant mass fit range is varied from the standard -8 -JHEP10(2020)141  The non-uniformity in the detector acceptance could lead to a residual effect in the calibration of the event flow vector Q n . The cross-term products of the event flow vector, Q x,A × Q y,B , are evaluated to verify that values are negligible compared to the linear products. In addition, possible impacts on the v n are checked by calculating the crossterm products between the components of the Q n vector and the unitary vector u n of the J/ψ candidates. No clear p T or centrality dependence is found for this contribution, and the corresponding systematic uncertainty is estimated to be less than 1%. Additional uncertainties related to the calculation of the reference flow vector are evaluated as the difference between the event flow factor R n obtained using MB events or dimuon-triggered events. For the dimuon analysis it amounts to 1% for R 2 and up to 3% for R 3 .
The variation of the J/ψ reconstruction efficiency with the local occupancy of the detector could bias the measured v n . At forward rapidity, this effect is evaluated using azimuthally isotropic simulated J/ψ → µ + µ − decays embedded into real Pb-Pb events. A maximum effect of 0.002 for v 2 and 0.001 for v 3 is observed in non-central collisions with no clear p T dependence. At midrapidity, the strongest dependence of reconstruction performance on the local detector occupancy is caused by the TPC particle identification (PID). A data driven study, using a clean electron sample from photon conversions, shows that the largest variation of the TPC electron PID response between the region along the event flow vector and the region orthogonal to it is approximately 2% of the dE/dx resolution. This leads to a decrease of the observed v 2 by less than 1% and is thus neglected.
The presence of a correlated background and its unknown flow contribution can affect the v n extraction. The contribution of the correlated background to the flow of the background can be introduced in eq. (3.3) by replacing the denominator N bkg +− /N mix +− by -9 -

JHEP10(2020)141
N bkg +− /(N mix +− + β(N bkg +− − N mix +− )), where β represents the relative strength of the correlated background flow with respect to the combinatorial background flow. The systematic uncertainty is defined as the difference between the default fit, equivalent to β = 0, and the modified fit with β left as a free parameter. This uncertainty is, as expected, negligible for central collisions and low p T but becomes significant for peripheral collisions and high-p T . The estimated systematic uncertainty for the v 2 and v 3 extraction reaches a maximum of about 0.01 for peripheral collisions and at high-p T .
In the dielectron analysis, the signal-to-background ratio can vary significantly depending on the TPC electron identification selection and centrality, which may impact the J/ψ v 2 fits. Thus, the v 2 was extracted for a set of nine electron PID cuts where both the electron selection and the hadron rejection were varied such that the J/ψ efficiency is changed by approximately 50%. The RMS of the v 2 obtained from all of these selections is assigned as a systematic uncertainty, which ranges between 0.010 and 0.023 depending on the centrality and p T interval, while the average value is taken as central value. In addition, the fit range of the v 2 (m ee ) is varied by either making it narrower or wider, but no significant systematic effects are observed.

Results and discussions
The J/ψ elliptic flow coefficient measured by ALICE in Pb-Pb collisions at √ s NN = 5.02 TeV at forward and central rapidity is shown in figure 2 as a function of p T , for the centrality intervals 0-10%, 10-30%, 30-50% and 0-50%. Systematic uncertainties, obtained as described in the previous section, are shown as boxes around the data points, while the statistical uncertainties are shown as error bars. Here, and in all figures as a function of p T , the J/ψ data points are located at the average p T of the reconstructed J/ψ. These results are compared with the midrapidity v 2 measurements for charged pions by ALICE [17] and prompt D mesons by ALICE [68] and CMS [43]. At forward rapidity and for all centrality intervals, the J/ψ v 2 values increase with p T , possibly reaching a maximum at intermediate values of p T , and decreasing or saturating towards high p T . Also, the J/ψ v 2 values increase when decreasing centrality from the 0-10% to 10-30%, then to 30-50%. This behavior is qualitatively similar to the one for light hadrons and D mesons. The J/ψ v 2 measurement at midrapidity is statistically compatible to the one at forward rapidity in both centrality intervals within uncertainties. Considering all the midrapidity data points as statistically independent measurements, it was found that the J/ψ v 2 is larger than zero with a significance of approximately 2.5 standard deviations in both centrality intervals. It is worth to remark that the ALICE apparatus is undergoing an ambitious upgrade programme in preparation of Runs 3 and 4 of the LHC that will enable the separation of the prompt and non-prompt J/ψ contributions to the measured flow coefficients, thus providing valuable information on both the charmonium and open beauty hadron production dynamics.
As also noted previously [48], a clear mass hierarchy of the v 2 values is seen in the lowp T region (p T < 6 GeV/c) for the light hadrons and D mesons measured at midrapidity and inclusive J/ψ, with the J/ψ exhibiting the lowest elliptic flow. Here, it is important to  note that in the considered η range, the η dependence of the v 2 at a given p T is expected to be negligible, as shown by the CMS measurement for charged particles [69], albeit in a somewhat narrower η range. At high-p T (p T > 8 GeV/c), the v 2 coefficients from all species converge into a single curve suggesting that, in this kinematic range, the anisotropy for all particles arises dominantly from path-length dependent energy-loss effects [70]. However, in the case of the much heavier J/ψ, one may also consider that the hydrodynamic flow, which arises from a common velocity field, still contributes significantly even at high p T , as can be expected from the particle mass dependence of the p T range where the flow reaches its maximum.
In figure 3, the p T -dependent inclusive J/ψ triangular flow coefficient measured at forward rapidity is shown in each of the considered centrality intervals. For most of the centrality and p T intervals, the measured inclusive J/ψ v 3 is positive and with no significant centrality dependence. In the 0-50% centrality range, the triangular flow coefficient is larger than zero (0.0250 ± 0.0045 (stat.) ± 0.0020 (syst.) in 2 < p T < 5 GeV/c) corresponding to a significance of 5.1σ, calculated adding quadratically the statistical and systematic uncertainties. The positive v 3 indicates that the initial state energy-density fluctuations, the dominant source of v 3 , are reflected also in the anisotropic flow of charm quarks. Also shown in figure 3 are similar measurements for charged pions [17] and D mesons [43,68]   midrapidity. The mass hierarchy observed for v 2 holds also in the case of v 3 . Together with the J/ψ v 2 , these observations provide a strong support for the hypothesis of charm quark being, at least partially, kinetically equilibrated in the dense and deconfined QGP medium.
The ratio of the triangular to elliptic flow coefficients, v 3 /v 2 , as a function of p T is shown in the left panel of figure 4 for the inclusive J/ψ at forward rapidity, D mesons and charged pions at midrapidity. In this ratio, the statistical uncertainties are considered to be uncorrelated due to the weak correlation between the orientation of the Q 2 and Q 3 flow vectors [71], while the systematic uncertainties related to α(m µµ ) and to the reconstruction efficiency discussed in section 4, cancel in the ratio. The same hierarchy observed for the individual v 2 and v 3 measurements is also observed in the v 3 /v 2 ratio, which suggests that higher harmonics are damped faster for heavy quarks than for the light ones. At RHIC [72,73] and LHC [74,75], it was observed that the flow coefficients of light particles from different harmonics follow a power-law scaling as v  Calculations from a transport model [39,40] are also shown.
In figure 5, the inclusive J/ψ v 2 as a function of p T in the 20-40% centrality interval is compared with the microscopic transport calculations by Du et al. [39,40]. In this model, the J/ψ are created both from the primordial hard partonic interactions but also from the recombination of thermalized charm quarks in the medium, which accounts for roughly 50% of all J/ψ at low p T . Non-prompt J/ψ mesons, created in the weak decays of beauty hadrons, are also included in the model. The amplitude of the inclusive J/ψ v 2 in the calculations is in good agreement with the experimental measurements for p T < 4 GeV/c. However, the overall trend of the model calculation does not describe the data well, especially in the intermediate p T range, 4 < p T < 10 GeV/c, where the J/ψ flow JHEP10(2020)141 is largely underestimated. The primordial J/ψ component, which is sensitive mainly to path length dependent effects, like survival probability, exhibits a monotonically increasing trend from low towards high p T , with this mechanism becoming the dominant source of the anisotropic flow for p T larger than 8 GeV/c. Path length dependent energy loss, widely seen as a major source of anisotropy at large p T , is not implemented for J/ψ mesons in this calculation. It is worth noting that this model provides a qualitative good description of the centrality and transverse momentum of the J/ψ nuclear modification factor [37,76]. Figure 6 shows the centrality dependence of the inclusive J/ψ v 2 (top panels) and v 3 (bottom panels) for a low-p T interval (0 < p T < 5 GeV/c) on the left, and a high-p T one (5 < p T < 20 GeV/c) on the right. Here, due to the large integrated p T range, the v n coefficients are corrected for the J/ψ acceptance and efficiency A × ε. Each dimuon pair is weighted using the inverse of the p T and y dependent A × ε factor before filling the invariant mass and v n (m µµ ) distributions. The p T and y dependent A×ε map was obtained from the embedded simulations discussed in section 4. Any possible systematic uncertainty related to the A × ε corrections is already included in the systematic uncertainty due to the dependence of the reconstruction efficiency on the local detector occupancy. The J/ψ results are compared with the flow coefficients of charged pions for a p T value similar to the corrected J/ψ p T , published by ALICE in ref. [17]. In addition, the ratio v π 2 /v J/ψ 2 is computed and shown in the bottom sub-panels. Both at low p T (1.75 < p T < 2 GeV/c) and high p T (6 < p T < 7 GeV/c), the v 2 of π ± increases from central to semi-central collisions, reaching a maximum at 40-50% centrality, and then decreases towards peripheral collisions. For the J/ψ at low p T , while the centrality trend is qualitatively similar, the maximum (or even saturation) of v 2 seems to be reached for more central collisions than for the pions. This is more clearly emphasized by the increasing trend of the ratio v π 2 /v J/ψ 2 , from central to peripheral collisions, which deviates from unity by a significance of 8.5σ. In the framework of transport models, this could be understood by the increasing fraction of regenerated J/ψ at low p T when moving from peripheral to central collisions. Alternatively, and independently of the regeneration scenario, the increase of the v π 2 /v J/ψ 2 from central to peripheral collisions, could also be understood in terms of partial or later thermalization of the charm quarks compared to light quarks. The decrease in energy density and lifetime of the system is counterbalanced by the increase of the initial spatial anisotropy towards peripheral collisions. The v 2 of the J/ψ will therefore reach its maximum at more central collisions compared to light particles because charm quarks require larger energy densities to develop flow [33,[77][78][79]. At high p T , J/ψ mesons and charged pions seem to exhibit the same centrality dependence, although the v 2 coefficients are systematically lower for the J/ψ mesons than for the pions. Such a similar centrality dependence could indicate a similar mechanism at the origin of the flow for both J/ψ mesons and pions at high p T .
The centrality dependence of the v 3 coefficient at low p T is less pronounced than that of the v 2 for both pions and J/ψ, as expected since initial state fluctuations only weakly depend on centrality. Also, the J/ψ v 3 is smaller relative to the one of charged pions, in both the p T intervals considered.
The flow of light and strange particles was shown to approximately scale with the number of constituent quarks (NCQ scaling) at both RHIC and LHC energies [80,81].
Systematic uncertainty (uncorrelated) Figure 6. The inclusive J/ψ v 2 and v 3 as function of the centrality of the collision, at forward rapidity, for the low-p T range 0 < p T < 5 GeV/c (left panel) and high-p T range 5 < p T < 20 GeV/c (right panel). The results are compared to the v n coefficients of midrapidity π ± [17] at low and high-p T corresponding to 1.75 < p T < 2 GeV/c and 6 < p T < 7 GeV/c, respectively. The ratio of midrapidity π ± v 2 to inclusive J/ψ v 2 is also shown.
This was typically interpreted to arise naturally in hadronization scenarios based on quark coalescence in which the flow of bound mesons and baryons depends solely on the collective flow of light and strange quarks (assumed to be identical) and the number of valence quarks [24,28]. In the case of charmed hadrons, the NCQ scaling assuming a flavor independent flow would obviously not work due to the large observed differences between the flow of light-flavor particles, D and J/ψ mesons. However, this scaling can be extended by assuming that the much heavier charm quark has a different flow magnitude [25] and that it can be derived from the flow of the J/ψ via the usual NCQ formula, v Then it is straightforward to show that the flow of the D meson can be constructed as the sum of the flow coefficients for light and charm quarks as where p q T and p c T are the p T of the light and charm quarks, respectively, corresponding to the D-meson p T , p D T . The light quark flow is obtained by interpolating the measured charged pions flow using v π n (p π T ) = 2 · v q n (p π T /2). Figure 7 shows a comparison of the Dmeson v 2 and v 3 as a function of p T , derived assuming the above described procedure, to the D-meson v n measured by the CMS Collaboration [43].
The red dashed curves show fits to the J/ψ v n employing an ad-hoc function, a third order polynomial at low-p T and a linear function at high-p T , used to extract the flow  of charm quarks needed to obtain the scaled D-meson flow according to eq. (5.1). The scaled D-meson flow is found to be very sensitive to the fraction of p T carried by each of the constituent quarks. In coalescence-like models, constituent quarks must have equal velocities which leads to a sharing of the D-meson p T proportional to the effective quark masses. This implies that by far the largest fraction of p T should be carried by the charm quark. Based on the simplistic and naive approach described here, a p T sharing between light and charm quarks [25,82] where the ratio p q T /p D T = 0.2 (black curve), is clearly disfavored by the data. Surprisingly, it was found that a good description of the D-meson flow measurement, as illustrated by the blue curves in figure 7, is obtained when the light quark carries a relatively large fraction of the D-meson p T (dark blue and green curves). The best agreement with the D-meson data of the CMS Collaboration [43] is obtained when the light-quark p T fraction has a value of p q T /p D T = 0.4 (dark blue curve), but a rather good description of the data is observed also when assuming that the light and charm quarks share equally the D-meson p T (green curve). Within uncertainties, the scaling seems to work well for both v 2 and v 3 over the entire covered p T range and in all centrality intervals. In summary, the inclusive J/ψ v 2 at forward and midrapidity and the J/ψ v 3 at forward rapidity were measured in Pb-Pb collisions at √ s NN = 5.02 TeV using the scalar product method. In non-central collisions, the J/ψ v 2 values are found to be positive up to the last interval corresponding to 12 < p T < 20 GeV/c and reach a maximum of approximately 0.1 around a p T of 5 GeV/c. The J/ψ v 3 values at forward rapidity reach 0.04 around a p T of 4 GeV/c and are positive in the 0-50% centrality interval for 2 < p T < 5 GeV/c with a significance of 5.1σ. The mass hierarchy observed for v 2 , v 2,π > v 2,D > v 2,J/ψ , seems to also hold in the case of v 3 and will be the subject of more detailed studies with the Run 3 and Run 4 data. At high p T , the v 2 for all particles converge to similar values, suggesting that path-length dependent effects become dominant there. The measured J/ψ v 3 /v 2 ratios exhibits the same hierarchy indicating that higher harmonics are damped faster for charmonia compared to lighter particles. The p T -integrated v 2 coefficient in a low and a high-p T region is in both cases dependent on centrality and reaches a maximum value of about 0.1, while the v 3 has no clear centrality dependence. Both J/ψ p T -integrated v 2 and v 3 coefficients, either at low-p T or at high-p T are found to be lower than the ones of charged pions at a p T similar to the J/ψ average p T . At low p T , the ratio of the charged pions v 2 to those of p T -integrated J/ψ increase from central to peripheral collisions, compatible with a scenario in which charm quarks thermalize later than the light ones. At high p T , this ratio is compatible with unity without any statistically significant centrality dependence.
Using an extension of the well known number of constituent quark scaling, the measured charged pion and J/ψ v n can be used as proxies in order to derive the D-meson v 2 and v 3 as a combination of the flow of light and charm quarks. Within this procedure, it is surprising to observe that the measured D-meson v 2 and v 3 can be described if one considers that the light and charm quarks share similar fractions of the D-meson p T , which is counterintuitive in a coalescence approach. The fact that such a simple scaling works suggests that the flow of charmonia and open charm mesons can be effectively explained assuming a common underlying charm quark flow in addition to the flow of light quarks.