Multiplicity dependence of light-flavor hadron production in pp collisions at $\sqrt{s}$ = 7 TeV

Comprehensive results on the production of unidentified charged particles, $\pi^{\pm}$, $\rm{K}^{\pm}$, $\rm{K}^{0}_{S}$, $\rm{K}$*(892)$^{0}$, $\rm{p}$, $\overline{\rm{p}}$, $\phi$(1020), $\Lambda$, $\overline{\Lambda}$, $\Xi^{-}$, $\overline{\Xi}^{+}$, $\Omega^{-}$ and $\overline{\Omega}^{+}$ hadrons in proton-proton (pp) collisions at $\sqrt{s}$ = 7 TeV at midrapidity ($|y|<0.5$) as a function of charged-particle multiplicity density are presented. In order to avoid auto-correlation biases, the actual transverse momentum ($p_{\rm{T}}$) spectra of the particles under study and the event activity are measured in different rapidity windows. In the highest multiplicity class, the charged-particle density reaches about 3.5 times the value measured in inelastic collisions. While the yield of protons normalized to pions remains approximately constant as a function of multiplicity, the corresponding ratios of strange hadrons to pions show a significant enhancement that increases with increasing strangeness content. Furthermore, all identified particle to pion ratios are shown to depend solely on charged-particle multiplicity density, regardless of system type and collision energy. The evolution of the spectral shapes with multiplicity and hadron mass shows patterns that are similar to those observed in p-Pb and Pb-Pb collisions at LHC energies. The obtained $p_{\rm{T}}$ distributions and yields are compared to expectations from QCD-based pp event generators as well as to predictions from thermal and hydrodynamic models. These comparisons indicate that traces of a collective, equilibrated system are already present in high-multiplicity pp collisions.


Introduction
Recently, several collective phenomena have been observed in high-multiplicity pp and p-Pb collisions that are reminiscent of observations attributed to the creation of a medium in thermal and kinematic equilibrium in Pb-Pb collisions. In p-Pb collisions, these include the observation of double-ridge structures on the near and away side in two-particle correlation studies [1], non-vanishing v 2 coefficients in multi-particle cumulant studies [2], mass dependent hardening of identified particle p T spectra [3][4][5] and consistency of integrated particle yield ratios with thermal model expectations at high multiplicities [6].
While double-ridge structures have already been observed in high-multiplicity pp collisions [7], a comprehensive study of identified-particle hadrochemistry as well as of the corresponding kinematics after hadronization has not yet been performed in these collisions: such a study is the main topic covered in this paper. The investigation of mass dependent effects as expected within a hydrodynamic evolution scenario requires the measurement of several particle species such as the ones presented here and relies on the excellent particle identification capabilities provided by the ALICE detector.
While similarities in the production of light-flavor hadrons between p-Pb and Pb-Pb collisions at comparable event multiplicities have been discussed previously [4,6], the measurements presented here allow a unique comparison of the observables with several QCD-inspired event generators such as PYTHIA [8] and EPOS [9]. Traditionally, bulk particle production in heavy-ion collisions is described on the basis of equilibrium many-body theories such as hydro-and thermodynamics (see for instance [10,11] and references therein). The continuous transition of light-flavor hadron measurements from pp to p-Pb and Pb-Pb collisions as a function of event multiplicity thus links the dynamic production of particles in individual 2→2 QCD parton-parton scattering processes and subsequent hadronization as an underlying equilibration mechanism to a thermodynamic description of the system.
In a recent letter [12], the ALICE Collaboration reported the multiplicity dependent enhancement of strange (K 0 S , Λ and Λ) and multi-strange (Ξ − , Ω − , Ξ + and Ω + ) particle production in pp collisions at √ s = 7 TeV. In this paper, those results are complemented by the measurement of π ± , K ± , p, p, K*(892) 0 , and φ (1020), as well as by an extended discussion on p T -differential and p T -integrated particle ratios and model comparisons. For the sake of brevity, in this work, (π + + π − ) and (K + + K − ) will be denoted by π ± and K ± , while p refers to (p+p) unless otherwise stated. In addition, (Ξ − + Ξ + ) and (Ω − + Ω + ) will be denoted by Ξ and Ω, while Λ refers solely to the particle and not the antiparticle unless otherwise stated. Finally, (K*(892) 0 +K * (892) 0 ) and φ (1020) will be denoted simply by K * 0 and φ throughout this document. The paper is organized as follows. In section 2, the details of the analysis techniques and of the event classification are described. The results are given in section 3 in which the transverse momentum spectra as well as the extraction of the p T -integrated yields and average transverse momenta are presented. Detailed model comparisons and an interpretation of the results are presented and discussed in section 4.

Analysis
For this analysis, data collected by ALICE in the LHC pp run of the year 2010 are used. In total, the analysis is based on up to 281 million minimum-bias events, corresponding to an integrated luminosity of 4.5 nb −1 . A detailed description of the ALICE apparatus and of its performance can be found in [13,14]. The main subdetectors used in this analysis are the Inner Tracking System (ITS) [15], the Time Projection Chamber (TPC) [16], the Time-Of-Flight detector (TOF) [17], and the V0 scintillator hodoscopes [18]. All tracking detectors are positioned inside a magnetic field B = 0.5 T.
The innermost barrel detector is the ITS, consisting of six cylindrical layers of high-resolution silicon tracking detectors using three different technologies. The two innermost layers are based on silicon pixel technology (SPD) with digital readout. The four outer layers, made of drift (SDD) and strip (SSD) detectors provide analogue readout and thus allow for particle identification via specific energy loss. The ITS, used as a standalone tracker, enables the reconstruction and identification of low momentum particles down to 100 MeV/c that do not reach the TPC.
The TPC is a large cylindrical drift detector of radial and longitudinal dimensions of approximately 85 cm < r < 250 cm and −250 cm < z < 250 cm, respectively. As the main tracking device, it thus provides full azimuthal acceptance for tracks in the pseudorapidity region |η| < 0.9. In addition, it provides particle identification via the measurement of the specific energy loss dE/dx. At low transverse momenta, the dE/dx resolution of 5.2% for a minimum ionizing particle allows a track-by-track identification, while at high transverse momenta the overlapping energy losses can still be statistically distinguished using a multi-gaussian fit to the dE/dx distributions.
Further outwards in radial direction from the beam-pipe and located at a radius of approximately 4 m, the TOF measures the time-of-flight of the particles, providing particle identification over a broad range at intermediate transverse momenta. It is a large-area array of multigap resistive plate chambers with an intrinsic time resolution of 50 ps. The total time resolution includes contributions from the start time determination and amounts to about 120 ps in pp collisions. As described in detail in [19], the start time contribution to the total time resolution improves with increasing number of hits in the TOF in a given event, thus leading to a slight dependence on the event multiplicity and results in a total time resolution of about 100 ps for the highest multiplicities under study.
The V0 detectors are two scintillator hodoscopes that are located on either side of the interaction region at z = 3.3 m and z = −0.9 m, respectively. They cover the pseudorapidity region 2.8 < η < 5.1 and −3.7 < η < −1.7 in full azimuth and are employed for triggering, background suppression, and eventclass determination.
Measurements of unidentified and identified primary particles are reported. Primary particles are defined as any hadron with a mean proper lifetime that is of at least 1 cm/c either produced directly in the interaction or emerging from decays of particles with lifetime shorter than 1 cm/c and excluding particles from interactions with the detector material [20]. The criteria for the selection of primary tracks for π ± , K ± , p and p as well as for the decay products of K* 0 and φ follow the procedures described in [21]. All measurements are corrected for detector acceptance and reconstruction efficiency using Monte Carlo events generated with PYTHIA6 Perugia 0 [22,23] and propagated through the full ALICE geometry with GEANT3 [24]. These events are then reconstructed using the same techniques employed in the case of real data. The corresponding detector acceptance and reconstruction efficiencies are found to be multiplicity-independent within 1% and thus the integrated values are used for all event classes to minimize statistical fluctuations.

