Multiplicity dependence of pi, K, and p production in pp collisions at root s=13 TeV

This paper presents the measurements of π ± , K ± , p and p transverse momentum ( p T ) spectra as a function of charged-particle multiplicity density in proton–proton (pp) collisions at √ s = 13 TeV with the ALICE detector at the LHC. Such study allows us to isolate the center-of-mass energy dependence of light-ﬂavour particle production. The measurements reported here cover a p T range from 0.1 to 20 GeV / c and are done in the rapidity interval | y | < 0 . 5. The p T -differential particle ratios exhibit an evolution with multiplicity, similar to that observed in pp collisions at √ s = 7 TeV, which is qualitatively described by some of the hydrodynamical and pQCD-inspired models discussed in this paper. Furthermore, the p T -integrated hadron-to-pion yield ratios measured in pp collisions at two different center-of-mass energies are consistent when compared at similar multiplicities. This also extends to strange and multi-strange hadrons, suggesting that, at LHC energies, particle hadrochemistry scales with particle multiplicity the same way under different collision energies and colliding systems.


Introduction
The unprecedented energies available at the Large Hadron Collider (LHC) provide unique opportunities to investigate the properties of strongly-interacting matter. Particle production at large transverse momenta ( p T ) is well-described by perturbative Quantum Chromodynamics (pQCD). The soft regime ( p T 2 GeV/c), in which several collective phenomena are observed in proton-proton (pp), proton-lead (p-Pb), and heavy-ion (A-A) collisions, is not calculable from first principles of QCD. Instead, in order to describe bulk particle production in A-A collisions, one usually relies on hydrodynamic and thermodynamic modelling, which assumes the system to be in kinetic and chemical equilibrium [1,2]. On See Appendix A for the list of collaboration members. e-mail: alice-publications@cern.ch the other hand, the description of lowp T particle spectra in smaller systems such as pp collisions is often based on phenomenological modelling of multi-partonic interactions (MPI) and color reconnection (CR) [3,4] or overlapping strings [5].
Recent reports on the enhancement of (multi)strange hadrons [6], double-ridge structure [7,8], non-zero v 2 coefficients [9], mass ordering in hadron p T spectra, and characteristic modifications of baryon-to-meson ratios [10] suggest that collective phenomena are present at the LHC energies also in p-Pb collisions. This is further extended to even smaller systems, such as pp collisions at √ s = 7 TeV, where similar observations have been reported in high multiplicity events, indicating that the collective effects are not characteristic of heavy-ion collisions only. Furthermore, a continuous transition of light-flavor hadron-to-pion ratios as a function of charged-particle multiplicity density d N ch /dη from pp to p-Pb and then to Pb-Pb collisions was found [11][12][13]. The observed similarities suggest the existence of a common underlying mechanism determining the chemical composition of particles produced in these three collision systems.
Results from pp [11] and p-Pb [10] collisions indicate that particle production scales with d N ch /dη independent of the colliding system. Measurements reported in previous multiplicity-dependent studies have considered different colliding systems, each at a different center-of-mass energy. In this work, we extend the existing observations by performing a detailed study of pp collisions at √ s = 13 TeV. A similar study has been reported by the CMS Collaboration, albeit in a limited p T range [14]. Thanks to the availability of Run 2 data from the LHC, for the first time, in pp collisions, we can disentangle the effect of center-of-mass energy from the multiplicity dependence of π ± , K ± and p (p) production in a wide p T range.
In this paper, we report on the multiplicity dependence of the production of primary π ± , K ± and p (p) at √ s = 13 TeV. Particles are considered as primary if their mean proper decay length cτ is larger than 1 cm and they are created in the collision (including products of strong and electro-magnetic decays), but not from a weak decay of other lightflavor hadrons or muons. An exception to this are products of weak decays, where cτ of the weakly decaying particle is less than 1 cm [15]. The reported particle spectra are measured in the rapidity region |y| < 0.5 with the ALICE detector [16], which offers excellent tracking and particle identification capabilities from p T = 0.1 GeV/c to several tens of GeV/c [17]. As particles and anti-particles are produced roughly in equal amounts at LHC energies [18], we adopt a notation where π , K, and p refer to (π + + π − ), (K + + K − ), and (p + p) unless stated otherwise. This paper is organized as follows. In Sect. 2, the details on particle identification techniques, systematic uncertainties, spectra corrections and normalization are provided. The results are presented and discussed in Sect. 3, together with comparisons to Monte Carlo model predictions. Finally, the most important findings are summarized in Sect. 4.

