Study of J/$\psi$ azimuthal anisotropy at forward rapidity in Pb-Pb collisions at $\sqrt{{\textit s}_{\rm NN}}$ = 5.02 TeV

The second ($v_2$) and third ($v_3$) flow harmonic coefficients of J/$\psi$ mesons are measured at forward rapidity (2.5 $<$ $y$ $<$ 4.0) in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV with the ALICE detector at the LHC. Results are obtained with the scalar product method and reported as a function of transverse momentum, $p_{\rm{T}}$, for various collision centralities. A positive value of J/$\psi$ $v_3$ is observed with 3.7$\sigma$ significance. The measurements, compared to those of prompt D$^0$ mesons and charged particles at mid-rapidity, indicate an ordering with $v_{\rm n}$(J/$\psi$) $<v_{\rm n}$(D$^0$) $<v_{\rm n}$(h$^\pm$) (n = 2, 3) at low and intermediate $p_{\rm{T}}$ up to 6 GeV/$c$ and a convergence with $v_2$(J/$\psi$) $\approx v_2$(D$^0$) $\approx v_2$(h$^\pm$) at high $p_{\rm{T}}$ above 6-8 GeV/$c$. In semi-central collisions (5-40% and 10-50% centrality intervals) at intermediate $p_{\rm{T}}$ between 2 and 6 GeV/$c$, the ratio $v_3/v_2$ of J/$\psi$ mesons is found to be significantly lower (4.6$\sigma$) with respect to that of charged particles. In addition, the comparison to the prompt D$^0$-meson ratio in the same $p_{\rm{T}}$ interval suggests an ordering similar to that of the $v_2$ and $v_3$ coefficients. The J/$\psi$ $v_2$ coefficient is further studied using the Event Shape Engineering technique. The obtained results are found to be compatible with the expected variations of the eccentricity of the initial-state geometry.