Event selection and classification
The data were collected using a minimum-bias trigger requiring a hit either in the A or C side of the V0 (denoted in what follows as V0A or V0C, respectively) or in the SPD, in coincidence with the arrival of proton bunches from both directions. Contamination from beam-gas events is removed offline by using timing information from the V0, which has a time resolution better than 1 ns. Events in which pile-up or beam-gas interactions occurred are also rejected by exploiting the correlation between the number of pixel hits and the number of SPD tracklets. Interactions used for the data analysis are further required to have a reconstructed primary vertex within |z| < 10 cm, where z is in the direction of the beam. Events containing more than one reconstructed vertex are tagged as pile-up occurring within the same bunch-crossing and discarded for the analysis, with up to 10% of all events being tagged in the highestmultiplicity event class considered for analysis. The pile-up tagging was estimated to be efficient enough so that the residual pile-up remaining in the analysed event sample is of no more than 10 −4 to 10 −2 for the lowest and highest multiplicity classes, respectively. The systematic uncertainty associated to pileup rejection was estimated to be smaller than 1% and is therefore not a dominant source of uncertainty for any of the analyses reported here.
The measurements shown here correspond to an event class (INEL > 0) in which at least one charged particle is produced in the pseudorapidity interval |η| < 1 with respect to the beam, corresponding to about 75% of the total inelastic cross-section. In order to study the multiplicity dependence of lightflavor hadron production, the sample is divided into event classes based on the total charge deposited in both of the V0 detectors (V0M amplitude). The V0M amplitude is found to be linearly proportional to the total number of charged particles produced in the pseudorapidity window corresponding to the acceptance of the V0 scintillators. Table 1 also lists the average charged-particle pseudorapidity densities dN ch /dη within |η| < 0.5 for the different event multiplicity classes. The relative standard deviations of the corresponding distributions range from 68% to 30% of the average dN ch /dη for the event class with the lowest average multiplicity to the one with the highest, respectively. These are obtained based on the reconstruction of SPD tracklets which have an acceptance of p T 50 MeV/c. The measurement has been fully corrected for acceptance, tracking and vertexing efficiency as well as for contamination from secondary particles and combinatorial background. Further details can be found in [25,26]. In addition, all quantities reported in this work are corrected for event detection efficiencies using a data-driven unfolding method. This correction is negligible for high-multiplicity event classes but is of up to ∼13% in multiplicity class X. The resulting percentages of the total INEL>0 cross-section, σ /σ INEL>0 , are also reported in Tab. 1. These values are reported after event detection efficiency corrections and do not match the integer boundaries that were used in analysis, e.g. high-multiplicity event classes such as I and II were selected as 0-1% and 1-5% for analysis but event losses at low-multiplicity compress these fractions into 0-0.95% and 0.95-4.7% of the true INEL > 0 cross-section. The analysis-level selection percentiles have been omitted as they are detector-dependent quantities. In previous studies, event classification was based on midrapidity charged-particle densities [26][27][28], as opposed to the forward and backward-pseudorapidity based selection utilized in this work. This choice is motivated by the fact that performing multiplicity selection and data analysis in the same pseudorapidity range may lead to auto-correlation biases and unphysical results. More specifically, hadrochemistry is significantly altered by selection biases, as exemplified by the progression of charged and neutral kaon abundances with multiplicity. If midrapidity-based selections were used, the integrated yields of K ± for high-multiplicity events would be higher than the ones for K 0 S because of the requirement of high charged-particle yields in the same pseudorapidity range. Conversely, if selection is performed with charged particle yields in a different pseudorapidity range than the one in which K ± and K 0 S production rates are measured, similar amounts of charged and neutral kaons would be found across multiplicity, as expected due to their similar masses. This can be readily tested in Monte Carlo simulations, as shown, for instance, in Fig. 1, where the charged and neutral kaon yields in pp collisions simulated with the PYTHIA8 event generator using the Monash 2013 tune are studied as a function of either midrapidity or forward/backward pseudorapidity charged-particle multiplicity. A significant bias towards charged Fig. 1: Multiplicity dependence of charged and neutral kaon yields obtained using mid-pseudorapidity charged particle multiplicities (|η| < 0.5, left) and the charged particle multiplicities within the pseudorapidity range corresponding to the V0A and V0C detectors (denoted by V0M, corresponding to −3.7 < η < −1.7 and 2.8 < η < 5.1, right) in PYTHIA8 simulations of inelastic pp collisions at √ s = 7 TeV.
kaons is observed in the former case, while the latter selection preserves the expected neutral-to-charged kaon ratio of approximately unity.
This discrepancy for charged and neutral kaons is understood to be a consequence of performing selections on charged-particle multiplicities whose probability distributions exhibit a rapid decrease and have low average values. Under such circumstances, any multiplicity selections are likely to isolate fluctuations of charged-particle yields in the reference region of phase space rather than uniformly affecting all particle species regardless of their charge. In the particular case of K ± and K 0 S production rates, residual differences in these kaon yields still arise from resonance decay products, given that φ mesons decay preferentially into charged kaons. However, Monte Carlo studies show that these different feed-down contributions introduce differences in K ± and K 0 S yields of no more than 1-2%, further corroborating the need to take into account the much larger selection bias effects shown in Fig. 1.
However, while multiplicity selections performed in different phase space regions will avoid selection biases, they are also naturally susceptible to the mid-to forward/backward pseudorapidity multiplicity correlation, which for small systems is not as strong as for nuclear collisions [29]. This has the consequence that the reach in midrapidity charged-particle densities is restricted in comparison to same-phase-space selections: when selecting high charged-particle multiplicities in forward/backward pseudorapidity detectors, midrapidity dN ch /dη will eventually saturate, while it will still increase if event selection is performed with detectors at midrapidity. Furthermore, the V0 scintillators that are used in this work for forward/backward charged-particle detection and event classification introduce an imperfect detector response into the analysis. In order to minimize potential biases coming from these factors, all observables are studied as a function of charged-particle density at midrapidity dN ch /dη . By doing so, both dN ch /dη and the variables under study are similarly folded with mid-to forward/backward pseudorapidity multiplicity correlations as well as the detector response within a given event class. This allows comparing results from this study with predictions from models by performing selections on chargedparticle production in the acceptance of the V0 and it has been verified that any residual effect because of the finite detector resolution is negligible.

