Observation of partonic flow in proton—proton and proton—nucleus collisions

TL;DR

Scientists observed quark-gluon plasma-like behavior in small collision systems, suggesting partonic flow occurs even in proton-proton and proton-nucleus collisions.

Contribution

First observation of quark-level anisotropic flow grouping in high-multiplicity proton–lead and proton–proton collisions.

Findings

Distinctive grouping of anisotropic flow for baryons and mesons observed in high-multiplicity proton–lead and proton–proton collisions.

Results are consistent with hydrodynamic flow followed by quark coalescence, indicating partonic flow in small systems.

Abstract

Quantum Chromodynamics predicts a phase transition from hadronic matter to quark–gluon plasma (QGP) at high temperatures and energy densities, where quarks and gluons (partons) are no longer confined within hadrons. The QGP forms in ultrarelativistic heavy-ion collisions. Anisotropic flow coefficients, quantifying the azimuthal expansion of produced matter, probe QGP properties. Flow measurements in high-energy heavy-ion collisions show a distinctive grouping of anisotropic flow for baryons and mesons at intermediate transverse momentum – a feature associated with flow imparted at the quark level, confirming QGP existence. The observation of QGP-like features in proton–proton and proton–ion collisions has sparked debate about QGP formation in smaller systems. For the first time, we demonstrate the distinctive grouping of anisotropic flow for baryons and mesons in high-multiplicity proton–lead and proton–proton collisions at the Large Hadron Collider (LHC). These results are described by a model including hydrodynamic flow followed by hadron formation via quark coalescence, consistent with the formation of partonic flowing systems in these collisions.

The authors from the ALICE collaboration identify multiple species of mesons and baryons and measure the anisotropic flow with non-flow removal techniques in pp and p-Pb collisions at the LHC, identifying the hallmark of quark flow associated with an expanding quark-gluon plasma.

Figures

Illustration of the overlap region^[67](#CR67)^.A schematic representation of the overlap region in a collision is shown in gray, along with overall particle emission patterns in the transverse (x-y) plane, represented by large arrows. a Non-flow sources: These are independent emissions, such as those from resonance decays or jets, where jets are collimated streams of hadrons created when a high-energy quark or gluon fragments after a collision. These effects lead to few-particle correlations but are not related to collective behavior in the system and have been subtracted from the final aniso

Full text

Introduction

Ultrarelativistic collisions of heavy ions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) create the quark–gluon plasma (QGP), a short-lived state of strongly interacting partonic matter, thought to have existed a few microseconds after the Big Bang^1^. The interactions among partons in the QGP, combined with the initial spatial anisotropy of the overlap region of colliding ions, create anisotropic pressure gradients in the transverse plane of the collision. These anisotropic pressure gradients result in momentum anisotropy of the emitted particles^2^. The anisotropic particle emission is quantified using the Fourier decomposition of the azimuthal distribution of the final state particles^3^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{{\rm{d}}}}N}{{{{\rm{d}}}}\varphi }\propto 1+\sum\limits_{n}2{v}_{n}({p}_{{{{\rm{T}}}}})\cos (n\varphi -n{\Psi }_{n}).$$\end{document}

Here, φ and pT denote the azimuthal angle and transverse momentum of the emitted particles, respectively, while Ψ_n_ is the azimuthal angle of the symmetry plane for the n-th harmonic. The largest contributions are the second and third Fourier coefficients, namely elliptic (v2) and triangular (v3) flow^4–6^, which result from the elliptic and triangular shapes in the initial overlap region of the colliding nuclei. The anisotropic flow extends along pseudorapidity ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta \equiv -{\rm ln}\left(\tan \frac{\theta }{2}\right)$$\end{document} ), where θ is the polar angle of the particle), forming an elongated structure known as the ridge^7^. In heavy-ion collisions, precise measurements of vn^4,6,8,9^ and detailed comparisons with models employing relativistic viscous hydrodynamics reveal that the QGP behaves as a liquid with a viscosity to entropy density ratio close to the lowest theoretical value allowed^10,11^.

