Measurements of azimuthal anisotropies at forward and backward rapidity with muons in high-multiplicity p-Pb collisions at $\sqrt{s_{\rm NN}} = 8.16$ TeV

The study of the azimuthal anisotropy of inclusive muons produced in p-Pb collisions at $\sqrt{s_{\rm NN}} = 8.16$ TeV, using the ALICE detector at the LHC is reported. The measurement of the second-order Fourier coefficient of the particle azimuthal distribution, $v_2$, is performed as a function of transverse momentum $p_{\rm T}$ in the 0-20% high-multiplicity interval at both forward ($2.03<y_{\rm CMS}<3.53$) and backward ($-4.46<y_{\rm CMS}<-2.96$) rapidities over a wide $p_{\rm T}$ range, $0.5<p_{\rm T}<10$ GeV/$c$, in which a dominant contribution of muons from heavy-flavour hadron decays is expected at $p_{\rm T}>2$ GeV/$c$. The $v_2$ coefficient of inclusive muons is extracted using two different techniques, namely two-particle cumulants, used for the first time for heavy-flavour measurements, and forward-central two-particle correlations. Both techniques give compatible results. A positive $v_2$ is measured at both forward and backward rapidities with a significance larger than $4.7\sigma$ and $7.6\sigma$, respectively, in the interval $2<p_{\rm T}<6$ GeV/$c$. Comparisons with previous measurements in p-Pb collisions at $\sqrt{s_{\rm NN}} = 5.02$ TeV, and with AMPT and CGC-based theoretical calculations are discussed. The findings impose new constraints on the theoretical interpretations of the origin of the collective behaviour in small collision systems.


