Strangeness Production from Proton–Proton Collisions at Different Energies by Using Monte Carlo Simulation

Abstract

Nuclear matter, at sufficiently energy density and high temperature, undergoes a transition to a state of strongly interacting QCD matter in which quarks and gluons are not confined known as the Quark–Gluon Plasma (QGP). QGP is usually produced in high-energy collisions of heavy nuclei in the laboratory, where an enhancement of strange hadrons’ production is observed. Many of the effects which are typical of heavy ion phenomenology have been observed in high-multiplicity proton–proton (pp) collisions. The enhancement of strange particles’ production in pp collisions was reported at $\sqrt{s}=7$ TeV and $\sqrt{s}=13$ TeV in 2017 and 2020, respectively, and it was found that the integrated yields of strange particles, relative to pions, increase notably with the charged-particle multiplicity of events. Here, we report the multiplicity dependence of strange particles at $\left|y\right|<0.5$ in pp collisions at $\sqrt{s}$ = 7 TeV, 13 TeV, 20 TeV, and 27 TeV from a Monte Carlo simulation using PYTHIA8, EPOS-LHC, and Herwig7.

1. Introduction

The study of strange hadrons is one technique used to investigate the properties of QCD, the theory of strongly interacting matter. The enhancement of strangeness production in heavy ion collisions has been proposed as a sign of the presence of a Quark–Gluon Plasma (QGP) [1,2,3]. QGP is typically produced during high-energy collisions of heavy nuclei in the laboratory, where enhanced strangeness production is observed [4,5,6,7,8]. There is a consensus now that QGP filled the universe in the first 20 $\mathsf{\mu }$s after the Big Bang. This state of matter impacting nuclei at high energies has been the subject of studies at the European Center for Particle Physics (CERN) and Brookhaven National Laboratory (BNL) for many years. The smoking gun of QGP creation was formerly thought to be heavy ion collisions. Because azimuthal correlations and mass-dependent hardening of $p_T$ distributions have been observed in high-multiplicity pp and proton–nucleus collisions at the LHC [9,10,11,12,13,14,15,16,17,18,19], and these phenomena in nuclear collisions are attributed to the formation of QGP, studying pp collisions at high-multiplicity is of significant interest as it aids in the microscopic understanding of the phenomena known from nuclear collisions.

There is a significant evidence that the universe was as a fireball at the beginning, which is called “Big Bang”, with extremely high energy density and temperature. At early stages, the temperature was high enough ($T>100$ GeV), such that all the known particles were extremely relativistic. Due to asymptotic freedom, even strong interactions weaken and an ideal plasma of quarks and gluons forms. So this was a QGP, a system of hot and weakly interacting colour-charged particles, in equilibrium with the other species. The discovery and thorough understanding of QGP is significant for QCD since it predicts the long-range behaviour where the theory is still poorly understood. Phase transitions in QCD at high temperatures are also relevant in cosmological investigations. If such a phase transition existed, it would have occurred in the universe within the first microsecond of the Big Bang [20]. Understanding the QGP phase transition can also provide essential information about the inner core of a neutron star with a large nuclear density [21]. During the expansion phase, the universe cooled and the quarks, anti-quarks, and gluons combined to make the hadron structures resulting in the baryonic matter that is observed today. Quantitatively, it has been challenging to address the transition from quarks and gluons to baryons. Substantial effort in theoretical physics is made to understand this transition by using lattice gauge theory, a high-level computational method. Because the generation of a huge number of particles in a finite volume of the collision shows that a big value of energy density occurs, ultra-relativistic nuclear collisions have aided us in our hunt for QGP formation in the laboratory [22]. The expectation of producing a QGP in ultra-relativistic heavy ion collisions led us to the research of QGP diagnostics. In the diagnostic investigations, we strive to find some particle spectra or ratios that should be significantly different if QGP production did not occur. There are several recommendations for QGP diagnostics or signatures. One way to identify these fingerprints is by measuring energy flow and density, particle spectra and their correlations, multiplicity fluctuations, direct photons and dileptons, increased creation of strange particles, and $j/\psi$ suppression.

A practical challenge arose when theories concerning QGP creation in relativistic heavy ion collisions developed: “How to distinguish the locally colour deconfined QGP state from a gas of confined hadrons?” The QGP strangeness signal was developed as a result between 1979 and 1986. The strange quark, the heaviest of the three light quark flavours, emerged as the critical signature for QGP in 1980–1982. When the colour bonds break, the deconfined state has a substantially higher abundance of strange quark pairs, according to [3,23]. It was observed that the gluon component of the QGP primarily creates strange quark pairs quickly and precisely on the necessary time scale [2]. If hadronization was primarily driven by the merging of pre-existing quarks and antiquarks, as suggested by [1], then the high strangeness density at the moment of QGP hadronization is the natural source of multi-strange hadrons [24]. The investigation of strangeness can reveal a wealth of information regarding the development and evolution of QGP fireballs. Since there is a plentiful supply of strange hadrons, the strangeness observable has gained popularity in experiments and also can be detected over a wide range of kinematics. So today, many experimental results are available, and all of these results are harmonious with the production of hadronic particles occurring from a dense source in which the deconfined strange quarks are already produced before the formation of hadrons. In February 2000, CERN reported the discovery of QGP based on the SPS findings that addressed strangeness and multi-strange anti-hyperon generation. Other observables served as the foundation for the RHIC community’s announcement of the QGP finding. The finding of strangeness enhancement in high-multiplicity proton–proton collisions for the first time at 7 TeV in 2017 marked the LHC’s continuous contribution to QGP physics.

The current study aims to investigate how the production of strangeness from various Monte Carlo event generators (MCEGs), relying on various hadronization mechanisms, depends on the event charged particle multiplicity and the center of mass energy at two different values $\sqrt{s}$ = 7 TeV and 13 TeV. The LHC has already been run at Those two energies, and to compare the results from the generators with those previously published by ALICE at the same energies [25,26]. We also examine how these generators behave at higher energies, $\sqrt{s}$ = 20 TeV and $\sqrt{s}$ = 27 TeV, where the LHC may operate in the upcoming runs.

2. Analysis Details