Data set and experimental setup
The dataset used for this study was recorded by the ALICE Experiment during the 2016 LHC pp data taking period. Overall ∼143M events have been analysed, corresponding to an integrated luminosity of 2.47 nb −1 considering the visible cross-section measured with the V0 detector [19]. A detailed description of the ALICE detector and its performance is provided in [16,17]. Measurements of identified particle spectra have been performed by using the central barrel detectors: the Inner Tracking System (ITS) (Sect. 3.1 of [16]), the Time Projection Chamber (TPC) [20] and the Time-of-Flight detector (TOF) [21]. The charged-particle multiplicity estimation is done by the V0 detector (Sect. 5.4 of [16]), which consists of two arrays of 32 scintillators each, positioned in the forward (V0A, 2.8 < η < 5.1) and backward (V0C, −3.7 < η < −1.7) rapidity regions. In addition, the V0 is also used for triggering purposes as well as background rejection. The determination of the event collision time [22] is performed by the T0 detector as well as the TOF detector. The former consists of two arrays of Cherenkov counters, positioned on both sides of the interaction region, and covering a pseudorapidity range of −3.3 < η < −2.9 (T0-C) and 4.5 < η < 5 (T0-A). The central barrel detectors are placed inside a solenoidal magnet, which provides a field strength of 0.5 T.
The ITS is the innermost detector and consists of six concentric cylindrical layers of high-resolution silicon detectors based on different technologies, covering pseudorapidity region |η| < 0.9. The two innermost layers form the Silicon Pixel Detector (SPD), which features binary readout and is also used as a trigger detector. The Silicon Drift Detector (SDD) and the Silicon Strip Detector (SSD), which form the four outer layers of the ITS, provide the amplitude of the charge signal, which is used for particle identification through the measurement of specific energy loss at low transverse momenta ( p T 100 MeV/c).
The TPC, which is the main tracking detector of the ALICE central barrel, is based on a cylindrical gaseous chamber with radial and longitudinal dimensions of 85 cm < r < 247 cm and −250 cm < z < 250 cm, respectively. The TPC is read out by multi-wire proportional chambers (MWPC) with cathode pad readout, located at its endplates. With the measurement of drift time, the TPC provides three-dimensional space-point information for each charged track in pseudorapidity range |η| < 0.8 with up to 159 samples per track. In the TPC, the identification of charged particles is based on the measurement of the specific energy loss, which in pp collisions is performed with a resolution of 5.2% [17].
The TOF is a large-area array of multigap resistive plate chambers (MRPC), formed into a ∼4 m radius cylinder around the interaction point and covering the pseudorapidity region |η| < 0.9 with full-azimuth coverage. The time-of-flight is measured as the difference between the particle arrival time and the event collision time, enabling particle identification at intermediate transverse momenta, 0.5 p T 4 GeV/c. The arrival time is measured by the MRPCs with an intrinsic time resolution of 50 ps, while the event collision time is determined by combining the T0 detector measurement with the estimate using the particle arrival times at the TOF [22].