In high-energy heavy-ion collisions, the v2 coefficient of identified hadrons exhibits a characteristic mass dependence at low pT, meaning that more massive particles show lower v2 values at a given pT^12–14^. Mass ordering arises from the interplay between average radial expansion velocity, anisotropic flow velocity, and thermal motion, which pushes heavier particles to higher pT^15,16^. This results in a mass-dependent reduction in v2 at low pT (pT < 3.0 GeV/c). In the intermediate pT region (3.0< pT < 8.0 GeV/c), a clear separation between the flow patterns of baryons (hadrons composed of three quarks or three antiquarks) and mesons (hadrons composed of quark–antiquark pairs) is observed with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${v}_{2}^{{{{\rm{baryons}}}}}$$\end{document} > \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${v}_{2}^{{{{\rm{mesons}}}}}$$\end{document} ^14,17,18^. A physical process that can explain this distinctive grouping of hadron v2 based on their valence quark number is hadron formation via quark coalescence^19,20^. In this process, a meson (baryon) is formed by combining two (three) quarks, and the meson (baryon) v2 is obtained by combining the v2 values of the two (three) quarks, as illustrated in Fig. 1. The experimental observation of baryon-meson grouping at intermediate pT is therefore interpreted as a consequence of a medium that includes a phase with collectively flowing partons.Fig. 1. Illustration of the overlap region^67^.A schematic representation of the overlap region in a collision is shown in gray, along with overall particle emission patterns in the transverse (x-y) plane, represented by large arrows. a Non-flow sources: These are independent emissions, such as those from resonance decays or jets, where jets are collimated streams of hadrons created when a high-energy quark or gluon fragments after a collision. These effects lead to few-particle correlations but are not related to collective behavior in the system and have been subtracted from the final anisotropic flow measurements (see Correlation function and template fit method in the methods subsection for details). b Anisotropic flow: This illustrates the development of anisotropic flow in a partonic system, propagated to the level of hadrons via the quark coalescence process, which describes the experimental measurements in the intermediate pT range (~3-8 GeV/c). In this process, two or three flowing partons coalesce to form mesons or baryons, which then interact with each other. The large arrows represent the overall anisotropy of particle emission in the transverse plane, with stronger expansion along the short (x) axis.

Proton–proton and proton–nucleus collisions were used as a baseline to study the QGP in heavy-ion collisions at RHIC and the LHC, as QGP formation was not expected in small collision systems. However, striking similarities have been observed between numerous observables in both small collision systems and heavy-ion collisions at RHIC and LHC energies. These observables include the ridge^21–24^, mass dependence of v2 at low pT^22,25,26^, azimuthal angle correlations carried by multiple particles^27^, and strangeness and baryon enhancement with increasing multiplicity^28^. These features are commonly considered indicators of QGP formation. The standard picture in heavy-ion collisions is that anisotropic flow is built up after the collision through final-state interactions among the partons combined with the initial spatial anisotropy of the overlap region of the colliding nuclei. In small collision systems, where the system evolution is shorter than in heavy-ion collisions and a QGP phase is not expected, a different scenario is proposed within the framework of Color Glass Condensate (CGC) effective theory^29^. According to this theory, the observed flow patterns come from the initial gluon momentum correlations in the colliding hadrons. These gluons scatter off specific regions, or color domains, during the collision and get a momentum boost in the same direction if they scatter from the same color domain. The current understanding is that initial-state momentum anisotropy alone cannot explain the existing data, and the measurements seem to favor the scenario of final-state effects driven by initial geometry^1,24^, following a similar scenario as in heavy-ion collisions. However, the impact of initial gluon momentum correlations on the development of anisotropic flow in small collision systems is not clear yet. At the same time, the precise mechanisms underlying the final-state effects remain unclear. The flow can develop during a partonic phase, transforming the initial spatial anisotropy into the measured flow^30–35^, or it can originate via other mechanisms without the need for a deconfined phase, such as rescatterings among hadrons^36^, via approaches involving initial state effects^37^, or via different string dynamics implemented in the PYTHIA 8 event generator^38^. None of the measurements performed so far have been able to give a clear answer to this question.

This work takes a step further in resolving the puzzle of the origin of collective flow in small collision systems by investigating the possibility of a partonic phase and its role in the system’s dynamic evolution. Utilizing the unique particle identification capabilities of the ALICE detector at the LHC^39,40^, the elliptic flow (v2) as a function of pT is presented for mesons (π^±^, K^±^, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} ) and baryons (p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ) in pp and p–Pb collisions at a nucleon–nucleon center-of-mass energy ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{{s}_{{{{\rm{NN}}}}}}$$\end{document} ) of 13 TeV (pp) and 5.02 TeV (p–Pb). In both pp and p–Pb systems, collisions are categorized into high-multiplicity (HM) and low-multiplicity (LM) events based on the charged particles detected within the pseudorapidity ranges 2.8 < η < 5.1 and  − 3.7 < η < − 1.7, respectively. Additionally, a selection criterion on the number of reconstructed and efficiency-corrected charged particles (Nch) with 0.2 < pT < 3.0 GeV/c at midrapidity (∣η∣ < 0.8) is applied, resulting in the same average Nch (〈Nch〉 ≈ 35) in high-multiplicity events for both p-Pb and pp collisions. The event, track selection, and particle identification are discussed in the Methods section.