Introduction
The study of collisions of ultra-relativistic heavy ions aims to characterize the Quark-Gluon Plasma (QGP), a strongly coupled state of matter comprising of deconfined quarks and gluons. One of the main features of heavy-ion collisions is the anisotropic particle flow [1,2]. It arises from initial collision geometry anisotropies being converted by the pressure gradients of the QGP medium to final-state particle momentum anisotropies. The anisotropic flow is described by the coefficients v n of a Fourier series decomposition of the azimuthal distribution of the produced particles [3] dN dϕ where ϕ is the azimuthal angle of the particle and Ψ n is the n-th harmonic symmetry plane angle. The dominant second-order flow coefficient (v 2 ) is called elliptic flow and mostly originates from the almondshaped overlap area between the colliding nuclei in non-central collisions. The third-order flow coefficient (v 3 ) is named triangular flow and is generated by fluctuations in the initial distribution of nucleons in the overlap region [4][5][6][7][8].
Heavy quarks, in particular their bound quark-antiquark states known as quarkonia, are important probes of the QGP. Heavy-quark pairs are created prior to the formation of the QGP through hard parton collisions and thus experience the full evolution of the system. Measurements of the J/ψ nuclear modification factor (R AA ) as a function of centrality in Pb-Pb collisions at the LHC [9][10][11] are reproduced by transport [12][13][14] and statistical hadronization [15,16] models including partial to full J/ψ (re)generation by recombination of thermalized charm quarks. Such (re)generation component is dominant at low transverse momentum (p T ) as shown by the comparison [11,17] of the R AA as function of p T with transport model calculations. In the case of the statistical hadronization model, the produced J/ψ reflects the dynamics of the charm quarks at the QGP phase boundary. The measured p T spectra seem to support this idea [18]. Measurements of the azimuthal anisotropies of J/ψ production in high-energy heavy-ion collisions can bring new important insights on the charm quark dynamics.
A recent measurement of the elliptic flow of J/ψ at forward rapidity in central and semi-central Pb-Pb collisions at the center of mass energy per nucleon pair of √ s NN = 5.02 TeV indicates a significant positive v 2 coefficient [19]. This result is compatible with the hypothesis of J/ψ production via recombination of thermalized c andc quarks from the QGP medium predominantly at low p T , but the magnitude and the transverse momentum dependence of the v 2 coefficient differ significantly from theoretical calculations [12][13][14]. Moreover, the v 2 coefficient is found to be quite significant at high p T , in contrast with the expectations of small azimuthal asymmetry originating mainly from path-length dependent J/ψ dissociation in the medium. Furthermore, a positive J/ψ v 2 coefficient at intermediate and high p T has been observed in p-Pb collisions [20,21], in which neither a significant contribution from charm-quark recombination nor sizable path-length effects are expected [22]. Recent measurements of D-meson azimuthal asymmetry in Pb-Pb collisions are interpreted as collective behavior of the charm quarks at low p T and path-length dependent charm-quark energy loss at high p T [23, 24].
Hydrodynamic calculations [25] show that v n ≈ κ n ε n for n = 2 and 3, where ε n is the eccentricity coefficient of the initial-state collision geometry. The parameters κ n encode the response of the QGP medium and depend on the particle type and mass as well as its transverse momentum. At low p T , the flow coefficients of light-flavoured particles increase with increasing p T [26,27]. This increase of v n coefficients as a function of p T depends of the particle mass and can be attributed to the radial expansion of the QGP medium. At 3-4 GeV/c, the flow coefficients reach a maximum. The position of the maximum, divided by the number of constituent quarks n q , does not dependent strongly on the particle mass as predicted by coalescence models [28]. Furthermore, the v n values at the maximum, divided by n q , are similar for all measured light-flavoured particles, with deviations of up to ±20% between mesons and baryons [27]. At high p T above 6-8 GeV/c, the observed azimuthal anisotropy of the final-state particles is believed to come from path-length dependent parton energy loss inside the QGP. Calculations [29] show that the corresponding v 2 and v 3 coefficients exhibit approximately linear dependence on ε 2 and ε 3 , respectively. Nevertheless, the correlation between the flow coefficients and the initial-state eccentricities is weaker with respect to the hydrodynamic case, especially between v 3 and ε 3 . Interestingly, the particle-mass dependence of v 2 and v 3 appears to be be strongly reduced in the ratio v 3 /v 2 in semi-central collisions for light-flavored particles [27]. Whether the above considerations also hold for heavy quarks and quarkonia is an open question whose answer could help to understand the origin of charm quark azimuthal anisotropies and characterize their interactions with the flowing medium.
In the present analysis, the J/ψ v 2 and v 3 coefficients as well as the ratio v 3 /v 2 as a function of the transverse momentum and the collision centrality are measured. Wherever possible, the data are compared to existing mid-rapidity charged-particle (predominantly π ± ) and prompt D 0 -meson results. In addition, the dependence of the J/ψ v 2 coefficient on the initial-state conditions is studied with the Event Shape Engineering (ESE) technique [30]. Fluctuations in the initial-state energy density distribution lead to event-by-event variations of the flow observed at a given centrality [31]. The ESE technique consists of selecting events with the same centrality but different flow and therefore initial-state geometry eccentricity [32,33]. Recently, the ESE technique has been applied to the measurement of mid-rapidity D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV [34]. The obtained results indicate a correlation between the D-meson azimuthal anisotropy and the flow of light-flavoured particles.
The J/ψ mesons are reconstructed at forward rapidity (2.5 < y < 4.0) via their µ + µ − decay channel. The measured J/ψ mesons originate from both prompt J/ψ (direct and from decays of higher-mass charmonium states) and non-prompt J/ψ (feed down from b-hadron decays) production.
This letter is organized as follows. A brief description of the ALICE apparatus and the data sample used is given in Sec. 2. Section 3 outlines the employed analysis technique. The evaluation of the systematic uncertainties is discussed in Sec. 4, while the results are reported in Sec. 5. Finally, conclusions are presented in Sec. 6.

