#### Abstract

We have studied the transverse momentum () spectra of the final-state strange particles, including , , , and , produced in high energy lead–lead (Pb–Pb), proton–lead (–Pb), xenon–xenon (Xe–Xe) collisions at the Large Hadron Collider (LHC). Taking into account the contribution of multiquark composition, whose probability density distribution is described by the modified Tsallis–Pareto-type function; we simulate the spectra of the final-state strange particles by a Monte Carlo method, which is shown to be in good agreement with the experimental data in most the cases. The kinetic freeze-out parameters are obtained. The present method provides a new tool for studying the spectra of various particles produced in high energy collisions, reflecting more realistically the collision process, which is of great significance to study the formation and properties of the produced particles.

#### 1. Introduction

Plenty of final-state particles are produced in relativistic heavy-ion collisions [1, 2] in various collision processes, resulting in different configurations of final-state particles. In high energy collisions, novel physical phenomena can appear with the most representative being formation of quark-gluon plasma (QGP) [3, 4]. In addition, multiparticle production provides abundant information on thermal and statistical properties of a system. Although quantum chromodynamics (QCD) can provide an important theoretical basis to study strong interactions among quarks and features of collision system [5, 6], thermal and statistical method is still a powerful tool for such analysis.

Investigation of the transverse momentum () (or transverse mass, ) spectra of the final-state particles is an effective and fast method [1, 2] to study the processes of multiparticle production and system evolution. This method can help to obtain the thermodynamic parameters of the final-state particles and collision system. The spectra of different final-state particles have been studied extensively. In addition to baryons, leptons, and other elementary particles have also attracted much attention.

In recent years, the production of single and multistrange particles [7–10] has attracted an increased interest, being extensively researched. Based on the exotic properties [11], the spectra of strange particles measured in high energy collisions by various collaborations have been analyzed and predicted utilizing Tsallis function [12–17], nonequilibrium chemical or kinetic freeze-out model [8], etc.

Relative to pions, enhanced production of multistrange hadrons [10–21] in high-multiplicity proton–proton collisions has been observed experimentally for the first time, and it is found that the integrated yields of strange particles increases significantly with the event charged-particle multiplicity. Besides, the azimuthal angular correlation and the mass dependence of due to the existence of QGP in high-energy collisions have been also reported. Further researches on strange particles can better explore and reveal the properties of QGP.

We are interested in the exploration of the properties of different strange particles. In contrast to our previous work [22–24], in this work we adopt a new algorithm to simulate and analyze the spectra of strange particles. We aim to restore the collision process more realistically and extract more accurate characteristic parameters. The multiquark composition of baryons is considered. A Monte Carlo method [25–27] is used to simulate different transverse momenta carried by various quarks. The modified Tsallis–Pareto-type function [22–24] is used to define the contribution of multiquarks, from which the kinetic freeze-out temperature () [28–31], average transverse flow velocity () [32, 33], and other related parameters can be directly extracted.

This paper is structured as follows. The formalism and Monte Carlo method are briefly introduced in Section 2. The simulated or fitted results, their comparison with the data, and discussion are given in Section 3. Finally, in Section 4 we summarize our main observations and conclusions.

#### 2. Formalism and Method

To deal with of strange particles produced in high energy collisions, we briefly introduce the analysis, based on the Monte Carlo method. Generally, the modified Tsallis–Pareto-type function [22–24]. is suitable to characterize the spectra of final-state particles in low- and intermediate- regions. Here, is the number of particles; is the normalization constant; is the effective temperature of a collision system; is an entropy-related index, which is used to describe the degree of nonequilibrium of the system; is the transverse mass of a particle; is the rest mass of the particle, and is the correction index, which makes the function to fit better the spectra in low- region.

It should be noted that when we set . Equation (1) is naturally converged to the Tsallis–Pareto-type function. That is to say, the introduction of in Equation (1) is a modification of Tsallis–Pareto-type function [17]. To extract characteristic parameters such as the thermal or kinetic freeze-out temperature of a collision system [28–31] and the average transverse flow velocity of the produced particles [32, 33] at the quark level, following References [34–37], we perform the Lorentz-like transformation on and in Equation (1) by using and , where is the Lorentz factor.

Thus, the transverse momentum of the -th quark contributed to the particle is assumed to obey the new modified probability density function