Event selection, classification and normalization
The analysed data were recorded using a minimum-bias trigger requiring signals in both V0A and V0C scintillators in coincidence with the arrival of the proton bunches from both directions. The background events produced outside the interaction region are rejected using the correlation between the SPD clusters and the tracklets reconstructed in SPD. The out-of-bunch pileup was rejected offline using the timing information from the V0 counter. The primary vertex was reconstructed either using global tracks (reconstructed using ITS and TPC information) or SPD tracklets (reconstructed using only the SPD information) with |z vtx | < 10 cm along the beam axis. Events with in-bunch pileup were removed if a second vertex was reconstructed within 8 mm of the primary vertex in the beam direction. The typical interaction rate of pp collisions in the 2016 data taking periods was around 120 kHz while beam-gas interactions occurred at a rate of 1.2 kHz.
In the analysis presented in this paper, we consider the event class INEL>0 with at least one charged particle produced in the pseudorapidity region |η| < 1, which corresponds to ∼75% of the total inelastic scattering crosssection [23]. To avoid auto-correlation biases [11,23], the events are classified using the total charge collected in the V0 detector (V0M amplitude), which scales linearly with the total number of the corresponding charged particles in its acceptance [24]. For each event class, the corresponding mean charged-particle multiplicity density d N ch /dη is measured at mid-rapidity (|η| < 0.5) as summarised in Table 1. 2.2 Identification of charged pions, kaons and protons In order to measure particle spectra in a wide p T range, several sub-analyses employing different detectors and particle identification (PID) techniques were performed and combined. As a result, the combined spectra cover transverse momenta ranges from 0.1/0.2/0.3 GeV/c to 20 GeV/c for π /K/p. The p T and (pseudo)rapidity ranges covered by each analysis for different particle species are summarized in Table 2.
At low p T , hadron spectra were measured by the ITS stand-alone (ITSsa) analysis. The dynamic range of the analogue readout of SDD and SSD allows for dE/dx measurements of highly ionizing particles, which otherwise do not reach the outer detectors. Hadron identification in the ITS is carried out by calculating the truncated mean of dE/dx and comparing it to the expected energy loss under different mass hypotheses. The difference between measured and expected dE/dx is then estimated in terms of the standard deviation σ and the particle mass hypothesis with the lowest score is assigned. This is feasible even for pp collisions with the highest multiplicities, as the number of charge clusters wrongly assigned to the reconstructed tracks is negligible. A detailed description of the method is provided in [11].
Hadrons at intermediate p T enter the fiducial volume of the TPC where they can be identified by measuring the charge generated in the gas. The truncated mean of dE/dx is calculated for the global tracks and compared to the expected energy loss under a given mass hypothesis. At low transverse momenta where the separation between different species is sufficiently large, tracks within three standard deviations from the expected dE/dx are assigned to a given hypothesis. In the regions where signals from several species overlap ( p T < 0.4 GeV/c for π , p T > 0.45 GeV/c for K, and p T > 0.6 GeV/c for p), dE/dx is fit with two Gaussian distributions, one to describe the signal and the other to describe the tail of the overlapping species. The fit of the overlapping species is then integrated in the signal region and subtracted from the signal [11].
In the p T region where the statistical unfolding of the TPC signal becomes unfeasible, particle identification is performed using the time-of-flight measurements. The results presented in this paper were obtained by combining the particle spectra estimated with two separate TOF analyses, taking into account the non-common part of the respective system-atic uncertainties. In the "TOF template fits", the PID is based on a statistical unfolding method, where the distribution of the difference between measured and expected time-of-flight (i.e. t) is fitted with templates for pions, kaons and protons in each p T and multiplicity bin [25]. An additional template is needed to take into account the background due to wrongly associated tracks with hits in the TOF detector. The template for each particle is built from data, considering the measured TOF time response function (Gaussian with an additional exponential tail for larger arrival times). The fits are repeated separately for each particle hypothesis in |y| < 0.5. In contrast to this, in the "TOF fits" analysis, the velocity β distribution is simultaneously fitted for all three particle types. For this purpose, four analytic functions, three for π , K and p, and one for mismatches, are employed. The analysis is performed in two narrow pseudorapidity slices (|η| < 0.2 and 0.2 < |η| < 0.4) and in momentum bins, which are then unfolded to transverse momenta. The corresponding rapidity interval is determined under the assumption of a flat d N ch /dη distribution in the aforementioned pseudorapidity bins [26].
Charged kaons can also be identified via the kink decay topology, where a charged particle decays into a charged and a neutral daughter (K ± → μ ± ν μ or K ± → π ± π 0 ). This secondary vertex where both decaying particle and the charged decay product have the same charge is reconstructed inside the ALICE TPC detector. This technique extends the charged kaon identification up to 6 GeV/c on a track-by-track basis. The algorithm for selecting kaons via their kink decay is used in a fiducial volume inside the TPC corresponding to a radial distance of 120 < R < 210 cm. This selection allows for an adequate number of TPC clusters to be associated with the decaying particle and its products. The track of the decaying particle is required to fulfil all the criteria of the global tracks except for the minimum number of clusters, which in this case is 30.
The topological selection of the kaon candidates and their separation from the pion decays (π ± → μ ± ν μ ) is based on the two-body decay kinematics. The transverse momentum of the decay product with respect to the decaying particle's direction (q T ) has an upper limit of 236 MeV/c for kaons and 30 MeV/c for pions for the two-body decay to μ ± ν μ . Similarly, for kaons decaying to pions, this limit is 205 MeV/c. Thus, a selection of q T < 120 MeV/c rejects the majority (85%) of pion decays. In addition, the angle between the mother and the daughter tracks is selected to be above the maximum allowed decay angle for pions and below the maximum allowed decay angle for kaons [27]. The invariant mass for the decay μ ± ν μ , M μν is calculated by assuming the daughter track to be a muon and the undetected track to be a neutrino. These selection criteria lead to a kaon sample with a purity of 97%.
The strategy employed to measure particle production in the region of the relativistic rise of the TPC was reported in [28]. The dE/dx signal in the relativistic rise (3 < βγ = p m < 1000) follows the functional form ln(βγ ). In addition to the logarithmic growth, the separation in number of standard deviations between pions and protons, pions and kaons, and kaons and protons as a function of momentum is nearly constant, which allows identification of charged pions, kaons, and (anti)protons with a statistical deconvolution approach from p T ≈ 2 − 3 GeV/c up to p T = 20 GeV/c. In order to describe the TPC response in the relativistic rise, clean external samples of secondary particles were used to parametrize the Bethe-Bloch and resolution curves. These correspond to pions (protons) from weak decays: K 0 S → π + + π − ( → p + π − ) and electrons from photon conversion. Moreover, primary pions measured with the TOF detector were used. The parametrization is done as a function of pseudorapidity. For short (long) tracks, i.e tracks within |η| < 0.2 (0.6 < |η| < 0.8), the resolution for protons is ≈ 6.2% (≈ 5.4%), while for pions it is ≈ 5.4% (≈ 5.0%). To extract the fraction of charged pions, kaons, and protons in the four different pseudorapidity intervals (|η| < 0.2, 0.2 < |η| < 0.4, 0.4 < |η| < 0.6, and 0.6 < |η| < 0.8) a 4-Gaussian fit (three for π , K, p and one to remove the unwanted electron contribution) to the dE/dx distribution in momentum bins is performed. The only free parameter in each of the Gaussian functions is the normalization, while the dE/dx and σ dE/dx are obtained and fixed using the Bethe-Bloch and resolution parametrizations, respectively. A weighted average of the four different measurements is calculated to obtain the particle fractions in |η| < 0.8. The yields are obtained by multiplying the particle fractions by the measured unidentified charged particle spectrum.

Corrections and normalization
The raw particle distributions are normalized to the total number of events analysed 1 in each multiplicity class. To obtain the p T distributions of primary π , K, and p, the raw particle distributions obtained from the different PID approaches need to be also corrected for the detector efficiency and acceptance, the ITS-TPC, and TPC-TOF matching efficiency, the PID efficiency, the trigger efficiency and the contamination from secondary particles.
Secondary particles are either produced in weak decays or from the interaction of particles with the detector material. The estimation of secondary particle contribution is based on the Monte Carlo (MC) templates of the distance of closest approach of the track to the primary vertex in the transverse plane with respect to the beam axis (DCA xy ), as carried out in previous works [11,25]. The DCA xy distributions of the tracks in data are fitted with three MC templates corresponding to the expected shapes of primary particles, secondaries from material and secondaries from weak decays to obtain the correct fraction of primary particles in the data. This procedure is repeated in each p T and multiplicity bin and thus takes into account the possible differences in the feed-down corrections due to the change in the abundances and spectral shapes of the weakly decaying particles. The contamination is different in each PID analysis due to different track selection criteria and PID techniques and hence it is estimated separately for each analysis. The contribution of secondary particles was found to be significant for π (up to 2%) and p (up to 15%) whereas the contribution for K is negligible.
The spectra are corrected for the detector acceptance and track reconstruction efficiencies based on a simulation using the Pythia8 (Monash-2013 tune) Monte Carlo event generator [29] and particle propagation through the full ALICE geometry using GEANT3 [30]. In this simulation, tracks are reconstructed using the same algorithms as for the data. The detector acceptance and reconstruction efficiencies are found to be independent of charged-particle multiplicity and thus the multiplicity-integrated values are used in all multiplicity classes. As GEANT3 does not fully describe the interaction of low-momentum p and K − with the detector material, an additional correction factor to the efficiency for these two particles is estimated with GEANT4 [31] and FLUKA [32], respectively, where the interaction processes are known to be better reproduced [25]. Additional corrections to the efficiency are applied when TPC or TOF information is used to take into account the track matching between ITS and TPC, and between TPC and TOF.
Signal losses due to the trigger selection are extracted from Pythia8 (Monash-2013 tune) MC simulation as performed in [23]. The correction is found to be 17-18% at low p T in the V0M class X (the lowest multiplicity), and reduces to ∼5%, ∼2% in classes IX and VIII, respectively. The correction is negligible in higher multiplicity pp collisions and for p T 4 GeV/c in all multiplicity bins except in class X. In the latter, the correction reaches ∼2% at p T = 7 GeV/c. Finally, an additional correction is applied to pass from triggered INEL>0 to true INEL>0 events, i.e. events with at least one primary charged particle in |η true | < 1 and with the primary vertex in the region |V true z | < 10 cm. The correction is independent of particle species and is found to be negligible from V0M I (the highest multiplicity) to V0M VI, while it ranges from 1% in class VII to 11% in class X. The correction is about 8% for multiplicity-integrated INEL>0 events.