Experimental setup and data sample
The ALICE detectors essential for the present analysis are briefly described below. A full overview of the ALICE apparatus and its performance can be found in Refs. [35,36]. The muon spectrometer, which covers the pseudorapidity range -4 < η < -2.5, is used to reconstruct muon tracks. The spectrometer consists of a front absorber followed by five tracking stations. The third station is placed inside a dipole magnet. The tracking stations are complemented by two trigger stations located downstream behind an iron wall. The Silicon Pixel Detector (SPD) [37] is employed to reconstruct the position of the primary vertex and to determine the flow direction. The SPD consists of two cylindrical layers covering |η| < 2.0 and |η| < 1.4, respectively. It is placed in the central barrel of ALICE. The central barrel is operated inside a solenoidal magnetic field parallel to the beam line. The SPD is also used to reconstruct the socalled tracklets, track segments formed by the clusters in the two SPD layers and the primary vertex [38]. The V0 detector [39] consists of two arrays of 32 scintillator counters each, covering 2.8 < η < 5.1 (V0A) and -3.7 < η < -1.7 (V0C), respectively. It provides the minimum-bias (MB) trigger and is used for event selection and determination of collision centrality [40]. In addition, two tungsten-quartz neutron Zero Degree Calorimeters (ZDCs), installed 112.5 meters from the interaction point along the beam line on each side, are used for event selection.
The present analysis is based on the data sample of Pb-Pb collisions collected by ALICE in 2015 at √ s NN = 5.02 TeV. The trigger required coincidence of MB and dimuon triggers. The MB trigger was provided by the V0 detector requesting signals in both V0A and V0C arrays. The dimuon unlike-sign trigger required at least a pair of opposite-sign track segments in the muon trigger stations. The transverse momentum threshold of the trigger algorithm was set such that the efficiency for muon tracks with p T = 1 GeV/c is 50%. The sample of single muons or like-sign dimuons were collected using the same trigger algorithm, but requiring at least one track segment or at least a pair of like-sign track segments, respectively. The integrated luminosity of the analyzed data sample is about 225 µb −1 .
The beam-induced background is filtered out offline by applying a selection based on the V0 and the ZDC timing information [41]. The interaction pile-up is 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. The primary vertex position is required to be within ±10 cm from the nominal interaction point along the beam direction. The data are split in intervals of collision centrality, which is obtained based on the total signal in the V0A and V0C detectors [40].
The muon selection is identical to that used in Ref. [20]. The dimuons are reconstructed in the acceptance of the muon spectrometer (2.5 < y < 4.0) and are required to have a transverse momentum between 0 and 12 GeV/c and an invariant mass M µ µ between 1.5 and 6.0 GeV/c 2 .