The temperature parameter at this time has been naturally converted from in Equation (1) to the kinetic freeze-out temperature . In Equation (2), is the constituent mass of the -th quark. As shown in the literature [38], we have GeV/ for up and down quarks and GeV/ for a strange quark. Thus, we have constructed the probability density function, which satisfies multiquark states.

Within the defined interval of of the -th quark, we have the normalization condition where the upper limit is usually quite a large value, but not infinite. In a Monte Carlo calculation, denotes the random number in , and we have the relation satisfied by to be where is a small shift from . According to Equation (4), one can obtain a series of discrete values of which satisfy Equation (2). In the calculation, we take GeV/. Considering the contributions of multiquarks, we have to 2 for mesons, and to 3 for baryons, because meson consists of two quarks, and baryon has three quarks.

In general condition, the particle transverse momentum is the vector superposition of of two or three quarks. In the right-handed Cartesian coordinate system -, let the axis be the beam direction, plane be the reaction plane, and plane be the transverse one. In the source rest frame, is assumed to be isotropic. The movement of the source along the axis constitutes the longitudinal flow, and the interactions among different sources causes the transverse flow. The longitudinal flow is not the focus of the present work, though it can be described by the rapidity distribution. The transverse flow is described by the average transverse flow velocity that is mentioned in the above discussion.

For the -th quark, we have the - and -components of to be where the isotropic azimuthal angle is distributed uniformly in . Using a random number distributed uniformly in , one has in the Monte Carlo calculation. According to the principle of vector superposition, we obtain the expressions of ’s components and of the final-state particle to be where the sums are considered due to the contributions of 2 or 3 quarks; the upper limit “2” for the sums corresponds to the meson, and the limit “3” for the sums corresponds to the baryon in the first two relations in Equation (6).

According to the above analysis, we may obtain a lot of discrete values of by the iterative calculations. Finally, we may perform a statistical analysis on the discrete values of and obtain suitable distributions such as and , which can be compared with the experimental data, where denotes the rapidity, and is the width of rapidity bin at midrapidity. In calculations, the statistical interval of is taken to be 0.1 GeV/.

Although we may simply use the modified Tsallis–Pareto-type function to fit the experimental data and adopt the -minimisation scheme to obtain the fit parameters, the results obtained by the simple method are at the particle level, which seem not to be deeper insight compared to the results at the quark and gluon level. To obtain the results at the partonic-level, we may consider the contribution of multiple partons to particle’s transverse momentum. However, the analytical expression obtained at the partonic-level is not available due to the complex calculations for the superposition of multiple transverse momenta with randomized azimuthal angles. Instead, we may use the Monte Carlo method to obtain numerical results. This is suitable for the case discussed in this paper.

It should be noted that several restrictions are used in the Monte Carlo calculations. *First*, the constituent masses of the considered quarks are determinate. These masses are taken from the literature [38]. *Second*, of a quark is restricted to obey a given function. This function is modified from the current Tsallis-like function. *Third*, an isotropic emission in the rest frame of emission source is assumed. The momentum in the rest frame of emission source is then obtained. *Fourth*, the movement of emission source is restricted by its rapidity which is assumed to be evenly distributed in the projectile or target thermalized region in the rapidity space. The rapidity of emission source is also related to longitudinal flow, which restricts the rapidity range. *Fifth*, the conservation of energy and momentum is obeyed. The combination of of 2- or 3-quarks with randomized azimuthal angles into of hadron can be conducted, in which any energy beyond hadron’s rest mass is converted into hadron’s kinetic energy. *Sixth*, the distribution of is restricted by the experimental data.

#### 3. Results and Discussion

##### 3.1. Comparison with Data

In the study, we analyze the spectra of strange particles generated in high energy lead–lead (Pb–Pb), proton–lead (–Pb), and xenon–xenon (Xe–Xe) collisions at the Large Hadron Collider (LHC). Figure 1 shows the spectra of (a) , (b) , (c) , (d) , (e) , (f) , (g) , (h) , and (i) produced in Pb–Pb collisions with different centrality intervals marked in the panels at the collision energy TeV per nucleon pair. The symbols represent the data measured by the ALICE Collaboration [39–41]. The mid-pseudorapidity range of panel (a) is , and those of other panels are . For better distinguishing and comparing, the data with different centralities in different panels are multiplied by 10 to the -th power. In addition, the gray rectangles represent the uncertainty of the experimental data, whose length is the uncertainty of the abscissa, whereas the width is the uncertainty of the ordinate, which denote the quadratic sum of the statistical and systematical uncertainties. The solid curves represent the Monte Carlo results calculated from the probability density function in Equation (2) and the expression in Equation (6).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

