#### Abstract

We study the transverse momentum () spectra of neutral pions and identified charged hadrons produced in proton–proton (pp), deuteron–gold (d–Au), and gold–gold (Au–Au) collisions at the center of mass energy GeV. The study is made in the framework of a multisource thermal model used in the partonic level. It is assumed that the contribution to the -value of any hadron comes from two or three partons with an isotropic distribution of the azimuthal angle. The contribution of each parton to the -value of a given hadron is assumed to obey any one of the standard (Maxwell-Boltzmann, Fermi-Dirac, and Bose-Einstein) distributions with the kinetic freeze-out temperature and average transverse flow velocity. The spectra of the final-state hadrons can be fitted by the superposition of two or three components. The results obtained from our Monte Carlo method are used to fit the experimental results of the PHENIX and STAR Collaborations. The results of the present work serve as a suitable reference baseline for other experiments and simulation studies.

#### 1. Introduction

The transverse momentum () distributions of the identified hadrons in the final state of high-energy collisions reflect the excitation degree of the particle emission source as well as the speed of collective motion of the particles. The excitation degree of particle emission source reflects the speed of thermal motion of the particles. Different distribution functions can be potential candidates for the description of the spectra measured in experiments. The distributions include, but are not necessarily limited to, the standard (Maxwell-Boltzmann, Fermi-Dirac, and Bose-Einstein) distributions obtained from the Boltzmann-Gibbs statistics, the Tsallis distribution obtained from the Tsallis statistics [1–6], the Tsallis form of the standard distribution (or the Tsallis standard distribution), the -dual distribution obtained from the -dual statistics [7], the -dual form of the standard distribution (or the -dual standard distribution), the Erlang distribution obtained from the multisource thermal model [8–12], and the Hagedorn function [13] (inverse power law [14–18]) obtained from the quantum chromodynamics (QCD) calculus and their superposition.

The flow effect may cause a red shift of the spectra. The influence of flow effect is not excluded in the above-mentioned distributions [1, 7, 8, 13, 14]. The temperature parameters used in these distributions are the effective temperature () of particle emission source which is larger than the real temperature one wants to extract from the spectra. Generally, contains the contributions of thermal motion and collective motion (flow effect) which are described by the kinetic freeze-out temperature () and average transverse flow velocity (), respectively. The contributions of both the thermal and collective motions to are not exactly separable, though the magnitudes of the two effects can be approximately calculated by some model methods.

One has at least four methods to separate the two types of motions. In method 1, one may use the blast-wave model with the Boltzmann-Gibbs statistics or Tsallis statistics, in which a determined velocity profile is assumed, and and can be obtained simultaneously [19–22]. The temperature given by the blast-wave model with the Boltzmann-Gibbs statistics is larger than that with the Tsallis statistics. In method 2, one may use the Lorentz-like transformation in a distribution, in which the collective motion is treated as the motion of the reference system and and can be also obtained simultaneously [23–28]. In method 3, one may use the intercept–slope method in which is the intercept in the linear relation of versus the rest mass () of particle and is the slope in the linear relation of the average () versus the average energy (i.e., the average mass () of moving particles in the source rest frame) [29–37]. Here, refers to a given kind of particle in all components in the data sets. In addition, in the linear relation of versus , is related to, but not equal to, the slope [38, 39]. In method 4, one may propose the contribution fractions of thermal and collective motions to to be determinate values which are from some models. For example, empirically, in the hydrodynamic simulations [40], and naturally, one obtains .

The above four methods were used in our previous work [27, 28, 33–37], though only few distributions were performed. The results from different methods are inconsistent in some cases. These inconsistent results appear not only for the magnitudes but also for the tendencies of and with increasing the collision energy and centrality. A robust method should be used to obtain a more reasonable result. The result of the blast-wave model is model dependent. Although the third and fourth methods seem to be model independent, they are hard to connect with the basic physics processes. Considering the fact that the standard distribution is the most basic one, which is from the relativistic ideal gas model in the thermodynamics, one prefers to use it with the Lorentz-like transformation to extract and . Here, the particles (or partons) emitted from the hot and dense system are assumed to obey the law of the relativistic ideal gas model, though the matter contained within the intermediate fireball is known to behave like a strongly interacting fluid of partons. To our best knowledge, this direct extraction method, using the standard distribution, is rarely reported in the literature.