Analysis
The flow coefficients v n of the selected dimuons are measured using the scalar product (SP) method [2,42], in which they are calculated from the expression where u n = exp(inϕ) is the unit flow vector of the dimuon, Q SPD n , Q V0A n and Q V0C n are the event flow vectors measured in the SPD, V0A and V0C detectors, respectively, and n is the harmonic number. The brackets · · · denote an average over all events, the double brackets · · · an average over all particles in all events, and * the complex conjugate. The SPD event flow vector Q SPD n is calculated from the azimuthal distribution of the reconstructed SPD tracklets. The V0A and V0C event flow vectors Q V0A n and Q V0C n are calculated from the azimuthal distribution of the signal in the V0 detector. The components of all three event flow vectors are corrected for non-uniform detector acceptance and efficiency using a recentering procedure (i.e. by subtracting of the Q n -vector averaged over many events from the Q n -vector of each event) [43]. The denominator R n in the above equation is called resolution and is obtained as a function of collision centrality. The gap in pseudorapidity between u n and Q SPD n (|∆η| > 1.0) suppresses short-range correlations ("non-flow"), which are unrelated to the azimuthal asymmetry in the initial geometry and come from jets and resonance decays [19]. In the following, the v n {SP} coefficients are denoted as v n .
The J/ψ flow coefficients are extracted by a fit of the superposition of the J/ψ signal and the background to the dimuon flow coefficients as a function of the dimuon invariant mass where v In previous measurements [19,20], the M µ µ dependence of the background flow coefficients was parameterized by an arbitrary function. This approach leads to an increase of the statistical uncertainty of the J/ψ flow coefficients, because the parameters of the function are not fixed. Moreover, an additional systematic uncertainty arises from the fact that the functional form of the background distribution is unknown. In the present analysis, we adopt a different approach. It is known that, in collisions of heavy ions, the dimuon background in the vicinity of the J/ψ is mostly combinatorial and can be described satisfactorily with the event-mixing technique [9,17]. This technique consists in forming dimuons by combining muons from two different events having similar collision centrality. The flow coefficients of the combinatorial background are fully determined by the flow coefficients of the single muons from which the background dimuons are formed. One can show that for any given kinematical configuration of the background dimuon, its flow coefficients can be expressed as where v (1) T , η 2 ) are the flow coefficients of the two muons as a function of their transverse momenta and pseudorapidities, ϕ 1 and ϕ 2 are the azimuthal angles of the two muons and ϕ is the azimuthal angle of the dimuon. The brackets · · · M µ µ denote an average over all dimuons (p T , η 1 , η 2 , ϕ 1 , ϕ 2 ) that belong to any given M µ µ interval. The details on the derivation of Eq.(4) are given in appendix A. In case of the event mixing, the numerator in Eq. (4) is calculated as where u (1) n and u (2) n are the unit vector of the two muons, Q n and Q (2) n are the SPD flow vectors for the events containing the two muons, and R (1) n and R (2) n are their resolutions. The brackets · · · M µ µ denote an average over all mixed-event dimuons belonging to any given M µ µ interval. The denominator in Eq. (4) reflects the modification of the dimuon yield due to the flow of single muons. Since the event flow vectors of the two mixed events are not correlated, the mixed-event dimuon yield is not modified by the single muon flow. Thus, the denominator is obtained directly as the ratio N B +− /N mix +− , where N mix +− is the number of mixed-event unlike-sign dimuons as a function of M µ µ . The ratio is calculated after a proper normalization of N mix +− using the like-sign dimuons from the same and mixed events. The normalization factor is obtained as [17] where N same ++ (N same −− ) and N mix ++ (N mix −− ) are the numbers of like-sign (positive and negative charges) sameevent and mixed-event dimuons, respectively. The integral is calculated in the invariant mass interval between 2.2 and 4.5 GeV/c 2 . Assuming a purely combinatorial background, the v B n (M µ µ ) coefficient, obtained with the event-mixing procedure described above, is used directly in order to fix the background term of the fit from Eq. (3). All the analysis steps discussed in this section are performed separately in each considered dimuon transverse momentum and centrality interval. The event mixing and the normalization of N mix +− are done in 5%-wide collision centrality intervals. Examples of the M µ µ fit and the mixed-event distribution N mix +− as a function of M µ µ in several centrality and p T intervals are shown in Fig. 1. At low and intermediate p T , the mixed-event distribution describes the dimuon background on a percent level with a residual difference presumably originating from the single muon flow. However, at high p T , this difference becomes much larger (up to ≈ 35% in the vicinity of the J/ψ peak in 8 < p T < 12 GeV/c and 30-50% centrality interval) and goes beyond a possible single muon flow contribution. This points to the presence of a correlated dimuon background. Such a background is believed to originate from production of heavy-flavor quark pairs and to become significant in semi-central and peripheral collisions at high p T [47,48].
Examples of the v 2 (M µ µ ) fit based on the analysis approach described above are presented in Fig. 2. As can be seen, the fit performs quite satisfactorily, with the mixed-event v 2 coefficient being able to describe the shape and amplitude of the background v 2 in the entire considered invariant mass interval from 1.5 to 4.5 GeV/c 2 . This is not surprising at low and intermediate p T , where the mixed-event dimuon distribution describes rather precisely the background dimuon distribution (top and middle panels in Figs. 1 and 2). Remarkably, however, the mixed-event approach performs satisfactorily also at high p T in semi-central collisions, where the contribution of the correlated background is significant (bottom right panels in Figs. 1 and 2). Given that the denominator in Eq. (4) is obtained as the ratio N B +− /N mix +− , this means that the flow coefficient of the correlated background is significantly lower than that of the combinatorial one. The systematic effect arising from the presence of the correlated background and the corresponding uncertainties are discussed in Section 4. The approach described above performs equally well also in case of the v 3 coefficient. This is illustrated in Fig.3, where the fits of the centrality and p T -integrated v 2 (M µ µ ) and v 3 (M µ µ ) distributions are compared.
The Event Shape Engineering (ESE) technique is performed following the procedure described in Ref. [33]. It is based on the magnitude of the second-order reduced V0A event flow vector defined as in Ref. [42] where |Q V0A 2 | is the magnitude of the second-order V0A event flow vector and S V0A is the total signal in the V0A detector. The large pseudorapidity gap between the V0A and the muon spectrometer (|∆η| > 5.3) greatly suppresses the non-flow contribution and guarantees a proper event-shape selection. Two event-shape classes with the lowest and highest q V0A 2 values corresponding to the 0-20% and 80-100% intervals, respectively, are investigated for the 5-40% centrality interval.