We use to quantify the deviation of the Monte Carlo results from the experimental data, where denotes the order number of the data. The smaller the mean of is, the closer the Monte Carlo results to the experimental data are. The free parameters, the entropy-related index , the kinetic freeze-out temperature , the correction index , and the average transverse flow velocity are extracted by the method of least squares. The values of free parameters, and the number of degree of freedom (ndof) are listed in Table 1. One can see that in most cases the Monte Carlo results agree well or approximately with the experimental spectra of strange particles. In few cases, the Monte Carlo results are in qualitative agreement with the experimental data, for which the values of /ndof are quite large.

The spectra of (a) , (b) , and (c) are displayed in Figure 2 for Pb–Pb collisions at TeV, quoted from the ALICE Collaboration [42, 43]. The mid-pseudorapidity is in panel (a) and in panels (b) and (c). Similarly to Figure 1, the data with different centralities are represented by different symbols, and the gray rectangles represent the uncertainty of the data. The solid curves show the Monte Carlo results. For clarity, we use a power index amplification from to . The free parameters , , , and , as well as /nd of are listed in Table 2. One can see that the Monte Carlo results show an approximate agreement with the data in some cases and a qualitative agreement with the data in the other cases.

**(a)**

**(b)**

**(c)**

In addition to Pb–Pb collisions, the experimental spectra in other collisions are also studied in order to better explore the properties of strange particles. Figure 3 shows the spectra of (a) , (b) , (c) , (d) , (e) , (f) , and (g) in –Pb collisions at TeV. The experimental data in the panels are all quoted from the ALICE Collaboration when the midrapidity is (a, c, e) and (b, d, f, g) [44, 45]. Different symbols indicate the different centrality intervals. The Monte Carlo results are given by the solid curves. The coefficients in parentheses are magnification factors for better distinguishing the spectra in different centrality intervals. The relevant parameters , , , and , as well as /nd of are listed in Table 3. One can see that the Monte Carlo results are in good agreement with the data in most cases and approximately in agreement with the data in few cases.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

Besides, we have also studied the spectra of (a) and (b) produced in Xe–Xe collisions at TeV in Figure 4. The symbols represent the data measured by the ALICE Collaboration, and the mid-pseudorapidity is [46]. Similarly, the data with different centralities are expressed as different symbols with gray rectangles, and the solid curves show the Monte Carlo results. Another, in legends are magnification factors for better distinguishing the spectra. The relevant parameters , , , and , as well as of are listed in Table 4. Once again, one can see that the Monte Carlo results are in good agreement with the experimental data.

**(a)**

**(b)**

##### 3.2. Tendencies of Parameters

To further understand the regularities shown by the spectra of strange particles, we analyze the tendencies of various parameters with changing centrality observed in the work. Figure 5 shows the dependence of entropy-related index (a, b) on centrality and in Pb–Pb, (c) –Pb, and (d) Xe–Xe collisions at the LHC energies. The symbols are quoted from Tables 1–4, where the results from different particles are represented by different symbols. As an entropy-related index, or describes the degree of nonequilibrium of the system, where () is the entropy index. Generally, a larger or a () closer to 1 corresponds to a higher degree of equilibrium. The present work shows that is large enough or () is close to 1, and the system stays in an approximate equilibrium. Meanwhile, the system in central collisions stays in higher degree of equilibrium. In addition, as the rest mass of strange particle increases, the value of increases. The multistrange particles correspond to larger than the single-strange particles. These results imply that the system stays in larger degree of equilibrium when it forms multistrange particles.

**(a)**

**(b)**

**(c)**

**(d)**

Similar to Figure 5, Figure 6 shows the dependence of kinetic freeze-out temperature on centrality in the mentioned collisions at the LHC. One can see that decreases slightly with the decrease of centrality from central to peripheral collisions in most cases. In few cases, the decrease is significant. In addition, extracted from the single-strange particle spectra is larger than that from the multistrange particle spectra. This finding suggests that the single-strange particles are formed earlier than the multistrange particles, though the latter may leave the system earlier than the former in the hydrodynamic evolution due to different masses.

**(a)**

**(b)**

**(c)**

**(d)**

As shown in Figure 7, we present the changing law of the correction index on centrality in the mentioned collisions at the LHC. One can see that with the decrease of centrality from central to peripheral collisions, is almost invariant, though in few cases shows very slight decrease. In most cases, the values of are far from 1, which means that the introduction of is necessary. We have also compared the present fits with with those by . An obvious difference appears in the low- region. This implies that the fits by are not suitable. In fact, in the fits of the present work, we use at the first, then we change to fit the spectra if is not satisfactory.

**(a)**

**(b)**

**(c)**

**(d)**

The dependence of the average transverse flow velocity on centrality in the mentioned collisions at the LHC is displayed in Figure 8. Similar conclusions to the dependences of on and on can be obtained. That is, decreases slightly with the decrease of centrality from central to peripheral collisions in most cases. In few cases, the decrease is significant. In addition, extracted from the single-strange particle spectra is larger than that from the multistrange particle spectra. This finding also suggests that the single-strange particles are formed earlier than the multistrange particles.

**(a)**

**(b)**

**(c)**

**(d)**

In order to further study the final states of the strange particles and the disorder degree of the system, we give the results of pseudoentropy [28], which is based on the probability density function of , inspired by the entropy [47–49], based on the probability density function of multiplicity , where is in the units of GeV/. As is done in our recent work [24], the width of bin is taken to be 0.1 GeV/, though the width is changeable. The unit of is neglected due to the fact that is dimensionless.

In Figure 9, we demonstrate the relation of the pseudoentropy versus the centrality . It is seen from the results that decreases slightly with the decrease of centrality from central to peripheral collisions in most cases. Only in few cases, decreases significantly with the decrease of . The dependence of on is similar to those of , , and on . Generally, the multistrange particles show larger than the single-strange particles, though both for multi and single-strange particles are negative. This also reflects the mass-dependence of , which is similar to those of the other parameters.

**(a)**

**(b)**

**(c)**

**(d)**

Analyzing Figures 5–9, one can see that the correlations between and , and , and , and , and , and and are positive; while the correlations between and , and , and , and and are very small or negligible. These results are understandable due to the fact that all , , and determine mainly, and is mostly defined by, the shape of the spectra in intermediate- and high- regions; while determines mainly the spectra in low- or even very low- region. Of course, because of the requirement of normalization, all parameters affect the spectra in the whole region, though the level of influence in different regions are different.

##### 3.3. Further Discussion

Although the extraction of and in present work is performed at the partonic-level, most extractions in the literature have been done at the particle level, showing inconsistent trends. For example, in terms of dependence of parameters on centrality from central to peripheral collisions, some works show an increase in and a decrease in , and in small system or peripheral collisions [50–52]. Other work shows a decrease in both and , and is considerable in small system or peripheral collisions [53]. Generally, the relative difference between two in central and peripheral collisions is small ().

We observe a slightly higher in central collisions at both the partonic- and particle level. If so, the centrality dependence of is similar to those of the chemical freeze-out temperature and effective temperature. However, it is still an open question which is larger. In any case, the explanation is understandable. A higher in central collisions means a higher excitation, and a lower in central collisions means a longer lifetime, of the hot and dense matter. Although we may extract the thermal parameters at the particle level, we aim to perform it at the partonic-level. The reason is that the similarity, commonality, and universality available in high energy collisions [54–61] should be related to partons, which are a deeper level in the structure of matter compared to the particle level.

The present work shows that central collisions correspond to a higher excitation degree and a larger blast due to the more energy being deposited. Meanwhile, central collisions correspond to higher degree of equilibrium due to many particles being produced. Compared with the parameters at the particle level, the parameters at the partonic-level are extracted at earlier time moment and correspond to larger and . However, the system at the particle level is at higher degree of equilibrium due to the fact that longer time is taken to approach the chemical and kinetic freeze-out. However, because the time span from parton phase to particle phase is very small, the difference between two sets of parameters at the partonic- and particle level is very small.

In some cases, the fitting quality is not so good due to large values. This is because only experimental uncertainty is considered. In fact, we have used the Monte Carlo method, which causes additional uncertainty. If we take the uncertainty induced by the Monte Carlo method to be approximately the same as the experimental uncertainty, values will be reduced by %. In the case of large values, we may regard the fits as the qualitative ones. In most cases, the fits are good or approximate, and in few cases the fits are qualitative.

Before summarizing and concluding, we would like to conduct some further discussions on the Monte Carlo results. As we know, most of baryons consist of up and down quarks, and the number of strange quarks is very small. Only when there is a formation of QGP, the system is likely to produce a lot of strange quarks and strange antiquarks, and then, strange quarks and strange antiquarks are able to combine with other neighbouring partons to form strange particles. Thus, the abundant yield of the final-state strange particles is an important signal of the existence of QGP matter [62–67]. Because of the formation of QGP with abundance of strange particles, the implementation of this work has become possible due to sufficient statistics.

In the above discussions, we have used different sets of the parameters to fit the spectra of different strange particles. This means that we have used the multiscenario of kinetic freeze-out. In some cases, the two-scenario of kinetic freeze-out is also applicable, if we consider the single-strange and multistrange particles, respectively. In our opinion, the single-scenario is a rough description, the two-scenario is a slight refined description, and the multiple-scenario is a more refined description, of the process of kinetic freeze-out. This situation is analogous to the atomic spectra and their fine structures.

#### 4. Summary and Conclusions

The transverse momentum spectra of strange particles, including , , , , and so on, in high energy collisions are analyzed by considering the contributions of constituent quarks. Each constituent quark contributes to the transverse momentum obeying the modified Tsallis–Pareto-type function with random azimuthal angle. The transverse momentum of strange particle is the vector superposition of the transverse momenta of the two or three constituent quarks. The results calculated by the Monte Carlo method are in good agreement with the experimental data in most cases in Pb–Pb, –Pb, and Xe–Xe collisions at a few TeV, measured by the ALICE Collaboration at the LHC. In few cases, the agreement is qualitative due to quite large values.

With the decrease of centrality from central to peripheral collisions, the free parameters , , and , as well as the derived parameter decrease slightly in most cases and decrease significantly in few cases. Meanwhile, the free parameter is almost invariant in most cases and decreases slightly in few cases. These results imply that central collisions stay in the state with larger degree of equilibrium, higher excitation, and larger blast than peripheral collisions, though both the central and peripheral collisions stay in approximate equilibrium when the system produces a lot of strange particles.

The correlations between and , and , and , and , and , and and are positive, because all , , and determine mainly, and is mostly defined by the shape of the spectra in intermediate- and high- regions. The correlations between and , and , and , and and are very small or negligible, because determines mainly the spectra in (very) low- region.

We have used the multiscenario of kinetic freeze-out, though in some cases the two-scenario of kinetic freeze-out may also be applicable if we consider a part of region. Meanwhile, if we consider the single-strange and multistrange particles, respectively, the single-strange particles are shown to form earlier than the multistrange ones. We may regard the relations of various scenarios as follows: the single-scenario is a rough description, the two-scenario is a slightly refined description, and the multiple-scenario is a more refined description of the process of kinetic freeze-out.

#### Data Availability

The data used to support the findings of this study are included within the article and are cited at relevant places within the text as references.

#### Ethical Approval

The authors declare that they are in compliance with ethical standards regarding the content of this paper.

#### Disclosure

The funding agencies have no role in the design of the study; in the collection, analysis, or interpretation of the data; in the writing of the manuscript or in the decision to publish the results. A preprint has previously been published [68].

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The work of X.-H.Z. was supported by the Innovative Foundation for Graduate Education in Shanxi University. The work of Shanxi Group was supported by the National Natural Science Foundation of China under Grant Nos. 12147215, 12047571, and 11575103; the Shanxi Provincial Natural Science Foundation under Grant Nos. 202103021224036 and 201901D111043; the Scientific and Technological Innovation Programs of Higher Education Institutions in Shanxi (STIP) under Grant No. 201802017, the Fund for Shanxi “1331 Project” Key Subjects Construction. The work of Kh.K.O. was supported by the Ministry of Innovative Development of the Republic of Uzbekistan within the fundamental project No. F3-20200929146 on analysis of open data on heavy-ion collisions at RHIC and LHC. The work of A.D. was partially supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq-Brazil), by Project INCT-FNA Proc. No. 464 898/2014-5, and by FAPESP under grant 2016/17612-7.