The results presented here are for primary strange hadrons [27]. We apply the same criteria and conditions applied on the sample of ALICE data so that the results obtained from the simulation process can be directly compared to ALICE results. We report the measurements for events that have one charged particle, at least, in the pseudorapidity interval $\left|\eta \right|<1$. To study the relation between the multiplicity and the production of strange and multi-strange hadrons, the sample is divided into event classes depending on the total ionization energy that is deposited in the forward detectors which cover the pseudorapidity regions $2.8<\eta <5.1$ and $-3.7<\eta <-1.7$. The mean pseudorapidity densities $⟨d{N}_{ch}/d\eta ⟩$ of primary charged particles are measured at midrapidity, $\left|\eta \right|<0.5$. The event multiplicity classes, their corresponding fraction of the cross-section ($\sigma /{\sigma }_{tot}$) and their corresponding $⟨d{N}_{ch}/d\eta ⟩$ at midrapidity ($\left|\eta \right|<0.5$) for PYTHIA8 (Monash tune), PYTHIA8 (Ropes model), EPOS-LHC, and Herwig7 are shown in Table 1 for $\sqrt{s}$ = 20 TeV and in Table 2 for $\sqrt{s}$ = 27 TeV. To fill the bins of $⟨d{N}_{ch}/d\eta ⟩$ at $\sqrt{s}$ = 20 TeV and $\sqrt{s}$ = 27 TeV, we used bin values generated from Herwig7 because it fits very well with ALICE data compared to the other models. Monte Carlo data sets were produced using the individual MCEGs and then the analysis procedure was after that processed using the Rivet [28] MCEG validation tool.

We chose the Rivet-specific definition of ALICE primary particles in order to directly compare the findings from the MCEGs with the data. The implementation of this definition in Rivet is described in [29]. The Rivet analysis is based on the analysis used for the ALICE collaboration’s published findings at $\sqrt{s}$ = 7 TeV [25,30].

In the generated plots, the x-axis represents the average charged particle multiplicity measured at $\left|\eta \right|<0.5$. The 10 bins were used for $K_S^0$, $\Lambda$ and $\Xi$ and 5 bins for $\Omega$ because its production is lower than that of the previous three because of its quark content; it is formed of three s-quarks, so generated events were divided into 5 classes for $\Omega$ rather than 10. Each bin represents a multiplicity class, so the events are divided into classes based on the charged particle multiplicity. Then, for each class the average charged particle multiplicity was measured as a representative for the whole multiplicity class. The division of events into classes and calculation of the average charged multiplicity is based upon the Rivet analysis described in [31].

3. Models

3.1. PYTHIA8

3.1.1. Monash Tune

PYTHIA [32] event generator is a tool that simulates collisions between different types of particles at high-energy. It comprises a consistent set of physics models to describe the evolution from a hard process of a few-body to a complex multiparticle final state. It includes a library of hard processes, as well as methods for matching and merging hard processes and parton showers, models for initial- and final-state parton showers, beam remnants, multiparton interactions, string fragmentation, and particle decays. Although exploring the experimental consequences of theoretical models is one of the main tasks performed by this event generator, there are several important tasks, such as interpreting experimental data, studying the detector performance, and developing search strategies that could be performed by it. Hadronization—the process by which the final outgoing coloured partons are converted into colourless particles—in the default tune of PYTHIA8, Monash tune, depends on Lund string model [33,34].

3.1.2. Ropes Model

We also present the results of another configuration for PYTHIA8 which uses the parameters of the Ropes model [35]. Models which are based on string hadronization usually treat the strings independently, so no interaction between the confined colour fields is allowed. The ropes model suggested that strings that are close in space can fuse to form “colour ropes” in nucleus collisions. The ropes are colour multiplets arising from strings overlapping with each other in densely populated events. Such ropes produce more strangeness and this increase in the strangeness production has been suggested as a marker for QGP production. In the Ropes model, since strings can be combined to form higher multiplets, this gives rise to production of baryons and strangeness compared to models that treat the strings independently. In PYTHIA8, the ropes are not the sole source of colour multiplets; there is also another source in the form of junctions made by QCD-CR [36]. Assuming they have all already formed, the rope implementation takes advantage of the junctions’ correspondence to lower multiplets. Due to this, we should activate both models to use the PYTHIA8 rope hadronization. Table 3 lists the PYTHIA8 settings configured to employ rope hadronization.

3.1.3. Description of the Lund String Model Underlying PYTHIA8

The process of hadronization or fragmentation converts the final outgoing coloured partons into colourless hadrons. Since this transition is non-perturbative, models must manage it. The Lund string model serves as the basis for PYTHIA [34,37]. For simplicity, we consider a model in one space dimension with only one quark flavour and a single mesonic state with mass m. In the Lund hadronization model, the probability, $\mathcal{P}$, for the production of a specific state with n mesons with momenta $p_i$ ($i=1,\dots ,n$) is given by the relation [38]:

$\mathcal{P}\propto \left\{{\displaystyle \prod _{i=1}^n\left[N{d}^2{p}_i\delta (p_i^2-m^2)\right]{\delta }^{\left(2\right)}(\sum p_i-P_{tot})}\right\}exp(-bA)$

A phase space factor is a term enclosed in curly parenthesis, where the weighting between states with various meson counts is determined by the dimensionless constant N. In the exponent, the word $bA$ denotes the imaginary portion of the massless string’s action, which causes the string to decay and have a finite lifetime. A measures the string’s space-time area before it splits apart, and b is a constant.

In a Monte Carlo simulation, the result of the above equation can be produced by iteratively creating mesons starting from one of the string ends, where each meson consumes a portion of the remaining energy, z. The probability distribution or splitting function provides each step’s pertinent z-value [35]:

$f\left(z\right)=N\frac{{(1-z)}^a}{z}exp(-\frac{b{m}^2}{z})$

where the constant a is related to N and b through the normalization constraint $\int f\left(z\right)dz=1$. The production points for the pairs will be located around a hyperbola in space-time, with a typical proper time determined by:

$⟨{\tau }^2⟩=\frac{1+a}{b{\kappa }^2}$

where $\kappa$ is the string tension. This timescale is related to the particle multiplicity by the relation [35]:

$dN/dy\sim \sqrt{⟨{\tau }^2⟩}\kappa /m=\sqrt{\frac{1+a}{b{m}^2}}$

3.1.4. Monte Carlo Samples from PYTHIA8

Producing the results from PYTHIA8 was done by applying the appropriate Rivet analysis code for the present study on the generated events from PYTHIA8. Rivet was used with PYTHIA8 through the Rivet interface to PYTHIA8 program [39].

3.2. EPOS-LHC

We present the results produced from EPOS-LHC [40], which is an MCEG for minimum bias hadronic interactions. EPOS-LHC is part of the CRMC (Cosmic Ray Monte Carlo) [41] package that provides access to various cosmic ray and non-cosmic ray event generators. EPOS-LHC serves to simulate heavy ion interactions, as well as cosmic ray air showers. EPOS is based on the Gribov Regge Theory [42], which was based on the VENUS model [43] for soft interactions and the QGSJET model [44] for the semi-hard scattering.

3.2.1. Elementary Interaction in the EPOS Model

The elementary interaction model in EPOS-LHC is the same as the previous version of EPOS 1.99 [45]. Nucleus–nucleus (proton–proton) scattering amounts to many elementary collisions happening in parallel. Such an elementary scattering is the so-called “parton ladder” shown in Figure 1. Parton evolutions from the projectile and target sides towards the center (small x) are represented by a parton ladder. In the most basic scenario, DGLAP states that an evolution equation controls the evolution. The initial scatterings cause strings to develop, fragmenting into segments typically referred to as hadrons. According to some critical values of density segments per unit volume, we discriminate between string segments in dense areas and those in low density areas at an appropriately early stage.The high density areas are referred to as core, and the low density areas as corona [46]. In the subsequent section, we focus on the core part, which is unique in EPOS and provides interesting effects not accounted for in other HEP models.