Systematic uncertainties
The systematic uncertainties are divided into two categories, those common to all analyses and those which are analysis specific. The common systematic uncertainties are those due to tracking, which includes track quality criteria and the p T -dependent ITS-TPC matching efficiency (except for the ITSsa analysis), the TPC-TOF matching efficiency (for TPC-TOF and TOF analyses), and the signal loss correction. In addition, the systematic uncertainty related to the effect of the material budget on the global tracking ( p T dependent) is also added. The uncertainties on global tracking and TPC-TOF matching due to material budget are calculated by varying the material budget in the simulation by ±5%. The uncertainty related to the hadronic interaction cross section in the detector material is estimated using GEANT4 [31] and FLUKA [32] transport codes. Finally, an additional systematic uncertainty of 2% is added to account for possible mul-tiplicity dependence of track reconstruction efficiency and signal loss correction calculated from a MC simulation. All common sources of systematic uncertainties are summarised in Table 3. In the same table, the individual analysis systematic uncertainties are also listed for each particle species.
The estimation of the systematic uncertainties for the ITSsa analysis is described in detail in [11,25]. The ITSsa tracking uncertainties are estimated by varying the main criteria for the track selection, namely those on the DCA xy , on the χ 2 of the track, and on the number of clusters required in the ITS layers. The uncertainty related to the particle identification is calculated by using a Bayesian technique and comparing the results obtained with the standard nσ method as already performed in [33]. Due to the Lorentz force, the positions of ITS clusters are shifted depending on the magnetic field polarity, giving rise to a 3% uncertainty. Finally, the energy-independent uncertainty related to the ITS material budget is estimated with a simulation of pp collisions at √ s = 900 GeV by varying the material budget of the ITS by ±7.5% [34]. For the TPC-TOF fits analysis at low p T (below 500, 600, and 800 MeV/c for π , K, and p, respectively), the systematic uncertainty associated with the PID technique is calculated by integrating the measured dE/dx of charged tracks in the ranges of ±3.5σ and ±2.5σ , where σ represents one standard deviation from the dE/dx under given mass hypothesis. At higher p T values, where only the time-offlight information is used, the associated uncertainties are calculated by simultaneously varying the width and tail parameters by 10%. An additional uncertainty is calculated by fixing the central values of the fit functions to the β calculated for each particle species in a given momentum range. This was found to be the dominant source of systematic uncertainty for π and K at the highest p T values ( 2.5 GeV/c). For the TOF template fits analysis, PID uncertainties are estimated by simultaneously varying the spread and tail slope of each t template by 10%. In addition to this, for both the TPC-TOF and TOF template fits analyses, systematic uncertainties associated with tracking are calculated by varying the track selection criteria: the number of crossed rows in the TPC, the distance of closest approach in beam and transverse directions, and the quality of the global track fit χ 2 . For the kink analysis the sources of systematic uncertainties are: the kink vertex finding efficiency (3% constant in p T ), the kink PID efficiency (calculated by taking into account the position of the kink vertex, the number of TPC clusters of the decaying particle track and the q T of the decay product), and the uncertainty related to the purity of the selected sample. The contamination due to the random association of tracks wrongly attributed to kaon decays is of the order of 2.3% at low transverse momenta and reaches the value of 3.4% above Eur. Phys. J. C (2020) 80:693 and resolution curves. To quantify this uncertainty, the variations of the Bethe-Bloch resolution parametrizations with respect to the measured dE/dx (σ dE/dx ) are used to vary the values of the mean and σ in the 4-Gaussian fit [28]. The largest relative deviation between the nominal particle ratios and the ones obtained after the variations are assigned as a systematic uncertainty.

