Multiplicity dependence of (multi-)strange hadron production in proton-proton collisions at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}$$\end{document} = 13 TeV

The production rates and the transverse momentum distribution of strange hadrons at mid-rapidity (y<0.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left| y\right| < 0.5$$\end{document}) are measured in proton-proton collisions at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}$$\end{document} = 13 TeV as a function of the charged particle multiplicity, using the ALICE detector at the LHC. The production rates of KS0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{K}^{0}_{S}$$\end{document}, Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document}, Ξ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Xi $$\end{document}, and Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} increase with the multiplicity faster than what is reported for inclusive charged particles. The increase is found to be more pronounced for hadrons with a larger strangeness content. Possible auto-correlations between the charged particles and the strange hadrons are evaluated by measuring the event-activity with charged particle multiplicity estimators covering different pseudorapidity regions. When comparing to lower energy results, the yields of strange hadrons are found to depend only on the mid-rapidity charged particle multiplicity. Several features of the data are reproduced qualitatively by general purpose QCD Monte Carlo models that take into account the effect of densely-packed QCD strings in high multiplicity collisions. However, none of the tested models reproduce the data quantitatively. This work corroborates and extends the ALICE findings on strangeness production in proton-proton collisions at 7 TeV.


Introduction
The production rates of strange and multi-strange hadrons in high-energy hadronic interactions are important observables for the study of the properties of Quantum Chromodynamics (QCD) in the non-perturbative regime. In the simplest case, strange quarks (s) in proton-proton (pp) collisions can be produced from the excitation of the sea partons. Indeed, in the past decades a significant effort has been dedicated to the study of the actual strangeness content of the nuclear wave function [1]. In QCD-inspired Monte Carlo generators based on Parton Showers (PS) [2] the hard (perturbative) interactions are typically described at the Leading Order (LO). In this picture the s quark can be produced in the hard pare-mail: alice-publications@cern.ch tonic scattering via flavour creation and flavour excitation processes as well as in the subsequent shower evolution via gluon splitting. At low transverse momentum, ss pairs can be produced via non-perturbative processes, as described for instance in string fragmentation models, where the production of strangeness is suppressed with respect to light quark production due to the larger strange quark mass [3]. However, these models fail to quantitatively describe strangeness production in hadronic collisions [4][5][6].
An enhanced production of strange hadrons in heavy-ion collisions was suggested as a signature for the creation of a Quark-Gluon Plasma (QGP) [7,8]. The main argument in these early studies was that the mass of the strange quark is of the order of the QCD deconfinement temperature, allowing for thermal production in the deconfined medium. The lifetime of the QGP was then estimated to be comparable to the strangeness relaxation time in the plasma, leading to full equilibration. Strangeness enhancement in heavy-ion collisions was indeed observed at the SPS [9] and higher energies [10,11]. However, strangeness enhancement is no longer considered an unambiguous signature for deconfinement (see e.g. [12]). Strange hadron production in heavy-ion collisions is currently usually described in the framework of statistical-hadronisation (or thermal) models [13,14]. In central heavy-ion collisions, the yields of strange hadrons turn out to be consistent with the expectation from a grandcanonical ensemble, i.e. the production of strange hadrons is compatible with thermal equilibrium, characterised by a common temperature. On the other hand, the strange hadron yields in elementary collisions are suppressed with respect to the predictions of the (grand-canonical) thermal models. The suppression of the relative abundance of strange hadrons with respect to lighter flavours was suggested to be, at least partially, a consequence of the finite volume, which makes the application of a grand-canonical ensemble not valid in hadron-hadron and hadron-nucleus interactions (canonical suppression) [15][16][17]. However, this approach cannot explain the observed particle abundances if the same volume is assumed for both strange and non-strange hadrons [18] and does not describe the system size dependence of the φ meson, a hidden-strange hadron [19,20].
The ALICE Collaboration recently reported an enhancement in the relative production of (multi-) strange hadrons as a function of multiplicity in pp collisions at √ s = 7 TeV [21] and in p-Pb collisions at √ s NN = 5.02 TeV [22,23]. In the case of p-Pb collisions, the yields of strange hadrons relative to pions reach values close to those observed in Pb-Pb collisions at full equilibrium. These are surprising observations, because thermal strangeness production was considered to be a defining feature of heavy-ion collisions, and because none of the commonly-used pp Monte Carlo models reproduced the existing data [3,21]. The mechanisms at the origin of this effect need to be understood, and then implemented in the state-of-the-art Monte Carlo generators [3].
In this paper, strangeness production in pp interactions is studied at the highest energy reached at the LHC, √ s = 13 TeV. We present the measurement of the yields and transverse momentum ( p T ) distributions of single-strange (K 0 S , , ) and multi-strange ( − , + , − , + ) particles at mid-rapidity, |y| < 0.5, with the ALICE detector [24]. The comparison of the present results with the former ones for pp and p-Pb interactions allows the investigation of the energy, multiplicity and system size dependence of strangeness production. Schematically, the multiplicity of a given pp event depends on (i) the number of Multiple Parton Interactions (MPI), (ii) the momentum transfer of those interactions, (iii) fluctuations in the fragmentation process. A systematic study of the biases induced by the choice of the multiplicity estimator along with the specific connections to the underlying MPI are also discussed in this paper. The paper is organised as follows. In Sect. 2 we discuss the data set and detectors used for the measurement; in Sect. 3 we describe the analysis techniques; in Sect. 4 we cover the evaluation of the systematic uncertainties; in Sect. 5 we present and discuss the results; in Sect. 6 we report our conclusions.