As Ψ_n_ used in Eq. (1) cannot be experimentally determined, the v2 of hadrons can be obtained using two-particle azimuthal angle correlations (2PC)^22^ via the three-subevent method^41^. This method employs reference particles selected in the forward and backward rapidity regions in addition to the identified hadrons selected at midrapidity, allowing a significant pseudorapidity separation between the two correlated particles, 1.1 < ∣Δη∣ < 7.8 in p–Pb, and 1.1 < ∣Δη∣ < 6.4 in pp collisions. The difference in the maximum available ∣Δη∣ separation between pp and p–Pb collisions is due to the different performance of the forward/backward detectors in different data-taking periods. However, for each collision system, the final results are estimated with varying ∣Δη∣ separations and found to be consistent with each other. This large Δη separation suppresses the non-flow contamination effects due to the azimuthal angle correlations between particles produced from resonance decays and jets^42^. To further minimize the remaining non-flow effects, a template fit method^23,43^ is employed to fit the two-particle correlation functions. The non-flow template is obtained from the analysis of LM pp and p–Pb data and is explained in detail in section Methods. This method reduces the non-flow contribution to the final v2 results for identified particles to 6% for pT  < 0.6 GeV/c and to less than 1% at higher pT . These estimates are based on applying the template-fit method to a pure non-flow model using the PYTHIA 8 event generator^44^.

Results

Figure 2 presents the pT -differential v2 measurement for mesons (π^±^, K^±^, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} ) and baryons (p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ) in semicentral Pb–Pb^45^, HM p–Pb and pp collisions. For Pb–Pb measurements, the two-particle correlation with ∣Δη∣ > 2.0 separation is used^45^, whereas for p–Pb and pp collisions, a ∣Δη∣ > 1.1 separation is applied. It was tested that varying the ∣Δη∣ separation does not significantly alter the results in Pb–Pb collisions^45^. Figure 2 shows a clear similarity in the characteristic features of v2 among the three collision systems. The difference in the magnitude of v2 among the three collision systems is consistent with previous measurements at the LHC^27^. For the low pT region, pT  < 2.0 GeV/c, a clear mass ordering of the v2 coefficients is observed, providing significant evidence of radial flow in small collision systems. The presence of radial flow in small collision systems is also supported by particle spectra measurements^46^. Around 2.0 < pT  < 3.0 GeV/c, the v2 coefficients of different particle species begin to cross. Beyond pT > 2.5 GeV/c, the v2 coefficients of baryons (p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ) are consistent with each other within 1 standard deviation ( ~ 1σ) up to 10 GeV/c in Pb–Pb and p–Pb collisions and up to 6 GeV/c in pp collisions. At the same time, the v2 of mesons (π^±^, K^±^, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} ) are compatible within  ~ 1σ at pT > 2(3) GeV/c for Pb–Pb and p–Pb (pp) collisions. Moreover, the v2 of baryons is larger than that of mesons by ~5σ at intermediate and higher pT (pT  > 3.0 GeV/c) in all three collision systems. In heavy-ion collisions, such distinctive baryon-meson v2 grouping at intermediate pT is explained by anisotropic flow development at the quark level, followed by particle production via the quark-coalescence mechanism^14,17,18^.Fig. 2pT-differential v2 measured with the two-particle correlation method^45^ for mesons (π^±^, K^±^, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} ) and baryons (p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ).Left: results for semicentral Pb–Pb collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{{s}_{{{{\rm{NN}}}}}}=5.02$$\end{document} TeV. Middle: results in high-multiplicity p–Pb collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{{s}_{{{{\rm{NN}}}}}}=5.02$$\end{document} TeV. Right: Same for pp collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13$$\end{document} TeV. 〈Nch〉 is the average number of reconstructed, efficiency-corrected charged particles with 0.2 < pT < 3.0 GeV/c at midrapidity (∣η∣ < 0.8). Horizontal bars (boxes) represent the statistical (systematic) uncertainties.

Existing measurements with identified particles in small collision systems presented either a single baryon or meson v2, limiting the opportunity to explore potential groupings among baryons and mesons^47–49^. Other similar measurements have not shown a clear grouping and splitting of baryon and meson v2 at intermediate pT in small systems, differing from similar measurements in heavy-ion collisions^25^. This difference may arise from difficulties in accounting for non-flow effects in small systems. For example, previous measurements either applied no non-flow removal technique^26,49^, or were based on the subtraction of the low-multiplicity correlation functions from high-multiplicity ones^22,25,47^, under the assumption that the former include only non-flow effects and no significant long-range correlations. However, recent ALICE measurements show that long-range correlations persist even in pp and p–Pb collisions with very low multiplicity (〈Nch〉 ≈ 10)^50^. This implies that the subtraction method also removes, along with non-flow effects, part of the real correlation signal, thus leading to over-subtraction. Such over-subtraction can vary by particle type and can potentially create a particle-type-dependent v2 pattern that does not originate from a true physics effect. As a result, subtraction-based methods are unreliable for studying the baryon–meson v2 splitting in small collision systems. The results presented in Fig. 2, after removal of non-flow effects (see Correlation function and template fit method in the methods subsection for details), show a distinctive baryon–meson v2 grouping (within 1σ) and a significant splitting (~5σ) at intermediate pT in both p–Pb and pp collisions at the LHC, similar to what is observed in heavy-ion collisions.