In the framework of multisource thermal model at the quark or gluon level [8–12], one may use the standard distribution with and to describe the behavior of partons. For the production of a given particle of any type, its contributors contain mainly two or three partons with isotropic azimuthal angles. Of course, one does not expect that the single-component nonanalytical description based on the superposition of two or three standard distributions with and by the Monte Carlo method is enough to fit the experimental data. Considering different violent degrees of binary nucleon-nucleon process in high-energy collisions, two or three or even more sets of parameters are possible, which results in a non-single-component distribution. Going from the binary process with the lowest intensity to the most violent one, the fractions of the corresponding components are getting smaller and smaller. This means that the fraction of the first component with smallest and is the largest.

In the case of using the two-component nonanalytical description, the spectrum in high- region is regarded as the result of hard scattering process which implies high- and , and the spectrum in low- region is regarded as the result of soft excitation process which implies low- and . In the case of using the three-component function, one needs one more component, i.e., the intermediate- and for the spectrum in intermediate- region. The multicomponent function corresponds to the multiregion fine structure of spectra [41–43], which is a natural result of the multisource thermal model [8–12] if the standard distributions with different values of parameter are used to describe different components.

In this article, the standard distribution with and will be used to describe the transverse momentum of partons. The transverse momentum of given particle is the sum of contributions of two or three partons. The related calculations are performed in the framework of multisource thermal model [8–12], where the two- or three-component description is available. The calculated results are fitted to the experimental data measured in high-energy collisions by the PHENIX [44] and STAR [45] Collaborations.

The remainder of this article is structured as follows. The picture and formalism of the multisource thermal model are described in Section 2. Results and discussion are given in Section 3. In Section 4, we give our summary and conclusions.

#### 2. Picture of Multisource and Formalism of Multicomponent

According to the multisource thermal model [8–12], one may assume that there are lots of energy sources to form in high-energy collisions. These energy sources can be quarks and/or gluons if one studies the production of particles. For a given particle of any type, its contributors may be generally two (for mesons) or three (for baryons) energy sources of contributor partons [8, 10]. The number of contributor partons is the same as that of constituent quarks of a given hadron. In most of the cases, the contributions of two or three partons are suitable to fit the hadronic spectra. If the two or three partons are not enough in the analysis, one may include the contributions from the fourth or more partons, which corresponds to the hadronic state of multiple quarks [8, 10]. Here, the contributor partons refer to the constituent quarks of identified hadrons. In the case of studying the spectra of leptons, one may consider two contributor partons as the energy sources, in which one is from the projectile and the other is from the target.

In the relativistic ideal gas model, the invariant particle momentum () distribution can be given by [1] where is the energy, is the transverse mass, is the rapidity, is the longitudinal momentum, is the particle number, is the degeneracy factor, is the volume, is the chemical potential, and , 1, and −1 correspond to the Maxwell-Boltzmann, Fermi-Dirac, and Bose-Einstein distributions, respectively.

The density function of momenta is obtained by

The unit-density function of transverse momentum and rapidity is written as [1]

The density function of transverse momentums is where and denote the minimum and maximum rapidities, respectively.

In the near midrapidity region, . In the case of having no collective flow, the transverse momentum of the -th parton contributing to the transverse momentum of the given particle is assumed to obey the probability density function of the standard distribution
where is the number of partons, is the normalization constant, is the transverse mass of the -th parton, is the constituent mass that is 0.31 GeV/c^{2} for up and down quarks [46], and is the chemical potential of the -th parton that is nearly 0 at high energy [47–50]. One can obtain easily from Eq. (5). In addition, if in Eq. (5) is replaced by , where denotes the number of events, one may obtain the mean multiplicity of particles in an event.

It should be noted that Eq. (5) is not the united probability density function of transverse momentum and rapidity (or longitudinal momentum), but only the probability density function of transverse momentum at the midrapidity which is the concerned major region in experiments. From a practical point of view, Eq. (5) is an approximate expression and easy to use. In addition, the constituent mass, but not the current mass, of a given quark is used in Eq. (5) due to the considered quarks being the constituents of the collision system and produced hadrons.

One may introduce the average transverse flow velocity and the Lorentz-like factor [23–28] at the parton level. The quantities and , as well as the transverse mass and transverse momentum containing the flow effect, can be transformed into each other. One has the Lorentz-like transformation where the absolute value is used due to being positive and being possibly negative in low- region. The Lorentz-like factor or transformation is called, instead of the Lorentz factor or transformation is called. The reason is that and , but not and , are used in the analysis.