Introduction
The study of high-energy heavy-ion collisions aims at investigating the properties of a state of stronglyinteracting matter characterised by high energy density and temperature, called the quark-gluon plasma (QGP) [1,2].One of the key observables to address the transport properties of the QGP is the azimuthal anisotropy of produced particles [3].In non-central collisions, the initial spatial anisotropy of the overlap region is converted into an anisotropy in momentum space via multiple interactions.The magnitude of the azimuthal anisotropies is usually quantified via a Fourier decomposition of the particle azimuthal distribution given by where ϕ and p T are the particle azimuthal angle and transverse momentum, respectively.The Fourier coefficients v n characterise the anisotropy of produced particles [4] and Ψ n is the azimuthal angle of the symmetry plane for the n th harmonic.The largest contribution to the asymmetry of non-central collisions is provided by the second Fourier coefficient v 2 referred to as elliptic flow, and is expressed as v 2 = ⟨cos[2(ϕ − Ψ 2 )]⟩ [3,4], where the brackets denote the average over all particles and all selected events.
Due to their large masses, heavy quarks (charm and beauty) are mostly produced via hard partonic scattering processes in the very early stage of the collisions before the formation of the QGP, and therefore they probe the properties and dynamics of the QGP through its full evolution [5].During their propagation through the medium, they lose energy via elastic and inelastic processes.Heavy-flavour hadrons and their decay products are thus sensitive probes of the QGP medium.The open heavy-flavour v 2 coefficient is expected to provide information on the collective expansion of the medium and thermalisation of heavy quarks at low p T [3,6], and is sensitive to the path-length dependence of the in-medium energy loss at high p T [7,8].These contributions are predicted to give positive v 2 values.Moreover, the elliptic flow is also expected to be sensitive to the hadronisation processes at low and intermediate p T [9,10].Extensive measurements of the elliptic flow of heavy flavours have been carried out in Pb-Pb collisions at a centreof-mass energy √ s NN = 2.76 and 5.02 TeV at the LHC, and in Au-Au collisions at √ s NN = 200 GeV at the RHIC [11] where a positive elliptic flow coefficient was observed, see Refs.[12][13][14][15][16][17][18][19][20][21][22][23][24] for results on open heavy flavours.The nuclear modification factor R AA , defined as the ratio of the particle yields in nucleus-nucleus (AA) collisions at the same centre-of-mass energy to that in binary-scaled pp collisions, is also an important observable to quantify in-medium effects.In addition to the non-zero v 2 coefficient, the R AA measurements in heavy-ion collisions at LHC energies show a strong suppression of the open heavy-flavour yields at intermediate and high p T [12,17,22,[25][26][27][28][29][30][31][32][33][34][35][36][37].These results reflect significant energy losses of heavy quarks due to strong interactions with the medium constituents.
The study of small collision systems such as p-Pb collisions at the LHC was initially proposed to address cold nuclear matter (CNM) effects relevant for the interpretation of the measurements in heavy-ion collisions, such as the nuclear modification of the parton distribution functions [38], k T broadening [39], and energy loss in cold nuclear matter [40].Surprisingly, long-range structures in two-particle correlations ("ridges") associated with a positive v 2 were first observed for light-flavour hadrons in high-multiplicity p-Pb collisions at √ s NN = 5.02 TeV in the midrapidity region by the ALICE, ATLAS, and CMS collaborations [41][42][43][44][45][46] and at forward rapidity by the LHCb collaboration [47], using the two-particle azimuthal correlation method.In the heavy-flavour sector, the ALICE and CMS collaborations revealed also hints of collectivity with the measurement of a positive v 2 of inclusive muons [48], inclusive and prompt J/ψ [49, 50], prompt and non-prompt D mesons [51-54], and electrons from heavy-flavour hadron decays [55], using two-particle correlations.A positive v 2 of muons from heavy-flavour hadron decays was measured by ATLAS in high-multiplicity pp collisions at √ s NN = 13 TeV [56] as well.Long-range correlations were also observed at lower beam energies, in d-Au [57-59] and 3 He-Au collisions [60], by the PHENIX and STAR collaborations at RHIC.Finally, the PHENIX collaboration studied the pseudorapidity (η) dependence of the charged-particle v 2 at high multiplicity in various asymmetric collision systems [61] as well.The v 2 signal was found to increase from the smallest to the largest system, and was also more pronounced in the backward rapidity region (−3 < η < −1) than at forward rapidity (1 < η < 3).Several strategies were developed to subtract the correlations not related to collectivity but rather due to jet correlations and resonance decays, referred to as nonflow effects.These nonflow contributions were usually suppressed by requiring a pseudorapidity separation between particles forming the pair and by subtracting correlations measured in low-multiplicity collisions [42,62].Later on, a standard template fit procedure [63] was implemented to isolate the long-range correlations, which was further corrected to consider the multiplicity dependence of the Fourier coefficients v n [64].It was also observed that multi-particle cumulants [45, [65][66][67][68][69][70] strongly suppress nonflow correlations.
On the other hand, the simultaneous observation of a positive v 2 and particle yields compatible with those in binary-scaled pp collisions, i.e. a nuclear modification factor R pPb ∼ 1, for high-p T charged particles measured in high-multiplicity p-Pb collisions at √ s NN = 5.02 TeV, is also puzzling [63,71,72].As reported in Ref. [63], a jet quenching calculation with two different initial geometries [73] slightly underestimates the measured v 2 and fails in describing the R pPb data which are in favour of no significant energy loss effects.Several other possible scenarios (see a review in Ref. [74]) relying either on final-state effects involving a hydrodynamic evolution of the produced particles [75][76][77] or initial-state effects such as gluon saturation in the framework of the colour glass condensate (CGC) [74,78,79], are investigated.Colour-charge exchanges in the final state [78] or the anisotropic escape of partons from the surface of the interaction region [80] might also contribute to the flow-like effects evidenced in small collision systems.A MultiPhase Transport (AMPT) model [81][82][83] addresses non-equilibrium dynamics and provides a microscopic evolution of parton interactions, hence including also the parton escape mechanism.The latter scenario is further investigated using the string-melting version of the AMPT model [83], employed here for the comparison with the data (see Section 4).In this version, the initial strings are converted into partons and the interactions are described via a parton cascade model [84].The partons are then combined into hadrons via a spatial quark coalescence model and the rescatterings between hadrons are described by a relativistic transport (ART) model [85].
In order to shed more light on the origin of the azimuthal anisotropies in small systems, this letter presents new results concerning the v 2 coefficient of inclusive muons in high-multiplicity p-Pb collisions at √ s NN = 8.16 TeV with the ALICE detector at the LHC.These measurements are performed at forward (2.03 < y CMS < 3.53) and backward (−4.46 < y CMS < −2.96) rapidities and cover the transverse momentum interval 0.5 < p T < 10 GeV/c.They are obtained in a significantly extended p T range dominated by decays of heavy-flavour hadrons at p T > 2 GeV/c.The total uncertainties are reduced by a factor up to about 2.1 and 1.3 at forward and backward rapidities, respectively, compared to the previous ALICE muon results in p-Pb collisions at √ s NN = 5.02 TeV in the interval 0.5 < p T < 4 GeV/c [48].Various analysis methods which exhibit different sensitivity to nonflow effects, are implemented to measure the p T -differential v 2 of inclusive muons.The two-particle correlation method and the two-particle cumulants with generic framework [86], the latter applied for the first time for heavy-flavour v 2 measurements with ALICE, are used.In order to reduce nonflow contributions, a novel procedure based on the subtraction of correlations measured in low-multiplicity events is developed.
The letter is organised as follows.Section 2 presents the ALICE apparatus with an emphasis on the detectors used in the analysis and the data taking conditions.Section 3 contains a description of the flow methods and a presentation of the analysis details. Section 4 presents the results, namely the p T -differential muon v 2 at forward and backward rapidities, measured using several multiplicity estimators with the two-particle correlation and two-particle cumulant methods.Comparisons with published measurements performed at √ s NN = 5.02 TeV in various kinematic regions are reported.Detailed comparisons with model calculations based on the colour glass condensate and AMPT predictions are also discussed.A summary and concluding remarks are given in Section 5.

Experimental apparatus and data samples
A detailed description of the ALICE apparatus and its performance can be found in Refs.[87,88].The main subdetectors used in the analysis are presented in the following.The analysis uses muons reconstructed in the muon spectrometer which covers the pseudorapidity interval −4.0 < η < −2.5.The muon spectrometer consists of a 10 nuclear interaction length (λ I ) absorber in front of five tracking stations, each composed of two planes of cathode pad chambers, with the central one inside a dipole magnet of 3 Tm integrated field.The tracking system is completed with two trigger stations, each equipped with two planes of resistive plate chambers, behind a 1.2 m thick iron wall (7.2 λ I ).The latter stops secondary hadrons escaping from the front absorber as well as low momentum muons from light-hadron decays.A conical absorber protects the muon spectrometer throughout its full length against secondary particles produced by the interaction with the beam pipe of primary particles at large η.Among the central barrel detectors, the two innermost layers of the Inner Tracking System composing the Silicon Pixel Detector (SPD) cover the pseudorapidity intervals |η| < 2 and |η| < 1.4.The SPD is employed for the determination of the position of the primary interaction vertex.The SPD tracklets, track segments joining hits in the two SPD layers, are also used in the flow analysis, either as associated particles in the two-particle correlation method or to determine the reference flow with two-particle cumulants (Section 3).The V0 detector composed of two scintillator arrays, covers the pseudorapidity intervals 2.8 < η < 5.1 (V0A) and −3.7 < η < −1.7 (V0C).It provides the minimum bias (MB) trigger defined by the coincidence of signals in the two sets of scintillators.The V0 is also used for the luminosity determination, and an independent measurement is obtained with the two T0 arrays of Cherenkov detectors located in the regions 4.6 < η < 4.9 and −3.3 < η < −3.0.The two sets of Zero Degree Calorimeters (ZDC), each including a neutron calorimeter (ZN) and a proton calorimeter (ZP), are located on both sides of the interaction point at z = ±112.5m, along the beam line.The timing information delivered by the V0 and ZDC is exploited offline to reject the beam-induced background.These two detectors are also used to estimate the event activity.
The results reported in this letter are obtained with the data samples recorded by ALICE during the 2016 p-Pb run at √ s NN = 8.16 TeV, with two different beam configurations obtained by reversing the rotation direction of the proton and lead beams.Due to the asymmetry of the beam energy per nucleon, the nucleon−nucleon centre of mass is shifted in rapidity with respect to the laboratory frame by ∆y = 0.465 in the direction of the proton beam.Therefore, inclusive muons are measured in the forward rapidity interval 2.03 < y CMS < 3.53 with the proton beam travelling in the direction of the muon spectrometer (p-going direction, p-Pb configuration) and in the backward rapidity interval −4.46 < y CMS < −2.96 (Pb-going direction, Pb-p configuration).The analysis is based on muon-triggered events, requiring the MB condition and at least a track registered in the muon trigger stations with a p T above a programmable threshold value.Data were collected with two programmable p T thresholds set to about 0.5 GeV/c and 4.2 GeV/c, referred to as MSL and MSH, respectively.These two thresholds are not sharp, and correspond to a ∼50% efficiency for muons [89].Based on statistical considerations, the measurement of the v 2 of inclusive muons is performed by combining MSL-and MSH-triggered events which are used up to p T = 2 GeV/c and above this value, respectively.It has been checked that in the overlapping region, both samples give same inclusive muon v 2 results within uncertainties.After applying an event selection which uses the information of the V0 and SPD together with an algorithm to tag events with multiple vertices, the pile-up contribution is found negligible for the data samples considered in this analysis.Moreover, only events with a primary vertex along the beam axis z vtx within ±7 cm are considered in order to reduce non-uniform acceptance effects.After the event selection, the integrated luminosity of the p-Pb and Pb-p data samples corresponds to about 0.22 (5.8) nb −1 and 0.22 (8.2) nb −1 for MSL-(MSH-)triggered events, respectively.The p-Pb and Pb-p data samples are further classified according to their activity using multiplicity-based estimators [90] such as the total charge deposited in the two V0 arrays (V0M), the number of clusters in the outer layer of the SPD (CL1), and the energy deposited by spectator neutrons in the ZN located in the direction of the Pb beam.The ZN estimator minimises biases on the binary scaling of hard processes [90].The multiplicity classes are defined as percentile intervals of the p-Pb and Pb-p hadronic cross section.In this analysis, the 0-20% and 60-90% multiplicity classes are studied.The evolution of the charged-particle pseudorapidity density at midrapidity in different multiplicity classes and with several multiplicity estimators is reported in Ref. [91].
The same selection criteria as in previous analyses [92,93] are applied to the muon candidates.Tracks in the muon spectrometer are reconstructed within the pseudorapidity region −4 < η < −2.5 to reject tracks at the edge of the muon spectrometer acceptance.The track polar angle at the end of the front absorber θ abs has to satisfy the condition 170 • < θ abs < 178 • in order to remove tracks crossing the high-density region of the absorber that undergo significant scattering.Tracks reconstructed in the muon tracking chambers are identified as muons by requiring their matching with corresponding track segments in the muon trigger chambers.Finally, a selection on the minimum distance of the track to the primary vertex in the transverse plane (DCA) weighted by its momentum (p) is applied to reject fake tracks and remaining beam-induced background tracks.The SPD tracklets are selected in the fiducial region |η| < 1.Moreover, a condition on the pseudorapidity dependence on the longitudinal position of the primary vertex z vtx is applied for the removal of edge effects in the region where the SPD acceptance is small.As in previous publications [48,49], a condition on the difference between the azimuthal angles of clusters in the two SPD layers with respect to the primary vertex, ∆ϕ SPD < 5 mrad, is applied in order to select in average reference tracks with larger p T , which exhibit a larger flow, and to reduce the contribution from fake and secondary tracklets.
3 Analysis methods and associated systematic uncertainties

Two-particle correlations
The method using two-particle correlations to extract the azimuthal anisotropy is extensively discussed in Refs.[41-45, 47-49, 55].The two-particle correlation between pairs of trigger particles, inclusive muons at forward or backward rapidities, and associated particles, SPD tracklets at midrapidity, is measured as a function of their azimuthal angle difference (∆ϕ) and pseudorapidity difference (∆η).The correlation is expressed in terms of Y , the associated yield per trigger particle defined as where N trig is the total number of trigger particles, i.e. the number of inclusive muons in a given multiplicity class, z vtx interval, and p T interval.The signal distribution S(∆η, ∆ϕ) given by 1 corresponds to the associated yield per trigger particle for particle pairs from the same event.The background distribution B(∆η, ∆ϕ) = α d 2 N mix d∆ηd∆ϕ is obtained by correlating trigger particles in an event with associated particles from other events in the same multiplicity class and z vtx interval.The parameter α is introduced to normalise the background distribution to unity in the region of maximum pair acceptance.Both the signal and background distributions are determined considering the same multiplicity class and same z vtx interval of 1 cm width.The final associated yield per trigger particle is obtained from an average over the z vtx intervals weighted by N trig .
The distribution of the associated yield per trigger particle (Eq.( 2)) measured in high-multiplicity collisions is usually composed of correlations arising from collective and nonflow effects, the latter consisting of near-side (|∆ϕ| < π/2) and away-side (π/2 < ∆ϕ < 3π/2) jet structures.These nonflow effects can be reduced by subtracting the per-trigger yield distribution measured in low-multiplicity collisions [42].A typical example of associated yield per trigger particle for muon-tracklet correlations as a function of ∆η and ∆ϕ after the subtraction of the per-trigger yield measured in low-multiplicity (60-90%) collisions [42], labelled as Y sub , is shown in Fig. 1 (top-left panel) for high-multiplicity (0-20%) p-Pb collisions at √ s NN = 8.16 TeV.The selected p T interval of trigger particles is 2 < p T < 2.5 GeV/c.A double-ridge structure is observed with a near-side ridge and a away-side ridge centred at ∆ϕ = 0 and ∆ϕ = π, respectively.This double-ridge structure indicates the presence of collective effects in p-Pb collisions.In order to quantify these remaining correlations, the resulting two-dimensional subtracted distribution in −5 < ∆η < −1.5 is projected onto ∆ϕ and is further fitted with a Fourier series up to the third order where, as shown in Fig. 1 (top-right panel), an azimuthal anisotropy dominated by the second-order coefficient a 2 is observed.
After the subtraction of the muon-tracklet distribution measured in low-multiplicity collisions, the awayside jet yield in high-multiplicity collisions has strongly decreased.However, since the jet correlation is observed to be diminished in low-multiplicity collisions [94], a residual jet contamination in the subtracted distribution, in particular in the away side, is still present.The subtraction of this potential remaining jet component in the away side relies on the scaling of the low-multiplicity event class as described hereafter (see also Refs. [45,48,94]).The per-trigger associated yield in the low-multiplicity event class is first fitted with a Gaussian function centered at ∆ϕ = π.The per-trigger yield in the high-multiplicity event class is further adjusted with Eq. ( 3) replacing the cos(∆ϕ) term by a Gaussian function whose width is fixed to the value obtained from the fit in low-multiplicity events.The per-trigger yield distribution in low-multiplicity events is scaled by a factor k defined as the ratio of the corresponding yields in the away side in high-multiplicity collisions to those in low-multiplicity collisions and the subtraction procedure is applied [48,49].The k factors are larger in the p-going than in the Pb-going direction.They reach maximum values up to 1.43 (p-going) and 1.23 (Pb-going) with the V0M estimator.The values increase up to 1.74 (1.43) and 1.32 (1.0) in the p-going (Pb-going) direction, with the CL1 and ZN multiplicity estimators.Figure 1 presents the resulting associated yield per trigger particle as a function of ∆η and ∆ϕ and the projected distribution onto ∆ϕ in bottom-left and bottom-right panels, respectively.One can notice that the amplitudes of the near-side and away-side ridges become comparable.
The subtraction of the scaled per-trigger yield in low-multiplicity collisions as well as the fit of the projected ∆ϕ distribution with Eq. ( 3) are repeated for each p T trigger-particle interval.The values of the reduced χ 2 are all smaller than 1.5, confirming that the distributions are well described by this Fourier series.The V can be factorised as the product of the anisotropies of single muons and SPD tracklets [62], the single muon second-order coefficient (v µ 2 {2PC}) can be expressed as In order to extract the V tracklet−tracklet 2∆ , the analysis is repeated by correlating only SPD tracklets, as also done in Refs.[42,44], and the region |∆η| < 1.2 is excluded to reduce a bias due to a possible remaining jet peak on the near side at (∆η ∼ 0, ∆ϕ ∼ 0) after the subtraction of the correlation distribution in low-multiplicity collisions.
The systematic uncertainties affecting the v µ 2 {2PC} coefficient originate from the measurement of the muon-tracklet and tracklet-tracklet correlation distributions.They arise from the selection criteria of SPD tracklets, the procedure to remove the jet contamination in the near side and away side, the v µ 2 {2PC} calculation including the fit stability and baseline determination, and the effect of the muon-track resolution.
The effect of the SPD acceptance is investigated by varying the range of the position of the reconstructed vertex along the beam axis, which is decreased down to |z vtx | < 5 cm.
The systematic effect related to the procedure implemented to subtract the jet contamination is investigated.A residual near-side jet peak in the tracklet-tracklet correlation is negligible since the V tracklet−tracklet 2∆ coefficient remains unchanged when varying the |∆η| gap in the range 0.8-1.2units.The effect of a possible incomplete away-side jet subtraction is studied with the scaling procedure of the correlations in low-multiplicity collisions prior to its subtraction, as previously discussed.The influence of such an effect is also investigated fitting the subtracted correlation distributions, (0-20%) − (60-90%), with Eq. ( 3) replacing the first order term by a Gaussian function with the width fixed to the one extracted from the fit of correlations in low-multiplicity collisions in the away side.The differences between the two procedures give an estimation of the corresponding systematic uncertainty.A potential bias could result from long-range correlations remaining in the 60-90% low-multiplicity class, which may lead to The stability of the v µ 2 {2PC} results is also investigated by excluding the third-order coefficient in the fit of the muon-tracklet and track-tracklet correlations with Fourier series.The quality of the fits is getting worse, although no significant change in the extracted v µ 2 {2PC} values is seen.A systematic uncertainty also arises from the procedure employed for the ∆ϕ projection.The baseline is a parameter which could influence the fits and is estimated following the same strategy as in Ref. [48].The latter is obtained by fitting the correlation distributions in low-multiplicity collisions using a Gaussian function in the away side and a constant for the baseline.Alternatively, the baseline can also be calculated in high-multiplicity collisions from the integral or from a second-order polynomial fit around the minimum at ∆ϕ ∼ π/2.Finally, the ∆ϕ projection is obtained from a constant fit instead of a first-order polynomial fit along ∆η for each ∆ϕ interval.
The systematic effect due to the angular and momentum resolution of reconstructed muon tracks is evaluated by means of a dedicated Monte Carlo simulation based on the DPMJET event generator [95], which uses the GEANT4 transport code [96,97] and the afterburner flow technique [98].
The various systematic uncertainties coming from each source are reported in Table 1 at forward (pgoing) and backward (Pb-going) rapidities for the V0M multiplicity estimator.They are added in quadrature to obtain the overall systematic uncertainty on v µ 2 {2PC}.Comparable uncertainty values are estimated in the 0-20% event class selected with CL1 and ZN multiplicity estimators.

Two-particle cumulants
Two-particle cumulants, which exhibit different sensitivities to nonflow compared to the two-particle correlation method, are also employed for the study of the second-order v 2 coefficient of inclusive muons [23].For the first time in the heavy-flavour sector, the muon v 2 is measured with the two-particle cumulants using the framework presented in Ref. [86] where a method to correct for non-uniform acceptance and inefficiencies is provided.The SPD tracklets (reference particles, RP) are selected under the same selection criteria as with two-particle correlations.Apart from that, a more restrictive selection of the primary vertex position along the beam direction of −5 < z vtx < 3 cm is considered to account for dead zones in the SPD detector in the plane z vtx − ϕ, which cannot be corrected with particle weights.The standard selections discussed in Section 3.1 are applied to identify muon tracks, the so-called particles of interest (POI).The two-particle second-order reference and p T -differential cumulants, c 2 {2} and d µ 2 {2}(p T ), are based on the construction of weighted Q and p vectors1 [86], defined on an event-by-event basis for a given multiplicity class and data sample (MSL-or MSH-triggered events) as where n is the n th harmonic number, l is an integer exponent of the k th -particle (RP or POI) weight w l k , M is the multiplicity of the SPD tracklets, m p is the muon multiplicity, and ϕ k is the SPD tracklet (muon) azimuthal angle needed for the computation of Q n,l (p n,l ).After applying the weights for both SPD tracklets and inclusive muons, the non-uniformities in the azimuthal acceptance and inefficiencies are found negligible.
The calculation of the p T -differential second-order coefficient of inclusive muons v µ 2 {2} in a given multiplicity interval is performed as c 2 {2} being related to the reference second-order coefficient V 2 {2} as The measurement of the muon p T -differential second-order coefficient with two-particle cumulants is influenced by nonflow effects which need to be isolated and subtracted.The v µ 2 {2} coefficient can be extracted from the equation As done in the two-particle correlation method [48], the long-range jet correlations are estimated in the low-multiplicity event class 60-90%.Then, the two-particle second-order reference and differential cumulants in these low-multiplicity collisions are scaled by a factor f defined as the ratio of the mean SPD-tracklet multiplicity in low-multiplicity collisions to that in high-multiplicity collisions and subtracted as shown in Eq. ( 8).Short-range jet correlations are usually suppressed applying a pseudorapidity gap between the correlated particles, i.e. for both muons and SPD tracklets.As a pseudorapidity gap between muons and SPD tracklets is naturally present, the procedure cannot be applied straightforwardly for SPD tracklets due to the dead zones in the SPD acceptance.Therefore, these nonflow effects are suppressed via the scaling of the muon v µ 2 {2} by another factor, f ∆η , estimated by means of AMPT simulations [81,99] and defined as the ratio of v µ 2 {2} extracted with |∆η| > 0.4 to that obtained without applying a pseudorapidity gap.The optimised condition |∆η| > 0.4 results from a compromise between statistical considerations and suppression of nonflow effects.The value of this scaling factor is about 1.2 and is found to be independent of p T and the multiplicity estimator.
As already discussed with the technique of two-particle correlations (see Section 3.1), remaining longrange jet correlations are not excluded even after the subtraction of correlations in low-multiplicity collisions.These remaining jet correlations can be suppressed by applying two additional correction factors.The factor g, applied to the p T -differential cumulants measured in low-multiplicity collisions, is estimated from the muon-tracklet correlation distribution as discussed in Section 3.1.The factor f RP , applied to V 2 {2}, is calculated from the tracklet-tracklet correlation function.This is the ratio of V 2 {2} from tracklet-tracklet correlations extracted with the scaling of the remaining jet contribution to that obtained without any scaling procedure.rapidities, using V0M for the multiplicity estimation.Large deviations appear at high p T , indicating the importance to subtract long-range jet correlations at both forward and backward rapidities.These trends could result from the larger long-range jet contamination in the p-fragmentation region (forward rapidity, left panel) compared to the Pb-fragmentation region (backward rapidity, right panel) at high p T , in particular.A similar behaviour is also seen with CL1 and ZN multiplicity estimators.
Several potential sources of systematic uncertainties affecting the measurement of the v µ 2 {2} coefficient are considered.
The non-uniform acceptance correction for SPD tracklets is expected to be independent of the analysed data sample.It is tested with the MB and muon-triggered samples.The difference of the results obtained with muon-triggered events with respect to the MB sample is taken as the corresponding systematic uncertainty.
The sensitivity of the results to the selection criteria of SPD tracklets is investigated by varying the η interval of SPD tracklets from |η| < 1 to |η| < 1.2.
The systematic uncertainty related to the procedure implemented for the subtraction of short-range jet correlations includes two components.A first contribution to the systematic uncertainty comes from the variation of the pseudorapidity gap up to |∆η| > 0.8.A second contribution is obtained conservatively from the comparison of the v µ 2 {2} computed with the scaling factor f ∆η estimated from AMPT simulations and the measured one, using |∆η| > 0 ignoring that the effects due to non-uniformities in the azimuthal acceptance depend on z vtx .In order to study the uncertainty from the away-side jet subtraction in the calculation of the p T -differential cumulant, the subtracted ] (60−90%) , i.e. the numerator of Eq. ( 8) with g set to unity, is scaled by a factor derived from the two-particle correlation method.This factor is estimated as the ratio of the V µ−tracklet 2∆ extracted fitting the subtracted correlation with Eq. (3) to that obtained replacing the first order term with the Gaussian function discussed in Section 3.1.
The systematic effect coming from the possible influence of long-range correlations in low-multiplicity events is assessed by changing the low-multiplicity interval from 60-90% to 70-90%.
Finally, the systematic uncertainty related to the angular and momentum resolution of the muon spectrometer is evaluated following the same strategy as with two-particle correlations (see Section 3.1).TeV.The results at forward and backward rapidities are depicted in the left and right panel, respectively.They are reported for the first time in a wide transverse momentum interval 0.5 < p T < 10 GeV/c.The event activity is determined with the V0M multiplicity estimator.Comparisons with the v µ 2 coefficient of inclusive muons measured in high-multiplicity events selected using CL1 and ZN multiplicity estimators will be discussed later.The measurements are carried out with two-particle correlations (open symbols) and two-particle cumulants (full symbols), which are denoted v µ 2 {2PC} and v µ 2 {2}, respectively.The statistical uncertainties (vertical bars) with the two-particle cumulants are estimated according to the Jackknife method [100].The latter is based on the resampling technique, where 6 out of 12 sub-samples are used to extract the RMS of the distribution for each p T interval which is further divided by √ 2 in order to obtain the statistical uncertainty.In the two-particle correlation method, the statistical uncertainty is obtained from the Fourier fit procedure.The systematic uncertainties are the empty and filled boxes.

Results and model comparisons
The two methods give compatible results within uncertainties at both forward and backward rapidities.One observes first an increase of the v µ 2 signal with increasing p T where it reaches a maximum value of about 0.07 (0.08) at p T ∼ 2 GeV/c in the forward (backward) rapidity region, followed by a decrease as p T increases.Although the p T dependence is similar in the two rapidity regions, there is a hint for a higher elliptic flow signal at backward rapidity than at forward rapidity as observed in Ref. [48].This is further investigated by the measurement of the p T -differential ratio of the measured v µ 2 {2PC} at backward rapidity to that at forward rapidity.This ratio exhibits a uniform behaviour within uncertainties and can be fitted with a constant function (not shown here).It amounts to 1.28 ± 0.07, the uncertainty being from the fit.This asymmetry could be a consequence of decorrelation effects of the flow vectors in different rapidity regions [101][102][103][104].
In the p T interval of interest, muons originate mainly from charged-pion and kaon decays at low p T (p T < 2 GeV/c), while muons from heavy-flavour hadron decays dominate over muons from lighthadron decays at higher p T [48].The contributions of these muon sources are estimated by means of simulations using the DMPJET event generator [95] and the GEANT4 transport package [96].The relative contribution of muons from primary charged-pion and kaon decays amounts to about 67% (68%) at 0.5 < p T < 1 GeV/c, and it decreases with increasing p T down to about 7.5% (14.5%) at 6 < p T < 10 GeV/c in the p-going (Pb-going) direction.The fraction of muons originating from charm and beauty decays represents about 60% (58%) of the total muon yield at p T = 2 GeV/c and it reaches about 87% (78%) in the p-going (Pb-going) direction at 6 < p T < 10 GeV/c.A similar fraction of muons from heavy-flavour hadron decays in 0.5 < p T < 4 GeV/c is reported for p-Pb collisions at √ s NN = 5.02 TeV [48].Moreover, it is worth to mention that based on fixed-order plus next-to-leading logarithms (FONLL) calculations [105,106] more than 60% of muons from heavy-hadron decays originate from beauty quarks in the highest p T interval (6 < p T < 10 GeV/c).The measured v µ 2 coefficient is positive with a significance which reaches values of 4.7σ -12σ (7.6σ -11.9σ ) at intermediate p T (2 < p T < 6 GeV/c) in the forward (backward) rapidity region, depending on the analysis technique.These results might suggest the existence of a collective behaviour of heavy quarks in high-multiplicity (0-20%) p-Pb collisions at √ s NN = 8.16 TeV at forward and backward rapidities.In the highest p T interval accessible in this analysis (6 < p T < 10 GeV/c), where muons from beauty-hadron decays take over charm as the dominant muon component, there is a hint for a positive v µ 2 although not significant within uncertainties (significance of 0.2σ (1.2σ ) at forward (backward) rapidity with two-particle correlations).
The analysis is repeated using the number of clusters in the outer layer of the SPD (CL1) and the energy deposited in the neutron ZDC (ZN) for the event activity selection.Figure 4 displays the measurements obtained with two-particle cumulants (top panels) and two-particle correlations (bottom panels) at forward (left panels) and backward (right panels) rapidities.It has been demonstrated that these multiplicity estimators select different event classes associated with different mean charged-particle multiplicity [71].The ZN is expected to be the least-biased estimator (see Section 2 and Ref. [90]) but the correlation between the energy measured with the ZDC and the charged-particle multiplicity at midrapidity is known to be weak [71].On the other hand, autocorrelation effects are present when using CL1 for the event class selection (Section 2) since the reference flow is also calculated with the SPD tracklets.Even in this context, the computed v µ 2 values are compatible within uncertainties, although there is a hint for smaller v µ 2 values with the ZN multiplicity estimator than with V0M and CL1 multiplicity estimators.The ZN estimator selects on the average smaller charged-particle multiplicity density in the high-multiplicity class than V0M and CL1 estimators, while the opposite trend is observed for the low-multiplicity class [90,107].Consequently, a smaller v µ 2 signal is expected after the subtraction of correlations from low-multiplicity events with the ZN estimator.
The present results are compared in Fig. 5 with previously published results of inclusive muons obtained   The results are also in good agreement within uncertainties with those obtained by ALICE for electrons from heavy-flavour hadrons measured in p-Pb collisions at √ s NN = 5.02 TeV at midrapidity (−0.8 < η < 0.8) and for 1.5 < p T < 6 GeV/c [55].It is interesting to point out that the magnitude of the v 2 of inclusive muons is comparable within uncertainties to the one of inclusive muons obtained at forward rapidity by ALICE in semicentral Pb-Pb collisions at √ s NN = 2.76 TeV for 2 < p T < 10 GeV/c [23].A positive v 2 also was reported by ALICE for J/ψ measured in 3 < p J/ψ T < 6 GeV/c and same rapidity intervals in p-Pb collisions at √ s NN = 8.16 TeV, where the contribution from recombination of thermalised charm quarks in the medium is expected to be negligible and path-length dependent effects are smaller with respect to Pb-Pb collisions [49].The present inclusive muon v 2 results also complement those obtained in high-multiplicity p-Pb collisions at √ s NN = 5.02 TeV for unidentified particles as well as pions, kaons and protons, although in a different kinematic region (different p T and y intervals) [108,109].
Further insights in the understanding of the observed azimuthal anisotropies in small collision systems can be gained by comparing the measurements with model predictions such as the AMPT and colour glass condensate (CGC) calculations.The v2-26t7b string-melting version of the AMPT model [81,82] which includes the improvements discussed in Ref. [83] is employed to compute the v 2 of heavy-flavour hadrons (D 0 and B mesons) and primary charged pions and kaons by means of the forward-central twoparticle corrections (Section 3.1).The event selection is performed by counting the charged particles in the acceptance of the V0 detector.The v 2 coefficient is then computed separately for muons from charm hadrons, beauty hadrons, charged pions and charged kaons by means of fast simulations which use the input p T and v 2 distributions of D 0 and B mesons, and primary charged pions and kaons, as well as PYTHIA 6.4 [110] for the decay kinematics.The D-meson species are assumed to have the same v 2 coefficient as D 0 mesons2 .The p T -differential inclusive muon v µ 2 is obtained from a weighted sum of the v 2 coefficient of muons from heavy-flavour hadron decays, v µ←b,c 2 , and the v 2 of muons from charged pion and kaon decays, , f being the relative abundance of muons from the decay of charged pions and kaons estimated from Monte Carlo simulations with the DPMJET event generator [95].Similarly, the v µ←b,c 2 coefficient is computed as a weighted sum of the v 2 coefficient of muons of charm-hadron and beauty-hadron decays, v , as The corresponding fractions of muons from charm-hadron and beauty-hadron decays with respect to the total yield of muons from heavy-flavour hadron decays, f c and f b , are obtained by means of the fixedorder plus next-to-leading logarithms (FONLL) approach [105,106].
Figure 6 presents a comparison of the p T -differential inclusive muon v µ 2 {2PC} with AMPT calculations, together with the different contributions of muons from charm-and beauty-hadron decays, separately, and muons from charged pion and kaon decays in the p-going (left) and Pb-going (right) direction.The AMPT model generates a positive v 2 for all particle species.Its magnitude increases significantly in the low-p T region up to about 2-3.5 GeV/c depending on the decay particle.At higher p T , the v 2 signal decreases smoothly with increasing p T or saturates, except in the Pb-going direction where a slightly increase with p T is seen for muons from beauty-hadron decays.In line with hydrodynamic calculations [111], the AMPT model predicts a larger v 2 for muons from charged-pion and kaon decays than for muons from heavy-flavour hadron decays at low p T in high-multiplicity p-Pb collisions.As observed with the data, the calculated v µ 2 values with the AMPT model are larger in the backward rapidity region compared to the forward rapidity region, which may be a consequence of rapidity-dependent flowvector fluctuations [101][102][103][104].The AMPT predictions are in fair agreement with the measured inclusive muon v µ 2 , although the model tends to slightly overestimate the data in the backward rapidity region.It is important to note that finite v µ←b,c 2 values are indeed needed to reach such agreement.These AMPT comparisons suggest that the azimuthal anisotropies are mainly driven by the anisotropic parton escape mechanism where partons have a higher probability to escape along the shorter axis of the interaction zone, as discussed in Ref. [80].TeV with AMPT calculations [81][82][83].The predictions are shown for muons from charm-hadron decays and beauty-hadron decays, separately, muons from charged-pion and kaon decays, and for muons from the combination of the various sources.The contribution of muons from both charm-and beauty-hadron decays is also displayed.[112,113] where interactions between partons from the proton projectile and dense gluons inside the target Pb nucleus at the early stage of the collision generate azimuthal anisotropies.The comparison is performed only for the p-going direction where the dilute-dense formalism is valid [114].The predictions of the v 2 coefficient have been provided separately for D 0 and B mesons from two-particle correlations.The v 2 coefficient of muons from heavy-flavour hadron decays is further obtained implementing the same strategy as with the AMPT predictions.One observes that in these CGC-based calculations, the correlations in the initial state generate a significant v 2 signal for muons from charm-hadron decays.Its magnitude increases significantly up to p T = 2 GeV/c where it reaches a maximum value of about 0.09, and it decreases smoothly with increasing p T .The predicted v 2 signal of muons from beauty-hadron decays is less pronounced, the maximum being of about 0.03 at p T = 3 GeV/c.In the highest p T region, the v 2 coefficient of muons from charm-hadron decays is similar to that of muons from beauty-hadron decays.The CGC-based calculations for muons from both charmand beauty-hadron decays reproduce qualitatively the measured v µ 2 coefficient of inclusive muons for p T > 2 GeV/c.It is worth pointing out that the contribution of muons from pion and kaon decays which contributes significantly to the measured inclusive muon yield at low p T (p T < 2 GeV/c) is not considered in the model predictions, preventing to draw any conclusion about the model comparisons in that kinematic region [48].A similar agreement with these CGC calculations is found for prompt D 0 , prompt J/ψ, and D 0 from beauty-hadron decays measured at midrapidity in p-Pb collisions at √ s NN = 8.16 TeV with the CMS detector [54].The qualitative agreement between data and model calculations may suggest possible contributions from initial-state effects to the measured positive v 2 in high-multiplicity p-Pb collisions.However, it should be mentioned that the measured inclusive muon v µ 2 is obtained using charged particles for the calculation of the reference v 2 coefficient, which is in general interpreted as a final-state effect [103].
The CGC and AMPT predictions for muons from charm-and beauty-hadron decays are displayed together in the right panel of Fig. 7 for the p-going direction, and compared to the measured inclusive muon v µ 2 .The CGC-based calculations provide a larger v 2 signal for muons originating from heavyflavour hadron decays at low p T , up to about p T = 3 GeV/c, compared to the AMPT calculations.Such behaviour indicates that the CGC calculations which do not incorporate muons from light-flavour hadron decays would have overestimated the data in this kinematic region.The two models provide compatible results and describe the data at high p T , where the muons from charged-pion and kaon decays do not affect the v µ 2 significantly.

Conclusion
The second-order coefficient v 2 coefficient is extracted using several event activity estimators with the well-known technique of two-particle correlations and, for the first time for open heavy-flavoured particles, with two-particle cumulants within the generic framework.Nonflow effects which cannot be neglected in the analyses involving correlations between two particles are subtracted using novel techniques.A positive v µ 2 signal is observed at forward and backward rapidities, with a significance larger than 4.7σ and 7.6σ , respectively, in the region 2 < p T < 6 GeV/c where muons from heavy-flavour hadron decays constitute the main source of muons.The results indicate that heavy quarks reveal a collective-like behaviour in high-multiplicity p-Pb collisions.The AMPT calculations which complement hydrodynamic models and address non-equilibrium dynamics provide a reasonable agreement with the data at both forward and backward rapidities over the whole p T interval.The measured v µ 2 coefficient at forward rapidity (p-going direction) for p T > 2 GeV/c is in qualitative agreement with CGC calculations within the dilute-dense formalism.Possible contributions from initial-state effects are not fully excluded at high p T .These comprehensive results on v µ 2 of inclusive muons spanning in a wide p T range, provide significant new insights for the understanding of the origin of the possible collective behaviour of heavy quarks in small systems such as p-Pb collisions and can be used to constrain the various approaches for modelling the azimuthal anisotropies in small collision systems.
µ−tracklet n∆ Fourier coefficients are further extracted from the fit parameters according to V µ−tracklet n∆ = a n /(a 0 + b), b corresponding to the baseline of the scaled 60-90% low-multiplicity class estimated from the integral in ∆ϕ of the correlation distribution around the minimum.By assuming that V µ−tracklet 2∆

Figure 2 Figure 2 :
Figure 2 presents the two-particle second-order p T -differential muon cumulant in high-multiplicity (0-20%) collisions without and with the subtraction of the scaled p T -differential muon cumulant in lowmultiplicity collisions and remaining long-range jet correlations at forward (left) and backward (right)

µ 2 {2} 2 vFigure 3 :
Figure 3: Inclusive muon v µ 2 as a function of p T at forward (left) and backward (right) rapidities in highmultiplicity p-Pb collisions at √ s NN = 8.16 TeV.The event activity is estimated with the V0M estimator.Open and full symbols refer to the measurements with two-particle correlations and two-particle cumulants, respectively.

Figure 3
Figure 3 presents the p T -differential second-order coefficient v µ 2 of inclusive muons after the subtraction of nonflow effects in the 0-20% high-multiplicity class for p-Pb collisions at √ s NN = 8.16 TeV.The results at forward and backward rapidities are depicted in the left and right panel, respectively.They are reported for the first time in a wide transverse momentum interval 0.5 < p T < 10 GeV/c.The event activity is determined with the V0M multiplicity estimator.Comparisons with the v

Figure 4 :
Figure 4: Inclusive muon v µ 2 as a function of p T at forward (left) and backward (right) rapidity in high-multiplicity p-Pb collisions at √ s NN = 8.16 TeV, extracted with two-particle cumulants (top) and two-particle correlations (bottom).The results are obtained with three different estimators of the event activity: V0M, CL1 and ZN.

Figure 5 :
Figure 5: Comparison of the p T -differential v 2 {2PC} of inclusive muons at forward (left) and backward (right) rapidity in high-multiplicity p-Pb collisions at √ s NN = 8.16 TeV, extracted with the two-particle correlation method with previous measurements performed in p-Pb collisions at √ s NN = 5.02 TeV for inclusive muons [48] and electrons from heavy-flavour hadron decays [55].

Figure 6 :
Figure 6: Comparison of the p T -differential v µ 2 {2PC} of inclusive muons at forward (left) and backward (right) rapidities in high-multiplicity p-Pb collisions at √ s NN = 8.16TeV with AMPT calculations[81][82][83].The predictions are shown for muons from charm-hadron decays and beauty-hadron decays, separately, muons from charged-pion and kaon decays, and for muons from the combination of the various sources.The contribution of muons from both charm-and beauty-hadron decays is also displayed.

Figure 7 (
Figure 7 (left) shows a comparison of the measured p T -differential muon v µ 2 {2PC} coefficient at forward rapidity with calculations based on the colour glass condensate (CGC) framework.The latter uses the dilute-dense formalism[112,113] where interactions between partons from the proton projectile and dense gluons inside the target Pb nucleus at the early stage of the collision generate azimuthal anisotropies.The comparison is performed only for the p-going direction where the dilute-dense formalism is valid[114].The predictions of the v 2 coefficient have been provided separately for D 0 and B mesons from two-particle correlations.The v 2 coefficient of muons from heavy-flavour hadron decays is further obtained implementing the same strategy as with the AMPT predictions.One observes that in these CGC-based calculations, the correlations in the initial state generate a significant v 2 signal for muons from charm-hadron decays.Its magnitude increases significantly up to p T = 2 GeV/c where it reaches a maximum value of about 0.09, and it decreases smoothly with increasing p T .The predicted v 2 signal of muons from beauty-hadron decays is less pronounced, the maximum being of about 0.03 at p T = 3 GeV/c.In the highest p T region, the v 2 coefficient of muons from charm-hadron decays is similar to that of muons from beauty-hadron decays.The CGC-based calculations for muons from both charmand beauty-hadron decays reproduce qualitatively the measured v

Figure 7 :
Figure 7: Left: comparison of the p T -differential v µ 2 {2PC} of inclusive muons at forward rapidity in highmultiplicity p-Pb collisions at √ s NN = 8.16 TeV with CGC-based calculations[112,113].The predictions are shown for muons from charm-hadron decays and beauty-hadron decays, separately, and for muons from the combination of the two sources.Right: comparison of the p T -differential v 2 of muons from heavy-flavour hadron decays as obtained with CGC and AMPT calculations with the measured inclusive muon v 2 .

µ 2 of
inclusive muons in high-multiplicity p-Pb collisions at √ s NN = 8.16 TeV is measured with the ALICE detector at the LHC at forward (2.03 < y CMS < 3.53) and backward (−4.46 < y CMS < −2.96) rapidities.These new measurements extend previous inclusive muon ALICE results in p-Pb collisions at √ s NN = 5.02 TeV to a significantly broader transverse momentum interval of 0.5 < p T < 10 GeV/c.The v µ

Table 1 :
Summary of absolute systematic uncertainties affecting the v NN = 8.16 TeV in p-going and Pb-going directions.The values are reported for the event class selected with the V0M multiplicity estimator.The systematic uncertainties vary within the indicated intervals depending on the inclusive muon p T .

Table 2 :
Summary of absolute systematic uncertainties affecting the v NN = 8.16 TeV in p-going and Pb-going directions.The values are reported for the event class selected with the V0M mutiplicity estimator.The systematic uncertainties vary within the indicated intervals depending on the inclusive muon p T .
µ 2 {2} determined with CL1 or ZN mutiplicity estimators follow the same trends as a function of p T and their values are compatible.