In refs. ^25,26^, the number-of-constituent-quark (NCQ) scaling of v2 has also been studied. This scaling was initially attributed to hadron production via the coalescence of thermal partons in heavy-ion collisions^19,20,51^. However, advanced coalescence models incorporate the recombination of thermal quarks with shower quarks originating from jet-medium interactions to describe the pT spectra and v2 of identified particles over a broad pT range^52^, differing from the coalescence mechanism^19,20,51^ associated with NCQ scaling of v2. In addition, contributions from radial flow and jet fragmentation at intermediate pT can also lead to deviations from NCQ scaling. Notably, the ALICE measurements exhibit deviations from NCQ scaling at the level of  ± 20% in Pb–Pb collisions^14,53^. This underscores the need for a better understanding of this scaling as evidence of partonic collectivity in relativistic collisions.

In Figs. 3 and 4, the v2 measurements are compared with the state-of-the-art calculations using the Hydro-Coal-Frag hybrid model^54,55^ for p–Pb and pp collisions, respectively. At low pT, this model incorporates the hydrodynamic evolution of a quark–gluon plasma with a partonic equation of state, followed by the formation of quarks before hadronization. At high pT, it accounts for interactions between high-energy partons and the medium using the linear Boltzmann transport (LBT) model, combined with hadronization via quark fragmentation. The intermediate pT hadrons are produced from the coalescence of quarks originating from both hydrodynamic evolution and jet-medium interactions. Finally, hadronic interactions occur after hadronization. A more detailed description of this model can be found in the methods subsection. This model provides a comprehensive explanation for both hadron production and anisotropic flow over a wide pT range in high-energy heavy-ion collisions^52^. It emphasizes the crucial role of partonic flow and particle production through quark coalescence in heavy-ion collisions, where the QGP is formed. In Fig. 3, the model parameters are tuned to describe the pT spectra of identified hadrons in high-multiplicity p–Pb collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{{s}_{{{{\rm{NN}}}}}}=5.02$$\end{document} TeV^54^. The figure demonstrates that the Hydro-Coal-Frag model successfully reproduces the baryon-meson v2 splitting and grouping features for pT < 8 GeV/c as observed in the experimental data. In contrast, the calculation from the Hydro-Frag model, which does not include the quark-coalescence process, strongly underestimates the v2 coefficients of all identified hadrons for pT > 4 GeV/c. Moreover, despite parameter adjustments, the Hydro-Frag model fails to even qualitatively reproduce the baryon-meson v2 splitting and grouping at intermediate pT ^55^.Fig. 3v2 in high-multiplicity p–Pb collisions.pT -differential v2 measured with two-particle correlation for mesons (π^±^, K^±^, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} ) and baryons (p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ) in high-multiplicity p–Pb collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{{s}_{{{{\rm{NN}}}}}}=5.02$$\end{document} TeV. Horizontal bars (boxes) represent the statistical (systematic) uncertainties. Comparisons with the calculations from the Hydro-Coal-Frag model (left) and the Hydro-Frag model (right) are also presented^54,55^. Only statistical uncertainties are shown for the calculations.

The comparison between the measurements and the model prediction for pp collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13$$\end{document} TeV is presented in Fig. 4. The model parameters are calibrated using pT spectra of identified hadrons from a different multiplicity interval than the one used in this paper. Still, the Hydro-Coal-Frag picture can explain the mass ordering of v2 for pT up to 3 GeV/c combined with the crossing between v2 of baryons and mesons at pT  ≈ 3 GeV/c, consistent with the data. Most importantly, the baryon-meson splitting and grouping of v2 can be qualitatively reproduced by the Hydro-Coal-Frag model up to approximately 6–7 GeV/c. In contrast, neither of these features is captured by the Hydro-Frag model. The crossing between v2 of different particles occurs at about 2 GeV/c, and only mass ordering is observed in this calculation^55^. The ordering is reversed for 2.5  < pT  < 5 GeV/c compared to pT < 2.5 GeV/c, similar to the Hydro-Frag calculations in p–Pb collisions shown in Fig. 3. Therefore, these results provide evidence of hadronization via coalescence of hydrodynamically flowing quarks in small collision systems at the LHC.Fig. 4v2 in high-multiplicity pp collisions.pT-differential v2 measured with two-particle correlation for mesons (π^±^, K^±^, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} ) and baryons (p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ) in high-multiplicity pp collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13$$\end{document} TeV. Horizontal bars (boxes) represent the statistical (systematic) uncertainties. Comparisons with the calculations from the Hydro-Coal-Frag model (left) and the Hydro-Frag model (right) are also presented^54,55^. Only statistical uncertainties are shown for the calculations.

