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 $\langle d{N}_{ch}/d\eta \rangle$ 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 $\langle d{N}_{ch}/d\eta \rangle$ at midrapidity ($\left|\eta \right|<0.5$) for PYTHIA8 (Monash tune), PYTHIA8 (Ropes model), EPOS-LHC, and Herwig7 are shown in

Table 1. Event multiplicity classes, their corresponding fraction of the cross-section ($\sigma /{\sigma }_{tot}$) and their corresponding $\langle d{N}_{ch}/d\eta \rangle$ 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}}$	${\langle {\mathit{dN}}_{\mathit{ch}}/\mathit{d}\mathit{\eta }\rangle }_{\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$

for $\sqrt{s}$ = 20 TeV and in

Table 2. Event multiplicity classes, their corresponding fraction of the cross-section ($\sigma /{\sigma }_{tot}$) and their corresponding $\langle d{N}_{ch}/d\eta \rangle$ 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}}$	${\langle {\mathit{dN}}_{\mathit{ch}}/\mathit{d}\mathit{\eta }\rangle }_{\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$

for $\sqrt{s}$ = 27 TeV. To fill the bins of $\langle d{N}_{ch}/d\eta \rangle$ 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. 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

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:

$$\langle {\tau }^2\rangle =\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{\langle {\tau }^2\rangle }\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 $\langle p_T\rangle$ 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\langle N^{pp}\rangle /N_{max}^{pp})}^2$, with $\langle N^{pp}\rangle$ 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, $\langle p_T\rangle$ as a function of $\langle d{N}_{c}h/d\eta \rangle$.

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 $\langle d{N}_{ch}/d\eta \rangle$ 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 $\langle p_T\rangle$

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 $\langle d{N}_{ch}/d\eta \rangle$ 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’ $\langle p_T\rangle$ exhibits a pattern that is somewhat consistent with the data and that slightly rises with increasing collision energy. We can observe that the $\langle p_T\rangle$ increases as a function of $\langle d{N}_{ch}/d\eta \rangle$ for all particles.

The present implementation of rope hadronization does not explain the $\langle p_T\rangle$, 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 $\langle p_T\rangle$ as well, the model’s parameters should be improved or new ones should be added. Although there is a tendency for $\langle p_T\rangle$ to rise with multiplicity, the values of $\langle p_T\rangle$ as a function of $\langle d{N}_{ch}/d\eta \rangle$ are lower than those of data.

The influence of the collective flow feature shown in the model may be seen in the $\langle p_T\rangle$ produced by EPOS-LHC. The results of $\langle p_T\rangle$ as a function of $\langle d{N}_{ch}/d\eta \rangle$ demonstrate that $\langle p_T\rangle$ 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 $\langle p_T\rangle$. It has lower values than the data for $K_S^0$. It fits badly for smaller multiplicity values for $\Omega$.

Herwig7 accurately depicts the $\langle p_T\rangle$ 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 $\langle p_T\rangle$ grows with increasing multiplicity, and when looking at its development with collision energy, the $\langle p_T\rangle$ 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 $\langle d{N}_{ch}/d\eta \rangle$, 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 $\langle d{N}_{ch}/d\eta \rangle$ 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, $\langle dN/dy\rangle$ of $K_S^0$ as a function of $\langle d{N}_{ch}/d\eta \rangle$ is described well compared to ALICE data. For $\Lambda$, $\Xi$, and $\Omega$ which have higher content of strangeness, $\langle dN/dy\rangle$ of $K_S^0$ as a function of $\langle d{N}_{ch}/d\eta \rangle$ 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 $\langle d{N}_{ch}/d\eta \rangle <15$ and for $\langle d{N}_{ch}/d\eta \rangle <11$ for $\Omega$ but for larger values of multiplicity, the values obtained from EPOS-LHC for $\langle dN/dy\rangle$ are higher than those of data.

The results of Herwig7 for $\langle dN/dy\rangle$, 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 $\langle d{N}_{ch}/d\eta \rangle$.

For all the event generators, the evolution of $\langle dN/dy\rangle$ as a function of $\langle d{N}_{ch}/d\eta \rangle$ 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 $\langle p_T\rangle$ 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 $\langle p_T\rangle$ 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 $\langle p_T\rangle$.

EPOS-LHC model gives a good description for the strange particles’ production as a function of multiplicity and $\langle p_T\rangle$ 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.