Results and discussion
The p T -differential spectra of π , K, and p measured as a function of the charged-particle multiplicity density in pp collisions at √ s = 13 TeV are shown in Fig. 1. For each V0M class, charged-particle multiplicity density has been measured in the central region (|η| < 0.5), as summarized in Table 1. The bottom panels in Fig. 1 show spectral ratios to the INEL>0 (sum of all V0M classes) class. We observe that the measured p T spectra become harder with increasing d N ch /dη , and the effect is more pronounced for protons. The hardening of the inclusive charged-hadron spectra with d N ch /dη has been also recently reported in [35], where different MPI models were shown to describe such effect. On the other hand, the mass dependence of spectral shape modifications is also observed in Pb-Pb collisions at √ s NN = 2.76 TeV [28], where it is usually associated with the hydrodynamical evolution of the system. At higher p T ( 8 GeV/c), we find that slopes of particle spectra become independent of the multiplicity class considered, as expected from pQCD calculations [36].
The p T -differential K/π and p/π ratios as a function of d N ch /dη measured at low, intermediate and high transverse momenta are shown in Fig. 2 together with those measured in pp collisions at √ s = 7 TeV [11] and predictions from several MC generators for pp collisions at √ s = 13 TeV. The measured K/π ratio shows no evident sign of evolution with multiplicity in all p T ranges considered, while the p/π ratio shows depletion at low p T , an increase at intermediate p T , and constant behavior at high p T . In addition, the measured K/π and p/π ratios are consistent between the two center-of-mass energies [11].
For MC predictions, the event classification is based on the number of charged tracks simulated at forward and backward pseudorapidities covered by the V0 detector, in a way similar to the way the event classification is done for the data. The mean charged-particle multiplicity density is then calculated in the central pseudorapidity region, |η| < 0.5. HERWIG 7 [37,38], where a clustering approach is used for hadronization, provides a good description of the evolution of the K/π and p/π ratios with d N ch /dη in the low and intermediate p T ranges and is consistent with the measured ratios within 1-2 standard deviations. Pythia8 [39] without color reconnection (CR) predicts no evolution of K/π and p/π ratios. The CR scheme, which has been shown to capture the modifications of the baryon-to-meson ratios [3], provides only a qualitative description of the evolution of the p/π ratio with d N ch /dη and overestimates the absolute values of the ratio at low and high p T . The implementation of color ropes [5,40,41] in Pythia8, which results in higher effective string tension and thus enhances strange-and di-quark production, provides a qualitative description K/π (p/π ) ratio only at low (intermediate) p T and overestimates the p/π ratio at low p T . This could be understood considering that larger effective string tension is mostly translated to hadronic mass and thus feeds down the low p T part of the spectrum.
In large collision systems such as Pb-Pb, multiplicitydependent modifications of hadron p T spectra can be interpreted as the hydrodynamical radial expansion of the system and studied in the context of the Boltzmann-Gibbs Blast-Wave model [42]. In this model, a thermalized medium expands radially and undergoes an instantaneous kinematic freeze-out. The average expansion velocity β T , the kinetic freeze-out temperature T kin , and the velocity profile exponent n can be extracted from simultaneous model fits to hadron spectra. As the trends observed in the evolution of particle spectra measured in pp collisions are highly reminiscent to those in p-Pb and Pb-Pb, it is interesting to check whether the Blast-Wave model can be extended to describe pp collisions. Such study has been previously reported in [11], where pp, p-Pb, and Pb-Pb collisions at √ s NN = 7, 5.02, and 2.76 TeV were considered. Now, for the first time, we can study the evolution of β T , T kin and n in pp collisions as a function of the collision energy.
At low transverse momenta ( p T 500 MeV/c), the dominant mechanism of π production is from resonance decays. To account for this in the Blast-Wave model fits, spectral measurements of all strongly decaying hadrons are required. Alternatively, one can choose to omit the lowp T pions. Noting that there is a strong dependence of Blast-Wave parameters on the fitting range [25], it is important to consider the same p T range in the fitting procedure in order to obtain a consistent comparison between different colliding systems. The comparison of the β T -T kin correlations measured in different systems and center-of-mass energies is shown in Fig. 3. In this paper we consider three different approaches to the Blast-Wave model fits to particle spectra measured in pp collisions at √ s = 13 TeV: (a) traditional fits as done in [10,11,25], where π , K, and p spectra are fitted and resonance feed-down is neglected (represented by full markers in Fig. 3), (b) simultaneously fitting K, p, and spectra [23] noting that are not significantly affected by resonance decays (represented by shaded ellipses in Fig. 3), and (c) a method proposed in [43,44], where the resonance feed-down is calculated before the Cooper-Frye freeze-out using a statistical hadronization model (represented by empty circles in Fig. 3). We find that the β T -T kin correlation in pp Fig. 1 Transverse momentum spectra of π , K, and p for different multiplicity event classes. Spectra are scaled by powers of 2 in order to improve visibility. The corresponding ratios to INEL>0 spectra are shown in the bottom panels Table 3 Sources of the relative systematic uncertainties of the p Tdifferential yields of π , K, and p. The uncertainties are split into two categories, the common and the individual-analysis specific for low, intermediate and high p T . Numbers in parenthesis in the p column refer to p uncertainties. In the last rows, the maximum (among multiplicity classes) total systematic uncertainty is reported Signal-loss correction 0. collisions at √ s = 13 TeV follows similar trends as seen at lower energies. When 's are considered instead of pions, the trends seen in β T -T kin correlation do not change significantly and only at highest multiplicities we find a larger T kin . On the other hand, when a proper treatment of resonance decays is used, we find a significantly lower T kin of around 135 MeV at the lowest multiplicities, which then grows with increasing d N ch /dη and approaches the pseudocritical QCD temperature T c = 156 ± 1.5 MeV [45] at the highest multiplicity pp collisions. In addition, the evolution of β T , T kin , and n with d N ch /dη is shown in Fig. 4 for different colliding systems. From the lowest multiplicities, T kin grows with β T until it saturates at around 180 MeV. At larger multiplicities ( d N ch /dη 16), T kin decreases and becomes similar to that measured in p-Pb collisions at √ s NN = 5.02 TeV, suggesting that the system decouples at lower temperature and thus is longer-lived. open, and shaded bands represent statistical, total systematic, and multiplicity uncorrelated systematic uncertainties, respectively. Numbers in the parenthesis in different panels represent different scale factors for data and MC predictions for better readability cating that small systems become more explosive at larger multiplicities. In contrast to this, β T measured in Pb-Pb collisions is lower than that in smaller systems for the common d N ch /dη range, see Fig. 4. This indicates that the size of the colliding system might have significant effects on the final state particle dynamics. This is also reflected in the expansion velocity profile power n shown in Fig. 4: in pp and p-Pb collisions, large n suggests high pressure gradients which lead to larger β T , while in Pb-Pb collisions, n∼1 could be interpreted as lower pressure gradient and thus smaller expansion velocity [46].
Previous studies on hadron production as a function of multiplicity have reported the factorization of p T -integrated particle yields with d N ch /dη [11], which extends across different colliding systems and collision energies. Now for the first time we can isolate the center-of-mass energy dependence of this scaling for π , K, and p in pp collisions. The p Tintegrated particle yields (d N/dy) and average transverse momenta ( p T ) are calculated by integrating the measured transverse momentum spectra and using the Lévy-Tsallis parametrization [47][48][49] to extrapolate to the low p T regions not covered by the measurements. The extrapolated fractions of the yields at low p T are 8% (10%) for π , 6% (13%) for K, and 7% (20%) for p for the highest (lowest) multiplicities. For systematic uncertainties on the extrapolation, Bylinkin, Bose-Einstein, Fermi-Dirac, m T -exponential and Hagedorn functions are used to fit particle spectra. The largest system- Fig. 3 Correlation of kinetic freeze-out temperature T kin and average expansion velocity β T , extracted from simultaneous Blast-Wave fits to π , K, and p spectra measured in pp, p-Pb, and Pb-Pb collisions. Contours represent 1σ uncertainty. The shaded ellipses represent the β T -T kin correlation calculated from Blast-Wave fit to K, p, and spectra [23] measured in pp collisions at √ s = 13 TeV. The empty circles represent Blast-Wave fits with resonance decays [44]. References from [10][11][12]44] atic uncertainties on d N/dy ( p T ) related to the extrapolation procedure are found to be 2% (2%), 2% (2%), and 3% (2%) for π , K, and p at low-multiplicity classes and become smaller at higher multiplicities.
The statistical uncertainties of d N/dy and p T are calculated by coherently shifting the central values of each spectra Eur. Phys. J. C (2020) 80:693 Fig. 4 Evolution of β T , T kin , and n with d N ch /dη . β T , T kin and n are extracted from Blast-Wave fits to π , K, and p p T spectra measured in pp, p-Pb, and Pb-Pb collisions at different √ s. The resonance feed-down contribution is neglected Fig. 5 Integrated K/π (top) and p/π (bottom) yield ratios as a function of charged-particle multiplicity density measured in pp, p-Pb, and Pb-Pb collisions at different center-of-mass energies. Empty (shaded) boxes represent total (multiplicity uncorrelated) systematic uncertainties. Black lines represent predictions from different MC generators for pp collisions at √ s = 13 TeV. References from [11,25,33,50] point by a fraction of its statistical uncertainty. The fraction is randomly drawn from Gaussian distribution and new values of integrated yields and mean transverse momenta are calculated. The procedure is repeated 1000 times to calculate the standard deviations of d N/dy and p T , which are then used as the statistical uncertainties. To estimate the systematic uncertainty on the integrated yields, the spectra points are moved to maximal/minimal values allowed by their respective systematic uncertainties before repeating the fit procedure. For p T , each point of the spectra is shifted to the upper/lower edge of the corresponding p T bin to obtain the hardest/softest particle distribution. The largest differences to the nominal yield and p T values are combined with the extrapolation uncertainties to calculate the total systematic uncertainties. The kaon-and proton-to-pion integrated yield ratios measured in pp collisions at √ s = 13 TeV are found to be in a good agreement within systematic uncertainties with those measured in pp, p-Pb, and Pb-Pb collisions at √ s NN = 7, 5.02, and 2.76 TeV, respectively, as shown in Fig. 5. In addition, with the availability of (multi)strange hadron yields [23] we can study the relative abundances of hyperons to pions, and the results are shown in Fig. 6. We find that the (multi)strange hadron-to-pion ratios measured in pp collisions at √ s = 13 TeV are in good agreement to those measured at √ s = 7 TeV and similar d N ch /dη . This indicates that hadrochemistry at LHC energies scales with charged-particle multiplicity density in a uniform way, despite the colliding system or collision energy.
The description of hadron-to-pion ratio factorization with multiplicity at lower center-of-mass energies in MC generators has been previously shown to be qualitative at best [11]. In fact, both Pythia8 with color reconnection and HER-WIG 7 [37,38] predict no evolution of the ratios with d N ch /dη . In this paper, we consider more recent versions of the two MC generators. In particular, the hadronization in Pythia8 now considers overlapping color strings, which form color ropes with a larger effective string tension and are then allowed to interact with each other [41]. On the other hand, hadronization in HERWIG 7 now includes baryonic ropes -a reconnection scheme that enhances the probability of partons forming a baryon [37]. We find that both  [6,11,12,23] Pythia8 and HERWIG 7 predict the enhancement of strange baryons which is more pronounced for hadrons with a larger strangeness content as shown in Fig. 6. The largest quantitative differences are seen for /π ratio at the lowest multiplicity in pp collisions. The /π ratios are in a better agreement with Pythia8 with color ropes, while HERWIG 7 shows a large deviation from the data at low d N ch /dη . Finally, /π ratios are well described by HERWIG 7, while Pythia8 with color ropes predicts an increasing trend in the whole multiplicity range available and overestimates the ratio at the highest multiplicities. Overall, the agreement between MC generators and measured hadron-to-pion ratios become worse for particles with a larger strangeness content. This might point to the need of a further refinement of MC generator tuning, as similar trends are already observed for e + e − data [51].
The integrated K/π yield ratio shown in Fig. 5 at high multiplicity pp collisions are captured by Pythia8 ropes and HERWIG 7, but the latter predicts a peak-like structure at low d N ch /dη which is not observed in the data. The predictions from Pythia8 Monash tune are inconsistent with the measured K/π ratios in pp collisions at √ s = 13 TeV, whether color reconnection is considered or not. The quantitative description of p/π ratio is given only by HERWIG 7, while all considered versions of Pythia8 overpredict the data. Moreover, Pythia8 with color ropes predicts an increase of the p/π ratio with d N ch /dη , which could be attributed to the enhanced production of strange-and di-quark in the rope fragmentation. Overall, we conclude that none of the models considered provide a consistent description of the data.
The average transverse momenta of identified particles are found to increase with multiplicity in pp collisions at √ s = 7 and 13 TeV as shown in Fig. 7. A clear mass ordering is observed among the particle species considered, where protons have the largest p T . Similar observations have been previously reported in pp [52] and p-Pb [10] collisions at lower energies and for strange hadrons in pp collisions at √ s = 13 TeV [23]. The solid red line in Fig. 7 represents a fit of the form a − b/(c − d N ch /dη ) to the √ s = 13 TeV data, which is then used for a better comparison of p T between the two center-of-mass energies, see lower panels of the same figure. We find a small hint of an increase with √ s for similar multiplicities for π , while the p T of protons is similar at the two center-of-mass energies. Note that similar observations have been already reported in [23], where spectra of K 0 s were found to become harder with √ s at similar multiplicities. In addition, we find that Pythia8 Monash tune with color reconnection, HERWIG 7, and Pythia8 with ropes give a very good description of π p T evolution with d N ch /dη . This is expected as pions are the most abundant particles produced in collisions, and the three generators are tuned to explicitly to describe the p T of charged-particles. On the other hand, we observe that the p T of K and p are well described only by HERWIG 7, while Pythia8 with rope implementation underestimates the p T in the whole d N ch /dη range considered. This could be understood considering that the additional energy available during the rope fragmentation predominantly enhances the production of heavier hadrons at low p T .