Other theoretical calculations that predict a flow pattern in high multiplicity events of small collision systems have also been studied. Since the presented measurements are based on long-range (i.e., large ∣Δη∣ separation) two-particle correlations, they are predominantly influenced by geometry-driven effects, as the initial momentum anisotropy in the CGC approach generates only short-range correlations^37^. The hadronic rescatterings in UrQMD can mimic the long-range two-particle correlation and mass dependence of v2 at low pT ^36^, but do not generate any baryon-meson v2 splitting and grouping. The PYTHIA 8 model, with color ropes enabled in the hadronization process, can describe the strangeness enhancement in pp collisions^28^ without the formation of a QGP^56^, but cannot generate long-range correlations. Interestingly, the string-string repulsion in the string-shoving version of PYTHIA 8 can generate long-range two-particle correlations^38,57^, but it produces negative flow coefficients after the template fit, in contrast to the positive flow coefficients observed in the data. In small systems, transport models like AMPT^35^, which generate only a few partonic interactions during system evolution and incorporate the quark-coalescence model of hadronization, can approximately describe the v2 at low pT. However, they fail to even qualitatively explain the baryon-meson v2 grouping and splitting feature observed at intermediate pT as shown in Fig. 5^58^. This indicates that anisotropic flow is developed in a dense partonic system and propagated to the level of hadrons via the quark coalescence process.Fig. 5v2 in high-multiplicity p–Pb collisions compared with AMPT calculations.pT -differential v2 measured with two-particle correlation for mesons (π^±^, K^±^, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} ) and baryons (p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ) in high-multiplicity p–Pb collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{{s}_{{{{\rm{NN}}}}}}=5.02$$\end{document} TeV. Horizontal bars (boxes) represent the statistical (systematic) uncertainties. Comparison with the calculations from the AMPT String-melting model^58^ is also presented. The AMPT curves are obtained by applying the same template fit method to the correlation distributions as used in the data analysis^58^. Only statistical uncertainties are shown for the AMPT calculations.

In summary, the v2 for the identified hadrons in high-multiplicity p–Pb collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{}}}}\sqrt{{s}_{{{{\rm{NN}}}}}}{{{\rm{}}}}=5.02$$\end{document} TeV and pp collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13$$\end{document} TeV has been presented as a function of pT and compared with measurements in semicentral Pb–Pb collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{}}}}\sqrt{{s}_{{{{\rm{NN}}}}}}{{{\rm{}}}}=5.02$$\end{document} TeV. A characteristic grouping (within  ~ 1σ) and splitting (with  ~ 5σ) of v2 for mesons (π^±^, K^±^, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} ) and baryons (p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ) at intermediate pT, similar to measurements in heavy-ion collisions, is observed in both p–Pb and pp collisions. The Hydro-Coal-Frag model, incorporating partonic flow and quark coalescence, provides the best possible description of the data to date in both heavy-ion and small collision systems. Alternative approaches fail to even qualitatively reproduce the baryon-meson v2 splitting and grouping at intermediate pT. The presented measurements and the comparisons with available theoretical model calculations provide evidence that the system created in high-multiplicity p–Pb and pp collisions includes a stage with hydrodynamically flowing partons, similar to the one observed in heavy-ion collisions.

Methods

Event selection

The analyzed data samples are from pp collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13$$\end{document} TeV and p–Pb collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{}}}}\sqrt{{s}_{{{{\rm{NN}}}}}}{{{\rm{}}}}=$$\end{document} 5.02 TeV, collected by the ALICE detector during the LHC Run 2 data-taking campaign between 2016 and 2018. An extensive description of all subdetectors of ALICE can be found in refs. ^39,40^. The collected data are classified based on the specific triggering conditions. Minimum bias (MB) events for both pp and p–Pb collisions are triggered using a coincidence signal in the two scintillator arrays of the V0 detector, which cover the pseudorapidity ranges 2.8 < η < 5.1 (V0A) and  − 3.7 < η < − 1.7 (V0C), respectively. To avoid the possibility of overlap between HM and LM event classes, additional requirements on the number of reconstructed and efficiency-corrected charged particles (Nch) within the acceptance of ∣η∣ < 0.8 and transverse momentum 0.2 < pT  < 3.0 GeV/c are introduced. For p–Pb collisions, the collisions of the 0–20% and 60–100% V0A-centrality are used as the HM and LM events, with an additional criterion of Nch < 20 applied to the LM sample. In pp collisions, HM events are selected from the 0.07% of events with the highest multiplicity. This selection uses a special HM trigger based on the amplitude of the V0M detector arrays (V0A + V0C) and requires Nch > 25. For the LM event class, MB events with Nch < 20 are selected. Only events with a reconstructed primary vertex (PV) within  ± 10 cm from the nominal interaction point along the beam line are selected. Background events due to interaction between the beam and the residual gas molecules in the beam pipe are removed using the information from the Silicon Pixel Detector (SPD) and V0 detectors. The out-of-bunch pileup events are rejected using a correlation between the multiplicities in the V0 and Forward-Multiplicity-Detector (FMD)^39^. The in-bunch pileup is reduced by rejecting events with multiple vertices. The event selections result in a sample of 226 × 10^6^ HM and 572 × 10^6^ LM pp collisions, corresponding to integrated luminosities of approximately 4 nb^−1^ and 10 nb^−1^, respectively^59^. For p–Pb, the analyzed data sample consists of 101 × 10^6^ HM collisions and 198 × 10^6^ LM collisions, corresponding to integrated luminosities of about 0.05 nb^−1^ and 0.09 nb^−1^, respectively^60^. In this analysis, the HM events are used for the flow measurements, whereas the LM events are used for the baseline non-flow estimation.