After the Lorentz-like transformation, the probability density function, , of is transformed into the probability density function, , of . The relation between the two probability density functions is

From Eqs. (5)–(7), one has in which is naturally rewritten as due to the introduction of . As the quantities at the parton level, and may show different tendencies with centrality. The reason is that multiple scattering of secondary particles may happen in the participants and spectators which are centrality dependent.

Equation (8) is obtained from Eq. (5) due to the conversion of probability densities of transverse momentums in which and are considered. Because of Eq. (5) being only for the case of midrapidity, Eq. (8) is deserved for the same case. Meanwhile, the application of and at the parton level, but not and of each parton, can avoid using too many parameters. In fact, the kinetic freeze-out temperature and transverse flow velocity extracted from the data fitting are usually the average quantities.

In the Monte Carlo calculations, let and denote random numbers distributed evenly in [0,1]. A concrete satisfies the relation where denotes the integral variable to differ from the integral upper limit and denotes a small amount shift from . The contributor partons are assumed to move isotropically in the transverse plane. To obtain a discrete azimuthal angle that satisfies the isotropic distribution, we have

In the transverse plane of the rectangular coordinate system, the - and -components of the vector of the parton transverse momentum are respectively. If partons contribute to , the - and -components of the vector of the particle transverse momentum are respectively. Then, we have

In particular, if meaning that two contributor partons taken part in the formation of a particle, we have

If meaning that three contributor partons taken part in the formation of a particle, we have

The spectra of particles can be divided into two or three regions. This means that one needs two- or three-component function to fit the spectra. If the first component describes the spectra in low- region, which corresponds to the contribution of soft excitation process, the last component describes the spectra in high- region which corresponds to the contribution of hard scattering process. Naturally, the intermediate component (in three-component function) describes the spectra in intermediate- region. Generally, at low energy or for narrow spectra, one or two components are needed. At high energy or for wide spectra, three or more components are needed.

In the Monte Carlo calculations, one may obtain the digitized probability density function, , for the contribution of the -th component, where and denote the kinetic freeze-out temperature and average transverse flow velocity corresponding to the -th component. Due to the multicomponent, one has the distribution measured in experiments to be where denote the number of components and denote the contribution fraction of the -th component. The normalization is naturally obeyed.

In the above discussions, from the physics point of view, the origin of multiple sources has two meanings. For a given kind of particle, the multiple sources originate from multiple mechanisms of interactions or different excitation degrees of the system, which results in the multicomponent distribution. If the particles in low-, intermediate-, and high- regions are produced by three different mechanisms or from three excitation degrees, the multiple sources become three sources. For a given particle, the multiple sources refer to multiple energy sources when the particle is formed. Generally, two or three energy sources are considered in the formation of the given particle.

It is noteworthy that the collective motion, which pertains to a common velocity or momentum of the particles, does not lead to the kinetic freeze-out temperature. In fact, the temperature is known to originate only from random thermal motion and reflects the degree of intensity of the thermal motion. In the present work, to obtain the kinetic freeze-out temperature, the transverse flow velocity is introduced to exclude the influence of collective motion. This treatment method is easy to use in the fit of experimental data. If the influence of collective motion is not excluded, i.e., if the transverse flow velocity is not considered, one will obtain the effective temperature which is larger than the kinetic freeze-out temperature.

In the present work, to fit the experimental invariant yield, , one needs to structure the relation, , where is the normalization constant that is generally the area under the data, . In the fit, is determined by the data itself and has no relation to the model. Although does not appear as a parameter in the present work to avoid triviality, one more “1” is subtracted when counting the number of degree of freedom (ndof).

#### 3. Results and Discussion

As an application of the above extraction method based on the description of spectra, the invariant yields, , of neutral pions () produced in mid-pseudorapidity () in gold–gold (Au–Au) collisions with different centrality percentages at the center-of-mass energy per nucleon pair GeV are presented in Figure 1, where is the rapidity as defined above, the pseudorapidity , and denotes the emission angle of the particle. The various symbols represent the experimental data measured by the PHENIX Collaboration [44]. The curves are our results fitted by the three-component function. For clarity, the samples for different centrality percentages are rescaled by different factors as mentioned in the legends in the figure. The values of the parameters , , , , and ndof are listed in Table 1, where , 2, and 3. Here, () is not listed due to the fact that it can be obtained from the normalization.

One can see that the PHENIX data are fitted satisfactorily by the three-component function. Although there is no data in the region of GeV/c in peripheral collisions, one may show an extension of the fitted curves based on the parameters extracted from the data in the region GeV/c. The extension of the fitted curves in peripheral collisions can be compared with those in central and semicentral collisions. The wavy structure in each case is caused by the low statistics in high- region.

Based on Table 1, the dependences of parameters on centrality percentage are also shown in Figure 2. From Figures 2(a)–2(c), the dependences are for , , and , respectively. It is observed that with the decrease in centrality (or with the increase in centrality percentage) from central to peripheral collisions, the parameters studied here do not have significant change and no obvious fluctuations are observed. If there are any fluctuations, they are insignificant. One can say that these parameters are centrality independent. The reason is that the spectra in Figure 1 are very similar, if not equal, in the shape in different centrality intervals in Au–Au collisions. This implies that the kinetic freeze-out parameters extracted from the spectra of are nearly independent of the centrality.

**(a)**

**(b)**

**(c)**

It should be noted that most of the parameters in Table 1 have the same uncertainties across the centrality, though one has seen some differences in the fourth decimal place. Only three decimal places are kept in the table. Because of parameter uncertainties being less than the symbol size, they are not visible in Figure 2.

Table 2 shows the values of correlation coefficients () between two parameters, where the derived parameter is also included. It is observed that the correlations are not present for the case of having no (). For example, there is no significant correlation between and , , and , as well as and , where both and are the sequence numbers of the component. Some correlations involved to are significant (). The reason is that the parameters and are restrained to the special regions and affect the local shapes of spectra, while the parameters affect the spectra in whole range due to the normalization. Concretely, only and (or due to ), as well as and (or ), are highly correlated. This does not give rise to the problem of multicollinearity.

To see the influence of azimuthal angular difference between the two contributor partons, Figure 3 displays a comparison of the results of the isotropic azimuthal angle (the solid curves), the parallel or identical case (, the dashed curves), the vertical case (, the dotted curves), and the special case of (the dot-dashed curves) for three centrality percentages as an example, where the tail parts with several sharp drops in the curves have been cut to avoid confusion. One can see that the shapes of the curves in low- ( GeV/c) and high- ( GeV/c) regions are significantly different. In the intermediate- region, the dashed curve is more similar to the solid curve if one renormalizes the former. This reflects that in the intermediate- region, one may use the result of to replace approximately that of the isotropic azimuthal angle. The cases of and need larger and to cater to that of the isotropic azimuthal angle.

To show the contribution of each component with the isotropic azimuthal angle, Figure 4 displays the contributions of the first (the dashed curves), the second (the dotted curves), and the third (the dot-dashed curves) components together with that of the total three components (the solid curves), where the tail parts with several sharp drops in the curves have been cut to avoid confusion. Naturally, the first component contributes mainly in the low- ( GeV/c) region. The second component contributes mainly in the low- and intermediate- ( GeV/c) regions. The third component contributes in the whole region. Although the second and third components contribute in wider region than the first one, the most important contributor is the first component with calculated from Table 1.

Due to very small amounts for (here and ), the second and third components do not affect the values of () and () significantly. One has and . In fact, including and causes the increase of and to be less than 0.5%. One may say that although there are 8 free parameters in the fit for the data sets in Figure 1, the main parameters are only and which represent reasonably and .

In the data analysis, for the purpose of the extraction of kinetic freeze-out parameters, one does not need too wide spectra. Generally, the range of GeV/c (even GeV/c) is enough, though the range is 0~20 GeV/c in some cases in Figures 1, 3, and 4. If the spectrum in low- region is contributed by the soft excitation process and that in high- region is contributed by the hard scattering process, the present work shows that only the contribution of soft process is enough to extract the kinetic freeze-out parameters. The contribution fraction of the hard process is smaller than 0.2%.

The value of is extracted together with at the parton level. Due to the introduction of , becomes smaller than the effective temperature because the flow effect is excluded. The introduction of also affects the tendency of with changing the centrality and other factors. At both the parton and particle levels, the obtained or may be inconsistent in the tendency and/or size due to different extraction methods. One should have a determined and uniform method to extract and . As a tentative work, the present work focuses on the parameters at the parton level. Meanwhile, the standard distribution used in the relativistic ideal gas model is applied for the system. In the aspect for extracting the parameters at the kinetic freeze-out, this work has an important significance in methodology.

Many data sets have been analyzed in our previous work [33–37] by the thermal-related models, though only the spectra of in Au–Au collisions at GeV are mainly analyzed in this work as an example to check the validity of the model in methodology. In our recent work [51], using the blast-wave model with fluctuations [52–54], the kinetic freeze-out parameters extracted from the spectra of charged pions () and kaons () produced in lead–lead (Pb–Pb) collisions at TeV [55], proton–lead (p–Pb) collisions at TeV [56], and xenon–xenon (Xe–Xe) collisions at TeV [57] show almost the independence of centrality, and the parameters extracted from the spectra of antiprotons and protons show the dependence of centrality.

The difference between the parameters from the spectra of mesons and is caused by different production mechanisms. Generally, mesons are newly produced and some protons already exist in the projectile and target nuclei before the collisions. In our opinion, the extracted parameter values for production are reduced by the influence of these pre-existing protons, which are the leading particles with low excitation in the collisions. They increase relatively the yield in low- region. The relative increase is more obvious in peripheral collisions due to lesser multiple scattering.

It is expected that for a narrow spectrum in low- region, a single- or two-component function is suitable; for a not too wide spectrum, the three-component function is suitable; and for a wider spectrum, the fourth component may be included if the three components are not enough. To check the methodology by the spectra in a narrow region, for which only one component is used in the fit, Figures 5(a)–7(a) show the spectra of , , and , and Figures 5(b)–7(b) show the spectra of , , and , produced in Au–Au, deuteron–gold (d–Au), and proton–proton (pp) collisions at GeV, where is simplified to for pp collisions. The symbols represent the experimental data measured by the STAR Collaboration [45]. The curves are our fitted results by a single-component function in which only and are free parameters which are listed in Table 3. For the productions of and , two contributor partons contribute to . For the production of , three contributor partons contribute to . One can see that the single-component function can fit the spectra in the narrow region. Note that some of the /ndof values in Table 3 are very small. This is because the combined errors in spectra are predominantly determined by large systematic uncertainties, which do not follow a normal distribution. The statistical uncertainties, which have a random nature, are very small and negligible.

The parameter values listed in Table 3 show that pp collisions are similar to peripheral d–Au and Au–Au collisions at the same . Central, semicentral, and peripheral d–Au collisions are similar to peripheral Au–Au collisions. Larger values of parameters in central Au–Au collisions imply that larger amounts of collision energy are deposited in the system and more multiple scatterings occur in the hot and dense matter.

Table 3 also shows that and decrease with the decrease in centrality from central to peripheral events and decrease with the decrease in particle mass for the considered three kinds of particles, in Au–Au and d–Au collisions at GeV. Figure 8 shows the dependence of and on centrality percentage in Au–Au collisions. One can see clearly the dependence of the parameters on centrality and particle kind. Comparing with that for , the dependence for charged massive particles ( and ) is significant and the distributions for are almost flat, also showing no dependence on centrality. The reason is that charged massive particles have more probability to interact with spectator nucleons, which causes a large red shift of the spectra in peripheral collisions.

**(a)**

**(b)**

Generally, and reflect the excitation and expansion degrees of the system, respectively. The larger these parameters are, the higher the excitation and expansion degrees of the system is. In the case of charged massive particles, central collisions correspond to higher excitation and expansion degrees due to larger amounts of collision energy being deposited. Meanwhile, charged massive particles have relative large elastic scattering cross-section due to their large sizes, while neutral and charged pions have very small elastic scattering cross-section in the system, which results their parameters showing almost independent of centrality.

The reason why one may apply the same distribution to the energy sources with different degrees of excitation and expansion in high-energy collision systems is because of the similarity, commonality, and universality, especially the universality, in the collisions [58–65]. In particular, in high-energy collisions, the underlying reason is the contributor partons appearing as the energy sources or influence factors. This also explains the consistency of some quantities in high-energy collision systems with different sizes and centralities. Some model independent dependencies, if available, are mainly caused by the effective energies of the contributor partons or energy sources.

In the above discussion, the classical concepts of temperature and equilibrium are tentatively used. However, the collision system is very small. In particular, only two or three contributor partons are considered in the formation of given particle. It seems that the mentioned concepts are not applicable. In fact, although few partons are considered to contribute directly and mainly to particle’s , lots of partons exist in the system of high-energy collisions. In addition, the experiments of high-energy collisions involve lots of events. One may use the grand canonical ensemble to study the particle production. Then, the concepts of temperature and equilibrium are applicable. At least, one may regard and as parameters that describe the average kinetic energies of thermal and collective motions of partons, respectively.

Before summary and conclusions, we would like to emphasize that although many parameters are used in this article, this is only for the wide spectra with multiregion structure. Indeed, mathematically, increasing the number of parameters increases the probability of a good fit as one has more free parameters to play around. Technically, one may not completely discard the usage because there are different particle production mechanisms in different regions. This work uses the standard distribution for a given parton source and the multicomponent distribution for the parton sources with different excitation degrees. In the case of extracting and , only two parameters and are enough. The distribution with and () describes wider spectra than that with only (where ). Because this work is based on the standard distribution which is widely used in statistical analysis in modern physics, the results are suitable to be the baseline for comparing with other experiments and simulation studies.

In addition, it has been established for almost 20 years that the high- spectra of hadrons in nuclear collisions are explained by jet quenching; i.e., high- partons lose energy and then fragment to hadrons. This work shows an alternative and uniform explanation for the statistical behavior of particle spectra in various regions. There is no contradiction between the two explanations. If the explanation of jet quenching is focused on the production mechanism, the present work is focused on the statistical law obeyed by the contributor partons and produced particles. Even if in the high- region, the standard distribution is applicable to extract the kinetic freeze-out parameters. Although this results in quite a large and (~0.365–0.435 GeV) in the fitting for all centralities, the value of weighted by is small. Here, the weighted factor is . is determined by that is extracted from the pion spectra to be ~0.132–0.175 GeV in central Au–Au collisions, ~0.123 GeV in central d–Au collisions, and 0.116 GeV in pp collisions at GeV. It is likely that the unexpected large and are obtained by the hard process due to violent head on collisions between partons, though this probability is very small.

Different methods or functions used in the extractions of temperature and flow velocity are different “thermometers” and “speedometers.” Although the tendencies of parameters based on different methods are almost the same or approximately the same, there are differences in concrete values. Obviously, before giving a comparison, these thermometers and speedometers should be uniformed according to the selected baseline. In our opinion, the standard distribution is a good candidate to be the baseline. The fact that and for pion production are almost independent of centrality is an indicator that there is something more universal than bulk medium flow that governs the physics in different regions. The underlying contributors are partons, but not nucleons, in various collisions at high energy. This also implies the similarity, commonality, and universality, in particular universality, in high-energy collisions [58–65].

#### 4. Summary and Conclusions

In the framework of multisource thermal model used in the parton level, the transverse momentum spectra of the final-state neutral pions and identified charged hadrons produced in mid-(pseudo)rapidity region in Au–Au and d–Au collisions with various centralities and in pp collisions at GeV have been studied. For a given particle of any type, its contributors may be two or three partons with isotropic azimuthal angle distribution. The contribution of each parton to transverse momentum of the hadron is assumed to obey the standard distribution with given kinetic freeze-out temperature and average transverse flow velocity. The transverse momentum spectra of the final-state hadrons can be fitted by the superposition of two or three components. The number of components is related to the width of transverse momentum spectra.

The results calculated by the Monte Carlo method fit satisfactorily the experimental data measured by the PHENIX and STAR Collaborations. With the decrease in centrality from central to peripheral collisions, the kinetic freeze-out temperature and average transverse flow velocity for each component in pion production do not change significantly. Due to the very small contribution fractions of the second and third components, the main parameters are determined by the first component in the low transverse momentum region. The result corresponding to the isotropic azimuthal angles is similar to that of the identical azimuthal angles. The kinetic freeze-out parameters decrease with the decrease in centrality for the production of the charged massive hadrons. The work based on the standard distribution is suitable to be the baseline in comparing with other experiments and simulation studies.

#### 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 in arXiv [66]. This paper applies the multisource thermal model, which has also been applied in our recent work [67, 68]. As a result, some similar statements are inevitably applied in this paper.

#### Conflicts of Interest

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

#### Acknowledgments

The work of Shanxi Group was supported by the National Natural Science Foundation of China under Grant No. 12147215, the Shanxi Provincial Natural Science Foundation under Grant No. 202103021224036, and the Fund for Shanxi “1331 Project” Key Subjects Construction. The work of K.K.O. was supported by the Agency of Innovative Development under the Ministry of Higher Education, Science and Innovations 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.