Experimental setup and data selection
A detailed description of the ALICE detector and its performance can be found in [24,25]. In this section, we briefly outline the main detectors used for the measurements presented in this paper. The ALICE apparatus comprises a central barrel used for vertex reconstruction, track reconstruction and charged-hadron identification, complemented by specialised forward detectors. The central barrel covers the pseudorapidity region |η| < 0.9. The main central-barrel tracking devices used for this analysis are the Inner Tracking System (ITS) and the Time-Projection Chamber (TPC), which are located inside a solenoidal magnet providing a 0.5 T magnetic field. The ITS is composed of six cylindrical layers of high-resolution silicon tracking detectors. The innermost layers consist of two arrays of hybrid Silicon Pixel Detectors (SPD), located at an average radial distance r of 3.9 and 7.6 cm from the beam axis and covering |η| < 2.0 and |η| < 1.4, respectively. The SPD is also used to reconstruct tracklets, short two-point track segments covering the pseudorapidity region |η| < 1.4. The outer layers of the ITS are composed of silicon strips and drift detectors, with the outermost layer having a radius r = 43 cm. The TPC is a large cylindrical drift detector of radial and longitudinal sizes of about 85 < r < 250 cm and −250 < z < 250 cm, respectively. It is segmented in radial "pad rows", providing up to 159 tracking points. It also provides charged-hadron identification information via specific ionisation energy loss in the gas filling the detector volume. The measurement of strange hadrons is based on "global tracks", reconstructed using information from the TPC as well as from the ITS, if the latter is available. Further outwards in radial direction from the beam-pipe and located at a radius of approximately 4 m, the Time of Flight (TOF) detector measures the time-of-flight of the particles. It is a large-area array of multigap resistive plate chambers with an intrinsic time resolution of 50 ps. The V0 detectors are two forward scintillator hodoscopes employed for triggering, background suppression and event-class determination. They are placed on either side of the interaction region at z = −0.9 m and z = 3.3 m, covering the pseudorapidity regions −3.7 < η < −1.7 and 2.8 < η < 5.1, respectively.
The data considered in the analysis presented in this paper were collected in 2015, at the beginning of Run 2 operations of the LHC at √ s = 13 TeV. The sample consists of 50 M events collected with a minimum bias trigger requiring a hit in both V0 scintillators in coincidence with the arrival of proton bunches from both directions. The interaction probability per single bunch crossing ranges between 2 and 14%.
The contamination from beam-induced background is removed offline by using the timing information in the V0 detectors and taking into account the correlation between tracklets and clusters in the SPD detector, as discussed in detail in [25]. The contamination from in-bunch pile-up events is removed offline excluding events with multiple vertices reconstructed in the SPD. Part of the data used in this paper were collected in periods in which the LHC collided "trains" of bunches each separated by 50 ns from its neighbours. In these beam conditions most of the ALICE detectors have a readout window wider than a single bunch spacing and are therefore sensitive to events produced in bunch crossings different from those triggering the collision. In particular, the SPD has a readout window of 300 ns. The drift speed in the TPC is about 2.5 cm/µs, which implies that events produced less than about 0.5 µs apart cannot be resolved. Pile-up events produced in different bunch crossings are removed exploiting multiplicity correlations in detectors having different readout windows.