Systematic uncertainties
The systematic effect related to the presence of correlated background is checked by modifying the definition of the background coefficient , where the parameter α represents the strength of the flow of the correlated background. The value of 0 corresponds to the default approach (e.g. assuming negligible flow of the correlated background), while the value of 1 corresponds to the assumption that the correlated background has the same flow coefficient as compared to the combinatorial background. The parameter α is left free in the fit of Eq. (3) and the differences in the resulting J/ψ v 2 with respect to the default approach are taken as systematic uncertainties. As expected, in central (0-10%) collisions and at low transverse momentum, the uncertainties are practically negligible. In semi-central (30-50% centrality interval) collisions and in the highest considered transverse momentum interval (8 < p T < 12 GeV/c), the uncertainty of the J/ψ v 2 reaches 0.013. The parameter α is found to be well below 1 in all centrality and p T intervals. The corresponding systematic uncertainty of the J/ψ v 3 coefficient is in general significantly smaller. No clear pattern is found as a function of collision centrality and p T . Conservatively, the parameter α is fixed to 1 and the difference in the results with respect to the ones obtained with default value of 0 is taken as systematic uncertainty. It is worth noting that even though the fraction of correlated background at high p T in semi-central collisions is significant, its effect on the J/ψ flow coefficients is suppressed by the high signal-to-background ratio N J/ψ /N B +− . As described in appendix A, a small additional v  term is present in v B 2 . Its estimated contribution is added to the fit to the v 2 (M µ µ ) distribution and the change in the J/ψ v 2 results with respect to the default approach is taken as systematic uncertainty. These uncertainties are found to be sizable only in 0 < p T < 2 GeV/c and 10-50% centrality interval, where they reach 0.002.   The effect of any residual non-uniform detector acceptance and efficiency in the calculation of the SPD event flow vector is checked via the imaginary part of the scalar product defined in Eq. (2) [49]. No systematic uncertainty is assigned as the terms are consistent with zero within statistical uncertainties. The resolution of the SPD event flow vector is calculated from the events containing at least one selected dimuon by default. Alternatively, it is calculated from all events recorded with the MB trigger and passing the offline event selection, as well as from the events containing at least one selected single muon. Differences up to 1% and 2% with respect to the default approach are observed for R 2 and R 3 , respectively, and are taken as systematic uncertainties. For the event-shape classes, a bias can arise from auto-correlations due to the usage of the V0A event flow vector for both q 2 and R 2 . This potential bias is assessed by replacing the ratio Q SPD n Q V0A * n / Q V0A n Q V0C * n in Eq. 2 with the one from the unbiased data sample. The resulting effect is smaller than 1% and is neglected.
The muon spectrometer occupancy affects the reconstruction efficiency and thus can bias (lower) the measured v n coefficients. The reconstruction efficiency as a function of centrality is taken from Ref. [11], where it is obtained by embedding simulated J/ψ → µ + µ − decays into real Pb-Pb events. It is found to decrease linearly with the signal in the V0C detector S V0C , which largely covers the geometrical acceptance of the muon spectrometer. Thus, the systematic deviations of the J/ψ v n are calculated as the product of the single muon v n , the first derivative of the reconstruction efficiency with respect to S V0C and the mean S V0C in the considered centrality interval. The single muon v n coefficients are obtained with the same SP approach as the one employed for J/ψ. Conservatively, the maximum of the single muon v n as a function of p T is used. The typical values of these systematic deviations are found to be up to 0.0025 and 0.0015 for the J/ψ v 2 and v 3 , respectively. Given the small magnitude of the effect, we do not correct the measured coefficients, but take the above deviations as systematic uncertainties. v n (h ± ) (n = 2, 3). At high p T , above 6-8 GeV/c, the v 2 results indicate a convergence between charged particles, prompt D 0 mesons and J/ψ. Such an observation suggests that, at high p T , the azimuthal asymmetry of the J/ψ mesons as well as that of charged particles and prompt D 0 mesons is possibly governed by in-medium path-length dependent energy-loss effects.