Track reconstruction

The charged particles at midrapidity are reconstructed using the Inner Tracking System (ITS) and the Time Projection Chamber (TPC). The selected tracks have at least 70 TPC space points (out of a maximum of 159) fitted for track reconstruction with χ^2^ per degree of freedom lower than four. The reconstructed tracks from the TPC and the ITS must coincide to ensure the consistency of the reconstruction method. Additionally, a minimum of two hits are required in the ITS to improve the momentum resolution. The pseudorapidity of the selected tracks is required to be within ∣η∣ < 0.8 to reject tracks with reduced reconstruction efficiency at the detector edges. To reduce the contamination from secondary particles, the distance-of-closest approach (DCA) of the selected tracks to the PV must be within 2 cm in the longitudinal direction. Furthermore, a pT-dependent DCA selection in the transverse plane, ranging from 0.2 cm at pT = 0.2 GeV/c to 0.02 cm at pT = 5.0 GeV/c, is applied. This criterion suppresses the residual contamination from secondary particles from weak decays and interactions in the detector material. The reference particles used for the construction of long-range di-hadron correlation functions are selected from the FMD detector segments located at forward (FMD1,2 with 1.7 < η < 5.1) and backward (FMD3 with  − 3.1 < η < − 1.7) rapidity regions.

Particle identification and reconstruction

The charged tracks are identified as π^±^, K^±^, and p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} based on the specific energy loss (dE/dx) information in the TPC and the velocity information from the Time-of-Flight (TOF) detector. A Bayesian approach^45^ is used to identify particle species at a given pT  using the correlation of the normalized differences between the measured and the expected signal in the TPC (n**σTPC) and the TOF (n**σTOF), respectively. In this method, the signals are converted into probabilities and folded with the expected abundances (priors) of each particle species. To ensure high purity of the selected sample, a minimal probability threshold of 0.95 for π^±^ and 0.85 for K^±^ and p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} is set. In addition, the tracks with proper TOF information are required to be within ∣n**σTPC∣ < 3 and ∣n**σTOF∣ < 3. The resulting purity, estimated using Monte Carlo (MC) simulations, is higher than 95% for π^±^ for 0.2 < pT  < 10 GeV/c, above 80% for K^±^ for 0.3 < pT  < 10 GeV/c, and reaches values larger than 90% for p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} for 0.5 < pT  < 10 GeV/c. The high purity of the studied sample reduces the uncertainties due to particle misidentification. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} and Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} are weakly decaying neutral particles, reconstructed by calculating the invariant mass of the daughter particles from the most probable decay channels of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}\,\to {\pi }^{+}\,+{\pi }^{-}$$\end{document} and Λ  → p + π^−^  ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }\,\to \overline{{{{\rm{p}}}}}\,+{\pi }^{+}$$\end{document} ) with branching ratios of 69.2% (±0.05%) and 64.1% (±0.5%)^61^, respectively. The combinatorial background is suppressed by using a set of selection criteria on the decay topology used in the previous \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} and Λ measurements in ALICE^28^. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} and Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} candidates are selected within the rapidity range ∣y∣ < 0.5 inside the TPC, and the daughter tracks are used to reconstruct the secondary decay vertex (SV) in the offline reconstruction. The SV is required to be more than 0.5 cm away from the PV, and the reconstructed proper lifetime, defined as m**L/p, (m being the particle mass, L the distance between the primary and secondary vertices, and p the particle momentum) should be smaller than 20 cm and 30 cm for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} and Λ ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ) candidates, respectively. The oppositely charged daughter tracks are combined only if they are identified as pions or protons based on the TPC dE/dx hypothesis (<3σ). A set of topological cuts, such as the distance of closest approach (DCA) of the daughter tracks to the primary vertex (> 0.06 cm), DCA between the daughter tracks (<1 cm), and cosine of the pointing angle, which is the angle between the momentum direction of the mother particle and the direction from the PV to the decay point (>0.97 for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} and >0.995 for Λ ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} )) are applied to reduce the combinatorial background contribution to the invariant mass spectrum.

Correlation function and template fit method