Analysis details
The results are presented for primary strange hadrons. 1 The measurements reported here have been obtained for events having at least one charged particle produced with p T > 0 in the pseudorapidity interval |η| < 1 (INEL > 0), corresponding to about 75% of the total inelastic cross-section. In order to study the multiplicity dependence of strange and multi-strange hadrons, for each multiplicity estimator the sample is divided into event classes based on the total charge deposited in the V0 detectors (V0M amplitude) or on the number of tracklets in two pseudorapidity regions: |η| < 0.8 and 0.8 < |η| < 1.5. The event classes are summarised in Table 1. Since the measurement of strange hadrons is performed in the region |y| < 0.5, the usage of these three estimators allows one to evaluate potential biases on particle production, arising from measuring the multiplicity in a pseudorapidity region partially overlapping with the one of the reconstructed strange hadrons (N |η|<0.8 tracklets ), or disjoint from it (V0M and N 0.8<|η|<1.5 tracklets ). In particular, the effect of fluctuations can be expected to be stronger if the multiplicity estimator and the observable of interest are measured in the same pseudorapidity region. The usage of two different non-overlapping estimators allows the study of the effect of a rapidity gap between the multiplicity estimator and the measurement of interest.
The events used for the data analysis are required to have a reconstructed vertex in the fiducial region |z| < 10 cm. As mentioned in the previous section, events containing more than one distinct vertex are tagged as pile-up and discarded. For each event class and each multiplicity estimator, the average pseudorapidity density of primary charged-particles dN ch /dη is measured at mid-rapidity (|η| < 0.5) using the technique described in [27]. The dN ch /dη values, corrected for acceptance and efficiency as well as for contamination from secondary particles and combinatorial background, are listed in Table 1. When multiplicity event classes are selected outside the |η| < 0.8 region, the corresponding charged particle multiplicity at mid-rapidity is characterized by a large variance. In the case of the V0M estimator, the variance ranges between 30 and 70% of the mean dN ch /dη for the highest and lowest multiplicity classes, respectively.
The following decay channels are studied [28]: 1 A primary particle [26] is defined as a particle with a mean proper decay length cτ larger than 1 cm, which is either (a) produced directly in the interaction, or (b) from decays of particles with cτ smaller than 1 cm, excluding particles produced in interactions with material.  Table 2 Track, topological and candidate selection criteria applied to K 0 S , and candidates. DCA stands for "distance of closest approach", PV represents the "primary event vertex" and θ is the angle between the momentum vector of the reconstructed V 0 and the displacement vector between the decay and primary vertices. The selection on DCA between V 0 daughter tracks takes into account the corresponding experimental resolution In the following, we refer to the sum of particles and antiparticles, + , − + + and − + + , simply as , and .
The details of the analysis have been discussed in earlier ALICE publications [5,6,18,22]. The tracks retained in the analysis are required to cross at least 70 TPC readout pads out of a maximum of 159. Tracks are also required not to have large gaps in the number of expected tracking points in the radial direction. This is achieved by checking that the number of clusters expected based on the reconstructed trajectory and the measurements in neighbouring TPC pad rows do not differ by more than 20%.
Each decay product arising from V 0 (K 0 S , , ) and cascade ( − , + , − , + ) candidates is verified to lie within the fiducial tracking region |η| < 0.8. The specific energy loss (dE/dx) measured in the TPC, used for the particle identification (PID) of the decay products, is also requested to be compatible within 5σ with the one expected for the corresponding particle species' hypothesis. The dE/dx is evaluated as a truncated mean using the lowest 60% of the values out of a possible total of 159. This leads to a resolution of about 6%. A set of "geometrical" selections is applied in order to identify specific decay topologies (topological selection), improving the signal/background ratio. The topological variables used for V 0 s and cascades are described in detail in [6]. In addition, in order to reject the residual out-of-bunch pile-up background on the measured yields, it is requested that at least one of the tracks from the decay products of the (multi-)strange hadron under study is matched in either the ITS or the TOF detector. The selections used in this paper are summarised in Table 2 for the V 0 s and in Table 3 for the cascades.
Strange hadron candidates are required to be in the rapidity window |y| < 0.5. K 0 S ( ) candidates compatible with the alternative V 0 hypothesis are rejected if they lie within ±5 MeV/c 2 (±10 MeV/c 2 ) of the nominal (K 0 S ) mass. A similar selection is applied to the , where candidates compatible within ±8 MeV/c 2 of the nominal mass are rejected. The width of the rejected region was determined according to the invariant mass resolution of the corresponding competing signal. 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 signal extraction is performed as a function of p T . A preliminary fit is performed on the invariant mass distribution using a Gaussian plus a linear function describing the background. This allows for the extraction of the mean (μ) and width (σ G ) of the peak. A "peak" region is defined within ±6(4)σ G for V 0 s (cascades) with respect to μ for any measured p T bin. Adjacent background bands, covering a combined mass interval as wide as the peak region, are defined on both sides of that central region. The signal is then extracted with a bin counting procedure, subtracting counts in the background region from those of the signal region. Alternatively, the signal is extracted by Identity assumption for cascades and for daughter tracks fitting the background with a linear function extrapolated under the signal region. This procedure is used to compute the systematic uncertainty due to the signal extraction. Examples of the invariant mass peaks for all particles are shown in Fig. 1.
Only the turns out to be affected by a significant contamination from secondary particles, coming from the decay of charged and neutral . In order to estimate this contribution we use the measured − and + spectra, folded with a p T -binned 2D matrix describing the decay kinematics and secondary reconstruction efficiencies. The → π decay matrix is extracted from Monte Carlo (see below for the details on the generator settings). The fraction of secondary particles in the measured spectrum varies between 10 and 20%, depending on p T and multiplicity. Further details on the uncertainties characterising the feed-down contributions are provided in the next section.
The raw p T distributions are corrected for acceptance and efficiency using Monte Carlo simulated data. Events are generated using the PYTHIA 6.425, (Tune Perugia 2011) [29,30] event generator, and transported through a GEANT 3 [31] (v2-01-1) model of the detector. With respect to previous GEANT 3 versions, the adopted one contains a more realistic description of (anti)proton interactions. The quality of this description was cross-checked comparing to the results obtained with the state-of-the-art transport codes FLUKA [32,33] and GEANT 4.9.5 [34]. It was found that a correction factor < 5% is needed for the efficiency of , + , + for p T < 1 GeV/c, while the effect is negligible at higher p T . Events generated using PYTHIA 8.210 (tune Monash 2013) [35,36] and EPOS-LHC (CRMC package 1.5.4) [37] and transported in the same way are used for systematic studies, namely to compute the systematic uncertainties arising from the normalisation and from the closure of the correction procedure (details provided in the next sections).
The acceptance-times-efficiency changes with p T , saturating at a value of about 40%, 30%, 30% and 20% at p T 2, 3, 3 and 4 GeV/c for K 0 S , , and , respectively. These values include the losses due to the branching ratio. They are found to be independent of the multiplicity class within 2%, limited by the available Monte Carlo simulated data. The dependence of the efficiency on the generated p T distributions was checked for all particle species. It is found to be relevant only in the case of the , where large p T bins are used. This effect is removed by reweighting the Monte Carlo p T distribution with the measured one using an iterative procedure.
In order to compute p T and the p T -integrated production yields, the spectra are fitted with a Tsallis-Lévy [38] distribution to extrapolate in the unmeasured p T region. The

Systematic uncertainties
Several sources of systematic effects on the evaluation of the p T distributions were investigated. The main contributions for three representative p T values are summarised in Table 4 for the INEL > 0 data sample.
The stability of the signal extraction method was checked by varying the widths used to define the "signal" and "background" regions, expressed in terms of number of σ G as defined in Sect. 3. The raw counts were also extracted with a fitting procedure and compared to the standard ones computed by the bin counting technique. An uncertainty ranging between 0.2 and 3.5% depending on p T is assigned to the signal extraction of the V 0 s and cascades based on these checks.
The stability of the acceptance and efficiency corrections was verified by varying all track, candidate and topological selection criteria within ranges leading to a maximum variation of ±10% in the raw signal yield. The results were compared to those obtained with the default selection criteria (Sect. 3). Variations not compatible with statistical fluctuations (following the prescription in [39] with a 2σ threshold) are added to the systematic uncertainty.
The resulting uncertainty from topological and track selections (except TPC dE/dx) depends on p T and amounts at most to 4%, 5%, 4% and 6% for K 0 S , , and , respectively.
The TPC dE/dx selection is used to reduce the combinatorial background in the strange baryon invariant mass distribution. The uncertainty was evaluated varying the TPC dE/dx selection requirements between 4 and 7σ and was found to be at most 1% (3%) for ( and ). For the K 0 S the uncertainty due to the TPC dE/dx usage was evaluated by comparing results obtained adopting the default loose PID requirements (5σ ) with those obtained without applying any PID selection. The difference is found to be negligible (< 1%).
The systematic uncertainty for the competing decay rejection was investigated by removing entirely this condition for and . It resulted in a deviation on the p T spectra of at most 4% and 6%, respectively. For the K 0 S the systematic uncertainty was evaluated by changing the width of the competing rejected mass window between 3 and 5.5 MeV/c 2 and the corresponding deviation was found to be at most 1%. For the strange baryons, the systematic uncertainty related to the proper lifetime was computed by varying the selection requirements between 2.5 and 5 cτ . The variation range for the K 0 S was set to 5-15 cτ . The statistically significant deviations were found to be at most 3% for the and negligible (< 1%) for all other particles.
An uncertainty related to the absorption in the detector material was assigned to the anti-baryons, mostly due to the interactions of anti-proton daughters. It was estimated on the basis of the comparison of the different transport codes mentioned in Sect. 3. The uncertainty on the absorption cross section for baryons and K 0 S was found to be negligible. Furthermore an additional 2% uncertainty is added to account for possible variations of the tracking efficiency with multiplicity (Sect. 3).
The uncertainty due to approximations in the description of the detector material was estimated with a Monte Carlo simulation where the material budget was varied within its uncertainty [25]. The assigned systematic uncertainty ranges between 8% at low p T to about 1% at high p T .
The p T spectrum is affected by an uncertainty coming from the feed-down correction, due to the uncertainties on the measured spectrum and on the multiplicity dependence of the feed-down fraction. Furthermore the contribution from neutral 0 was taken into account by assuming ± / 0 = 1 or using the ratio provided by the Monte Carlo (using the reference PYTHIA 6 sample described in the previous section). The difference between these two estimates was taken into account in the calculation of the total uncertainty due to the feed-down correction, which ranges from 2% to 4% depending on p T and multiplicity.
The systematic uncertainty due to the out-of-bunch pileup rejection was evaluated by changing the matching scheme with the relevant detectors, considering the following configurations: matching of at least one decay track with the ITS (TOF) detector below (above) 2 GeV/c of the reconstructed (multi-) strange hadron; ITS matching of at least one decay track in the full p T range. Half of the maximum difference between these configurations and the standard selection was taken as the systematic uncertainty, which was found to increase with transverse momentum and to saturate at high p T , reaching a maximum value of 2.4% (3%) for V 0 s (cascades).
The effect of a possible residual contamination from inbunch pile-up events was estimated varying the pile-up rejection criteria and dividing the data sample in three groups with an average interaction probability per bunch crossing of 3%, 6%, 13%, respectively. The resulting uncertainty is larger at low multiplicity, and ranges between 1% (3%) for the K 0 S to 2% (3%) for the baryons in high-(low-)multiplicity events.
Several additional consistency checks have been performed which are described in the following. The analysis has been repeated separately for events with positive and negative z-vertex position, as well as considering candidates reconstructed in positive and negative rapidity windows. The    resulting p T -spectra were found to be statistically compatible with the standard analysis.
In order to ensure that the estimated systematic uncertainties are not affected by statistical fluctuations, two checks were performed. First of all, the threshold used to consider a variation statistically significant was varied between one and three standard deviations. Then, the analysis was repeated with wider p T bins. These checks showed that statistical fluc-      Transverse momentum distribution of K 0 S , , , and , for multiplicity classes selected using the tracklets in 0.8 < |η| < 1.5. Statistical and total systematic uncertainties are shown by error bars and boxes, respectively. In the bottom panels ratios of multiplicity dependent spectra to INEL > 0 are shown. The systematic uncertainties on the ratios are obtained by considering only contributions uncorrelated across multiplicity. The spectra are scaled by different factors to improve the visibility. The dashed curves represent Tsallis-Lévy fits to the measured spectra p T spectra well for all the examined strange hadrons over the measured p T range.
The uncertainty on the extrapolated fraction was estimated repeating the fit to the spectra with five alternative functions (Blast-Wave, Boltzmann, Bose-Einstein, m T -exponential, Fermi-Dirac). Since these alternative functions do not, in general, describe the full p T -distribution, the fit range was limited to a small number of data points in order to obtain a S , , , and in multiplicity event classes selected according to different estimators (see text for details). Statistical and systematic uncertainties are shown by error bars and empty boxes, respectively. Shadowed boxes represent uncertainties uncorrelated across multiplicity good description of the fitted part of the spectrum. This procedure was repeated separately for the low and for the high p T extrapolation, with the final uncertainties being dominated by the low p T contribution. The reliability of the extrapolation uncertainty estimate was checked using a simple linear extrapolation to p T = 0 as an extreme case and in a full Monte Carlo closure test where the EPOS model was used with data and PYTHIA was used to estimate the corrections. The resulting uncertainty on the integrated yields is around 2.5% for the and ranges between 3% (4%) at high multiplicity to 19% (12%) at low multiplicity for the ( ). The extrapolation uncertainty on p T is ∼2% for the and ranges between 2% (3%) at high multiplicity to 12% (7%) at low multiplicity for the ( ).
The main focus of this paper is the study of the multiplicity dependence of strangeness production. In this light, the different systematic uncertainties can be categorised in the following way: 1. Fully uncorrelated uncertainties: the change in the data is completely uncorrelated across multiplicity classes. These sources are rare, as most systematic effects have a smooth evolution with multiplicity. 2. Fully correlated uncertainties: lead to a correlated shift of the data in the same direction, independently of the multiplicity class being studied. The common part of this shift has to be considered separately, since it does not affect the shape of the multiplicity-dependent observable.
3. Anti-correlated uncertainties: the effect is opposite in low and high multiplicity events.
In the following we quote separately uncertainties that affect the trends as a function of multiplicity and those that lead to a constant fractional shift across all multiplicity classes. The effect of every systematic variation was evaluated simultaneously in each multiplicity class and in minimum bias events to separate the fully correlated part of the uncertainty.
Sources leading to a global, fully-correlated shift in the spectra were subtracted from the total uncertainty while the remaining contribution could in principle belong to any of the three aforementioned categories. However, since different (independent) sources are combined in the final uncertainties, it is a reasonable assumption to consider them as uncorrelated. For the figures shown in Sect. 5, total uncertainties are shown as boxes, while shadowed boxes represent uncertainties uncorrelated across different multiplicity classes. Statistical uncertainties are depicted as error bars.

Results and discussion
Particles and anti-particles turn out to have compatible transverse momentum distributions within uncertainties, consistently with previous results at the LHC. In the following, unless specified otherwise, we present results for their sum. The p T distributions of strange hadrons, measured in |y| < 0.5, are shown in Figs. 2, 3, 4 for the different multiplicity classes, selected using the estimators V0M, tracklets in |η| < 0.8 and tracklets in 0.8 < |η| < 1.5, as summarised in Table 1. The Lévy-Tsallis fit to the p T distributions, used for the extrapolation, are also displayed. The bottom panels depict the ratio to the minimum bias (INEL > 0) p T distribution. The p T spectra become harder as the multiplicity increases, as also shown in Fig. 5, which shows the average p T ( p T ) as a function of the mid-rapidity charged particle multiplicity. While the V0M and N 0.8<|η|<1.5 tracklets results show the same trend, the spectra obtained with the N |η|<0.8 tracklets estimator are systematically softer for comparable dN ch /dη values (though still compatible within uncertainty in the case of the strange baryons).
The increase of the p T as a function of the chargedparticle multiplicity (Fig. 5) is compatible, within uncertainties, for all particle species. The hardening of p T spectra with the charged-particle multiplicity was already reported for pp [40] and p-Pb collisions [22] at lower energies.
It is interesting to notice that the ratio of the measured p T spectra to the minimum bias one, shown in the bottom panel of Figs. 2, 3, 4, reaches a plateau for p T 4 GeV/c. This applies to all particle species and to all multiplicity syst. syst. uncorr.    Fig. 8 dN /dy (integrated over the full p T range) as a function of multiplicity for different strange particle species reported for different multiplicity classes (see text for details). Error bars and boxes represent statistical and total systematic uncertainties, respectively. The bin-tobin systematic uncertainties are shown by shadowed boxes estimators. The trend at highp T is highlighted in Fig. 6, which shows the integrated yields for p T > 4 GeV/c as a function of the mid-rapidity multiplicity. Both the yields of strange hadrons and the charged particle multiplicity are selfnormalised, i.e. they are divided by their average quoted on the INEL > 0 sample. The highp T yields of strange hadrons increase faster than the charged particle multiplicity. Despite the large uncertainties, the data also hint at the increase being non-linear. The self-normalised yields of strange hadrons reach, at high multiplicity, larger values for the N |η|<0.8 tracklets multiplicity selection as compared to the other estimators. For all multiplicity selections the self-normalised yields of baryons   5  10  15  20  25  30  35  5  10  15  20  25  30  35   5  10  15  20  25  30  35  5  10  15  20  25  30 Fig. 9 dN /dy (integrated over the measured p T ranges 0-12, 0.4-8, 0.6-6.5 and 0.9-5.5 GeV/c for K 0 S , , − and − , respectively) as a function of multiplicity for different strange particle species reported for different multiplicity classes (see text for details). Error bars and boxes represent statistical and total systematic uncertainties, respectively. The bin-to-bin systematic uncertainties are shown by shadowed boxes are higher than those of K 0 S mesons. The Monte Carlo models EPOS-LHC, PYTHIA 8 and PYTHIA 6, also shown in Fig. 6, reproduce the overall trend of strange hadrons seen in the data, with EPOS-LHC showing a clear difference between the K 0 S and the baryons, as observed in the data. Indeed, a non-linear increase can be explained by the combined effect of multiplicity fluctuations and interactions between different MPIs, which produces a collective boost [41]. Both color reconnection effects implemented in PYTHIA and a collective hydrodynamic expansion can account for the non-linear increase pattern.

PYTHIA6-Perugia2011
The left (right) panel of Fig. 7 shows the ratio of the p T spectra obtained with different estimators for multiplicity classes with comparable average dN ch /dη value tracklets estimators shows that the bias introduced by the latter estimator does not depend on the particle species and is more pronounced at low and intermediate p T , although the uncertainties are large.
The p T -integrated yields of strange hadrons are shown in Fig. 8 for the three estimators considered in this paper. For reference, Fig. 9 shows the results integrated in the measured p T range with no extrapolation. These are characterised by a smaller uncertainty and can be useful for Monte Carlo tuning. In the rest of this section, we focus on the discussion of the fully extrapolated yields. The results obtained with the V0M and N 0.8<|η|<1.5 tracklets estimators follow a similar linear trend, while the results obtained with N |η|<0.8 tracklets tend to saturate, showing a lower dN /dy for a similar high-multiplicity class. In order to gain insight on this difference, a generatorlevel PYTHIA (tune Perugia 2011) simulation study was performed. The abundance of strange hadrons in |y| < 0.5 was studied as a function of the estimated mid-rapidity charged shown by error bars and empty boxes, respectively. Shadowed boxes represent uncertainties uncorrelated across multiplicity  Fig. 11 Integrated yields of K 0 S , , , and as a function of dN ch /dη in V0M multiplicity event classes at √ s = 7 and 13 TeV. Statistical and systematic uncertainties are shown by error bars and empty boxes, respectively. Shadowed boxes represent uncertainties uncorrelated across multiplicity. The corresponding results obtained for INEL > 0 event class are also shown particle multiplicity for several event classes, selected using charged primary particles measured in the η ranges corresponding to the experimental estimators presented in this paper: this generator level study confirms the trend observed in the data, which can be understood in terms of a selection bias sensitive to fluctuations in particles yields. Indeed, an estimator based on charged tracks enhances charged primary particles over neutral particles or secondaries. If the multiplicity estimator and the observable under study are measured in the same η region one expects primary charged-particles to be enhanced with respect to strange hadrons, which explains the saturation of yields observed in Fig. 8 for the N |η|<0.8 tracklets estimator. As a further check of this interpretation, Fig. 10 shows the correlations between the yields of different strange hadrons for the multiplicity estimators studied in this paper. The trend is linear, and similar for all estimators, confirming that the selection bias on primary charged-particles is stronger than that on strange hadrons.
The energy dependence of the strangeness yields and p T versus the charged particle multiplicity at mid-rapidity is studied in Figs. 11 and 12, where our results are compared to the previous pp measurements at √ s = 7 TeV [21]. The minimum bias results in the INEL > 0 event class at √ s = 7 and 13 TeV [5] are also shown.
As can be seen in Fig. 11, the yields of strange hadrons increase with the charged particle multiplicity following a power law behaviour, and the trend is the same at √ s = 7 and 13 TeV. The INEL > 0 results also follow the same trend at all the tested centre-of-mass energies. This result indicates that the abundance of strange hadrons depends on the local charged particle density and turns out to be invariant with the collision energy, i.e. an energy scaling property applies for the multiplicity-dependent yields of strange hadrons. It should also be noted that the yields of particles with larger strange quark content increase faster as a function of multiplicity as already reported in [21]. Figure 13 shows the /K 0 S (no contribution considered here), /K 0 S and /K 0 S ratios compared to calculations from grandcanonical thermal models [13,14], which were found to satisfactorily describe central Pb-Pb data at √ s NN = 2.76 TeV. In the context of a canonical thermal model for a gas of hadrons, an increase of the relative strangeness abundance depending on the strange quark content can be interpreted as a consequence of a change in the system volume, called canonical suppression. Indeed, it was recently shown in [18] that the existing ALICE data can be described within this framework, introducing an additional parameter to quantify the rapidity window over which strangeness is effectively correlated. It was also suggested that strangeness follows a universal scaling behaviour in all colliding systems, when the transverse energy density [42] or the multiplicity per transverse area [43] are used as a scaling variable. While there are some caveats in the observation reported in these papers (due to the uncertainties on the transverse size of the system) this is a very intriguing observation, that could help to clarify the origin of strangeness enhancement.
The p T is seen in Fig. 12 to be harder at 13 TeV than at 7 TeV for event classes with a similar dN ch /dη. All the tested Monte Carlo models describe in a qualitative way the observed smooth rise of the p T with dN ch /dη; from a quantitative point of view, EPOS-LHC provides a slightly better description of the p T multiplicity evolution, especially for what concerns the strange baryons.
In Fig. 14 the results on the strange hadron yields as a function of the charged particle multiplicity at mid-rapidity are compared with some commonly-used general purpose QCD inspired models, focusing on the multiplicity classes defined by the V0M estimator. The yields of all measured strange  S results, however it shows a less pronounced increase of the strange baryons yields versus the charged particle multiplicity than what is observed in the data. The description worsens with increasing strangeness content. Recent attempts to improve the PYTHIA color reconnection scheme to account for the observed strangeness enhancement are discussed in [3], how-ever they cannot explain the reported experimental observations.
The DIPSY Monte Carlo [44,45] is based on a BFKL inspired initial state dipole evolution model [46] interfaced to the Ariadne model [47] for final state dipole evolution and to the PYTHIA hadronisation scheme, where the latter has been modified to take into account the high density of strings that occurs in events with several MPI. More specifically, depending on their transverse density, the strings are allowed to form "colour-ropes" leading to increased string tension, which leads to an increase of strangeness production stronger than in default PYTHIA [48]. However, despite a qualitative improvement in the description of the data, the discrepancy for the distributions remains large. Moreover, as discussed in [21], DIPSY also predicts an increase of protons normalised to pions as a function of multiplicity, which is not observed in our √ s = 7 TeV data. The Monte Carlo generators based on EPOS [49] rely on the Parton Based Gribov Regge Theory (PBGRT) [50]. A common feature of all the EPOS versions is a collective evolution of matter in the secondary scattering stage in all reactions, from pp to AA, with a core-corona separation mechanism [51] which defines the initial conditions of the secondary interactions. In EPOS, the initial parton scatterings create "flux tubes" that either escape the medium and hadronise as jets or contribute to the core, described in terms of hydrodynamics. The core is then hadronised in terms of a grand-canonical statistical model. The relative amount of multi-hadron production arising from the core grows with the number of MPIs. The EPOS-LHC model adopted in this analysis (distributed with the CRMC package 1.5.4) shows a pronounced increase of the strangeness yields as a function of the charged particle multiplicity at mid-rapidity. Despite the differences in the underlying physics, the comparison to data leads to a similar conclusion for EPOS-LHC and DIPSY.
The comparison of the Monte Carlo model predictions to the data indicates that the origin of the strangeness enhancement in hadronic collisions and its relation to the QCD deconfinement phase transition are still open problems for models. While the current versions of these models do not reproduce the data quantitatively, it is possible that the agreement can be improved with further tuning or with new implementations. A quantitative description of strangeness production and enhancement in a microscopic Monte Carlo model would represent a major step to understand flavour production in hadronic collisions at high energy.