Results
Discussing the above observations, should be noted the different rapidity interval of the J/ψ measurement. The effect of the decorrelation of the symmetry plane angles Ψ n (n = 2, 3) between mid and forward pseudorapidity has been estimated to be less than 1% and 3% for v 2 and v 3 , respectively [51, 52]. An η dependence of the p T -integrated v n coefficients for charged particles has been observed in Pb-Pb collisions at √ s NN = 2.76 TeV [53]. However, the ratio v 3 /v 2 has shown no significant dependence on η. Furthermore, the p T -differential v 2 was found to be independent of η (up to |η| < 2.4) [54], thus indicating that the η dependence of the p T -integrated v 2 arises mainly from changes in the transverse momentum spectra. The present analysis of the J/ψ v 2 coefficient, performed in the centrality intervals used in Ref.
[19], yields consistent results. It is worth noting that the statistical uncertainties are reduced by up to 15% due to the event-mixing approach described in Section 3.
In Fig. 4, the J/ψ v 3 is positive in most of the intervals, although it is also compatible with zero given the large uncertainties. A positive value of v 3 is found integrating the data over the centrality intervals, as seen in Fig. 5. The Fisher's combined probability test [58] is used to quantify the probability that J/ψ v 3 is zero. The data in all p T intervals are treated as independent measurements. The statistical and  systematic uncertainties are added in quadrature. The total combined probability of the zero hypothesis is found to be 1.23×10 −4 , which corresponds to about 3.7σ significance of the measured positive J/ψ v 3 coefficient.
The flow coefficients of the J/ψ, prompt D 0 mesons and charged particles are further compared in Fig. 6, where the ratio v 3 /v 2 is shown as a function of p T . In order to increase the significance of the ratio, the central collisions (0-5% and 0-10% centrality intervals), where v 2 has small magnitude, are excluded. The uncertainties of v 2 and v 3 coefficients are considered uncorrelated due to the weak correlation between the Ψ 2 and Ψ 3 angles [59]. Taking into account all p T intervals, the obtained J/ψ v 3 /v 2 ratio is found to be significantly lower (4.6σ ) with respect to that of charged particles. Moreover, at intermediate p T between 2 and 6 GeV/c, the prompt D 0 -mesons v 3 /v 2 ratio is 2.3σ below that of charged particles and 3.4σ above that of the J/ψ mesons. Thus, the data seem to suggest an ordering similar to the one observed for the v 2 and v 3 coefficients in semi-central collisions. It is interesting to note that, in con- The left panel of Fig. 7 presents the J/ψ v 2 as a function of p T for event-shape selected and unbiased events in the 5-40% centrality interval. The systematic uncertainties of the results from the event-shape selected and unbiased events are considered fully correlated and therefore cancel out in the ratios shown in the right panel of Fig. 7. The values of the J/ψ v 2 coefficient in low (high) q V0A 2 event classes are found to be lower (higher) with respect to those in the unbiased events. The v 2 coefficient of single muons is also measured in the same event-shape selected and unbiased samples. The corresponding ratios between the results in the event-shape selected and unbiased events show no p T dependence up to 10 GeV/c (Fig. 7, right panel). This behavior demonstrates that the applied ESE technique based on q V0A 2 allows the selection of a global property of the collisions, most likely linked to the eccentricity ε 2 of the initial-state geometry [33]. The mean values of the ratios for single muons v 2 {low-q V0A 2 }/v 2 {unbiased} and v 2 {high-q V0A 2 }/v 2 {unbiased} are estimated from a fit with constant and are found to be 0.87 and 1.15, respectively. These values reflect the sensitivity of the V0A-based event-shape selection. The corresponding mean values of the J/ψ ratios, 0.79±0.14 and 1.35±0.14, are consistent with the muon ratios. This implies that the J/ψ v 2 results are compatible with the expected variations of the eccentricity of the initial-state geometry within the uncertainties.

Conclusions
In summary, the elliptic and triangular flow coefficients of inclusive J/ψ mesons at forward rapidity have been measured in Pb-Pb collisions at √ s NN = 5.02 TeV over a broad range of transverse momentum and in various centrality intervals. This is the first measurement of the v 3 coefficient for inclusive J/ψ production, indicating a positive value with 3.7σ significance for 0 < p T < 12 GeV/c. The obtained inclusive J/ψ v 2 and v 3 coefficients as well as the ratio v 3 /v 2 are compared to the results for charged particles and prompt D 0 mesons at mid-rapidity. At low and intermediate p T , the v 2 and v 3 results exhibit an ordering with the charged particles having largest values, followed by the prompt D 0 mesons and finally the J/ψ having the smallest values. In semi-central collisions at intermediate p T , the J/ψ v 3 /v 2 ratio is found to be significantly lower compared to that of charged particles. Despite the large uncertainties, the values of the prompt D 0 ratio are somewhat lower than the charged particles and higher than the J/ψ mesons, hinting at a possible ordering similar to that observed for the v 2 and v 3 coefficients.
At high p T , the v 2 of the charged particles, the prompt D 0 mesons and the J/ψ seem to converge to similar values. The uncertainties of the v 3 coefficients do not allow one to draw firm conclusions about their convergence, although the centrality-and p T -integrated J/ψ v 3 is compatible with that of high-p T charged particles.
The analysis using Event Shape Engineering technique shows that the J/ψ v 2 coefficients increase (decrease) for classes of events with high (low) reduced event flow vector. Compared to single muons reconstructed in the same rapidity interval, the J/ψ results are found compatible with the expected variations of the eccentricity of the initial-state geometry.

Acknowledgements
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centres and the Worldwide LHC Computing Grid (WLCG) collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the ALICE detector:       [38] ALICE Collaboration, J. Adam et al., "Charged-particle multiplicities in proton-proton collisions at √ s = 0.9 to 8 TeV," Eur. Phys. J. C77 no.

A Appendix
The azimuthal distribution of the combinatorial background dN B /dϕ is a product of the azimuthal distributions of the single muons from which the background dimuons are formed. Thus, using Eq. (1)  T , η 2 ) are the flow coefficients of the two muons as a function of their transverse momenta and pseudorapidities, ϕ 1 and ϕ 2 are the azimuthal angles of the two muons, ϕ is the azimuthal angle of the dimuon and ∆ϕ 1,2 = ϕ 1,2 − ϕ. . (A.9) Finally, the v B n as a function of M µ µ is obtained by averaging the numerator and denominator in Eq.(A.9) over all dimuons, which belong to a given M µ µ interval: where the brackets · · · denote an average over all events. The contribution is estimated as described in the following. First, the v 2 and v 4 coefficients of single muons are measured with the SP method, averaged over pseudorapidity and parameterized as a function of p T . The obtained parameterizations v 2,4 (p T ) are then combined with of opposite-sign dimuons (p (1) T , p T , η 1 , η 2 , ϕ 1 , ϕ 2 ) in the data outside the J/ψ mass peak. The values of cos[4(Ψ 4 − Ψ 2 )] , which ranges from 0 in central collisions to about 0.8 in peripheral collisions, are taken from Ref. [59]. Finally, the magnitude of the effect is calculated via interpolation of the results at the J/ψ mass peak. In general, the magnitude is found to be at the order of 10 −4 , reaching at most 7 × 10 −4 for 0 < p T < 2 GeV/c and the 30-50% centrality interval.
A similar effect is present in the numerator of Eq.(A.10) for v B 3 , due to the correlation of the Ψ 3 and Ψ 6 angles. In practice, however, this contribution can be certainly neglected, because of the small magnitude of the v 6 coefficient.