3.2.2. New Features in EPOS-LHC

In EPOS 1.99, the maximal radial rapidity, $y_{rad}$ was parameterized as function of the system energy (s) and size as [40]:

$y_{rad}=y_{rad}^{max}+y_{rad}^{mi}·log(1+\sqrt{\frac{s}{N_{pair}}})$

where $N_{pair}$ is the number of possibly interacting pairs of nucleons, $y_{rad}^{max}$ and $y_{rad}^{mi}$ are parameters. The evolution with $N_{pairs}$ was safe and easy to test with heavy ion data as SPS and RHIC, but the evolution with energy especially for $N_{pair}=1$ in case of p-p could lead to wrong extrapolation at high energy. Based on the results of identified particle spectra from CMS [47], results show that the increase in the $⟨p_T⟩$ as a function of the multiplicity does not depend on the center-of-mass energy but it increases with increasing of multiplicity. So it was natural to parameterize all flows as a function of the total mass $M_{core}$, which is directly linked to the final multiplicity. So, as a consequence of that, the final formula for the maximal radial rapidity in case of p-p collisions, $y_{rad}^{pp}$, can be written as the following equation [40]:

$y_{rad}^{pp}=y_{rad}^{px}·F_{pp}·log(\frac{M_{pp}}{M_{min}})$

where $F_{pp}=min{(1,2⟨N^{pp}⟩/N_{max}^{pp})}^2$, with $⟨N^{pp}⟩$ being the average number of segments going to core per participating pair of nucleons, is a normalization factor going to 0 in case of heavy ion collisions. $M_{pp}$ is the mass of the part of the core coming from the pair of nucleons with the largest multiple scattering, in this case, p-p. $F_{pp}=1$ and $M_{pp}=M_{core}$ in case of p-p scattering.

We present the strangeness production using the EPOS-LHC model to examine the effects of the collective flow feature on the results of strangeness production as a function of multiplicity, particularly the study of the average transverse momentum, $⟨p_T⟩$ as a function of $⟨d{N}_{c}h/d\eta ⟩$.

3.2.3. Monte Carlo Samples from EPOS-LHC

For the generation of Monte Carlo samples from EPOS-LHC, we used the EPOS-LHC model distributed under the CRMC package [48]. Further analysis of the data was done using the Rivet analysis, which applies the appropriate final state conditions so that the results can be compared directly to ALICE data.

3.3. Herwig7

Results from Herwig7 [49] are also presented. Herwig7 is a multi-purpose high-energy physics event generator for simulating hard hadron–hadron collision. Herwig7 has benefited from the experience gained with the HERWIG [50] and the Herwig++ [51] event generators. Herwig7 provides highly improved and extended physics capabilities compared to both of its predecessors. Herwig7 uses the cluster hadronization model [52] to describe the formation of hadrons from the quarks and gluons produced in the parton shower.

3.3.1. Description of the Cluster Hadronization Model Underlying Herwig7 and Cluster Decays

After the elementary hard subprocess of the event, all partons engaged in hard scatterings, further scattering, and partonic decays down to low scales have evolved in the parton showers, the final state typically consists of coloured partons close in momentum space to partons with which they share a colour index, called their colour partner. Herwig7 projects these colour–anticolour pairings onto singlet states called clusters, which decay into hadrons and hadron resonances using the cluster hadronization model. Herwig7 uses spin correlations between individual decays and a matrix element description of the decay product distributions to simulate hadron decays.

The final step of the cluster hadronization model is the decay of the cluster into a pair of hadrons. A quark–antiquark or diquark–antidiquark pair (q, $\overline{q}$) is extracted from the vacuum for a cluster of a given flavour ($q_1$, ${\overline{q}}_2$), and a pair of hadrons with flavours ($q_1$,$\overline{q}$) and (q,${\overline{q}}_2$) is formed. Based on the available phase space, spin, and hadrons’ flavour, the relevant hadrons are chosen from among all conceivable hadrons with the appropriate flavour. All cluster models follow the same general methodology; however, minor differences exist. In Herwig7, the original model of ref. [52], the approach of ref. [53], and a new variant that addresses the issue of the low rate of baryon production in the approach of ref. [53] are implemented. As a result the weight for the production of the hadrons $a(q_1,\overline{q})$ and $b(q,{\overline{q}}_2)$ is given as in [51] by:

$W(a_{(q_1,\overline{q})},b_{(q,{\overline{q}}_2)}|q_1,{\overline{q}}_2)=P_q{w}_a{w}_b{s}_a{s}_b{p}_{a,b}^{*}$

where $P_q$ is the weight for the production of the given quark–antiquark or diquark–antidiquark pair, $w_{a,b}$ are the weights for the production of individual hadrons, $s_{a,b}$ are the suppression factors for the hadrons, which allow the production rates of individual meson multiplets, and singlet and decuplet baryons to be adjusted, and $p_{a,b}^{*}$ is the momentum of the hadrons in the rest frame of the decaying cluster.

3.3.2. Monte Carlo Samples from Herwig7

To generate the Monte Carlo samples with Herwig7, Herwig7 program which is based on ThePEG program was used, both available here [54]. In the simulation process, we implemented the same Rivet analysis mentioned above to obtain our results.

4. Results and Discussion

4.1. Ratio of Yields to the Pion Yield

In Figure 2, Figure 3, Figure 4 and Figure 5, the ratios of the yields of $K_S^0$, $\Lambda$, $\Xi$, and $\Omega$ to the pion (${\pi }^{+}+{\pi }^{-}$) yield as a function of $⟨d{N}_{ch}/d\eta ⟩$ are shown as produced from PYTHIA8 (Monash tune), PYTHIA8 (Ropes model), EPOS-LHC, and Herwig7, respectively, at $\sqrt{s}=7$ TeV, $\sqrt{s}=13$ TeV, $\sqrt{s}=20$ TeV and $\sqrt{s}=27$ TeV. The results are compared to the ALICE results at $\sqrt{s}=7$ TeV.

The Monash tune of PYTHIA8 shows an apparent disagreement with ALICE data. In contrast to the other event generators, strange hadrons’ generation is suppressed with increasing multiplicity. For this reason, it was necessary to set the model’s settings so that strings can overlap in highly crowded events, where QGP is most likely to be formed. PYTHIA8 event generator has to not only regard the strings as independent but also consider some strings interaction to allow the formation of colour multiplets to account for the increased generation of strange particles in events with high energy density. The PYTHIA8 Monash tune does not differentiate particle strangeness content; the evolution of the particle to pion ratio as a function of multiplicity is almost the same for all particles.