Summary
We studied the production of primary strange and multistrange hadrons at mid-rapidity in pp collisions at √ s = 13 TeV, focusing on the multiplicity dependence. The  boxes, respectively. Shadowed boxes represent uncertainties uncorrelated across multiplicity. The results are compared to predictions from several Monte Carlo models main feature of this analysis is the usage of multiplicity estimators defined in different pseudorapidity regions, providing a different sensitivity to the fragmentation and MPI components of particle production and allowing for a detailed study of the selection biases due to fluctuations.
Hardening of the p T spectra with the increase of the multiplicity is observed, as already reported for pp [40] and p-Pb collisions [22] at lower energies.
The p T -integrated yields of strange hadrons increase as a function of multiplicity faster than the ones of unidentified charged-particles. This behaviour is even more pronounced for particles with higher strangeness content, confirming the earlier observations in pp collisions at √ s = 7 TeV [21]. This leads to a multiplicity-dependent increase of the ratio of the strange baryons and over K 0 S , while the over K 0 S ratio turns out to be constant within uncertainties. In the context of a canonical thermal model, an increase of the relative strangeness abundance depending on the strange quark content can be understood as a consequence of an increase in the system volume leading to a progressive removal of canonical suppression.
Comparing the 13 TeV results to the 7 TeV ones, the data exhibit an interesting scaling property: for multiplicity esti-mators selected in the forward region, the strange hadron yields turn out to be independent of the centre of mass energy, the √ s increase resulting just in harder p T spectra. The use of high-multiplicity triggered data collected during the full Run 2 period will allow us to test these scaling ansatzes extending the measurement to higher multiplicities, comparable to peripheral Pb-Pb collisions at the LHC (dN ch /dη ∼ 40-50).
The color reconnection effects implemented in PYTHIA 8 and DIPSY as well as the collective hydrodynamic expansion implemented in EPOS-LHC could not account quantitatively for all the reported results. Some of the qualitative features of the data are described by these models, including the indication of a non-linear increase of the yields of strange hadrons at high p T , similar to what was previously reported in [52] for J/ψ and D meson production. Although none of the tested models can reproduce the data from a quantitative point of view, EPOS-LHC and DIPSY do provide a better qualitative description of the strangeness enhancement. 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:

Data Availability Statement
This manuscript has associated data in a data repository. [Authors' comment: The numerical values of the data points will be uploaded to HEPData.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .