Search for collectivity with azimuthal J/$\psi$-hadron correlations in high multiplicity p-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 and 8.16 TeV

We present a measurement of azimuthal correlations between inclusive J/$\psi$ and charged hadrons in p-Pb collisions recorded with the ALICE detector at the CERN LHC. The J/$\psi$ are reconstructed at forward (p-going, 2.03 $<$ y $<$ 3.53) and backward (Pb-going, $-$4.46 $<$ y $<$ $-$2.96) rapidity via their $\mu^+\mu^-$ decay channel, while the charged hadrons are reconstructed at mid-rapidity ($|\eta|$ $<$ 1.8). The correlations are expressed in terms of associated charged-hadron yields per J/$\psi$ trigger. A rapidity gap of at least 1.5 units is required between the trigger J/$\psi$ and the associated charged hadrons. Possible correlations due to collective effects are assessed by subtracting the associated per-trigger yields in the low-multiplicity collisions from those in the high-multiplicity collisions. After the subtraction, we observe a strong indication of remaining symmetric structures at $\Delta\varphi$ $\approx$ 0 and $\Delta\varphi$ $\approx$ $\pi$, similar to those previously found in two-particle correlations at middle and forward rapidity. The corresponding second-order Fourier coefficient ($v_2$) in the transverse momentum interval between 3 and 6 GeV/$c$ is found to be positive with a significance of about 5$\sigma$. The obtained results are similar to the J/$\psi$ $v_2$ coefficients measured in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV, suggesting a common mechanism at the origin of the J/$\psi$ $v_2$.


Introduction
The measurement of angular correlations between particles produced in hadron and nucleus collisions is a powerful tool to study the particle production mechanisms. Usually the two-particle correlation function is expressed in terms of differences in the azimuthal angle (∆ϕ) and pseudorapidity (∆η) of the emitted particles. In minimum-bias proton-proton (pp) collisions, the dominant structures in the correlation function are a near-side peak at (∆ϕ, ∆η) ≈ (0, 0) and an away-side ridge located at ∆ϕ ≈ π and elongated in ∆η [1]. The near-side peak originates from jet fragmentation, resonance decays and femtoscopic correlations. The away-side ridge results from fragmentation of recoil jets. In collisions of heavy ions, the two-particle correlation function exhibits additional long-range structures elongated in ∆η [2]. These structures are usually interpreted as signatures of collective particle flow produced during the hydrodynamic evolution of the fireball. They are analyzed in terms of the Fourier coefficients of the relative angle distributions. Assuming factorization, these coefficients are then related to the Fourier coefficients (v n ) of the particle azimuthal distribution relative to the common symmetry plane of the colliding nuclei's overlap area.
The discovery of a near-side ridge in high-multiplicity pp [3] and p-Pb [4] collisions has increased the interest in two-particle angular correlations in small collision systems. These discoveries were followed by the observation that the near-side ridge in p-Pb collisions is accompanied by an away-side one [5,6]. Long-range structures have also been reported in two-particle correlations in d-Au collisions at RHIC [7,8]. Further studies using multi-particle correlations have proven that the observed long-range correlations are of a collective origin [9][10][11]. Moreover, the transverse-momentum and particle-mass dependencies of the v n coefficients in p-Pb collisions have been found to be similar to those measured in A-A collisions, suggesting a common hydrodynamic origin of the observed correlations [12,13]. Alternative interpretations, including Color-Glass Condensate based models [14] and final-state partonparton scattering [15], have also been proposed. Long-range correlations of forward and backward muons with mid-rapidity hadrons have also been found in p-Pb collisions at a centre-of-mass energy per nucleon pair √ s NN = 5.02 TeV [16]. The results show that these correlations persist across wide rapidity ranges and extend into the high muon transverse-momentum interval, which is dominated by decays of heavy flavours.
In pp collisions, the J/ψ resonance is formed mainly from pairs of c andc quarks produced in hard scattering reactions during the initial stage of the collision. The theoretical models describing the J/ψ production combine calculations of the production of cc pairs within a perturbative Quantum Chromodynamics approach with the subsequent non-perturbative formation of the cc bound state [17]. In p-Pb collisions, the production is affected by the modification of parton distribution functions inside the nucleus [18] as well as possible energy loss and inelastic scattering inside nuclear matter [19,20]. In A-A collisions, there are two additional competing phenomena that influence the J/ψ production. First is the suppressed production due to the dissociation of the cc pairs in the quark-gluon plasma [21]. Second is the J/ψ enhancement via recombination of charm quarks thermalized in the medium [22,23]. The recombination is expected to become prevalent in central collisions at the LHC energies.
Recently, the ALICE Collaboration has published a precise measurement of the second-order Fourier coefficient, v 2 , of the azimuthal distribution of the J/ψ production in Pb-Pb collisions at √ s NN = 5.02 TeV [24]. The results show significant v 2 in central and semi-central collisions. The measured J/ψ v 2 at low and intermediate transverse momentum can be qualitatively described by a transport model in which the J/ψ azimuthal anisotropy is inherited from that of recombined charm quarks [25,26]. However, at higher transverse momentum the data still indicates significant v 2 while the transport model predicts significantly smaller values coming mostly from path-length dependent suppression in the almond-shaped interaction region of the colliding nuclei and from non-prompt J/ψ produced from b-hadron decays assuming thermalized b quarks. Given these results in Pb-Pb collisions, it is of interest to study the J/ψhadron azimuthal correlations also in the smaller p-Pb system. The recombination of charm quarks, if any, should have much smaller impact, due to the smaller number of initially produced charm quarks with respect to Pb-Pb collisions. The small system size should not lead to a sizeable path-length dependent suppression. Nevertheless, the study of the J/ψ-hadron azimuthal correlations could allow to determine whenever J/ψ production is affected by the medium possibly created in these collisions [27][28][29].
In this Letter, we present results for long-range correlations between forward (p-going, 2.03 < y < 3.53) and backward (Pb-going, −4.46 < y < −2.96) inclusive J/ψ and mid-rapidity charged hadrons in p-Pb collisions at √ s NN = 5.02 and 8.16 TeV. Inclusive J/ψ refers to both prompt J/ψ (direct and decays from higher mass charmonium states) and non-prompt J/ψ (feed down from b-hadron decays).

Experimental setup and data samples
A detailed description of the ALICE apparatus can be found in Ref. [30]. Below, we briefly describe the detector systems essential for the present analysis.
In the following, η and y lab will denote the pseudorapidity and rapidity in the ALICE laboratory system. The muons are reconstructed in the muon spectrometer covering the range of −4 < η < −2.5. The spectrometer contains a front absorber located between 0.9 and 5 m from the nominal interaction point. The absorber is followed by five tracking stations, each made of two planes of Cathode Pad Chambers. The third station is placed inside a dipole magnet with 3 Tm field integral. The tracking stations are followed by an iron wall with a thickness of 7.2 interaction lengths and two trigger stations, each one consisting of two planes of Resistive Plate Chambers.
The position of the interaction point is obtained using the clusters reconstructed in the Silicon Pixel Detector (SPD) [31,32]. The SPD is located in the central barrel of the ALICE apparatus and operated inside a large solenoidal magnet providing a uniform 0.5 T magnetic field parallel to the beam line. The SPD consists of two cylindrical layers which cover |η| < 2.0 and |η| < 1.4 with respect to the nominal interaction-point, for the inner and outer layer, respectively. The associated charged hadrons at midrapidity are reconstructed via the so called SPD tracklets, short track segments formed from the clusters in the two layers of the SPD and the primary vertex [32].
The data samples presented here were collected during the 2013 and 2016 p-Pb LHC runs. The collision energy was √ s NN = 5.02 and 8.16 TeV for the 2013 and 2016 data samples, respectively. Part of the 5.02 TeV data were collected during the 2016 p-Pb run. Data with both beam configurations, namely Pb-nucleus momentum (denoted as Pb-p collisions) or proton momentum (denoted as p-Pb collisions) oriented towards the muon spectrometer, have been analyzed. The asymmetric beam energies, imposed by the two-in-one LHC magnet design, resulted in collisions whose nucleon-nucleon centre-of-mass reference system is shifted in rapidity by 0.465 in the direction of the proton beam with respect to the ALICE laboratory system. The data were taken with a trigger that required coincidence of minimum-bias (MB) and dimuon triggers. The MB trigger was provided by the V0 detector requesting a signal in both V0-A and V0-C rings. Its efficiency is found to be about 98% [34]. The dimuon trigger required at least a pair of opposite-sign track segments in the muon trigger system, each with a transverse momentum (p T ) above the threshold of the online trigger algorithm. This threshold was set to provide 50% efficiency for muon tracks with p T = 0.5 GeV/c.
The collected data samples of p-Pb and Pb-p collisions at 5.02 TeV (8.16 TeV) correspond to integrated luminosities of 8.1 and 5.8 (8.7 and 12.9) nb −1 , respectively. The maximum interaction pile-up probability ranged up to 3% and 8% during 2013 and 2016 data taking, respectively.

Event, track and dimuon selection
The beam-induced background is rejected by requiring that the timing signals from both rings of the V0 detector are compatible with particles coming from collision events. Events containing multiple collisions (pile-up) are rejected by requiring one single interaction vertex reconstructed in the SPD and by exploiting the correlation between the number of clusters in the two layers of the SPD and the number of the reconstructed SPD tracklets.
The longitudinal position of the reconstructed primary vertex (z vtx ) is required to be within ±10 cm from the nominal interaction point. The reconstructed SPD tracklets are selected by applying a z vtx -dependent pseudorapidity cut. The cut is adjusted to exclude the contribution from the edges of the SPD where the detector acceptance is low. For example, we select tracklets within −1.8 < η < 0.5, −1.3 < η < 1.3 and −0.5 < η < 1.8 for events with z vtx = 10, 0 and −10 cm, respectively. The contribution from fake and secondary tracklets is reduced by applying a |∆Φ| < 5 mrad cut on the difference between the azimuthal angles of the clusters in the two layers of the SPD with respect to the primary vertex. With this cut, the mean p T of the selected charged hadrons is found to be approximately 0.75 GeV/c [16].
The tracks reconstructed in the muon spectrometer are required to emerge at a radial transverse position between 17.6 and 89.5 cm from the end of the front absorber in order to avoid regions with higher material budget. The tracks reconstructed in the tracking chambers are identified as muons by requiring their matching with corresponding track segments in the trigger chambers. Background tracks are removed with a selection on the product of the total track momentum and the distance of closest approach to the primary vertex in the transverse plane [35]. The selected dimuons are defined as pairs of oppositesign muon tracks having −4 < y µ µ lab < −2.5, transverse momentum p µ µ T between 0 and 12 GeV/c and invariant mass M µ µ between 1 and 5 GeV/c 2 . Only events with at least one dimuon satisfying these selection criteria are considered.
The data samples are split into multiplicity classes based on the total charge deposited in the two rings (V0-A and V0-C) of the V0 detector (V0M) [34]. The high-multiplicity (low-multiplicity) event class is defined as 0-20% (40-100%) of the MB trigger event sample.

Analysis
The M µ µ distribution in each event-multiplicity class and p µ µ T bin is fit with the combination of an extended Crystal Ball (CB2) function for the J/ψ signal and a Variable-Width Gaussian (VWG) function for the background [36]. The tail parameters of the CB2 function were fixed to the values used in [37, 38]. The J/ψ peak position and width were obtained from the fit in the 0-100% event class and fixed to these values in the other event-multiplicity classes. Examples of the M µ µ fit in the 0-20% and the 40-100% event classes in the 3 < p µ µ T < 6 GeV/c interval are shown in Fig.1. The angular correlations between J/ψ and charged hadrons are obtained from the associated-particle (SPD tracklets) yields per dimuon trigger. The yields are defined as where is the number of associated SPD tracklets corrected for acceptance and combinatorial effects (as shown in the second line of the is the yield of associated SPD tracklets from the same event. The distribution is constructed using the event-mixing technique, i.e. combining dimuons from one event with SPD tracklets from other events selected in the same event-multiplicity class and z vtx interval. It serves both to correct for detector acceptance and efficiency and to take into account the combinatorial background. Within each event-multiplicity class and bin of M µ µ , p µ µ T , ∆ϕ and ∆η, the yields Y i averaged over z vtx are obtained by fitting the distribution Y i N trig (z vtx ) i ME i (z vtx ) to the distribution SE i (z vtx ). A Poisson likelihood fit is used in order to properly deal with the cases of low number of tracklets. Then, the average yields are projected on the ∆ϕ axis in the range of 1.5 < |∆η| < 5 using the method described in [16].
In order to extract the yields per J/ψ trigger, the yields per dimuon trigger in each event-multiplicity class, p µ µ T and ∆ϕ bins are fit as a function of M µ µ using the following superposition where S and B are the number of J/ψ and the background dimuons in each bin of M µ µ obtained from the invariant mass fit (using a CB2 function for the J/ψ signal and a VWG function for the background) described above, Y J/ψ is the associated yield corresponding to the J/ψ trigger and Y B (M µ µ ) is a secondorder polynomial function aimed to describe the associated yields corresponding to the background. The fit range is chosen between 1.5 and 4.5 GeV/c 2 . Examples of fits in high-multiplicity and low-multiplicity event classes are shown in Fig. 2.   Figure 3 shows the obtained associated tracklet yields per J/ψ trigger for p-Pb and Pb-p collisions at √ s NN = 5.02 and 8.16 TeV. As expected, in low-multiplicity collisions we observe a significant correlation structure on the away side (Fig. 3, top panels), presumably originating from the fragmentation of recoil jets. In high-multiplicity collisions (Fig. 3, middle panels), a possible enhancement on both near (∆ϕ ≈ 0) and away (∆ϕ ≈ π) side can be spotted on top of the away-side structure. In order to isolate possible correlations due to collective effects between the J/ψ and the associated tracklets, we apply the same subtraction method as in previous measurements [5, 6, 12, 16], namely subtracting the Y J/ψ yields in low-multiplicity collisions from those in high-multiplicity collisions (Fig. 3, bottom panels). The subtraction method relies on the assumptions that the jet correlations on the away side remain unmodified as a function of the event multiplicity and that there are no significant correlations due to collective effects in low-multiplicity collisions (see discussion in Section 6).
In order to quantify the remaining correlation structures, the subtracted yields Y sub J/ψ (∆ϕ) are fit with The second-order Fourier coefficient V 2 {J/ψ − tracklet, sub} of the azimuthal correlation between the J/ψ and the associated charged hadrons is finally calculated as a 2 /b high 0 . The denominator b high 0 = a 0 + b low 0 corresponds to the combinatorial baseline of the high-multiplicity collisions, where the parameter b low 0 is the combinatorial baseline of the low-multiplicity collisions obtained at the minimum of the pertrigger yields, namely in ∆ϕ < π/6. The parameter b low 0 is the normalization factor used in Fig. 3. The parameter a 1 , which describes the strength of the remaining away-side correlation structure, is found to be compatible with zero in practically all p J/ψ T intervals, in both p-Pb and Pb-p collisions at both 5.02 and 8.16 TeV.
As an alternative extraction method, the calculation of b low 0 , the subtraction of low-multiplicity from high-multiplicity collision yields and the fit to Eq. (3) is done in each bin of M µ µ separately. Then the V 2 {J/ψ − tracklet, sub} coefficient is extracted by fitting V 2 {µ µ − tracklet, sub}(M µ µ ) with a superposition similar to the one defined in Eq. (2)    with 1-2% relative statistical uncertainty and 5-6.5% relative systematic uncertainty.

Systematic uncertainties
The combined statistical and systematic uncertainties of the measured v tracklet A possible inaccurate correction for the SPD acceptance is assessed by varying the z vtx range between ±8 and ±12 cm. Systematic uncertainties are assigned only in the cases of a significant change of the results. The significance is defined according to the procedure described in Ref. [39].
The systematic effect related to the uncertainty of the shape of the dimuon background yields Y B (M µ µ ) is estimated by performing the fit with Eq. (2) using a linear function for the background term and varying the fit range. The systematic effect coming from the uncertainty of the signal-to-background ratio S/B is checked by employing various invariant mass fit functions, both for the background and for the J/ψ signal. The maximal difference of the results obtained with the above checks with respect to the default approach is taken as the corresponding systematic uncertainty.
The uncertainty arising from the employed analysis approach is obtained as the difference between the two extraction methods described in Section 4.
As described in Section 4, by default the mixed-event distribution ME(∆ϕ, ∆η) is normalized to unity in the ∆η region corresponding to the maximal acceptance. As an alternative approach, normalizing the integral of ME(∆ϕ, ∆η) to unity is used. No significant effect on the obtained results is observed and thus no systematic uncertainty is assigned.
The used event-mixing technique can introduce systematic biases. The event multiplicity distribution of the selected dimuons (1 < M µ µ < 5 GeV/c 2 ) differs from that of the J/ψ signal. Since the charged-hadron spectra and the charged-hadron density as a function of η change with event multiplicity [34], the non- uniform (both in the azimuthal and longitudinal directions) SPD acceptance can introduce a bias. The corresponding systematic uncertainty is evaluated by doing the event mixing in finer event-multiplicity bins.
The non-uniform acceptance of the muon spectrometer coupled to sizeable correlations between the dimuons and SPD tracklets can bias azimuthally the sample of SPD tracklets used for event mixing. In order to check for possible effects on our measurement, the event mixing is performed in intervals of azimuthal angle of the selected dimuons. We observe no significant systematic effect as the obtained results show negligible deviations with respect to the results using the default event-mixing technique.
The effect of a possible residual near-side peak is checked by varying the rapidity gap between the trigger dimuons and associated charged-hadrons from 1.0 to 2.0 units. We observe no indication of increasing v 2 with reduced gap and thus consider the default gap of 1.5 units sufficient to eliminate any significant residual near-side peak contribution.
As shown in Section 4, the recoil-jet away-side correlation structure in the high-multiplicity event class is greatly diminished after the subtraction of the low-multiplicity event class. By default, any remaining away-side structure is supposed to be taken into account by the cos ∆ϕ term in Eq. (3). In order to check for residual effects we proceed in the following way. First, the correlation function in the low-multiplicity event class is fit with a Gaussian function centered at ∆ϕ = π. Then, the correlation function in the highmultiplicity event class is fit with the function from Eq. (3), where the cos ∆ϕ term is replaced by a Gaussian function with a width fixed to the value obtained from the fit in the low-multiplicity collisions. No clear signature of systematic change of the results is seen, except some hints of a possible effect in the highest p J/ψ T interval. Conservatively, we assign systematic uncertainty as the difference with respect to the default analysis approach. Since the typical values of the Gaussian width are around 1 rad, one-sided (negative) systematic uncertainty is assigned.
In Table 1 we present a summary of the assigned systematic uncertainties of the v As additional cross-checks the analysis is done using alternative event-multiplicity estimators, varying the tracklet |∆Φ| cut, applying a cut on the asymmetry of transverse momentum of the two muon tracks, removing the pile-up cuts and excluding the SPD regions with non-uniform acceptance in pseudorapidity. The corresponding results are found to be compatible with those obtained with the default analysis approach and therefore no further systematic uncertainties are assigned.

Results
In Fig. 5 we report the measured v significance. In the calculation of the above probabilities, both statistical and systematic uncertainties of the measured values are taken into account. The global systematic uncertainty is not taken into account as it is irrelevant in the case of the zero hypothesis.
The analysis method presented in this Letter relies on the assumption that there are no significant correlations due to collective effects in the low-multiplicity event class. In case of a presence of such correlations, the measured V 2 {J/ψ − tracklet, sub} is equal to where V 2 {J/ψ − tracklet, high} and V 2 {J/ψ − tracklet, low} are the second-order Fourier coefficients of the azimuthal correlation between the J/ψ and the associated charged hadrons in the high-multiplicity and the low-multiplicity collisions, respectively, and b low 0 /b high 0 ≈1/3 is the ratio of the combinatorial baseline in the low-multiplicity and high-multiplicity collisions (see Fig. 3). As is demonstrated in Ref. [44], the assumption of no significant collective correlations in the low-multiplicity collisions is certainly questionable for light-flavour hadrons. Our data indicates the same, as we observe a statistically significant increase of the measured values of v tracklet 2 {2, sub} when subtracting a lower event-multiplicity, e.g. 60-100%, class. Ultimately, the value of the v tracklet 2 coefficient is found to be about 17% higher in case no subtraction is applied. Therefore, replacing the subtracted v tracklet TeV also appear to be consistent with each other. The largest absolute difference between the results at the two collision energies is observed in Pb-p collisions in the 3 < p J/ψ T < 6 GeV/c interval. The significance of this difference is rather low (below 1.5σ ), because of the large uncertainties of the measurement at √ s NN = 5.02 TeV. Hence, the data for the two collision energies are combined as a weighted average taking into account both statistical and systematic uncertainties. In Fig. 6, we present these combined results for p-Pb and Pb-p collisions together with measurements and model calculations for Pb-Pb collisions at √ s NN = 5.02 TeV [25].
In Pb-Pb collisions, the positive v J/ψ 2 coefficients at p J/ψ T below 3-4 GeV/c are believed to originate from the recombination of charm quarks thermalized in the medium and are described fairly well by the transport model [25] (see Fig. 6). In p-Pb collisions, the amount of produced charm quarks is small and therefore the contribution from recombination should be negligible. Our measured values at p J/ψ T < 3 GeV/c are compatible with zero, in line with this expectation. There is one publication [28] which suggests that even in p-Pb collisions a sizeable contribution from recombination could occur due to canonical enhancement effects. The uncertainties of our results do not allow to confirm or to rule out this scenario.