When compared to the Monash tuning, the Ropes model gives a better representation of the strange particle ratio to pions as a function of multiplicity, suggesting that the colour multiplets presented in the model are responsible for the excess in strange particles’ rates which may simulate effects of plasma generation.

In EPOS-LHC, in addition to the new features with respect to the prior EPOS 1.99, the parton ladder description shows the behaviour of strangeness enhancement with increasing multiplicity. Although the results demonstrate some enhancement, they do not match the ALICE results very well, especially in the high multiplicity regions. The EPOS-LHC model does not well represent the saturation behaviour of strange particle ratios to pions at high multiplicity.

The results demonstrate that Herwig7 well describes the increase in strangeness generation observed in high-multiplicity proton–proton collisions. Compared to the other Monte Carlo models in the present study, it has the best agreement with the data. It describes well the evolution of particles’ production at low and high multiplicity, especially the saturation region of particles’ production ratio relative to pions at high multiplicity.

Overall, we note that the behaviour of increasing ratio to pions with increasing multiplicity is more clear for particles with more strangeness content. Additionally, for all event generators, we notice that the results are almost the same for each individual generator at different energies, which means that strangeness production in proton–proton collisions is determined by the characteristics of the final state rather than by the collision energy.

4.2. Average Transverse Momentum $⟨p_T⟩$

In Figure 6, Figure 7, Figure 8 and Figure 9, the average transverse momentum of $K_S^0$, $\Lambda$, $\Xi$, and $\Omega$ as a function of $⟨d{N}_{ch}/d\eta ⟩$ are shown as produced from the same Monte Carlo models at $\sqrt{s}=7$ TeV, $\sqrt{s}=13$ TeV, $\sqrt{s}=20$ TeV and $\sqrt{s}=27$ TeV. The results are compared to the ALICE results at $\sqrt{s}=13$ TeV.

The Monash tuning of PYTHIA8 deviates from the data in terms of the strangeness production ratio to pions, but the generated particles’ $⟨p_T⟩$ exhibits a pattern that is somewhat consistent with the data and that slightly rises with increasing collision energy. We can observe that the $⟨p_T⟩$ increases as a function of $⟨d{N}_{ch}/d\eta ⟩$ for all particles.

The present implementation of rope hadronization does not explain the $⟨p_T⟩$, despite the fact that the Ropes model provides better rates for the creation of strange particles. In order to provide a good description for the $⟨p_T⟩$ as well, the model’s parameters should be improved or new ones should be added. Although there is a tendency for $⟨p_T⟩$ to rise with multiplicity, the values of $⟨p_T⟩$ as a function of $⟨d{N}_{ch}/d\eta ⟩$ are lower than those of data.

The influence of the collective flow feature shown in the model may be seen in the $⟨p_T⟩$ produced by EPOS-LHC. The results of $⟨p_T⟩$ as a function of $⟨d{N}_{ch}/d\eta ⟩$ demonstrate that $⟨p_T⟩$ depends on the multiplicity—it rises with increasing multiplicity—but only very slightly on the collision energy. This illustrates the effect of parameterizing the flows as a function of the core mass as previously discussed. For $\Lambda$ and $\Xi$, EPOS-LHC provides a good description of the particle’s $⟨p_T⟩$. It has lower values than the data for $K_S^0$. It fits badly for smaller multiplicity values for $\Omega$.

Herwig7 accurately depicts the $⟨p_T⟩$ of particles for $\Lambda$ and $\Xi$. It has lower values than data for $K_S^0$. For $\Omega$, it fits poorly for smaller multiplicity values. In general, the $⟨p_T⟩$ grows with increasing multiplicity, and when looking at its development with collision energy, the $⟨p_T⟩$ increases only minimally with rising energy.

It is clear from the results that for all models of the study that the average transverse momentum increases as the multiplicity increases. Additionally, as a function of $⟨d{N}_{ch}/d\eta ⟩$, the average transverse momentum increases very slightly with increasing collision energy.

4.3. Integrated Yields

In Figure 10, Figure 11, Figure 12 and Figure 13, the results on $K_S^0$, $\Lambda$, $\Xi$, and $\Omega$ yields as a function of $⟨d{N}_{ch}/d\eta ⟩$ are shown as produced from the generators at $\sqrt{s}=7$ TeV, $\sqrt{s}=13$ TeV, $\sqrt{s}=20$ TeV, and $\sqrt{s}=27$ TeV. The results are compared to the ALICE results at $\sqrt{s}=7$ TeV and 13 TeV.

For the Monash tune of PYTHIA8, $⟨dN/dy⟩$ of $K_S^0$ as a function of $⟨d{N}_{ch}/d\eta ⟩$ is described well compared to ALICE data. For $\Lambda$, $\Xi$, and $\Omega$ which have higher content of strangeness, $⟨dN/dy⟩$ of $K_S^0$ as a function of $⟨d{N}_{ch}/d\eta ⟩$ has lower values than those of data due to the suppression of strange particles’ production resulting from treating strings as independent.

The results from the Ropes model have a very good agreement with data for all the strange particles except for $\Lambda$; the agreement is not like the other particles, but as a whole it is much better than the results from the Monash tune of PYTHIA8.

For certain particles, the EPOS-LHC results matched the ALICE data well, whereas the matching was less good for others. For $K_S^0$ and $\Lambda$, the results fit well with data in nearly the whole range of multiplicity. For $\Xi$, it fits well for $⟨d{N}_{ch}/d\eta ⟩<15$ and for $⟨d{N}_{ch}/d\eta ⟩<11$ for $\Omega$ but for larger values of multiplicity, the values obtained from EPOS-LHC for $⟨dN/dy⟩$ are higher than those of data.

The results of Herwig7 for $⟨dN/dy⟩$, compared to the other event generators of the study, has the best matches with ALICE data for all the particles in the whole range of $⟨d{N}_{ch}/d\eta ⟩$.

For all the event generators, the evolution of $⟨dN/dy⟩$ as a function of $⟨d{N}_{ch}/d\eta ⟩$ for all the particles is nearly independent on the collision energy.

5. Conclusions

The production of strangeness in proton–proton collisions for high multiplicity events is investigated using a variety of event generators, including PYTHIA8 (Monash tuning), PYTHIA8 (Ropes model), EPOS-LHC, and Herwig7. The level of agreement between the Monte Carlo event generators and the data varies; some models match the ALICE data, while others have a poor agreement or none. That gives an idea of how well the underlying physics models of the generators can explain the strangeness enhancement phenomenon as a signal for QGP.

The PYTHIA8 with Monash tune cannot be used to describe the strange particles’ enhancement as a function of multiplicity; however, it gives a good description for the average transverse momentum $⟨p_T⟩$ of the particles. Using the Ropes model parameters for PYTHIA8 did well for the strange particles production as a function of multiplicity but the Ropes model poorly describes $⟨p_T⟩$ of the particles. The ropes model is considered a good step in PYTHIA8 to study strange particles’ production, but it needs some modification to get good agreement with data, when investigating $⟨p_T⟩$.

EPOS-LHC model gives a good description for the strange particles’ production as a function of multiplicity and $⟨p_T⟩$ in general, although it miss fits with the data in some ranges. Some adjustments also should be considered for the EPOS-LHC model concerning the issue of strangeness production.

Herwig7 gives a very good description for the strangeness enhancement issue. So, if we use Herwig7 to study the other signals of QGP and find that it gives a good description for these signals too, then this would give a better understanding of QGP characteristics.

This investigation of strangeness production from several Monte Carlo event generators should provide guidance on how to adjust the models of various event generators and their settings to match with real data. However, this should be paired with the studies of the other QGP signals utilizing these generators in order to have a clearer description of the QGP phase transition as a whole and as a further step to map the QCD phase diagram within the scope of the existing real data.

Author Contributions

Formal analysis, Y.M.; Funding acquisition, M.A.M.; Project administration, M.A.M. and Y.M.; Resources, Y.M.; Software, A.H.; Supervision, M.A.M.; Validation, A.A.A. and M.N.E.-H.; Writing—original draft, A.H.; Writing—review & editing, A.A.A. and M.N.E.-H. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

In this work we used data published by the the ALICE collaboration at $\sqrt{s}$ = 7 TeV [55] and $\sqrt{s}$ = 13 TeV [56] for comparison.

Acknowledgments

This paper is based upon work supported by Science, Technology and Innovation Funding Authority (STDF) under Project ID: 30163. The authors would like to thank the Rivet authors for their permanent existence to provide solutions to the problems faced us while using the Rivet framework. In addition, we would like to thank the PYTHIA authors, EPOS authors, and Herwig authors for their guidance through the analysis process and their help to settle issues faced us during the work. The authors would like to thank Christian Bierlich for his suggestions on rope hadronization and for providing the parameters for the PYTHIA8 Ropes model.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

QCD	Quantum Chromodynamics
QGP	Quark–Gluon Plasma
MCEG	Monte Carlo event generator
CRMC	Cosmic Ray Monte Carlo

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Elementary interaction in the EPOS model [40].

Figure 2. Integrated yield ratios to pions ([Forumla omitted. See PDF.]) as a function of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using PYTHIA8 (Monash tune) [25].

Figure 3. Integrated yield ratios to pions ([Forumla omitted. See PDF.]) as a function of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using PYTHIA8 (Ropes model) [25].

Figure 4. Integrated yield ratios to pions ([Forumla omitted. See PDF.]) as a function of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using EPOS-LHC [25].

Figure 5. Integrated yield ratios to pions ([Forumla omitted. See PDF.]) as a function of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using Herwig7 [25].

Figure 6. [Forumla omitted. See PDF.] of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using PYTHIA8 (Monash tune) [25].

Figure 7. [Forumla omitted. See PDF.] of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using PYTHIA8 (Ropes model) [25].

Figure 8. [Forumla omitted. See PDF.] of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using EPOS-LHC [25].

Figure 9. [Forumla omitted. See PDF.] of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using Herwig7 [25].

Figure 10. Integrated yields of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] as a function of of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using PYTHIA8 (Monash tune) [25,26].

Figure 11. Integrated yields of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] as a function of of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using PYTHIA8 (Ropes model) [25,26].

Figure 12. Integrated yields of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] as a function of of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using EPOS-LHC [25,26].