Summary
We have studied π , K, and p production as a function of multiplicity in pp collisions at √ s = 13 TeV. To avoid autocorrelation biases, the event classification has been based on multiplicity measurements at forward (backward) pseudorapidity, while event activity d N ch /dη has been correspondingly estimated at central pseudorapidity region, |η| < 0.5. We find that hadron p T spectra become harder with multiplicity, and the effect is more pronounced for heavier particles. The hardening of the spectra is predicted by Pythia8 with rope hadronization, Pythia8 Monash with color recon- nection, and HERWIG7 MC generators. In addition, all three generators provide a quantitative description of π p T , while K and p are described qualitatively only by HERWIG7. At high p T ( 8 GeV/c) we find that spectral shapes become independent of d N ch /dη as predicted by pQCD calculations [36]. The measured p T -differential K/π ratios show no evolution with multiplicity in the p T range considered. In contrast to this, a depletion (enhancement, saturation) is visible for the p/π ratios at low (intermediate, high) p T . In addition, we find that the ratios measured in pp collisions at √ s = 13 TeV are consistent with those measured at √ s = 7 TeV. The saturation at high p T is captured by Pythia8 Monash tunes, while HERWIG 7 and Pythia8 with color ropes show signs of enhancement. While some of the most common MC generators capture the trends seen in the p T -differential K/π and p/π ratios, it is interesting to see that none of them provides a consistent description of the data and predict the absolute values of the ratios at high p T .
The study of hadron p T spectra in the context of the Blast-Wave model reveals that the kinetic freeze-out temperature T kin , average expansion velocity β T , and the velocity profile exponent n show little or no dependence on the centerof-mass energy and are consistent within uncertainties with those extracted from particle spectra measured in pp colli- . On the other hand, we observe a strong dependence of the extracted parameters on d N ch /dη .
The p T -integrated hadron-to-pion ratios as a function of multiplicity show no center-of-mass energy dependence and the measurement in pp collisions at √ s = 13 TeV are compatible to those in pp, p-Pb, and Pb-Pb collisions at √ s NN = 7, 5.02, and 2.76 TeV, respectively. This suggests that, at the LHC energies, the chemical composition of primary hadrons scales with charged-particle multiplicity density in a uniform way, despite the colliding system and collision energy. Comparisons of the integrated hadron-topion ratios to the predictions from MC generators show that Pythia8 with color ropes provides the best description of (multi)strange hadrons, but overestimates the measured p/π ratio. HERWIG7 also captures the evolution of the ratios with d N ch /dη , but underestimates the absolute values of /π and /π . Overall, none of the generators are able to provide a consistent quantitative description of the measured hadron-to-pion ratios.
Acknowledgements The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centres and the Worldwide LHC Computing Grid (WLCG) collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the ALICE detector: