Multiplicity Dependence of Pion, Kaon, Proton and Lambda Production in p-Pb Collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV

In this Letter, comprehensive results on ${\rm\pi}^\pm$, K$^\pm$, K$^0_S$, p, $\rm\bar{p}$, $\rm \Lambda$ and $\rm \bar{\Lambda}$ production at mid-rapidity ($0<y_{\rm cms}<0.5$) in p-Pb collisions at $\sqrt{s_{\rm NN}} = 5.02$ TeV, measured by the ALICE detector at the LHC, are reported. The transverse momentum distributions exhibit a hardening as a function of event multiplicity, which is stronger for heavier particles. This behavior is similar to what has been observed in pp and Pb-Pb collisions at the LHC. The measured $p_{\rm T}$ distributions are compared to results at lower energy and with predictions based on QCD-inspired and hydrodynamic models.


Introduction
High-energy heavy-ion (AA) collisions offer a unique possibility to study nuclear matter under extreme conditions, in particular the deconfined quark-gluon plasma which has been predicted by quantum chromodynamics (QCD) [1,2,3,4].The interpretation of heavy-ion results depends crucially on the comparison with results from smaller collision systems such as proton-proton (pp) or proton-nucleus (pA).
The bulk matter created in high-energy nuclear reactions can be quantitatively described in terms of hydrodynamic and statistical models.The initial hot and dense partonic matter rapidly expands and cools down, ultimately undergoing a transition to a hadron gas phase [5].The observed ratios of particle abundances can be described in terms of statistical models [6,7], which are governed mainly by two parameters, the chemical freeze-out temperature T ch and the baryochemical potential µ B which describes the net baryon content of the system.These models provide an accurate description of the data over a large range of center-of-mass energies (see e.g.[8]), but a surprisingly large deviation (about 50%) was found for the proton production yield at the LHC [9,10].During the expansion phase, collective hydrodynamic flow develops from the initially generated pressure gradients in the strongly interacting system.This results in a characteristic dependence of the shape of the transverse momentum (p T ) distribution on the particle mass, which can be described with a common kinetic freeze-out temperature parameter T kin and a collective average expansion velocity β T [11].
Proton-nucleus (pA) collisions are intermediate between proton-proton (pp) and nucleus-nucleus (AA) collisions in terms of system size and number of produced particles.Comparing particle production in pp, pA, and AA reactions has frequently been used to separate initial state effects, linked to the use of nuclear beams or targets, from final state effects, linked to the presence of hot and dense matter.At the LHC, however, the pseudorapidity density of final state particles in pA collisions reaches values which can become comparable to semi-peripheral Au-Au (∼60% most central) and Cu-Cu (∼30% most central) collisions at top RHIC energy [12].Therefore the assumption that final state dense matter effects can be neglected in pA may no longer be valid.In addition, pA collisions allow for the investigation of fundamental properties of QCD: the relevant part of the initial state nuclear wave function extends to very low fractional parton momentum x and very high gluon densities, where parton shadowing and novel phenomena like saturation, e.g. as implemented in the Color Glass Condensate model (CGC), may become apparent [13,14].
Recently, measurements at the LHC in high multiplicity pp and p-Pb collisions have revealed a near-side long-range "ridge" structure in the two-particle correlations [15,16].The observation of an unexpected "double-ridge" structure in the two-particle correlations in high-multiplicity p-Pb collisions has also been reported [17,18,19,20].This is flat and long-range in pseudo-rapidity ∆η and modulated in azimuth approximately like cos(2∆φ ), where ∆η and ∆φ are the differences in pseudo-rapidity η and azimuthal angle φ between the two particles.Various mechanisms have been proposed to explain the origin of this double-ridge like structure.Both a CGC description [21], based on initial state nonlinear gluon interactions, as well as a model based on hydrodynamic flow [22,23], assuming strong interactions between final state partons or hadrons, can give a satisfactory description of the p-Pb correlation data.However, the modeling of small systems such as p-Pb is complicated because uncertainties related to initial state geometrical fluctuations play a large role and because viscous corrections may be too large for hydrodynamics to be a reliable framework [24].Additional experimental information is therefore required to reveal the origin of these correlations.The p T distributions and yields of particles of different mass at low and intermediate momenta of less than a few GeV/c (where the vast majority of particles is produced), can provide important information about the system created in high-energy hadron reactions.Previous results on identified particle production in pp [25,26,27,28,29] and Pb-Pb [9,10] collisions at the LHC have been reported.In this paper we report on the measurement of π ± , K ± , K 0 S , p(p) and Λ( Λ) production as a function of the event multiplicity in p-Pb collisions at a nucleon-nucleon center-of-mass energy √ s NN = 5.02 TeV.The results are presented over the following p T ranges: 0.1-3, 0.2-2.5, 0-8, 0.3-4 and 0.6-8 GeV/c for π ± , K ± , K 0 S , p(p) and Λ( Λ), respectively.Results on π, K, p production in p-Pb collisions have been recently reported by the CMS collaboration [30].

Sample and Data analysis
The results presented in this letter are obtained from a sample of the data collected during the LHC p-Pb run at √ s NN = 5.02 TeV in the beginning of 2013.Because of the 2-in-1 magnet design of the LHC [31], the energy of the two beams cannot be adjusted independently and is 4 ZTeV, leading to different energies due to the different Z/A.The nucleon-nucleon center-of-mass system, therefore, was moving in the laboratory frame with a rapidity of y NN = −0.465 in the direction of the proton beam.The number of colliding bunches was varied from 8 to 288.The total number of protons and Pb ions in the beams ranged from 0.2×10 12 to 6.5×10 12 and from 0.1×10 12 to 4.4×10 12 , respectively.The maximum luminosity at the ALICE interaction point was for the data used in this paper 5 × 10 27 cm −2 s −1 resulting in a hadronic interaction rate of 10 kHz.The interaction region had an r.m.s. of 6.3 cm along the beam direction and of about 60 µm in the direction transverse to the beam.For the results presented in this letter, a low-luminosity data sample has been analyzed where the event pile-up rate has been estimated to have negligible effects on the results.The integrated luminosity corresponding to the used data sample was about 14 µb −1 (7 µb −1 ) for the neutral (charged) hadron analysis.The LHC configuration was such that the lead beam circulated in the "counter-clockwise" direction, corresponding to the ALICE A direction or positive rapidity as per the convention used in this paper.
A detailed description of the ALICE apparatus can be found in [32].The minimum-bias trigger signal was provided by the VZERO counters, two arrays of 32 scintillator tiles each covering the full azimuth within 2.8 < η lab < 5.1 (VZERO-A, Pb beam direction) and −3.7 < η lab < −1.7 (VZERO-C, p beam direction).The signal amplitude and arrival time collected in each tile were recorded.A coincidence of signals in both VZERO-A and VZERO-C detectors was required to remove contamination from single diffractive and electromagnetic events [33].The time resolution is better than 1 ns, allowing discrimination of beam-beam collisions from background events produced outside of the interaction region.In the offline analysis, background was further suppressed by the time information recorded in two neutron Zero Degree Calorimeters (ZDCs), which are located at +112.5 m (ZNA) and −112.5 m (ZNC) from the interaction point.A dedicated quartz radiator Cherenkov detector (T0) provided a measurement of the event time of the collision.
The ALICE central-barrel tracking detectors cover the full azimuth within |η lab | < 0.9.They are located inside a solenoidal magnet providing a magnetic field of 0. the Inner Tracking System (ITS).It consists of six layers of silicon devices grouped in three individual detector systems which employ different technologies (from the innermost outwards): the Silicon Pixel Detector (SPD), the Silicon Drift Detector (SDD) and the Silicon Strip Detector (SSD).The Time Projection Chamber (TPC), the main central-barrel tracking device, follows outwards.Finally the Transition Radiation Detector (TRD) extends the tracking farther away from the beam axis.The primary vertex position was determined separately in the SPD [33] and from tracks reconstructed in the whole central barrel (global tracks).The events were further selected by requiring that the longitudinal position of the primary vertex was within 10 cm of the nominal interaction point and that the vertices reconstructed from SPD tracklets and from global tracks are compatible.In total from a sample of 29.8 (15.3) million triggered events about 24.7 (12.5) million events passing the selection criteria were used in the neutral (charged) hadron analysis.
In order to study the multiplicity dependence, the selected event sample was divided into seven event classes, based on cuts on the total charge deposited in the VZERO-A detector (V0A).The corresponding fractions of the data sample in each class are summarized in Tab. 1.The mean charged-particle multiplicity densities ( dN ch /dη ) within |η lab | < 0.5 corresponding to the different centrality bins are also listed in the table.These are obtained using the method presented in [33] and are corrected for acceptance and tracking efficiency as well as for contamination by secondary particles.The relative standard deviation of the track multiplicity distribution for the event classes defined in Table 1 ranges from 78% to 29% for the 80-100% and 0-5% classes, respectively.It should be noted that the average multiplicity in the 80-100% bin is well below the corresponding multiplicity in pp minimum-bias collisions [34] and therefore likely to be subject to a strong selection bias.Contrary to our earlier measurement of dN ch /dη [33], the values in Tab. 1 are not corrected for trigger and vertex-reconstruction efficiency, which is of the order of 2% for NSD events [33].The same holds true for the p T distributions, which are presented in the next section.
Charged-hadron identification in the central barrel was performed with the ITS, TPC [35] and Time-Of-Flight (TOF) [36] detectors.The drift and strip layers of the ITS provide a measurement of the specific energy loss with a resolution of about 10%.In a standalone tracking mode, the identification of pions, kaons, and protons is thus extended down to respectively 0.1, 0.2, 0.3 GeV/c in p T .The TPC provides particle identification at low momenta via specific energy loss dE/dx in the fill gas by measuring up to 159 samples per track with a resolution of about 6%.The separation power achieved in p-Pb collisions is identical to that in pp collisions [37].Further outwards at about 3.7 m from the beam line, the TOF  array allows identification at higher p T measuring the particle speed with the time-of-flight technique.The total time resolution is about 85 ps for events in the multiplicity classes from 0% to ∼ 80%.In more peripheral collisions, where multiplicities are similar to pp, it decreases to about 120 ps due to a worse start-time (collision-time) resolution [37].The start-time of the event was determined by combining the time estimated using the particle arrival times at the TOF and the time measured by the T0 detector [36].
Since the p-Pb center-of-mass system moved in the laboratory frame with a rapidity of y NN = −0.465, the nominal acceptance of the central barrel of the ALICE detector was asymmetric with respect to y CMS = 0.In order to ensure good detector acceptance and optimal particle identification performance, tracks were selected in the rapidity interval 0 < y CMS < 0.5 in the nucleon-nucleon center-of-mass system.Event generator studies and repeating the analysis in |y CMS | < 0.2 indicate differences between the two rapidity selections smaller than 2% in the normalization and 3% in the shape of the transverse momentum distributions.
In this paper we present results for primary particles, defined as all particles produced in the collision,  including decay products, but excluding weak decays of strange particles.The analysis technique is described in detail in [9,10,38].Here we briefly review the most relevant points.
Three approaches were used for the identification of π ± , K ± , and p(p), called "ITS standalone", "TPC/TOF" and "TOF fits" [9,10] in the following.In the "ITS standalone" method, a probability for each particle species is calculated in each layer based on the measured energy loss signal and the known response function.The information from all layers is combined in a bayesian approach with iteratively determined priors.Finally, the type with the highest probability is assigned to the track.This method is used in the p T ranges 0.1 < p T < 0.7 GeV/c, 0.2 < p T < 0.6 GeV/c and 0.3 < p T < 0.65 GeV/c for π ± , K ± , and p(p), respectively.In contrast to the analysis in the high multiplicity environment of central heavy-ion collisions, the contribution of tracks with wrongly associated clusters is negligible in p-Pb collisions.In the "TPC/TOF" method, the particle is identified by requiring that its measured dE/dx and time-of-flight are within ±3σ from the expected values in the TPC and/or TOF.This method is used in the p T ranges 0.2 < p T < 1.5 GeV/c, 0.3 < p T < 1.3 GeV/c and 0.5 < p T < 2.0 GeV/c for π ± , K ± , and p(p), respectively.In the third method the TOF time distribution is fitted to extract the yields, with the expected shapes based on the knowledge of the TOF response function for different particle species.This method is used in the p T range starting from 0.5 GeV/c up to 3, 2.5 and 4 GeV/c for π ± , K ± , and p(p), respectively.Contamination from secondary particles was subtracted with a data-driven approach, based on the fit of the transverse distance-of-closest approach to the primary vertex (DCA xy ) distribution with the expected shapes for primary and secondary particles [9,10].The results of the three analyses were combined using the (largely independent) systematic uncertainties as weights in the overlapping ranges, after checking for their compatibility.
The K 0 S and Λ( Λ) particles were identified exploiting their "V 0 " weak decay topology in the channels K 0 S → π + π − and Λ( Λ) → pπ − (pπ + ), which have branching ratios of 69.2% and 63.9%, respectively [39].The selection criteria used to define two tracks as V 0 decay candidates are listed in Tab. 2 (see [26] for details).Since the cosine of pointing angle (the angle between the particle momentum associated with the V0 candidate and a vector connecting the primary vertex and the V0 position [26]) resolution changes significantly with momentum, the value used in the selection is p T dependent and such that no more than 1% of the primary particle signal is removed.

The ALICE Collaboration
GeV/c), increasing to about 70% for K 0 S and 55% for Λ( Λ) at higher momenta (p T > 3 GeV/c).The signal is extracted from the reconstructed invariant mass distribution subtracting the background from the peak region with a bin counting method.The background and signal regions are defined on the basis of the mass resolution as the windows in [−12σ , −6σ ], [6σ , 12σ ] and [−6σ , 6σ ], respectively.The value of σ changes with p T to account for the actual mass resolution and ranges from about 3 MeV/c 2 to 7 MeV/c 2 for K 0 S and from about 1.4 MeV/c 2 to 2.5 MeV/c 2 for Λ( Λ).More details on V 0 reconstruction can be found in [26,38].The contribution from weak decays of the charged and neutral Ξ to the Λ( Λ) yield has been corrected following a data-driven approach.The measured Ξ − ( Ξ+ ) spectrum is used as input in a simulation of the decay kinematics to evaluate the fraction of reconstructed Λ( Λ) coming from Ξ − ( Ξ+ ) decays.The contribution from the decays of Ξ 0 is taken into account in the same way by assuming the ratio Ξ − ( Ξ+ )/Ξ 0 = 1, as supported by statistical models and Pythia or DMPJET Monte Carlo simulations [40,41].The raw transverse momentum distributions have been corrected for acceptance and reconstruction efficiency using a Monte Carlo simulation, based on the DPMJET 3.05 event generator [40] and a GEANT3.21[42] model of the detector.As compared to the version used in [9,10], GEANT3.21 was improved by implementing a more realistic parameterization of the antiproton inelastic cross-section [43].A correction factor based on FLUKA [44] estimates was applied to negative kaons as in [9,10].
The study of systematic uncertainties follows the analysis described in [9,10] for π ± , K ± and p(p).The main sources are the correction for secondary particles (4% for protons, 1% for pions, negligible for kaons), knowledge of the material budget (3% related to energy loss), hadronic interactions with the detector material (from 1% to 6%, more important at low p T and for protons), tracking efficiency (4%), TOF matching efficiency (from 3 to 6%, depending on the particle) and PID (from 2% to 25%, depending on the particle and the p T range).For the neutral Λ and K 0 S particles, the main sources are the level of knowledge of detector materials (resulting in a 4% uncertainty), track selections (up to 5%) and the feed-down correction for the Λ and Λ (5%), while topological selections contribute 2-4% depending on transverse momentum.The main sources of systematic uncertainties for the analysis of charged and neutral particles are summarized in Tables 3 and 4, respectively.The study of systematic uncertainties was repeated for the different multiplicity bins in order to separate the sources of uncertainty which are dependent on multiplicity and uncorrelated across different bins (depicted as shaded boxes in the figures).

Results
The p T distributions of π ± , K ± , K 0 S , p(p) and Λ( Λ) in 0 < y CMS < 0.5 are shown in Fig. 1 for different multiplicity intervals, as defined in Tab. 1. Particle/antiparticle as well as charged/neutral kaon transverse momentum distributions are identical within systematic uncertainties.
The p T distributions show a clear evolution, becoming harder as the multiplicity increases.The change is most pronounced for protons and lambdas.They show an increase of the slope at low p T , similar to the one observed in heavy-ion collisions [9,10].The stronger multiplicity dependence of the spectral shapes of heavier particles is evident when looking at the ratios K/π = (K + + K − )/(π + + π − ), p/π = (p + p)/(π + + π − ) and Λ/K 0 S as functions of p T , shown in Fig. 2 for the 0-5% and 60-80% event classes.The ratios p/π and Λ/K 0 S show a significant enhancement at intermediate p T ∼ 3 GeV/c, qualitatively reminiscent of that measured in Pb-Pb collisions [9,10,38].The latter are generally discussed in terms of collective flow or quark recombination [45,46,47].However, the magnitude of the observed effects differs significantly between p-Pb and in Pb-Pb.The maximum of the p/π (Λ/K 0 S ) ratio reaches ∼ 0.8 (1.5) in central Pb-Pb collisions, but only 0.4 (0.8) in the highest multiplicity p-Pb events.The highest multiplicity bin in p-Pb collisions exhibits ratios of p/π and Λ/K 0 S which have maxima close to the corresponding ratios in the 60-70% bin in Pb-Pb collisions but differ somewhat in shape at lower p T .The value of dN ch /dη in central p-Pb collisions (45 ± 1) is a factor ∼ 1.7 lower than the one in the 60-70% Pb-Pb bin.A similar enhancement of the p/π ratio in high-multiplicity d-Au collisions has also    been reported for RHIC energies [48].
It is worth noticing that the ratio p/π as a function of dN ch /dη in a given p T -bin follows a power-law behavior: As shown in Fig. 3 (top), the same trend is also observed in Pb-Pb collisions.The exponent of the power-law function exhibits the same value in both collision systems (Fig. 3, middle).The same feature is also observed in the Λ/K 0 S ratio (Fig. 3, bottom).The p T -integrated yields and p T are computed using the data in the measured range and extrapolating them down to zero and to high p T (up to 10 GeV/c).The fraction of extrapolated yield for high (low) multiplicity events is about 8% (9%), 10% (12%), 7% (13%), 17% (30%) for π ± , K ± , p and p, Λ and Λ respectively and is negligible for K 0 S .Several parametrizations have been tested, among which the blast-wave function [11] (see below) gives the best description of the data over the full p T range (Fig. 1).Other fit functions [49] (Boltzmann, m T -exponential, p T -exponential, Tsallis-Levy, Fermi-Dirac, Bose-Einstein) have been used to estimate the systematic uncertainty on the extrapolation, restricting the range to low p T for those functions not giving a satisfactory description of the data over the full range.The uncertainty on the extrapolation amounts to about 2% for π ± , K ± , p(p), 3% (8% in low multiplicity events) for Λ( Λ), and it is negligible for K 0 S (since the p T coverage ranges down to 0).The p T increases with multiplicity, at a rate which is stronger for heavier particles, as shown in Fig. 4. A similar mass ordering is also observed in pp [28] and Pb-Pb [10] collisions as a function of multiplicity.
In Fig. 5, the ratios to the pion yields are compared to Pb-Pb results at the LHC and Au-Au and d-Au results at RHIC [50,49,51,48,52,53].While the p/π ratio shows no evolution from peripheral to central events, a small increase is observed in the K/π and Λ/π ratios, accounting for the bin-to-bin correlations of the uncertainties.A similar rise is observed in Pb-Pb, Au-Au and d-Au collisions.This is typically attributed to a reduced canonical suppression of strangeness production in larger freeze-out volumes [54] or to an enhanced strangeness production in a quark-gluon plasma [55].
The observations reported here are not strongly dependent on the actual variable used to select multiplicity classes.Alternative approaches, such as using the total charge in both VZERO-A and VZERO-C detectors, the energy deposited in the ZNA (which originates from neutrons of the Pb nucleus) and   the number of clusters in the first ITS layers reveal very similar trends.In the cases where the largest deviation is observed, the p/π ratio is essentially the same in 0-5% events and it is ∼ 15% higher at p T ∼ 3 GeV/c in the 60-80% class.Part of this difference is due to the mild correlation of events at forward and central rapidity: the lowest multiplicity class selected with ZNA leads to a larger multiplicity at midrapidity than the corresponding class selected with the VZERO-A.

Discussion
In heavy-ion collisions, the flattening of transverse momentum distribution and its mass ordering find their natural explanation in the collective radial expansion of the system [56].This picture can be tested in a blast-wave framework with a simultaneous fit to all particles for each multiplicity bin.This parameterization assumes a locally thermalized medium, expanding collectively with a common velocity field and undergoing an instantaneous common freeze-out.The blast-wave functional form is given by [11] where the velocity profile ρ is described by Here, m T = p 2 T + m 2 is the transverse mass, I 0 and K 1 are the modified Bessel functions, r is the radial distance from the center of the fireball in the transverse plane, R is the radius of the fireball, β T (r) is the transverse expansion velocity, β s is the transverse expansion velocity at the surface, n is the exponent of the velocity profile and T kin is the kinetic freeze-out temperature.The free parameters in the fit are T kin , β s , n and a normalization parameter.
In contrast with the individual fits discussed above, the simultaneous fit to all particle species under consideration can provide insight on the (common) kinetic freeze-out properties of the system.It has to be kept in mind, however, that the actual values of the fit parameters depend substantially on the fit range [10].In spite of this limitations, the blast-wave model still provides a handy way to compare the transverse momentum distributions and their evolution in different collision systems.
The fit presented in this Letter is performed in the same range as in [9,10], also including K 0 S and Λ( Λ).The ranges 0.5-1 GeV/c, 0.2-1.5 GeV/c, 0-1.5 GeV/c, 0.3-3 GeV/c and 0.6-3 GeV/c have been used for π ± , K ± , K 0 S , p(p) and Λ( Λ) respectively.They have been defined according to the available data at low p T and based on the agreement with the data at high p T , justified considering that the assumptions underlying the blast-wave model are not expected to be valid at high p T .Excluding the K 0 S and Λ( Λ) from the fit causes a negligible difference in the fit parameters.
The results are reported in Tab. 5 and Fig. 6.Variations of the fit range lead to large shifts (∼ 10%) of the fit results (correlated across centralities), as discussed for Pb-Pb data in [9,10] As can be seen in Fig. 6, the parameters show a similar trend as the ones obtained in Pb-Pb.Within the limitations of the blast-wave model, this observation is consistent with the presence of radial flow in p-Pb collisions.A detailed comparison of the resulting fit parameters between Pb-Pb [9,10] and p-Pb (Tab.5) collisions shows that at similar dN ch /dη the values of parameters for T kin are similar for the two systems, whereas the β T values are significantly higher in p-Pb collisions.

0.84
Table 5: Blast-wave parameters for simultaneous p-Pb fit of π ± , K ± , K 0 S , p(p) and Λ( Λ) in the fit ranges 0.5-1 GeV/c, 0.2-1.5 GeV/c, 0-1.5 GeV/c, 0.3-3 GeV/c and 0.6-3 GeV/c, respectively.Positive and negative variations of the parameters using the different fit ranges as done in [9,10] are also reported.collisions high multiplicity events are obtained through multiple soft interactions, in p-Pb collisions the high multiplicity selection biases the sample towards harder collisions [57].This could lead to the larger β T parameter obtained from the blast-wave fits.Under the assumptions of a collective hydrodynamic expansion, a larger radial velocity in p-Pb collisions has been suggested as a consequence of stronger radial gradients in [58].
In a hydrodynamically expanding system, the flow coefficients v n are also expected to exhibit a characteristic mass-dependent ordering depending on the transverse expansion velocity.To probe this picture, the p T distributions are fitted simultaneously with the elliptic flow coefficient extracted from two particle correlations v 2 of π ± , K ± , p(p) measured in [59], with the extension of the blast-wave model of [60].This global fit is found to describe the v 2 of pions, kaons and protons relatively well, even if the quality of the fit is slightly worse than that of similar fits in Pb-Pb collisions, in particular for the proton v 2 .Compared to the case where only the particle p T -differential yields are used, the fit results of T kin and β T differ by about 2% only.
Other processes not related to hydrodynamic collectivity could also be responsible for the observed results.This is illustrated in Fig. 6, which shows the results obtained by applying the same fitting procedure to transverse momentum distributions from the simulation of pp collisions at √ s = 7 TeV with the PYTHIA8 event generator (tune 4C) [61], a model not including any collective system expansion.PYTHIA8 events are divided into several classes according to the charged-particle multiplicity at midrapidity |η lab | < 0.3, namely N ch < 5, 5 ≤ N ch < 10, 10 ≤ N ch < 15, 15 ≤ N ch < 20 and N ch ≥ 20.The fit results are shown for PYTHIA8 simulations performed both with and without the color reconnection mechanism [62,63].This mechanism is necessary in PYTHIA tunes to describe the evolution of p T with multiplicity in pp collisions [57].With color reconnection the evolution of PYTHIA8 transverse momentum distributions follows a similar trend as the one observed for p-Pb and Pb-Pb collisions at the LHC, while without color reconnection it is not as strong.This generator study shows that other final state mechanisms, such as color reconnection, can mimic the effects of radial flow [64].
The p T distributions in the 5-10% bin are compared in Fig. 7 with calculations from the DPMJET, Kraków [65] and EPOS LHC 1.99 v3400 [66] models.The QCD-inspired DPMJET [40] generator, which is based on the Gribov-Glauber approach, treats soft and hard scattering processes in an unified way.It has been found to successfully reproduce the pseudorapidity distribution of charged particles in NSD p-Pb collisions at the LHC as reported in [33].On the other hand, it cannot reproduce the p T distribution [67] and the p T of charged particles [57].In the Kraków hydrodynamic model, fluctuating initial conditions are implemented based on a Glauber model using a Monte Carlo simulation.The expansion of the system is calculated event-by-event in a 3+1 dimensional viscous hydrodynamic approach and the S , p(p) and Λ( Λ) in p-Pb Collisions at . . .The ALICE Collaboration freeze-out follows statistical hadronization in a Cooper-Frye formalism.In the EPOS model, founded on "parton-based Gribov Regge theory", the initial hard and soft scattering creates "flux tubes" which either escape the medium and hadronize as jets or contribute to the bulk matter, described in terms of hydrodynamics.The version of the model used here implements a simplified treatment of the collective expansion [66].EPOS predictions including the full hydrodynamic calculation [68] are not available at the time of writing.
The transverse momentum distributions in the 5-10% multiplicity class are compared to the predictions by Kraków for 11 ≤ N part ≤ 17, since the dN ch /dη from the model matches best with the measured value in this class.DPMJET and EPOS events have been selected according to the charged particle multiplicity in the VZERO-A acceptance in order to match the experimental selection.DPMJET distributions are softer than the measured ones and the model overpredicts the production of all particles for p T lower than about 0.5-0.7 GeV/c and underpredicts it at higher momenta.At high-p T , the p T spectra shapes of pions and kaons are rather well reproduced for momenta above 1 and 1.5 GeV/c respectively.Final state effects may be needed in order to reproduce the data.In fact, The Kraków model reproduces reasonably well the spectral shapes of pions and kaons below transverse momenta of 1 GeV/c where hydrodynamic effects are expected to dominate.For higher momenta, the observed deviations for pions and kaons could be explained in a hydrodynamic framework as due to the onset of a non-thermal component.EPOS can reproduce the pion and proton distributions within 20% over the full measured range, while larger deviations are seen for kaons and lambdas.The yield and the shape of the p T distributions of protons are rather well described by both models.In contrast to a similar comparison for Pb-Pb collisions [9,10], in the Kraków calculation the yield of pions and kaons seems to be overestimated.It is interesting to notice that when final state interactions are disabled in EPOS, the description of many pp and p-Pb observables worsens significantly [66].

Conclusions
In summary, we presented a comprehensive measurement of π ± , K ± , K 0 S , p(p) and Λ( Λ) in p-Pb collisions at √ s NN = 5.02 TeV at the LHC.These data represent a crucial set of constraints for the modeling of proton-lead collisions at the LHC.The transverse momentum distributions show a clear evolution with multiplicity, similar to the pattern observed in high-energy pp and heavy-ion collisions, where in the latter case the effect is usually attributed to collective radial expansion.Models incorporating final state effects give a better description of the data.

Fig. 1 :
Fig. 1: (color online) Invariant p T -differential yields of π ± , K ± , K 0 S , p(p) and Λ( Λ) in different V0A multiplicity classes (sum of particle and antiparticle states where relevant) measured in the rapidity interval 0 < y CMS < 0.5.Top to bottom: central to peripheral; data scaled by 2 n factors for better visibility.Statistical (bars) and full systematic (boxes) uncertainties are plotted.Dashed curves: blast-wave fits to each individual distribution.

Fig. 3 :Fig. 4 :
Fig.3: (color online) p/π ratio as a function of the charged-particle density dN ch /dη in three p T intervals in p-Pb (measured in the rapidity interval 0 < y CMS < 0.5) and Pb-Pb collisions (measured at midrapidity).The dashed lines show the corresponding power-law fit (top).Exponent of the p/π (middle) and Λ/K 0 S (bottom) power-law fit as a function of p T in p-Pb and Pb-Pb collisions.The empty boxes show the total systematic uncertainty; the shaded boxes indicate the contribution uncorrelated across multiplicity bins (not estimated in Pb-Pb).

Fig. 5 :
Fig. 5: (color online) Particle yields dN/dy of kaons, protons, and lambdas normalized to pions as a function of dN ch /dη in each V0A multiplicity class (see text for details) measured in the rapidity interval 0 < y CMS < 0.5.The values are compared to results obtained from Pb-Pb collisions at the LHC and Au-Au and d-Au collisions at RHIC measured at midrapidity.The empty boxes show the total systematic uncertainty; the shaded boxes indicate the contribution uncorrelated across multiplicity bins (not estimated in Pb-Pb).

Fig. 6 :
Fig. 6: (color online) Results of blast-wave fits, compared to Pb-Pb data and MC simulations from PYTHIA8 with and without color reconnection.Charged-particle multiplicity increases from left to right.Uncertainties from the global fit are shown as correlation ellipses.

Fig. 7 :
Fig. 7: (color online) Pion, kaon, and proton transverse momentum distributions in the 5-10% V0A multiplicity class measured in the rapidity interval 0 < y CMS < 0.5 compared to the several models (see text for details).

Table 1 :
S , p(p) and Λ( Λ) in p-Pb Collisions at . . .Definition of the event classes as fractions of the analyzed event sample and their corresponding dN ch /dη within |η lab | < 0.5 (systematic uncertainties only, statistical uncertainties are negligible).

Table 4 :
Main sources of systematic uncertainty for the K 0 S and Λ( Λ).
While in Pb-Pb