Most hadrons produced in high energy collisions have a relatively low (transverse) momentum perpendicular to the beam axis. The production of soft hadrons is mainly driven by the soft QCD processes, and in particular the non-perturbative hadronization. Experimental and theoretical studies of soft hadron production are important for understanding the properties of the soft parton system created in collisions, and for testing and/or developing phenomenological models. Heavy-ion physics at the SPS, RHIC and LHC energies shows that quark-gluon plasma (QGP) is created in the early stage of collisions. In and/or collisions, it is usually assumed that QGP is not created, at least not up to the RHIC energies. However, recent measurements at the LHC energies show a series of new hadron production phenomena in collisions, such as the ridge and collectivity behavior [1-3], increased baryon-to-meson ratio, and increased strangeness [4-6]. Theoretical studies of these new phenomena mainly focus on how the features of a small partonic system are related to these observations by considering different mechanisms, such as the color re-connection, string overlap and/or color rope [7-10], or by considering the creation of mini-QGP or the phase transition [11-16].

In our latest works [17-21], by studying the available data for the hadronic spectra and yields, we proposed a new understanding of the novel features of hadron production in small quark/parton systems created in and/or -Pb collisions at the LHC energies, i.e. a change of the hadronization mechanism from the traditional fragmentation to quark (re-)combination. In the quark (re-) combination mechanism (QCM), some typical features of the identified hadron production appear, such as the enhanced baryon-to-meson ratio and quark number scaling of hadron elliptical flow at intermediate . These features were observed in the relativistic heavy-ion collisions [22-24], and recently also in and -Pb collisions at the LHC energies in the high multiplicity classes [3, 4, 6, 25]. In particular, a quark number scaling property of the hadron transverse momentum spectra was observed in -Pb collisions at TeV [17].

Recently, the ALICE collaboration reported spectra of the identified hadrons in different multiplicity classes in collisions at 7 TeV [26], and preliminary data for the inelastic events in collisions at 13 TeV [27]. For the first time, a clear signal of the quark number scaling property in the hadronic spectra in collisions is seen. Considering the production of baryons and mesons , the momentum distribution functions in QCM with the equal velocity combination approximation read as

Here, and are coefficients independent of momentum. is the ( ) quark distribution at hadronization, and we assume in the center rapidity region at the LHC energies. With the above two formulas, we get a correlation between the production of and in QCM

where is independent of momentum. In order to check this scaling property, we apply the following operation on the data for and at midrapidity [26]: (i) we divide the bin for ( ) by 3 (2), (ii) take the ( power of the measured for ( ), and (iii) multiply by a constant factor , so that the data points for small ( 0.5 GeV/c) are in coincidence with the scaled data for as much as possible. We show in Fig. 1 the scaled data for and in different multiplicity classes in collisions at 7 TeV. The relative statistical uncertainties of the scaled data are only a few percent, and are shown as rectangles with filled colors in the figure. We see that in the high multiplicity classes, e.g. Fig.1(a) and (b), the scaled data for are consistent with those for , and therefore the quark number scaling property holds. This verifies our argument in the recent work [21], and is a clear signal of quark combination hadronization in collisions at the LHC energies. In the low multiplicity classes, Fig. 1(c) and (d), the scaled data for are somewhat flatter than for at GeV/c, and the quark number scaling property seems to be broken to a certain extent. We note that this is probably due to the threshold effects of strange quark production [21].

In Fig. 2, we show the scaled data for and in collisions at 7 and 13 TeV [26-29] as a guide for energy dependence. We see that the quark number scaling property in the inelastic events in collisions at 7 TeV is broken to a certain extent, but it holds well in the inelastic events in collisions at 13 TeV. This is an indication of the quark combination hadronization at higher collision energies.

On the other hand, we performed a run of event generators Pythia8 [30, 31] and Herwig6.5 as a naive test of the string and cluster fragmentation mechanisms in collisions at TeV. Fig. 3 shows the results for the scaled spectra of and at mid-rapidity given by the two event generators. Here, we adopted the Pythia version 8240 and Herwig version 6521. We chose two event classes, the inelastic non-diffractive events (INEL) and the high multiplicity events with 15, to check the multiplicity dependence of the predictions. In the Pythia8 simulations, we further checked the predictions with the default string fragmentation tune (marked as Pythia8 in Fig. 3), and with the rope hadronization mechanism (marked as Pythia8 rope in Fig. 3). Panels (a)-(c) show the scaled spectra of and where the coefficient is chosen so that the two spectra are coincident at small . Panel (d) shows the ratio of the two scaled spectra. We see that the constituent quark number scaling property given by the two event generators with the current tunes is violated by more than 20% at 1.5 GeV/c.

In this paper, we apply the quark combination model proposed in our recent works [17, 21] to systematically study the production of identified hadrons in collisions at 13 TeV. We calculate the distributions and yields of identified hadrons and focus on various ratios or correlations of the hadronic yields and spectra. We compare our results with the available experimental data to test the quark combination hadronization in collisions at the LHC energies. Predictions are made for future tests.

The paper is organized as follows: Sec. 2 briefly introduces the model of quark (re)combination mechanism with the equal velocity combination approximation. Sec. 3 and Sec. 4 present our results and relevant discussions for the inelastic events and different multiplicity classes. A summary and discussion is given in Sec. 5.

Quark (re-)combination/coalescence mechanism was proposed in the 1970s [32] and has many applications in high energy , and heavy-ion collisions [33-39]. In particular, ultra-relativistic heavy-ion collisions create deconfined hot quark matter in a large volume whose microscopic hadronization process can be naturally described by QCM [40-45]. In this section, we briefly introduce the quark combination model, proposed in previous works [17, 21] in the framework of QCM with the equal velocity combination approximation. We take the constituent quarks and antiquarks as the effective degrees of freedom of the soft parton system created at hadronization. Combinations of constituent quarks and antiquarks with equal velocity result in formation of identified baryons and/or mesons.

The momentum distributions of the identified baryons and mesons are denoted as

Here, and are the momenta of baryons and mesons , and and are the momentum-integrated multiplicities of and , respectively. The superscript means that the distribution function is normalized to one. In the equal velocity combination approximation, also called co-moving approximation, the momentum distributions of baryons and mesons can be simply obtained as a product of the distributions for the constituent quarks and/or antiquarks. We have, for

and for

where is the momentum distribution of quarks normalized to one. and are the normalization coefficients for baryons and mesons , respectively. The momentum fraction is given by, recalling that the momentum ,

where indexes are for baryons, and for mesons. The quark masses are taken to be the constituent masses MeV and MeV.

The multiplicities of baryons and mesons are

Here, is the number of all possible combinations of three quarks related to the formation, and is taken as , and for the case of three different flavors, two identical flavors and three identical flavors, respectively. The factors 6 and 3 are the number of permutations related to the different quark flavors. is the number of all possible pairs related to the formation.

Considering the flavor independence of the strong interaction, we assume that the probability of forming a baryon and the probability of forming a meson are flavor independent. The combination probability can then be written as

Here, denotes the average (or flavor blinding) probability of three quarks combining into a baryon. is the average number of total baryons and is the number of all possible three-quark combinations, with the total quark number. is the probability of selecting the correct discrete quantum number, such as spin, in the formation of as is destined to form a baryon. Similarly, denotes the average probability for a quark and antiquark to combine into a meson, and is the branching ratio for the formation of as is destined to from a meson. is the total meson number, and is the number of all possible quark-antiquark pairs for meson formation.

In this paper, we only consider the ground state mesons and baryons in the flavor SU(3) group. For mesons

where we introduce a parameter , which represents the relative production weight of the vector mesons to the pseudo-scalar mesons with the same flavor composition. For baryons

except that , , . The parameter stands for the relative production weight of the decuplet to the octet baryons of the same flavor content. Here, is taken as 0.45 by fitting the data for the ratio in collisions at 7 TeV and -Pb collisions at 5.02 TeV [46], and is taken as 0.5 by fitting the data for and [47]. The fraction of baryons relative to mesons is for vanishing net-quarks [18, 21, 45]. Using the unitarity constraint for hadronization, , and can be calculated using the above formulas for the given quark numbers at hadronization.

We summarize the main underlying dynamics of the model. The constituent quarks and antiquarks are assumed to be the effective degrees of freedom of the soft parton system at hadronization. Combinations of constituent quarks and antiquarks with equal velocity result in the formation of baryons and mesons. This is similar to the constituent quark model, i.e. the sum of the masses (and momenta) of the constituent quarks is used to construct the mass (and momentum) of a hadron. The model parameters and contain the non-perturbative dynamics and are obtained by fitting the relevant experimental data. They are assumed to be relatively stable in/at different collision systems/energies. Also, the normalization of the hadronization process is a prerequisite for quark combination. The quark number conservation is not only globally satisfied via and , but is also satisfied for each quark flavor via . Here, runs over all hadron species, and . is the number of constituent quarks in a hadron . Therefore, this is a statistical model based on the constituent quark degrees of freedom, and is different from the popular parton recombination/coalescence models [40, 41] which adopt the Wigner wave function method with the instantaneous hadronization approximation.

As the quark number at hadronization is small, the identified hadron production will suffer some threshold effects. For example, baryon production is forbidden for events with . For events with , baryon production is forbidden. In collisions at the LHC energies, the event-averaged number of strange quarks is in the midrapidity region ( ) for the inelastic events and not-too high multiplicity event classes. Therefore, the yield of is no longer completely determined by the average number of strange quarks but is also strongly influenced by the distribution of the strange quark number. The case of , which needs two strange quarks, is similar. We use to denote the distribution of the quark number around the event average, and obtain the averaged multiplicity of identified hadrons by

where is given by Eqs. (9) and (10), and is a function of .

For simplicity, we assume the flavor-independent quark number distribution

where runs over , , flavors. We neglect the fluctuation of net-charges and take in all events. The distribution of and quarks is based on the Poisson distribution . As discussed above, we in particular tune the strange quark distribution. Since for the minimum bias events and small multiplicity classes in collisions , and the Poisson distribution in this case has a long tail for which may over-weight the events with , we distort the Poisson distribution by the suppression factor , i.e. we take , where is the Heaviside step function and is the normalization constant. is taken as 0.8 for the inelastic (INEL > 0) events and the various multiplicity classes.

There are other possible effects of the small number of quarks. For example, in the events with , as and are most likely created from the same vacuum excitation and therefore are not likely to directly constitute a color singlet, the production is suppressed. In addition, as the distributions of quark momenta are dependent on the number of quarks (i.e. on the system size), we neglect such dependence in the given multiplicity classes. Their potential effects will be studied in a future work.

We use the above quark combination model to describe the transverse production of hadrons at midrapidity in collisions.The approximation of equal velocity combination in the model is reduced to the equal transverse-velocity combination. Here, we only study one dimensional distribution of hadrons by further integrating over the azimuthal angle. The distribution functions of quarks at hadronization and at midrapidity are input for the model and are denoted as with . is the number of in the rapidity interval , and is the quark spectrum normalized to one. We assume the iso-spin symmetry between up and down quarks, and also assume the charge conjugation symmetry between a quark and antiquark. Finally, we have only two input functions and , which can be determined by fitting the data for identified hadrons.

Using the scaling property in Eq. (3) and the experimental data shown in Fig. 2, we can directly obtain the normalized distribution of strange quarks at hadronization, which can be parametrized in the form

where is the normalization constant, and the parameters are GeV/c, , , GeV and . By fitting the data for the other hadrons such as protons and , we also obtain the distribution of up/down quarks at hadronization. Taking the parametrization in Eq. (17), the parameters of up/down quarks are GeV/c, , , GeV and . In Fig. 4, we plot and as a function of , and their ratio in the inelastic events in collisions at 13 TeV.

We emphasize that by taking advantage of the quark number scaling property, we can conveniently extract the momentum distributions of soft quarks at hadronization from the experimental data for the hadronic spectra. The extracted quark spectra carry important information about the soft parton system created in collisions at the LHC energies. First of all, since the parameters and in the quark distribution function in Eq. (17) are obviously smaller than one, the extracted and deviate from the Boltzmann distribution in the low range. This indicates that thermalization may be not reached in the small partonic system created in collisions at the LHC energies. Secondly, we see that the ratio , Fig. 4 (b), increases for small and then saturates (or only slightly decreases) with . This property is similar to that observed in collisions at 7 TeV [21] and in -Pb collisions at 5.02 TeV [17], and in heavy-ion collisions at RHIC and the LHC [48-50]. This information about the constituent quarks provides an important constraint for developing more sophisticated theoretical models of soft parton system created in high energy collisions.

Among the hadrons that are often measured by experiments, pions and kaons are the most abundant particles. However, because the pion and kaon masses are significantly smaller than the sum of the masses of their constituent (anti-)quarks, the pion and kaon momenta can not be calculated by a simple combination of the constituent (anti-)quark momenta at hadronization [21]. Therefore, the momentum spectra of pions and kaons are not the most direct probes of the quark combination model, and these results are not shown here. On the other hand, protons, , , , and can be constructed from the constituent quarks and antiquarks. These hadrons can be used to effectively test the quark combination model.

In Fig. 5, we show the calculated results for the spectra of protons, , , , and in the inelastic (INEL > 0) events in collisions at 13 TeV using the quark spectra in Fig. 4 and quark numbers and . The quark numbers are fixed by globally fitting the data for the -integrated yield densities of these hadrons [27]. The solid lines are the QCM results which include the contribution of the strong and electromagnetic decay of the resonances. The symbols are the preliminary data for the hadronic spectra at midrapidity measured by the ALICE collaboration [27]. We see that the data can in general be fitted well by QCM. Our results for are slightly below the data. If we multiply the spectrum by a constant factor, we see that the shape is in good agreement with the data.

Besides the scaling property between the spectra of and shown in the Introduction, the ratio as a function of can also give an intuitive picture of the microscopic mechanism of hadron production. Protons and have similar masses but totally different quark content. In the central (0%-10% centrality) Pb-Pb collisions at TeV, the data for the ratio [51], black squares in Fig. 6, are almost flat with respect to . The flat ratio is related to the similar masses of protons and , and is usually attributed to the strong radial flow and statistical hadronization in the chemical/thermal equilibrium in relativistic heavy-ion collisions. However, the data for in the inelastic events in collisions at 13 TeV [52], solid circles in Fig. 6, show a rapid decrease with increasing . This is an indication of the out-of-thermal equilibrium in collisions. In QCM, the distributions of identified hadrons are determined by the spectra of (anti-)quarks at hadronization. The ratio in QCM reflects the ratio or the correlation between the third power of the quark spectrum and the square of the quark spectrum. With the quark spectra in Fig. 4 , which self-consistently describe the data for the hadronic spectra in Fig. 5, the calculated ratio in QCM, solid line in Fig. 6, shows a decreasing behavior with and is in good agreement with the data for collisions [52].

Hyperons , and contain one, two and three constituent quarks, respectively. Therefore, the ratios and reflect the difference in momentum distributions of ( ) quark and quark at hadronization. Fig. 7 shows our prediction for the ratios and as a function of in the inelastic events in collisions at 13 TeV. We see that the two ratios initally increase with , and then tend to saturate at intermediate GeV/c, which is due to the difference between the spectra of strange quarks and of up/down quarks in the range GeV/c, see the ratio shown in Fig. 4(b).

Using the preliminary data for the spectra of protons, and in different multiplicity classes [52, 53], we can determine the corresponding spectra of constituent quarks at hadronization and predict the spectra of the other identified hadrons. Fig. 8 shows the extracted quark spectra (using the parametrization in Eq. (17)) at midrapidity in different multiplicity classes. Since the parameter in the quark spectrum in high multiplicity classes tends to one, the distribution function in Eq. (17) asymptotically tends to the Boltzmann distribution in the low range, and therefore we see a thermal behavior of the quark spectrum. This is related to the increasing multiple parton interactions in these event classes. In small multiplicity classes, the parameter is relatively small and the quark spectrum deviates from the thermal behavior.

In Fig. 9, we show the -integrated yields of identified hadrons (including kaons)^① in different multiplicity classes and compare them with the preliminary data in collisions at 13 TeV [52, 54]. In general, the results of QCM, solid lines, are in good agreement with the data (with the maximum deviation of about 10%).

The yield ratios of different hadrons can considerably cancel the dependence on the model parameters and/or model input. Therefore, they are a more direct test of the basic physics of the model when confronted with the experimental data. In Fig. 10, we show the yield ratios of hyperons , and to pions divided by the number of inclusive INEL>0 events. The data for collisions at 7 [5] and 13 TeV [52], and for -Pb collisions at 5.02 TeV [55, 56], are presented in order to get a clear tendency with respect to the multiplicity of charged particles at midrapidity. The solid lines are the numerical results of QCM, which are found to be in agreement with the data. We emphasize that the strangeness-related hierarchy is closely related to the strange quark content of these hyperons at hadronization, which can be easily understood using the analytical relation in QCM. Taking the yield formulas Eqs. (9) and (11) and considering the strong and electromagnetic decays, we have

where we neglect the effects of small quark numbers, and adopt the strangeness suppression factor . Due to the contribution of complex decays, the pion yield has a complex expression [57], and here we write with the coefficient almost constant. The double ratios in Fig. 10 then have simple approximate expressions

where is the strangeness suppression factor in the INEL >0 events in collisions. Here, we see a clear hierarchy structure among the three double ratios in terms of . The dotted lines in Fig. 10 are the results of the above analytic formulas with a naively tuned strangeness suppression with . They fit well the experimental data for the double ratios for . In the small multiplicity classes , the small quark number effects are not negligible so that the analytic approximations give larger values than the experimental data. Our numerical results include the small quark number effects and are found to be closer to the data.

The yield ratio is also influenced by the small quark number effects. If we neglect them, we have

and using Eq. (19), we get the ratio

which slightly decreases with increasing , and therefore slightly decreases with an increase of multiplicity because increases with . This is in contradiction with the experimental data. However, taking into account the small quark number effects in QCM, we get the correct behavior of the ratio , as shown by the solid line in Fig. 11. The production of needs not only two quarks but also a quark, which is different from which needs only an and . Therefore, in small multiplicity events or for small quark numbers, the formation of is suppressed to a certain extent (or occasionally forbidden) in comparison with . We see that the calculated ratio using QCM increases with system multiplicity , and is consistent with the experimental data for collisions at 13 TeV and for -Pb collisions at 5.02 TeV [52, 58].

Protons and have similar masses but different quark content. The yield ratio can test the flavor-dominated features of hadron production in QCM. Neglecting the small quark number effects, the proton yield, after taking into account the decay of resonances, has a simple expression

Using Eq. (24), we get the yield ratio

which shows a significant dependence on the strangeness suppression factor . We get with in the low multiplicity classes, and with in the high multiplicity classes. The short dashed lines in Fig. 12 are for the above two values in the analytical approximation. They are slightly lower than the experimental data [52, 58], symbols in the figure. The small quark number effects increase the ratio to a certain extent as they suppress the proton yield. We show the numerical results of our model including the small quark number effects, solid line, and see a good agreement with the data.

In Fig. 13, we show the fit of the data for the spectra of protons, and [52, 53] using QCM, and the predictions for the other identified hadrons in different multiplicity classes for collisions at 13 TeV. Note that the classes IV and V are combined for the data and for our results. Beside comparing the predictions of single hadron spectra with the data, we emphasize that QCM can be more effectively tested by spectrum ratios and/or scaling. The first test is whether the constituent quark number scaling holds in the spectra for and in different multiplicity classes. The second is to study the ratio as a function of . The ratio in QCM is solely determined by the strange quark spectrum at hadronization, and this ratio usually exhibits a nontrivial dependence, as shown in Fig. 14(a), which is a typical behavior of the baryon-to-meson ratio in QCM and is absent or unapparent in the traditional fragmentation picture. We also see that the ratio in higher multiplicity classes can reach higher peak values. The peak position of in high multiplicity classes is also enlarged compared with the low multiplicity classes. The third test is to study the ratio as a function of so as to clarify if its dependence is flavor dominated or mass dominated. The results of QCM are shown in Fig. 14(b) , which decrease with and show a relatively weak multiplicity dependence.

Taking advantage of the available experimental data for the hadronic spectra and yields at midrapidity, we have systematically studied the production of soft hadrons in collisions at 13 TeV in the framework of the quark combination mechanism at hadronization. We applied the quark combination model which assumes that the constituent quarks and antiquarks are the effective degrees of freedom of the parton system at hadronization, and used the equal velocity combination approximation in hadron formation. We applied the model to systematically calculate the spectra and yields of soft strange hadrons in the inelastic events (INEL>0) and in different multiplicity classes.

We found several interesting results which are sensitive to the hadronization mechanism. (1) The spectra of and in the inelastic events (INEL>0) in collisions at 13 TeV exhibit a constituent quark number scaling property. The high multiplicity classes in collisions at = 7 TeV also show this scaling property. This is the first time that such a scaling property is observed in high energy collisions, and is an obvious experimental signal of the quark combination mechanism at hadronization in high energy collisions. (2) The ratio in the inelastic events (INEL > 0) in collisions at 13 TeV shows a clear decrease with increasing , which indicates that the statistical hadronization model is not the reason for this behavior. We demonstrated that the data are naturally explained by the quark combination model. (3) The yield ratios , and divided by the number of inelastic events as a function of system multiplicity at midrapidity show a strangeness-related hierarchy structure. We demonstrated that the hierarchy structure is closely related to the strange quark content of these hyperons as the result of the combination of strange quarks and up/down quarks during their production.

Using the quark number scaling property, the spectrum of strange quarks was extracted from the data for and . The spectra of up/down quarks was extracted from the data for the other hadrons containing up/down constituent quarks. The extracted quark momentum distribution functions are important results which describe the properties of the strongly interacting partonic system at hadronization in the language of constituent quarks.

To confirm the new features of hadronization dynamics in high energy collisions, we should carefully study all related experimental data. We have made predictions of the spectra and spectrum ratios of strange hadrons in collisions at 13 TeV to further test our model with the experimental data that will be available in the future. On the other hand, compared to our previous study of collisions at 7 TeV [21] which gave a first indication, the current study provides a stronger signal of quark combination hadronization in high energy collisions. We still need further systematical studies of collisions at the other LHC energies to test the universality of the new hadronization features, and to study their relation with the possible creation of mini-QGP in small collision systems.