Unidentified Charged Particles
Spectra of positively and negatively charged particles were obtained separately and summed afterwards. The differences between the final spectra for particles and antiparticles were found to be around 1.5%. The unidentified charged particles were reconstructed using the combined information from ITS and TPC. The p T range of the spectra in all multiplicity classes based on the V0M amplitudes is 0.  GeV/c and the pseudorapidity was limited to |η| < 0.5. This pseudorapidity limit allows a comparison of the charged hadron spectrum with the sum of pions, kaons and protons analyzed in the rapidity range |y| < 0.5 by transforming them to the corresponding pseudorapidity window with the appropriate Jaco- dηdp T for each p T -interval. This cross-check showed a difference of less than 5% which is consistent within the non-common systematics with the expected contributions from electrons, muons, and heavier charged baryons that are counted in addition to pions, kaons, and protons in the charged hadron spectrum.
The contribution from secondary particles was calculated in the same manner as described in detail in Section 2.3. The additional corrections based on FLUKA [30,31] for kaons and anti-protons, which are needed for those specific identified particle measurements in order to account for an imperfect description of absorption cross-section in the detector material, were found to have negligible impact on the unidentified charged particle spectra and were therefore not applied.
The systematic uncertainties are summarized in Table 2. The multiplicity dependence of the tracking efficiency and the feed-down correction were found to be less than 1% and were included in the final systematic uncertainty. The total systematic uncertainty is p T dependent, with values around 6-7% up to 20 GeV/c. It reaches 9.6% for the highest p T bin. The main contributions to the total systematic uncertainty come from the global tracking efficiency (5%) and the parameter variation for the track selection criteria (3-7%). The other sources have a p T dependent contribution of less than 2% each. The systematic uncertainty related to the dependence of the reconstruction efficiency on the MC event generator and the particle composition have been studied as described in [32]. All sources of uncertainty are assumed to be uncorrelated and the total uncertainty was calculated as the quadratic sum of the different contributions. The systematic uncertainty contribution that is uncorrelated across multiplicities was estimated to be 2.1% for all the V0M multiplicity bins over the entire p T range.

Charged pions, kaons, and protons
For the measurement of charged pions, kaons, and protons, several sub-analyses are combined for the comprehensive exploitation of the available PID techniques in ALICE. The spectra cover a range from 0.1/0.2/0.3 GeV/c to 20 GeV/c for π ± / K ± / p(p), respectively, with the exact ranges reported in Tab. 3. Similar approaches were followed in earlier analyses in pp, p-Pb and Pb-Pb collisions [21,33]. An overview of the individual analyses is presented in Tab. 3. Here, we briefly review the most relevant aspects of previously employed techniques: ITS standalone and TPC-TOF. Additionally, we describe methods which are used for the first time for the measurement of p T -spectra of charged kaons, pions, and protons: Bayesian PID, TPC-TOF fits, and TPC template fits.
In the "ITS stand-alone" technique, the average energy loss in the four outer ITS layers is calculated as a truncated mean. For each particle mass hypothesis, the distance between the measured and the expected value is calculated in multiples of the standard deviation σ of the measured energy loss distribution and the particle mass hypothesis with the smallest value is assigned. In contrast to the analysis in the high track density environment of central heavy-ion collisions, the contribution of tracks with wrongly assigned signal clusters is negligible even in the highest pp event multiplicity class. In the intermediate p T -range where a track-by-track identification is feasible, the "TPC-TOF" analysis identifies particles by requiring that the measured energy loss signals in the TPC and time-of-flight in the TOF are within 3σ of the expected value assuming a specific mass hypothesis. This approach finds its natural limitation towards higher momenta, as the expected energy losses and flight times for different species are insufficiently different to allow for a clear separation. The p T -ranges in which this procedure is applicable are given in Tab. 3. Two alternative methods, namely Bayesian PID and TPC-TOF fits, were employed in order to unfold the measured dE/dx and TOF distributions. The Bayesian method of particle identification for the extraction of the minimum-bias spectra of pions, kaons, and protons is described in detail in [34]. The a priori probabilities used in the Bayesian-approach analysis were extracted from the experimental data for the minimum bias event sample using an iterative procedure. The influence of different sets of a priori probabilities, determined for the lowest and highest event multiplicity bins, was evaluated and included in the systematic uncertainties. The actual identification of particles is based on the maximum probability method in which the most likely particle type is assigned to the track.
In the "TPC-TOF Fits" method, the energy loss distribution in the TPC is simultaneously fitted by three gaussian distributions corresponding to charged pions, kaons and protons in each p T and multiplicity bin. Similarly, the velocity distribution of the TOF is fitted for all three species simultaneously. In order to guarantee a sufficient separation of the particle species by minimizing the difference between total and transverse momentum p T , the TOF fits were performed in a narrow η-window (|η| < 0.2) and afterwards transformed to the common rapidity window of |y| < 0.5 assuming a flat distribution in y. Above a p T of ∼ 2 GeV/c, particle identification can still be achieved statistically, rather than on a trackby-track basis, by fitting the specific energy loss in the relativistic rise region with a multi-component fit function, as done in the "TPC Template fits" approach. In this method, the measured dE/dx distribution, in which the distributions of several particle species are overlapping, is fitted with a sum of templates (one for each particle type). The templates are extracted in a data-driven procedure from a pure sample of tagged particles of a given type. This pure sample is obtained from weak decay daughter tracks (p and p from Λ and Λ as well as π ± from K 0 S ) and tracks identified with the TOF (π ± , K ± , p, and p). After a further strict selection of primary-particle-like topologies, the expected dE/dx response is determined in fine bins of momentum and pseudorapidity. The template for each particle species in a given transverse momentum bin in the rapidity window |y| < 0.5 is then obtained by sampling the measured momenta and pseudorapidity values of the tracks in this bin. The individual particle yields are the only free parameters in the fit of the templates to the measured dE/dx distribution.
For all particle species and sub-analyses, contamination from secondary particles at low transverse momenta was subtracted in a data-driven approach on the basis of the measured distance of closest approach of the track to the primary vertex in the transverse plane (DCA xy ), as done in previous work [21]. The DCA xy distribution of the selected tracks was fitted with three Monte Carlo templates corresponding to the expected shapes of primary particles, of secondaries from material (including electrons from photon conversions) and of secondaries from weak decays. The procedure was repeated for each p T and event multiplicity bin and thus takes into account possible differences in the feed-down correction due to a change of the abundances and spectral shapes of weakly decaying strange particles.
The efficiencies obtained for anti-protons and kaons have been additionally corrected based on a comparison of the absorption cross-section used in GEANT3 and the more realistic description of hadronic cross-sections in FLUKA, as in [21,35].
The determination of systematic uncertainties follows the procedures established in previous analyses [21,33]. All the considered contributions are summarized in Tab. 4. Corrections for secondary particles lead to uncertainties of up to 4% for protons and 1% for pions while they are negligible for kaons. The uncertainty in the material budget is of 5% at very low momenta and is related to the energy loss of the particles in the detector material. In addition, inelastic and elastic hadronic scattering processes inside the detector material are described by the transport codes only with limited precision and lead to uncertainties of up to 6% for p (for which the respective cross-section is largest) at low transverse momenta. The track quality selection criteria and the matching of TPC tracks with ITS hits give rise to a systematic uncertainty of the global tracking efficiency that amounts to 4%, independent of p T and

Source
Uncertainty (%) p T (GeV/c) 0. 16 Table 3: Overview of p T ranges used for the combination of the various techniques used for identifying pions, kaons and protons. Since the true rapidity is not known at reconstruction level, fit based analyses ("TPC template fits" and "TPC-TOF fits"), which determine the yield of pions, kaons and protons simultaneously, require an additional η-cut.
particle species. The Lorentz force causes shifts of the cluster position in the ITS, pushing the charge in opposite directions depending on the polarity of the magnetic field of the experiment (E × B effect). In the ITS stand-alone analysis, the uncertainty related to this effect is estimated by analysing data samples collected with opposite magnetic field polarities, for which a difference of 3% is observed. For those sub-analyses (Bayesian PID, TPC-TOF, TPC-TOF fits) that require in addition that the track under study is matched to a hit in the TOF, an additional uncertainty of 3%/6%/4% is taken into account for pions, kaons, and protons, respectively. Following the approach presented in [33], this matching efficiency uncertainty was estimated by repeating the analysis separately for those regions in azimuth in which modules of the Transition Radiation Detector were already present in 2010 and for those in which they were not yet installed.  Specific source 8.1 4.4 6.9 6.7 6.1 12.2 10.5 13.5 11.5 Table 4: Main sources and values of the relative systematic uncertainties of the p T -differential yields of π ± , K ± and p(p). The values are reported for low, intermediate and high p T . The contributions that act differently in the various event classes are removed from the Total (quadratic sum of all contributions), defining the N chindependent ones, which are correlated across different multiplicity intervals. The contribution from the global tracking efficiency is common to all analyses except for the ITS stand-alone (ITSsa).
All sub-analyses were found to be in agreement in the overlapping p T ranges within the uncorrelated part of their respective systematic uncertainties. The final combined spectrum for each particle species was then obtained by calculating the average over all sub-analyses using the uncorrelated part of their systematic errors as weights [? ]. The uncertainties originating from common sources were then added in quadrature to each other and to the uncertainty attributed to the specific particle identification methods.

Weakly-decaying Strange Hadrons
The strange hadrons K 0 S , Λ, Λ, Ξ − , Ξ + , Ω − and Ω + are reconstructed at mid-rapidity (|y| < 0.5) via their characteristic weak decay topology in the channels [36] K 0 Charged particle tracks are selected on the basis of compatibility of their energy loss in the TPC with the expected losses under the pion, kaon and proton mass hypotheses. They are then combined into weak decay candidates following the topology of a V-shaped decay for K 0 S , Λ and Λ (denoted 'V0' decays) and a combination of a V0 decay and one additional charged track for Ξ − , Ξ + , Ω − and Ω + (denoted 'cascade' decays). In addition to several geometrical criteria on the arrangement of decay daughter tracks, K 0 S , Λ and Ω − candidates are required to have a calculated mass that is incompatible with other species that decay in a similar topological arrangement, which are Λ, K 0 S and Ξ − , respectively. This selection is commonly denoted 'competing decay rejection' and the exact numerical value depends on the invariant mass resolution for the competing particle species. Furthermore, candidates whose proper lifetimes are unusually large for their expected species are also rejected to avoid combinatorial background from interactions with the detector material. The selection criteria used to define V0 and cascade decay candidates are listed in Tab. 5 and Tab. 6, respectively.

V0 selection criterion
Value Table 5: Selection criteria parameters utilized in the K 0 S , Λ and Λ analyses presented in this work. If a criterion for Λ and K 0 S finding differs, the criterion for the Λ hypothesis is given in parentheses. The acronym DCA stands for "distance of closest approach" and PV for "primary event vertex". The pointing angle θ is the angle between the momentum vector of the reconstructed V0, and the line segment bound by the decay and primary vertices and R 2D denotes the transverse distance from the detector center.
Particle yields are calculated in p T and event multiplicity intervals by extracting the relevant signals from invariant mass distributions as done in previous work [4,6,37].

V0 from cascade selection criterion
Value > 0.6 (0.5) cm Cascade pointing angle cos(θ casc ) > 0.97 Proper Lifetime < 3cτ Competing Cascade Rejection Window (Ω ± only) ±8 MeV/c 2 Table 6: Selection criteria for V0 (Λ) from cascades, and cascades (Ξ ± and Ω ± ) presented in this work. If a criterion for Ξ ± and Ω ± finding differs, the criterion for Ω ± hypothesis is given in parentheses. DCA stands for "distance of closest approach" and PV for "primary event vertex.". The pointing angle θ is the angle between the momentum vector of the reconstructed V0 or cascade and the line segment bound by the decay and primary vertices and R 2D denotes the transverse distance from the detector center. The cascade track curvature is neglected, and τ refers to the average lifetime for the two different cascade species. Approximately 20% of the measured Λ (Λ) signals are from Ξ − (Ξ + ) and Ξ 0 (Ξ 0 ) decays. These feeddown contributions were subtracted using a data-driven approach in which the measured Ξ − (Ξ + ) spectra are used as input and a simulation is used to evaluate the fraction of reconstructed Λ (Λ) coming from Ξ − (Ξ + ) decays. Since production rates of Ξ 0 and Ξ 0 have not been measured, their contribution is estimated by assuming that they are as abundant as their charged counterparts and that their momentum distributions are identical.
Because in the specific case of the cascade analysis the measurement is performed in large momentum intervals because of the limited amount of data, efficiencies are reweighted in each p T bin to take into account differences between generated and real data spectral shapes.
Systematic uncertainties for K 0 S , Λ, Λ, Ξ − , Ξ + , Ω − and Ω + are estimated following the procedure described in [4,6]. The main sources of systematic uncertainty in these measurements are track selections (up to ∼6%), knowledge of detector materials (4%), feed-down from Ξ − (Ξ + ) and Ξ 0 (Ξ 0 ) for the Λ (Λ) (up to ∼4%) and topological selections, which contribute with a ∼1-8% uncertainty. The contributions to systematic uncertainties are summarized in Tab. 7. As in previous work, the study of systematic uncertainties was repeated for all event classes to determine differences in how each contribution affects results from each of these classes.

Resonances
The K * 0 and φ -mesons are reconstructed at mid rapidity |y| < 0.5 via their hadronic decay channels into charged particles, Both the TPC and TOF information are used to identify charged particles as pions or kaons from K * 0 decays whereas only TPC information is used to identify charged particles as kaons from decays of φ -mesons, as in the latter case the combinatorial background is significantly smaller.
Pairs of pions and kaons (pairs of kaons) of opposite charge are considered to obtain the invariant mass distribution of K * 0 (φ ) decay candidates. An event mixing technique is used to estimate the combinatorial background. The mixed-event distribution is normalized in the mass region outside of the mass peak, i.e at 1.1 < M πK (GeV/c 2 ) < 1.15 and 1.035 < M KK (GeV/c 2 ) < 1.045 for K * 0 and φ -mesons, respectively. The normalized mixed-event distribution is subtracted from the same event unlike-sign distribution to isolate the relevant signals. After mixed-event background subtraction, each invariant mass distribution is fitted with a Breit-Wigner function (Voigtian function) for the signal and a 2nd-order polynomial for any residual background. The parameterizations for the signal are given in Eq. 1 for the K * 0 and Eq. 2 for the φ -meson: where M πK and M KK are the reconstructed invariant masses of K * 0 and φ -meson candidates, M 0 , Γ and Y are the mass, width and raw yield of the resonances, respectively. The parameter σ represents the mass resolution. Figure 3 shows the invariant mass of πK (KK) in the left (right) panel for 2 < p T < 2.5 GeV/c in the V0M event multiplicity class I.   The raw yields are extracted in each p T bin and event multiplicity interval as done in previous works [5,38,39]. In this analysis, detector acceptance and reconstruction efficiency are re-weighted in each p T bin to take into account the differences between generated and real data spectral shapes.
The sources of systematic uncertainties for K * 0 and φ -meson production in pp collisions are the TPC-ITS matching efficiency, track selection criteria, PID, yield extraction method, hadronic interaction and material budget, and were evaluated following the same prescription used in previous works [38,39]. The main source of uncertainty for K * 0 and φ comes from the determination of the TPC-ITS track matching efficiency. This contribution has been estimated to be a p T -independent effect of 3% for charged particles [40], which results in a 6% effect when any two primary tracks are combined in the invariant mass analysis of K * 0 and φ . For both K * 0 and φ , the uncertainties due to various track selection cuts from low to high-p T are found to be 3.0-0.9% and 1.6-2.4%, respectively. The systematic uncertainty due to the signal extraction includes variations in the fit range, fit function and normalization range and is of ∼5-10% (∼3-9%) from low to high p T for K * 0 (φ ). The uncertainty due to different PID selection methods is estimated to be ∼2-4% (∼1-2%) for K * 0 (φ ). The knowledge of the material budget for both K * 0 and φ contributes to ∼4% and ∼6% at low p T and is negligible at high p T . The contribution from the estimate of the hadronic interaction cross section in the detector material at low-p T is ∼4% (∼6%) for K * 0 (φ ) and negligible at high-p T . The total systematic uncertainty for K * 0 and φ is estimated to be about 12% and 10%, respectively. The maximum value of the multiplicity-independent systematic uncertainty is found to be ∼8% (∼5%) for K * 0 (φ ). The main contributions to the systematic uncertainties are summarized in Tab. 8.
The systematic uncertainties were studied independently for all event classes, in order to separate the sources that are multiplicity-dependent and uncorrelated across multiplicity bins. In particular, signal extraction and PID are fully uncorrelated sources, whereas global tracking, track cuts, material budget and hadronic cross sections are correlated among different event multiplicity classes.    Table 9: Overview of the systematic uncertainties associated to the low-p T extrapolation used to calculate dN/dy and p T for the various particle species. Values for the highest and the lowest multiplicity classes are given. If the dependence with multiplicity is negligible, only a single value is given.

Transverse momentum distributions
The transverse momentum distributions measured at midrapidity for the event classes defined in Tab. 1 are shown in Fig. 4 for unidentified charged particles (|η| < 0.5) and Fig. 5 for π ± , K ± , K 0 S , K * 0 , p, p, φ , Λ, Λ, Ξ − , Ξ + , Ω − and Ω + (|y| < 0.5). In the particular case of the φ , K * 0 and Ω − measurements, some event classes were merged to allow for sufficient statistics. Particle and antiparticle as well as charged and neutral kaon production rates are compatible within uncertainties.
Transverse momentum spectra are observed to become harder with increasing charged-particle multiplicity, with absolute changes in the spectrum shapes being more pronounced for particles with larger mass. The evolution of the p T distributions with respect to the spectra in the INEL>0 event class for the various particle species is shown in Fig. 6 and is observed to be identical for the two π ± and K ± mesons as well as for the p, Λ, Ξ − baryons and their corresponding anti-particles. The spectra modification of φ and K * 0 resonances follows the trend observed for baryons at p T < 2 GeV/c while for larger momenta  the modification is similar to the one observed for other mesons. Given that these mesonic resonances have a significantly higher mass than that of π ± and K ± , this suggests that the spectra evolution with multiplicity is driven by the hadronic mass at low p T and by the number of constituent quarks at higher p T . It is also interesting to note that such behavior is not unique to high multiplicity, but is present even for the lowest multiplicity class, where mass-dependent mechanisms such as radial flow are not expected to play a significant role.
A comparison of the multiplicity dependence of the transverse momentum distributions can be performed by studying the ratios p/φ = (p + p)/φ , K/π = (K + + K − )/(π + + π − ), p/π = (p + p)/(π + + π − ) and Λ/K 0 S as a function of p T , as shown in Fig. 7 for the lowest and highest multiplicity classes in this work.
Results are compared with measurements in p-Pb collisions at √ s NN = 5.02 TeV [4,6,41] as well as in Pb-Pb at √ s NN = 2.76 TeV [42,43]. The p/φ ratio is observed to decrease with p T in pp and low-multiplicity p-Pb collisions but is seen to become progressively flatter for high-multiplicity p-Pb and Pb-Pb collisions. Given the similar mass of the particles involved in this ratio, it is possible that this flattening is a signature of significant radial flow in the larger systems. Furthermore, the baryon-tomeson ratios p/π and Λ/K 0 S exhibit a characteristic depletion at p T ∼ 0.7 GeV/c and an enhancement at intermediate p T (∼ 3 GeV/c), which is qualitatively similar to that observed in p-Pb and Pb-Pb collisions. Finally, the K/π ratio is observed to be relatively multiplicity-independent, except for central Pb-Pb collisions, where a weak depletion (enhancement) at low (intermediate) p T is visible.
While the observed changes in these particle ratios are quantitatively different in the various collision systems, it is worth noting that the final state charged particle multiplicities also cover very distinct ranges. If considered as a function of dN ch /dη, the ratios measured in specific low-, mid-and highp T intervals shown in Fig. 8 are seen to depend on multiplicity in a remarkably similar manner for all collision systems, despite differences in energy and collision geometry.       Transverse momentum dependence of p/φ = (p + p)/φ , K/π = (K + + K − )/(π + + π − ), p/π = (p + p)/(π + + π − ) and Λ/K 0 S (from top to bottom row) yield ratios in pp (left), p-Pb (middle), Pb-Pb (right) collisions for highand low-multiplicity classes, respectively. Cancellation of common systematic uncertainties in the numerator and denominator was carried out only for the Λ/K 0 S , as in other cases the cancellation is non-trivial because of the use of various combined identification techniques or, in the case of resonances, of significantly different analysis strategy. Reference p-Pb and Pb-Pb data from [4,6,[41][42][43].   Fig. 8: Transverse momentum dependence of p/φ = (p + p)/φ , K/π = (K + + K − )/(π + + π − ), p/π = (p + p)/(π + + π − ) and Λ/K 0 S (from top to bottom row) yield ratios in low-, mid-and high-p T intervals in pp, p-Pb, Pb-Pb collisions as a function of dN ch /dη . Reference p-Pb and Pb-Pb data from [4,6,[41][42][43].

Integrated yields and average momenta
The p T -integrated yields dN/dy and mean transverse momenta p T have been calculated using data in the measured momentum range and a Lévy-Tsallis parametrization to extrapolate the spectra down to zero p T , similarly to what was done in previous studies [21,37,38]. Several functional forms, such as Boltzmann, m T -exponential, p T -exponential, Fermi-Dirac and Bose-Einstein functions, were used to estimate the systematic uncertainties associated to this procedure. For those functions that are unable to describe the p T distributions in the full measured range, the fitting was restricted only to low-p T . The uncertainty on dN/dy and p T associated to this procedure is shown in Tab  The p T values are observed to increase with multiplicity for all measured particle species, with the increase being more pronounced for heavier particles. This observation resembles that of previous measurements in p-Pb [4] and Pb-Pb [33] in which p T values exhibit mass-ordering for π ± , K ± , p, p, Λ, Λ, Ξ − , Ξ + , Ω − and Ω + . However, while in central Pb-Pb collisions particles with similar mass such as φ and p had similar p T for central collisions, this is not the case for high-multiplicity pp events, where the p T of the φ is significantly higher compared to that of the p, as can be seen in Fig. 9. This has also been observed in inelastic pp and in p-Pb [39] and is further indication that radial flow is the dominant mechanism determining spectra shapes only in very high multiplicity Pb-Pb. It is also interesting to note that at similar multiplicities p T values for identified particles in pp and p-Pb are compatible within uncertainties despite differences in initial state and collision energy, pointing to a common mechanism at play in these two systems.
The multiplicity dependence of identified particle yields relative to pions is compared to p-Pb and Pb-Pb results in Fig. 10. Particles with a larger strangeness content are observed to be produced more abundantly with multiplicity, as can be seen in the Λ/π, Ξ/π and Ω/π ratios [12]. The p/π ratio is observed to be constant within uncertainties except for the lowest multiplicity event class. This indicates that the increase of hyperon production with respect to pions is a phenomenon that does not originate from mass differences but is connected to strangeness content. Furthermore, the relative production of the φ increases with dN ch /dη by approximately 20%, similarly to the Λ, a single-strange baryon. This suggests that φ production cannot be described solely by considering net strangeness or number of strange quark constituents. The K * 0 /π ratio, on the other hand, exhibits a hint of a decrease with multiplicity. In nuclear collisions, this decrease is more pronounced and is typically considered a consequence of the rescattering of K * 0 decay daughters during a hadronic phase of the system evolution [5,44,45]. All these observations are consistent with previous measurements in p-Pb [4,6,39] and indicate that relative         particle abundances can be described as a universal function of charged-particle density in pp and p-Pb collisions.

Comparison to Monte Carlo models
The multiplicity dependence of the p T -differential K/π, p/π and Λ/K 0 S ratios is compared to several Monte Carlo (MC) event generators, see Fig. 11. All predictions are obtained by performing selections on charged-particle multiplicities in the V0M acceptance and are compared to data parametrically as a function of dN ch /dη , as discussed in section 2.1. The PYTHIA8 event generator in its Monash tune [8,46] is only successful in describing the qualitative features of the evolution of the baryon-tomeson ratios if color reconnection (CR) is allowed to occur, as observed already in previous work [47]. In contrast to the string-model based PYTHIA, the HERWIG code implements hadronization in a clustering approach [48]. As shown in Fig. 11, the abilities of HERWIG7 to describe particle production at low and intermediate p T are still limited, but are currently being improved [49]. DIPSY, a model in which fragmenting strings are allowed to form color ropes which then hadronize with a higher effective string tension [50], is also able to reproduce the decreasing (increasing) trend of the baryon-to-meson ratios at low (high) p T , but fails in describing the absolute values of these ratios. Furthermore, the EPOS-LHC event generator, which relies on parton-based Gribov-Regge theory and includes elements from hydrodynamics [9], predicts increased baryon-to-meson ratios at intermediate p T as a consequence of radial flow, but overestimates the multiplicity dependence of these ratios and thus fails to quantitatively reproduce the measurements reported here. Finally, essentially all models are able to reproduce the fact that the K/π ratio is multiplicity-independent, while not necessarily describing the absolute value well. These findings suggest that there is more than one physical mechanism that would lead to the dynamics observed in p T -differential identified particle ratios, and a systematic study of other observables such as the flow coefficients v n is required in order to discriminate among the various possibilities.
The p T -integrated ratios are compared to predictions from the same Monte Carlo event generators in Fig. 12. The PYTHIA8 and HERWIG generators incorrectly predict no multiplicity dependence of relative strangeness production and therefore fail especially in the description of multi-strange baryon production, while DIPSY and EPOS-LHC exhibit increased strangeness in high-multiplicity pp collisions but fail to predict a flat p/π ratio. None of the tested generators correctly reproduces the multiplicity dependence of the K * 0 /π ratio, which is observed to decrease slightly. This comprehensive set of measurements provides essential input for all models aiming to describe flavor generation in pp collisions.

Comparison to Blast-Wave model predictions
The evolution of the p T distribution with multiplicity in pp collisions is remarkably similar to the evolution observed in larger colliding systems, as underlined in Section 3. In larger systems, this evolution can be interpreted as originating from the hydrodynamical radial expansion of the produced medium [4,33] that can be studied by means of the Boltzmann-Gibbs Blast-Wave model (BG-BW) [51]. This model assumes a locally thermalized medium which expands with a common velocity field and then undergoes an instantaneous freeze-out phase. The average expansion velocity β T and the kinematic freeze-out temperature T kin can be extracted with a simultaneous fit to the p T distribution of several particles for each multiplicity bin and the result can be used to predict the p T spectra of particles with different masses. The result of the simultaneous BG-BW fit to π, K and p for the combined I+II V0M multiplicity class in pp collisions at √ s = 7 TeV is shown as solid lines in Fig. 13 (top) for the p T ranges 0.5-1 GeV/c, 0.2-1.5 GeV/c and 0.3-3 GeV/c respectively. The predicted spectra for K 0 S , Λ, Ξ, Ω, K * 0 and φ are shown as dashed lines. The ratios between the data points and the fits or predictions are plotted in the bottom panels. Strange and multi-strange hadron spectra are reasonably well predicted by the BG-BW model in a p T range which gets larger as the mass of the particle increases. This indicates that strange particles may follow a common motion together with lighter hadrons and suggests the presence of radial flow even in pp collisions. Resonances seem not to follow a similar expansion pattern, as there is no p T region where the ratio data/prediction is flat. This discrepancy between BG-BW predictions and resonance  The fit to pp spectra of π, K and p has been performed for all the analyzed multiplicity bins and the values of β T and T kin are compared to those obtained for p-Pb and Pb-Pb collisions in Fig. 14. The fitting ranges are the same for all the three systems and a common color palette is used to highlight the average multiplicity corresponding to each point. Ellipses correspond to the 1-sigma contour, estimated by fitting the p T -differential spectra with total uncertainties, i.e. after adding statistical and systematic uncertainties in quadrature.
Spectra in pp and p-Pb lead to very similar β T and T kin values when considering similar multiplicities, while in Pb-Pb at similar multiplicities, lower β T are observed with respect to the other two systems. This behavior has been interpreted to be a consequence of stronger radial flow gradients in the smaller collision systems [52].
The CMS Collaboration has recently reported a similar study in [53], where a BG-BW fit to the spectra   of K 0 S and Λ has been performed for the three colliding systems, selecting the multiplicity with a central rapidity estimator. The β T -T kin progression is found to be different for pp and p-Pb collisions at high multiplicity, but the numerous differences with respect to the analysis discussed here (multiplicity estimator, particles included and treatment of systematic errors in the fits) do not allow for a quantitative comparison.

Comparison to statistical hadronization models
The measured abundances of hadrons produced in heavy-ion collisions have been successfully described by statistical hadronization models over a wide range of energies [54][55][56]. Statistical model calculations for central ultra-relativistic heavy-ion collisions are typically carried out in the grand-canonical ensemble. However, a grand-canonical description of particle production is only applicable if the volume V = R 3 of the system is sufficiently large and as a rule of thumb one needs V T 3 > 1 for a grand-canonical description to hold [57,58]. This condition must be fulfilled for each conserved charge separately. Several attempts were made to extend the picture of statistical hadronization to smaller systems such as pp Kinematic freeze-out temperature parameter T kin versus average expansion velocity β T from a simultaneous BG-BW fit to π, K and p spectra measured in pp, p-Pb and Pb-Pb collisions. The color of the datapoints indicates the corresponding average charged-particle multiplicity density.
or even e + e − collisions [59][60][61]. While particle ratios of non-strange particles are observed to be similar in small and large systems, the relative production of strange particles appears to be significantly lower in smaller systems. The data presented in [12] show for the first time that there is a continuous increase of strangeness production with increasing event multiplicity across various collision systems. In the strangeness canonical approach, it is assumed that the total amount of strange hadrons in the volume is small with respect to non-strange hadrons. Thus the conservation of strangeness is guaranteed locally and not only globally while the bulk of the particles is still described in the grand-canonical ensemble. For further details on this approach, we refer for instance to [62][63][64]. The study presented here utilizes the implementation in the THERMUS 3.0 code [65]. Alternative implementations of the statistical model are for instance adopted in the codes of the GSI-Heidelberg [11,56] and SHARE [66] groups.

Correlation volume for strangeness production
Previous studies based on THERMUS were targeted to describe the evolution of multi-strange particle production in p-Pb collisions as a function of event multiplicity [6]. In this case, the volume for particle production was chosen such that the charged pion multiplicity dN π /dy at midrapidity corresponding to a window of y = ±0.5 units was correctly described by the model. This approach is equivalent to calculating strangeness suppression for a system whose total extension only corresponds to one unit of rapidity. Consequently, the model describes qualitatively the suppression pattern, but overestimates the reduction of strangeness at small multiplicities. The discrepancy increases even further if this approach is extended to the smaller multiplicities which are covered in pp collisions.
In the study presented here, a similar approach as in [6] is followed: the strangeness saturation parameter is fixed to γ S = 1, the chemical potentials associated to baryon and electric charge quantum numbers are set to zero. The chemical freeze-out temperature T ch is varied from 146 to 166 MeV as in [6] following recent results from Lattice QCD calculations and their uncertainties [67] as well as by fits to experimental data from central Pb-Pb collisions [11,68]. Ratios of the production yields to pions are investigated for several particle species. In order to cancel the influence of the freeze-out temperature and to isolate the volume dependence, all ratios except for K * 0 are normalized to the high multiplicity limit, i.e. the grandcanonical saturation value for the model and the mean ratio in the 0-60% most central Pb-Pb collisions  Fig. 15: Ratios of the production yields to pions for several particle species as a function of the midrapidity pion yields for pp, p-Pb and Pb-Pb colliding systems compared to the THERMUS model prediction for the strangeness canonical suppression (black line), in which only the system size is varied. All values except for the K * 0 are normalized to the high multiplicity limit (see text for details). The upper axis shows the radius R of the correlation volume V = R 3 which corresponds to the predicted particle ratios. The width of the model prediction line corresponds to a variation of the chemical freeze-out temperature between 146 MeV and 166 MeV. Reference p-Pb and Pb-Pb data from [4,6,[41][42][43].
for the data. As the production rates of short-lived resonances in central heavy-ion collisions might be reduced by re-scattering effects in the hadronic phase [5,39], the values for K 0 * were normalized to the most peripheral bin in Pb-Pb collisions. The resulting double ratios are shown in Fig. 15.
In contrast to [6], the total volume V of the fireball was determined differently. While the strangeness conservation volume is also imposed to be of the size of the fireball (V = V C ), its absolute magnitude is not fixed to reproduce the pion multiplicity in a window of y = ±0.5 at midrapidity. Instead, it is fixed to reproduce a pion multiplicity which is larger by a factor k and thus corresponds to a larger rapidity window assuming a flat dependence of particle production as a function of rapidity. The same factor k was used for all particles and multiplicities. In practice, k corresponds to a constant scaling factor of dN π /dy (the x-axis in Fig. 15), and this feature is used for its numerical determination: the exact value of k was optimized in a one-parameter fit of the functions describing the evolution of the double ratios versus dN π /dy to the experimental data. For the determination of the systematic uncertainty on the value of k, an alternative normalization scheme for the data was applied (normalization to the highest available centrality bin) and the procedure for the determination of k was repeated. A value of k = 1.35 ± 0.28 is obtained corresponding to a rapidity window of y = ±0.67. The results thus indicate that the total correlation volume for strangeness production extends over about 1.35 units in rapidity. In other words, strangeness as a conserved quantity in QCD can be effectively equilibrated over this distance in the system. Similar values can be obtained from purely theoretical considerations on causality constraints [69]. Furthermore, the size of the correlation window is also compatible with similar estimates for charm production [70,71]. We also note that this value is smaller than the plateau in the rapidity distribution at LHC energies which extends typically over three to four units [72,73] and is thus meaningful from a physical point of view.

Comparison to experimental data
As shown in Fig. 15 this approach allows for a qualitative description of particle ratios as a function of event multiplicity. They can be naturally described within the framework of strangeness canonical suppression. The deviations observed for the K * 0 meson in central Pb-Pb collisions can be ascribed to the aforementioned re-scattering effects in the hadronic phase [74]. Furthermore, differences between the model and the data in the most peripheral Pb-Pb collisions in case of Ω and Ξ can be potentially reduced with core-corona corrections [75].
From a quantitative point of view, essentially all data points can be described within 1-2 standard deviations. A potentially different trend is only observed for the φ -meson for which -as a strangeness-neutral particle -a flat dependence as a function of event multiplicity is expected from the model, but which shows a rising and falling trend in data. Future experimental data will be needed in order to clarify the significance of this deviation. It must be noted, though, that the φ -meson also deviates from a common blast-wave fit to other light-flavor hadrons in peripheral Pb-Pb collisions indicating an out-of-equilibrium production except for most central Pb-Pb collisions [5].
Independently of the experimental precision and possible higher order effects in the particle production mechanisms, we find that the strangeness canonical approach can reproduce the multiplicity dependence of all measured light-flavor hadrons across various collision systems to within 10-20%.

Summary and conclusions
A comprehensive set of unidentified and identified particle spectra has been measured in pp collisions at √ s = 7 TeV as a function of charged-particle multiplicity, complementing the existing measurements in p-Pb and Pb-Pb collisions and allowing for a comparison of these different collision systems. In pp collisions, all transverse momentum spectra are observed to become harder with progressively larger charged-particle multiplicity density, with the effect being more pronounced for heavier particles. In addition, the modification of spectra with respect to the inclusive measurement follows a different pattern for mesons and baryons, except for resonances, which follow baryons at low-p T and tend to be modified similarly to mesons at mid-to high-p T . Furthermore, it has been demonstrated that the evolution of the baryon/meson ratios as a function of dN ch /dη exhibits a universal pattern for all collision systems. This behavior might indicate a common mechanism at work that depends solely on final-state multiplicity density. A similar statement can also be made for integrated particle ratios, which are observed to depend on dN ch /dη in approximately the same way for any colliding system despite crucial differences in the initial states as well as colliding energies. The simplest interpretation of this similarity is that in both cases the final state is a thermalized system whose volume at hadronization is proportional to the chargedparticle multiplicity density.
In order to test the assertion of equilibration more quantitatively, the Blast-Wave model was employed to check for kinematic equilibrium, and a statistical hadronization model employing a strangenesscanonical approach was used to check for chemical equilibrium. In all cases, all particle species except for resonances are described within 10-20% by these models. It is also interesting to note that, within the statistical hadronization model employed here, the correlation volume over which strangeness production is equilibrated extends over approximately 1.35 units of rapidity.
The results are also compared to predictions from event generators, which are only able to describe the evolution of p T -differential particle spectra with dN ch /dη if mechanisms such as color reconnection, color ropes or radial flow, as is the case for PYTHIA, DIPSY and EPOS-LHC, respectively, were present. To distinguish between these different mechanisms requires more studies, both experimental and theoretical.
The multiplicity dependence of relative abundances of identified particles was also compared to several event generators and it was found that no generator is able to fully describe the whole observed dynamics satisfactorily. This comprehensive set of results therefore provides a challenge to the theory community and represents an opportunity to study not only pp collisions specifically, but also hadronization at high energies in general.