The correlation function is obtained between two sets of particles classified as trigger and associated. Trigger particles are used as a reference, and the angular distribution of associated particles is measured relative to the trigger particles^22^. In this analysis, the 2D correlation function is constructed as a function of the difference in azimuthal angle Δφ = φtrigger -φassociated and pseudorapidity Δη = ηtrigger -ηassociated with trigger and associated particles from different detectors. Three sets of correlation functions are constructed to estimate the v2 of identified particles (π^±^, K^±^, p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} , and Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ). The identified particles in the TPC are correlated with unidentified reference particles reconstructed within the FMD acceptance in positive (FMD1,2) and negative (FMD3) rapidity regions to construct two sets (TPC–FMD1,2 and TPC–FMD3) of correlation functions. The third correlation function is constructed using two reference particles from the FMD1,2 and FMD3 detector segments. The pair acceptance effect due to the finite size of the detectors is corrected by dividing the same-event correlation functions with mixed-event correlation functions. The mixed-event correlation function is constructed by correlating the trigger particles in one event with the associated particles from other events belonging to the same multiplicity event class and with PV within a given 2 cm wide interval. The mixed event is normalized using a constant estimated by averaging over all Δφ bins at the Δη value where the mixed event correlation function reaches its maximum. The corrected correlation function is obtained as a ratio of the same and mixed event correlation functions for each PV position. The final correlation function for an event class (HM or LM events) is calculated after averaging the correlation functions over all PV positions. For each of the 2D correlation functions (TPC–FMD1,2, TPC–FMD3, and FMD1,2–FMD3), the projection along the Δφ axis is calculated for both HM and LM cases. The Δφ projections from LM collisions serve as a template for subsequent fitting of the Δφ projections from HM collisions to extract the v2 coefficients using the template fit method. The template fit method assumes that high-multiplicity (HM) collisions are a superposition of low-multiplicity (LM) collisions, which primarily contain non-flow effects with some residual flow, along with an additional flow modulation, i.e.,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\rm{Y}}}}{(\Delta \varphi )}^{{{{\rm{HM}}}}}=F{{{\rm{Y}}}}{(\Delta \varphi )}^{{{{\rm{LM}}}}}+G\left[1+\sum\limits_{n=2}^{\infty }2{{{{\rm{V}}}}}_{n\Delta }\cos (n\Delta \varphi )\right],$$\end{document}

where Y(Δφ)^HM^ and Y(Δφ)^LM^ are the one dimensional Δφ projections of the 2D correlation functions obtained in HM and LM collisions with F and G being the scaling factors. The V_nΔ_ coefficients are estimated by fitting the correlation function with the equation (2). The scaling factors F and G are free parameters in this template fit procedure. The final v2 of the identified particles in the TPC is calculated by combining the V_2Δ_ estimated from the TPC–FMD1,2, TPC–FMD3, and FMD1,2–FMD3 correlation functions

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${v}_{2}^{{{{\rm{PID}}}}}({p}_{{{{\rm{T}}}}})=\sqrt{\frac{{{{{\rm{V}}}}}_{2\Delta }^{{{{\rm{TPC-FMD1,2}}}}}\,{{{{\rm{V}}}}}_{2\Delta }^{{{{\rm{TPC-FMD3}}}}}}{{{{{\rm{V}}}}}_{2\Delta }^{{{{\rm{FMD1,2-FMD3}}}}}}}.$$\end{document}

In this work, all available non-flow suppression methods (low-multiplicity subtraction^22,25,47^, template fit^23^, and improved template fit^62^) have been tested, and the residual non flow has been estimated using PYTHIA8 for each method. Among these, the template fit provides the most effective non flow subtraction, yielding the lowest residual non flow (~5-7%) across the considered kinematic range. This residual non flow has been included in the systematic uncertainties. The inclusion of the remaining non flow enables better comparisons with theoretical models and supports robust, data driven physics conclusions.

Systematic uncertainty

The systematic uncertainties are evaluated by varying the event, track, and PID selection criteria with respect to the default ones, one at a time. For each variation, the difference between the default and varied result is estimated using the Barlow criterion^63^, and a difference higher than 1σ is considered as a possible source of systematic uncertainty in the measurement. The Barlow difference is calculated for each particle species and for each pT interval for which the final v2{2PC} results are presented in this paper. The Barlow ratio is calculated as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=\frac{| {v}_{2}^{{{{\rm{default}}}}}-{v}_{2}^{{{{\rm{syst}}}}}| }{\sqrt{| {\sigma }_{{{{\rm{default}}}}}^{2}-{\sigma }_{{{{\rm{syst}}}}}^{2}| }}.$$\end{document}

If the Barlow difference is higher than 1σ for more than 1/3 of the total pT intervals for any species, the contribution of that particular systematic source is included in the uncertainty of the final result. Otherwise, the contribution from that systematic source is considered negligible and does not contribute to the final systematic uncertainty. The minimum and maximum values of the relative systematic uncertainties from individual sources are presented in Tables 1 and 2 for p–Pb and pp collisions, respectively. The systematic sources listed in the tables from top to bottom include different PV intervals used for event selection, correlation between multiplicity from V0 and FMD detectors to reduce contamination in the FMD, track selection criteria, and particle identification criteria affecting the purity of the π^±^, K^±^ and p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} samples. Other factors are topological reconstruction criteria, invariant mass reconstruction and fitting requirements impacting the signal-to-background ratios for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} and Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} candidates, the definition of the low-multiplicity template used for non-flow removal, and the estimation of residual secondary contamination in the FMD. The latter is done using a Monte Carlo event generator, by transporting the generated particles through GEANT3-simulated detector response and performing track reconstruction in the ALICE framework. The contributions from the different sources are added in quadrature to estimate the total systematic uncertainty.Table 1. Systematic uncertainties in p–Pb collisionsv2{2PC, 1.1 < ∣Δη∣ < 7.8}Uncertainty sourceπ^±^K^±^p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} Primary vertex position1–2%1–3%1–2%negl.negl.FMD–V0 correlation0–4%0–2%0–2%1–4%2–4%Primary track qualitynegl.negl.negl.––Bayesian thresholdnegl.negl.negl.––Topological criteria–––1–2%1–4%Invariant mass acceptance–––negl.3–4%Invariant mass fit–––negl.1–3%Template variationnegl.negl.1–2%negl.negl.Secondary contamination1%1%1%1%1%Total2–4%2–3%1–3%2–5%5–7%The minimum and maximum values of the relative systematic uncertainties from individual sources for π^±^, K^±^, p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} , and Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} in p–Pb collisions. Percentage ranges are given to account for variations with pT. The fields marked as “negl." (negligible) denotes that the uncertainties have been tested but are not statistically significant.Table 2. Systematic uncertainties in pp collisionsv2{2PC, 1.1 < ∣Δη∣ < 6.4}Uncertainty sourceπ^±^K^±^p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} Primary vertex position0–1%negl.negl.2–3%1–3%FMD–V0 correlation1–2%negl.0–2%1–4%1–3%Primary track quality0–2%negl.negl.––Bayesian thresholdnegl.0–2%negl.––Topological criteria–––2–4%1–5%Invariant mass acceptance–––negl.1–3%Invariant mass fit–––1–2%3–4%Template variation0–4%0–2%1–3%1–3%2–4%Secondary contamination1%1%1%1%1%Total2–5%2%2–4%4–6%4–8%The minimum and maximum values of the relative systematic uncertainties from individual sources for π^±^, K^±^, p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\rm{K}}}}}_{{{{\rm{S}}}}}^{0}$$\end{document} , and Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} in pp collisions. Percentage ranges are given to account for variations with pT. The fields marked as “negl." (negligible) denotes that the uncertainties have been tested but are not statistically significant.

Hydro+Coal+Frag model description

The Hydro-Coal-Frag model^52,54,55^ provides a unified theoretical framework for hadron production in high-energy nuclear collisions, bridging soft and hard processes across transverse momentum (pT) regimes. It describes low-pT hadrons via viscous hydrodynamics, intermediate-pT hadrons through quark coalescence, and high-pT hadrons via string fragmentation. The modeling sequence begins with the TRENTo model^64^, which generates event-by-event initial entropy profiles based on the nuclear geometry. These profiles serve as initial conditions for the (2+1)-dimensional viscous hydrodynamics model VISH2+1^65^, which governs the space-time evolution of the quark-gluon plasma (QGP). As the system cools toward the hydrodynamic freeze-out temperature, thermal hadrons are emitted according to the Cooper-Frye prescription^65^, and thermal partons are sampled at low transverse momentum (pT) for subsequent hadronization. High-pT partons (hard partons), generated using PYTHIA8, traverse the quark-gluon plasma and undergo medium-induced interactions, which are modeled using the Linear Boltzmann Transport (LBT) framework^66^.

For intermediate-pT hadrons, the quark coalescence mechanism is used to recombine thermal–thermal, thermal–hard, and hard–hard partons produced by the hydrodynamics and LBT processes. Meson and baryon momentum distributions are derived from Wigner functions, which encode the spatial and momentum proximity of coalescing partons. Excited hadronic states, formed according to the invariant masses of parton pairs, subsequently decay into ground states, with binding energy differences and conservation of energy and momentum explicitly taken into account. Remaining hard partons without coalescence partners generate high-pT hadrons via string fragmentation. The transverse momentum (pT) cut-off values for thermal parton sampling at freeze-out, along with the criteria governing whether partons undergo coalescence or fragmentation following the LBT stage, and the gluon virtuality parameters, are tuned to reproduce the pT spectra of pions, kaons, and protons, as well as the (p( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} )/π^±^) ratio in the intermediate-pT region of high-multiplicity p–p and p–Pb collisions at the LHC. The final hadronic evolution, including scatterings and resonance decays, is simulated using the Ultrarelativistic Quantum Molecular Dynamics (UrQMD) model.

AMPT calculations

Figure 5 presents the pT-differential v2 measured from two-particle correlations for mesons (π^±^, K^±^) and baryons (p+ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{{{\rm{p}}}}}$$\end{document} , Λ + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\Lambda }$$\end{document} ) in HM p–Pb collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{{s}_{{{{\rm{NN}}}}}}=5.02$$\end{document} TeV, compared with estimations from the AMPT string melting model^58^. The AMPT curves are obtained by applying the same template fit method to the correlation distributions as used in the analysis of the data. Both the data and AMPT calculations^58^ select particles within similar rapidity regions, allowing a pseudorapidity separation between the two correlated particles of 1.1 < ∣Δη∣ < 7.8 in the data and 1.7 < ∣Δη∣ < 8.0 in the AMPT simulations.

Supplementary information

Transparent Peer Review file