Figure 13. Integrated yields of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] as a function of of [Forumla omitted. See PDF.] measured in [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 7 TeV, 13 TeV, 20 TeV, and 27 TeV using Herwig7 [25,26].

Table 1

Event multiplicity classes, their corresponding fraction of the cross-section ($\sigma /{\sigma }_{tot}$) and their corresponding $⟨d{N}_{ch}/d\eta ⟩$ at midrapidity ($\left|\eta \right|<0.5$) for PYTHIA8 (Monash), PYTHIA8 (Ropes), EPOS-LHC, and Herwig7 at $\sqrt{s}$ = 20 TeV.

Class Name	$\mathit{\sigma }/{\mathit{\sigma }}_{\mathit{tot}}$	${⟨{\mathit{dN}}_{\mathit{ch}}/\mathit{d}\mathit{\eta }⟩}_{\left|\mathit{\eta }\right|<0.5}$
		PYTHIA8 (Monash)	PYTHIA8 (Ropes)	EPOS-LHC	Herwig7
I	0–1%	$30.22\pm 0.042$	$29.13\pm 0.042$	$33.33\pm 0.047$	$30.53\pm 0.053$
II	1–5%	$24.19\pm 0.018$	$23.64\pm 0.019$	$25.98\pm 0.020$	$22.87\pm 0.020$
III	5–10%	$19.82\pm 0.015$	$19.41\pm 0.015$	$20.59\pm 0.015$	$18.48\pm 0.015$
IV	10–15%	$16.67\pm 0.014$	$16.28\pm 0.014$	$16.85\pm 0.014$	$15.46\pm 0.015$
V	15–20%	$14.22\pm 0.013$	$13.83\pm 0.012$	$14.02\pm 0.012$	$13.05\pm 0.012$
VI	20–30%	$11.39\pm 0.008$	$11.09\pm 0.008$	$10.82\pm 0.008$	$10.34\pm 0.008$
VII	30–40%	$8.44\pm 0.007$	$8.51\pm 0.007$	$7.71\pm 0.006$	$7.85\pm 0.007$
VIII	40–50%	$6.32\pm 0.006$	$6.65\pm 0.006$	$5.62\pm 0.005$	$5.92\pm 0.006$
IX	50–70%	$4.16\pm 0.003$	$4.69\pm 0.003$	$3.81\pm 0.003$	$3.86\pm 0.003$
X	70–100%	$2.37\pm 0.002$	$2.64\pm 0.002$	$2.18\pm 0.002$	$2.16\pm 0.002$

Table 2

Event multiplicity classes, their corresponding fraction of the cross-section ($\sigma /{\sigma }_{tot}$) and their corresponding $⟨d{N}_{ch}/d\eta ⟩$ at midrapidity ($\left|\eta \right|<0.5$) for PYTHIA8 (Monash), PYTHIA8 (Ropes), EPOS-LHC, and Herwig7 at $\sqrt{s}$ = 27 TeV.

Class Name	$\mathit{\sigma }/{\mathit{\sigma }}_{\mathit{tot}}$	${⟨{\mathit{dN}}_{\mathit{ch}}/\mathit{d}\mathit{\eta }⟩}_{\left|\mathit{\eta }\right|<0.5}$
		PYTHIA8 (Monash)	PYTHIA8 (Ropes)	EPOS-LHC	Herwig7
I	0–1%	$33.08\pm 0.044$	$32.44\pm 0.052$	$37.19\pm 0.050$	$34.25\pm 0.093$
II	1–5%	$26.56\pm 0.020$	$26.29\pm 0.024$	$28.83\pm 0.021$	$25.37\pm 0.034$
III	5–10%	$21.75\pm 0.016$	$21.51\pm 0.019$	$22.80\pm 0.016$	$20.41\pm 0.026$
IV	10–15%	$18.21\pm 0.014$	$17.85\pm 0.017$	$18.48\pm 0.014$	$16.94\pm 0.024$
V	15–20%	$15.39\pm 0.013$	$15.12\pm 0.016$	$15.11\pm 0.013$	$14.26\pm 0.022$
VI	20–30%	$12.0.6\pm 0.008$	$11.99\pm 0.010$	$11.37\pm 0.008$	$11.26\pm 0.014$
VII	30–40%	$8.91\pm 0.007$	$9.05\pm 0.008$	$8.00\pm 0.006$	$8.27\pm 0.012$
VIII	40–50%	$6.71\pm 0.006$	$7.09\pm 0.007$	$5.88\pm 0.006$	$6.17\pm 0.010$
IX	50–70%	$4.29\pm 0.003$	$4.96\pm 0.004$	$3.94\pm 0.003$	$4.05\pm 0.005$
X	70–100%	$2.38\pm 0.002$	$2.70\pm 0.002$	$2.23\pm 0.002$	$2.22\pm 0.003$

Table 3

Configuration parameters for PYTHIA8 rope hadronization.

Parameter	Value	Parameter Description
MultiPartonInteractions:pT0Ref	2.15	Parameter of the MPI model to keep total multiplicity reasonable
BeamRemnants:remnantMode	1
BeamRemnants:saturation	5
ColourReconnection:mode	1
ColourReconnection:allowDoubleJunRem	off
ColourReconnection:m0	0.3	Parameters related to Junction formation/QCD based CR
ColourReconnection:allowJunctions	on
ColourReconnection:junctionCorrection	1.2
ColourReconnection:timeDilationMode	2
ColourReconnection:timeDilationPar	0.18
Ropewalk:RopeHadronization	on	Parameters of the rope model
Ropewalk:doShoving	on
Ropewalk:tInit	1.5
Ropewalk:deltat	0.05
Ropewalk:tShove	0.1
Ropewalk:gAmplitude	0
Ropewalk:doFlavour	on
Ropewalk:r0	0.5
Ropewalk:m0	0.2
Ropewalk:beta	0.1
PartonVertex:setVertex	on
PartonVertex:protonRadius	0.7	Settings of vertex information
PartonVertex:emissionWidth	0.1

References

1. Koch, P.; Muller, B.; Rafelski, J. Strangeness in Relativistic Heavy Ion Collisions. Phys. Rep.; 1986; 142, pp. 167-262. [DOI: https://dx.doi.org/10.1016/0370-1573(86)90096-7]

2. Rafelski, J.; Muller, B. Strangeness Production in the Quark–Gluon Plasma. Phys. Rev. Lett.; 1982; 48, 1066.Erratum in Phys. Rev. Lett. 1986, 56, 2334 [DOI: https://dx.doi.org/10.1103/PhysRevLett.48.1066]

3. Rafelski, J.; Hagedorn, R. From Hadron Gas to Quark Matter. 2. Proceedings of the International Symposium on Statistical Mechanics of Quarks and Hadrons; Bielefeld, Germany, 24–31 August 1980.

4. Andersen, E.; Antinori, F.; Armenise, N.; Bakke, H.; Ban, J.; Barberis, D.; Beker, H.; Beusch, W.; Bloodworth, I.J.; Bohm, J. et al. Strangeness enhancement at mid-rapidity in Pb Pb collisions at 158-A-GeV/c. Phys. Lett. B; 1999; 449, pp. 401-406. [DOI: https://dx.doi.org/10.1016/S0370-2693(99)00140-9]

5. Anticic, T.; Afanasiev, S.V.; Barna, D.; Bartke, J.; Barton, R.A.; Betev, L.; Bialkowska, H.; Billmeier, A.; Blume, C.; Blyth, C.O. et al. Ξ⁻ and Ξ¯⁺ production in central Pb + Pb collisions at 158-GeV/c per nucleon. Phys. Lett. B; 2002; 538, pp. 275-281. [DOI: https://dx.doi.org/10.1016/S0370-2693(02)01970-6]

6. Antinori, F.; Bacon, P.A.; Badala, A.; Barbera, R.; Belogianni, A.; Bhasin, A.; Bloodworth, I.J.; Bombara, M.; Bruno, G.E.; Bull, S.A. et al. Energy dependence of hyperon production in nucleus nucleus collisions at SPS. Phys. Lett. B; 2004; 595, pp. 68-74. [DOI: https://dx.doi.org/10.1016/j.physletb.2004.05.025]

7. Abelev, B.I.; Aggarwal, M.M.; Ahammed, Z.; Anderson, B.D.; Arkhipkin, D.; Averichev, G.S.; Bai, Y.; Balewski, J.; Barannikova, O.; Barnby, L.S. et al. Enhanced strange baryon production in Au + Au collisions compared to p + p at s(NN)**(1/2) = 200-GeV. Phys. Rev. C; 2008; 77, 044908. [DOI: https://dx.doi.org/10.1103/PhysRevC.77.044908]

8. Alme, J.; Erdal, H.A.; Helstrup, H.; Hetl, K.F.; Kileng, B.; Altinpinar, S.; Djuvsl, O.; Haal, O.S.; Huang, M.; Langoy, R. et al. Multi-strange baryon production at mid-rapidity in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV. Phys. Lett. B; 2014; 728, pp. 216-227. Erratum in Phys. Lett. B 2014, 734, 409–410 [DOI: https://dx.doi.org/10.1016/j.physletb.2014.05.052]

9. Khachatryan, V.; Sirunyan, A.M.; Tumasyan, A.; Adam, W.; Bergauer, T.; Dragicevic, M.; Ero, J.; Fabjan, C.; Friedl, M.; Fruhwirth, R. et al. Observation of Long-Range Near-Side Angular Correlations in Proton-Proton Collisions at the LHC. JHEP; 2010; 9, 91. [DOI: https://dx.doi.org/10.1007/JHEP09(2010)091]

10. Khachatryan, V.; Sirunyan, A.M.; Tumasyan, A.; Adam, W.; Asilar, E.; Bergauer, T.; Brandstetter, J.; Brondolin, E.; Dragicevic, M.; Ero, J. et al. Evidence for collectivity in pp collisions at the LHC. Phys. Lett. B; 2017; 765, pp. 193-220. [DOI: https://dx.doi.org/10.1016/j.physletb.2016.12.009]

11. Abelev, B.B. ALICE Collaboration. Multiplicity Dependence of Pion, Kaon, Proton and Lambda Production in p-Pb Collisions at $\sqrt{s_{NN}}$ = 5.02 TeV. Phys. Lett. B; 2014; 728, pp. 25-38. [DOI: https://dx.doi.org/10.1016/j.physletb.2013.11.020]

12. Adamova, D.; Adam, J.; Aggarwal, M.M.; Aglieri Rinella, G.; Agnello, M.; Agrawal, N.; Ahammed, Z.; Ahmad, S.F.; Ahn, S.U.; Aiola, S. et al. Multi-strange baryon production in p-Pb collisions at $\sqrt{s_{\mathbf{N}\mathbf{N}}}$ = 5.02 TeV. Phys. Lett. B; 2016; 758, pp. 389-401. [DOI: https://dx.doi.org/10.1016/j.physletb.2016.05.027]

13. Chatrchyan, S.; Khachatryan, V.; Sirunyan, A.M.; Tumasyan, A.; Adam, W.; Aguilo, E.; Bergauer, T.; Dragicevic, M.; Ero, J.; Fabjan, C. et al. Observation of Long-Range Near-Side Angular Correlations in Proton-Lead Collisions at the LHC. Phys. Lett. B; 2013; 718, pp. 795-814. [DOI: https://dx.doi.org/10.1016/j.physletb.2012.11.025]

14. Alme, J.; Erdal, H.A.; Helstrup, H.; Hetl, K.F.; Kileng, B.; Altinpinar, S.; Djuvsl, O.; Fehlker, D.; Haal, O.S.; Huang, M. et al. Long-range angular correlations on the near and away side in p-Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV. Phys. Lett. B; 2013; 719, pp. 29-41. [DOI: https://dx.doi.org/10.1016/j.physletb.2013.01.012]

15. ATLAS Collaboration. Observation of Associated Near-Side and Away-Side Long-Range Correlations in $\sqrt{s_{NN}}$ = 5.02 TeV Proton-Lead Collisions with the ATLAS Detector. Phys. Rev. Lett.; 2013; 110, 182302. [DOI: https://dx.doi.org/10.1103/PhysRevLett.110.182302]

16. Aguilar-Saavedra, J.A.; Amor Dos Santos, S.P.; Anjos, N.; Cantrill, R.; Carvalho, J.; Castro, N.F.; Conde Muino, P.; Da Cunha Sargedas de Sousa, M.J.; Do Valle Wemans, A.; Fiolhais, M. et al. Measurement with the ATLAS detector of multi-particle azimuthal correlations in p+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV. Phys. Lett. B; 2013; 725, pp. 60-78. [DOI: https://dx.doi.org/10.1016/j.physletb.2013.06.057]

17. Chatrchyan, S.; Khachatryan, V.; Sirunyan, A.M.; Tumasyan, A.; Adam, W.; Bergauer, T.; Dragicevic, M.; Ero, J.; Fabjan, C.; Friedl, M. et al. Multiplicity and Transverse Momentum Dependence of Two- and Four-Particle Correlations in pPb and PbPb Collisions. Phys. Lett. B; 2013; 724, pp. 213-240. [DOI: https://dx.doi.org/10.1016/j.physletb.2013.06.028]

18. ALICE Collaboration. Long-range angular correlations of ß, K and p in p-Pb collisions at $\sqrt{s_{\mathrm{N}\mathrm{N}}}$ = 5.02 TeV. Phys. Lett. B; 2013; 726, pp. 164-177. [DOI: https://dx.doi.org/10.1016/j.physletb.2013.08.024]

19. Khachatryan, V.; Sirunyan, A.M.; Tumasyan, A.; Adam, W.; Asilar, E.; Bergauer, T.; Brandstetter, J.; Brondolin, E.; Dragicevic, M.; Ero, J. et al. Multiplicity and rapidity dependence of strange hadron production in pp, pPb, and PbPb collisions at the LHC. Phys. Lett. B; 2017; 768, pp. 103-129. [DOI: https://dx.doi.org/10.1016/j.physletb.2017.01.075]

20. Weinberg, S. Beyond the First Three Minutes. Phys. Scr.; 1980; 21, pp. 773-781. [DOI: https://dx.doi.org/10.1088/0031-8949/21/5/033]

21. Iwamoto, N. Chemical kinetics of beta decay reactions in degenerate quark matter. Phys. Rev. D; 1983; 28, pp. 2353-2362. [DOI: https://dx.doi.org/10.1103/PhysRevD.28.2353]

22. Kajantie, K.; McLerran, L. Probes of the Quark Gluon Plasma in High Energy Collisions. Annu. Rev. Nucl. Part. Sci.; 1987; 37, pp. 293-323. [DOI: https://dx.doi.org/10.1146/annurev.ns.37.120187.001453]

23. Rafelski, J. Extreme States of Nuclear Matter—1980. Eur. Phys. J.; 2015; 51, pp. 282-324. [DOI: https://dx.doi.org/10.1140/epja/i2015-15115-y]

24. Rafelski, J. Formation and Observables of the Quark-Gluon Plasma. Phys. Rep.; 1982; 88, 331.

25. Adam, J.; Adamova, D.; Aggarwal, M.M.; Rinella, G.A.; Agnello, M.; Agrawal, N.; Ahammed, Z.; Ahmad, S.; Ahn, S.U.; Aiola, S. et al. Enhanced production of multi-strange hadrons in high-multiplicity proton-proton collisions. Nature Phys.; 2017; 13, pp. 535-539. [DOI: https://dx.doi.org/10.1038/nphys4111]

26. Acharya, S.; Adamova, D.; Adhya, S.P.; Adler, A.; Adolfsson, J.; Aggarwal, M.M.; Aglieri Rinella, G.; Agnello, M.; Agrawal, N.; Ahammed, Z. et al. Multiplicity dependence of (multi-)strange hadron production in proton-proton collisions at $\sqrt{s}$ = 13 TeV. Eur. Phys. J. C; 2020; 80, 167. [DOI: https://dx.doi.org/10.1140/epjc/s10052-020-7673-8]

27. ALICE Collaboration. The ALICE Definition of Primary Particles; Alice Publications: Farmington, ME, USA, 2017.

28. Bierlich, C.; Buckley, A.; Butterworth, J.; Christensen, C.H.; Corpe, L.; Grellscheid, D.; Grosse-Oetringhaus, J.F.; Gutschow, C.; Karczmarczyk, P.; Klein, J. et al. Robust Independent Validation of Experiment and Theory: Rivet version 3. SciPost Phys.; 2020; 8, 26. [DOI: https://dx.doi.org/10.21468/SciPostPhys.8.2.026]

29. Rivet Team. Available online: https://rivet.hepforge.org/code/dev/classRivet_1_1ALICE_1_1PrimaryParticles.html (accessed on 23 September 2021).

30. Available online: https://rivet.hepforge.org/analyses/ALICE_2016_I1471838 (accessed on 23 September 2021).

31. Available online: https://rivet.hepforge.org/analyses/ALICE_2015_PPCentrality.html (accessed on 23 September 2021).

32. Bierlich, C.; Chakraborty, S.; Desai, N.; Gellersen, L.; Helenius, I.; Ilten, P.; Lonnblad, L.; Mrenna, S.; Prestel, S.; Preuss, C.T. et al. A comprehensive guide to the physics and usage of PYTHIA 8.3. arXiv; 2022; arXiv: 2203.11601

33. Andersson, B.; Gustafson, G.; Ingelman, G.; Sjostrand, T. Parton Fragmentation and String Dynamics. Phys. Rep.; 1983; 97, pp. 31-145. [DOI: https://dx.doi.org/10.1016/0370-1573(83)90080-7]

34. Sjöstrand, T. Jet fragmentation of multiparton configurations in a string framework. Nucl. Phys. B; 1984; 248, pp. 469-502. [DOI: https://dx.doi.org/10.1016/0550-3213(84)90607-2]

35. Bierlich, C.; Gustafson, G.; Lönnblad, L.; Tarasov, A. Effects of Overlapping Strings in pp Collisions. JHEP; 2015; 3, 148. [DOI: https://dx.doi.org/10.1007/JHEP03(2015)148]

36. Christiansen, J.R.; Skands, P.Z. String Formation Beyond Leading Colour. JHEP; 2015; 8, 3. [DOI: https://dx.doi.org/10.1007/JHEP08(2015)003]

37. Andersson, B.; Gustafson, G.; Soderberg, B. A General Model for Jet Fragmentation. Z. Phys. C; 1983; 20, 317. [DOI: https://dx.doi.org/10.1007/BF01407824]

38. Andersson, B.; Gustafson, G.; Söderberg, B. A probability measure on parton and string states. Nucl. Phys. B; 1986; 264, pp. 29-59. [DOI: https://dx.doi.org/10.1016/0550-3213(86)90471-2]

39. Available online: https://www.pythia.org/ (accessed on 6 May 2021).

40. Pierog, T.; Karpenko, I.; Katzy, J.M.; Yatsenko, E.; Werner, K. EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider. Phys. Rev. C; 2015; 92, 034906. [DOI: https://dx.doi.org/10.1103/PhysRevC.92.034906]

41. Ulrich, R.; Pierog, T.; Baus, C. Cosmic Ray Monte Carlo Package, CRMC. 2021; Available online: https://zenodo.org/record/5270381#.Y2h1InZByUk (accessed on 6 May 2022). [DOI: https://dx.doi.org/10.5281/zenodo.5270381]

42. Drescher, H.J.; Hladik, M.; Ostapchenko, S.; Pierog, T.; Werner, K. Parton based Gribov-Regge theory. Phys. Rep.; 2001; 350, pp. 93-289. [DOI: https://dx.doi.org/10.1016/S0370-1573(00)00122-8]

43. Werner, K. Strings, pomerons, and the venus model of hadronic interactions at ultrarelativistic energies. Phys. Rep.; 1993; 232, pp. 87-299. [DOI: https://dx.doi.org/10.1016/0370-1573(93)90078-R]

44. Kalmykov, N.N.; Ostapchenko, S.S.; Pavlov, A.I. Quark-Gluon String Model and EAS Simulation Problems at Ultra-High Energies. Nucl. Phys. B Proc. Suppl.; 1997; 52, pp. 17-28. [DOI: https://dx.doi.org/10.1016/S0920-5632(96)00846-8]

45. Pierog, T.; Werner, K. EPOS Model and Ultra High Energy Cosmic Rays. Nucl. Phys. B Proc. Suppl.; 2009; 196, pp. 102-105. [DOI: https://dx.doi.org/10.1016/j.nuclphysbps.2009.09.017]

46. Werner, K. Core-corona separation in ultra-relativistic heavy ion collisions. Phys. Rev. Lett.; 2007; 98, 152301. [DOI: https://dx.doi.org/10.1103/PhysRevLett.98.152301]

47. Chatrchyan, S.; Khachatryan, V.; Sirunyan, A.M.; Tumasyan, A.; Adam, W.; Aguilo, E.; Bergauer, T.; Dragicevic, M.; Ero, J.; Fabjan, C. et al. Study of the Inclusive Production of Charged Pions, Kaons, and Protons in pp Collisions at $\sqrt{s}$ = 0.9, 2.76, and 7 TeV. Eur. Phys. J. C; 2012; 72, 2164. [DOI: https://dx.doi.org/10.1140/epjc/s10052-012-2164-1]

48. Available online: https://gitlab.iap.kit.edu/AirShowerPhysics/crmc (accessed on 6 May 2021).

49. Bellm, J.; Gieseke, S.; Grellscheid, D.; Platzer, S.; Rauch, M.; Reuschle, C.; Richardson, P.; Schichtel, P.; Seymour, M.H.; Siodmok, A. et al. Herwig 7.0/Herwig++ 3.0 release note. Eur. Phys. J. C; 2016; 76, 196. [DOI: https://dx.doi.org/10.1140/epjc/s10052-016-4018-8]

50. Corcella, G.; Knowles, I.G.; Marchesini, G.; Moretti, S.; Odagiri, K.; Richardson, P.; Seymour, M.H.; Webber, B.R. HERWIG 6: An Event generator for hadron emission reactions with interfering gluons (including supersymmetric processes). JHEP; 2001; 1, 010. [DOI: https://dx.doi.org/10.1088/1126-6708/2001/01/010]

51. Bahr, M.; Gieseke, S.; Gigg, M.A.; Grellscheid, D.; Hamilton, K.; Latunde-Dada, O.; Platzer, S.; Richardson, P.; Seymour, M.H.; Sherstnev, A. et al. Herwig++ Physics and Manual. Eur. Phys. J. C; 2008; 58, pp. 639-707. [DOI: https://dx.doi.org/10.1140/epjc/s10052-008-0798-9]

52. Webber, B.R. A QCD Model for Jet Fragmentation Including Soft Gluon Interference. Nucl. Phys. B; 1984; 238, pp. 492-528. [DOI: https://dx.doi.org/10.1016/0550-3213(84)90333-X]

53. Kupco, A. Cluster hadronization in HERWIG 5.9. Proceedings of the Workshop on Monte Carlo Generators for HERA Physics (Plenary Starting Meeting); Hamburg, Germany, 27–30 April 1998; pp. 292-300.

54. Available online: https://herwig.hepforge.org/ (accessed on 6 May 2021).

55. ALICE Collaboration. Enhanced Production of Multi-Strange Hadrons in High-Multiplicity Proton-Proton Collisions; HEPData (Collection) Alice Publications: Farmington, ME, USA, 2017; [DOI: https://dx.doi.org/10.17182/hepdata.77284]

56. ALICE Collaboration. Multiplicity Dependence of (Multi-)Strange Hadron Production in Proton-Proton Collisions at $\sqrt{s}$ = 13 TeV; HEPData (Collection) Alice Publications: Farmington, ME, USA, 2021; [DOI: https://dx.doi.org/10.17182/hepdata.93535]