Constituent quark number scaling from strange hadron spectra in $pp$ collisions at $\sqrt{s}=$ 13 TeV

We show that the data of $p_{T}$ spectra of $\Omega^{-}$ and $\phi$ at midrapidity in inelastic events in $pp$ collisions at $\sqrt{s}=$ 13 TeV exhibit a constituent quark number scaling property, which is a clear signal of quark combination mechanism at hadronization. We use a quark combination model under equal velocity combination approximation to systematically study the production of identified hadrons in $pp$ collisions at $\sqrt{s}$= 13 TeV. The midrapidity data of $p_{T}$ spectra of proton, $\Lambda$, $\Xi^{-}$, $\Omega^{-}$, $\phi$ and $K^{*}$ in inelastic events are simultaneously well fitted by the model. The data of multiplicity dependency of yields of these hadrons are also well understood. The strong $p_{T}$ dependence for data of $p/\phi$ ratio is well explained by the model, which further suggests that the production of two hadrons with similar masses is determined by their quark contents at hadronization. $p_{T}$ spectra of strange hadrons at midrapidity in different multiplicity classes in $pp$ collisions at $\sqrt{s}=$ 13 TeV are predicted to further test the model in the future. The midrapidity $p_{T}$ spectra of soft ($p_T<2$ GeV/c) strange quark and up/down quark at hadronization in $pp$ collisions at $\sqrt{s}=$ 13 TeV are extracted.


I. INTRODUCTION
Most of hadrons produced in high energy collisions are of relatively low (transverse) momentum perpendicular to beam axis. Production of soft hadrons is mainly driven by soft QCD processes and, in particular, nonperturbative hadronization. Experimental and theoretical studies of soft hadron production are important to understand the property of the soft parton system created in collisions and to test and/or develop existing phenomenological models. Heavy-ion physics at SPS, RHIC and LHC energies shows the creation of quark-gluon plasma (QGP) in early stage of collisions. In elementary pp and/or pp collisions, it is usually presumed that QGP is not created, at least up to RHIC energies. However, recent measurements at LHC energies show a series of highlights of hadron production in pp collisions such as ridge and collectivity behaviors [1][2][3], the increased baryon-to-meson ratio and the increased strangeness [4][5][6]. Theoretical studies of these new phenomena mainly focus on how to build the new feature of small partonic system relating to these observations by considering different mechanisms such as the color re-connection, string overlap and/or color rope [7][8][9][10], or by considering the creation of mini-QGP or phase transition [11][12][13][14][15][16], etc.
In our latest works [17][18][19][20][21] by studying the available data of hadronic p T spectrum and yield, we proposed a new understanding for novel features of the hadron production in the small quark/parton systems created in pp and/or p-Pb collisions at LHC energies, that is, the change of hadronization mechanism from the traditional fragmentation to the quark (re-)combination. In quark (re-)combination mechanism (QCM), there exists some typical behaviors for the production of identified * shaofl@mail.sdu.edu.cn † songjun2011@jnxy.edu.cn hadrons such as the enhanced baryon-to-meson ratio and quark number scaling of hadron elliptical flow at intermediate p T . These behaviors are already observed in relativistic heavy-ion collisions [22][23][24] and, recently, are also observed in pp and p-Pb collisions at LHC energies in high multiplicity classes [3,4,6,25]. In particular, a quark number scaling property for hadron transverse momentum spectra is first observed in p-Pb collisions at √ s N N = 5.02 TeV [17].  Figure 1. The scaling property for the dN/dpT dy data of Ω − and φ at midrapidity in different multiplicity classes in pp collisions at √ s = 7 TeV. The coefficient κ φ,Ω in four multiplicity classes is taken to be (1.76, 1.82, 1.83, 1.93), respectively. Data of Ω − and φ are taken from Ref. [26].
Recently, ALICE collaboration reported the data of p T spectra of identified hadrons in different multiplicity classes in pp collisions at √ s = 7 TeV [26] and the preliminary data in inelastic events in pp collisions at √ s = 13 TeV [27]. Here, we find, for the first time, a clear signal of the quark number scaling property for hadronic p T spectra in pp collisions. Considering the production of baryon Ω − (sss) and meson φ (ss), their momentum distribution functions f (p T ) ≡ dN/dp T in QCM under equal velocity combination approximation read as Here, κ φ and κ Ω are coefficients independent of momentum. f s,s (p T ) is s (s) quark distribution at hadronization and we assume f s (p T ) = fs (p T ) in the center rapidity region at LHC energies. With above two formulas, we get a production correlation between Ω − and φ in QCM where κ φ,Ω = κ is independent of momentum. In order to check this scaling property, we take following operation on the dN/dp T dy data of Ω − and φ at midrapidity [26]: (i) divide the p T bin of Ω − (φ) by 3 (2) and (ii) take the 1/3 (1/2) power of the measured dN/dp T dy for Ω − (φ), and (iii) multiply (dN Ω /dp T dy) 1/3 by a constant factor κ φ,Ω so that data points at small p T (p T 0.5 GeV/c) are in coincidence with the scaled data of φ as much as possible. We show in Fig.1 the scaled data of Ω − and φ in different multiplicity classes in pp collisions at √ s = 7 TeV. The relative statistical uncertainties of the scaled data are only a few percentages 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 of Ω − are consistent very well with those of φ, and therefore the quark number scaling property holds well. This verifies our argument in the recent work [21] and is a clear signal of quark combination hadronization in pp collisions at LHC. In the low multiplicity classes, Fig. 1(c) and (d), the scaled data of Ω − are somewhat flatter than those of φ as p T 1 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 of Ω − and φ in pp collisions at √ s = 7 and 13 TeV [26][27][28][29] as a guide of energy dependence. We see that the quark number scaling property in inelastic events in pp collisions at √ s = 7 TeV is broken to a certain extent but it is well fulfilled in inelastic events in pp collisions at √ s = 13 TeV. This is an indication of the quark combination hadronization at higher collision energies.  Figure 2. The scaling property for the dN/dpT dy data of Ω − and φ at midrapidity in inelastic events in pp collisions at √ s = 7 and 13 TeV. The coefficient κ φ,Ω is taken to be (2.0, 1.5), respectively. Data of Ω − and φ are taken from Refs. [26,27].
On the other hand, we run event generators Pythia8 [30,31] and Herwig6.5 as the naive test of the prediction of string and cluster fragmentation mechanism in pp collisions at √ s = 13 TeV. Fig. 3 shows results of the scaled p T spectra of Ω and φ at mid-rapidity in two event generators. Here we adopt Pythia version 8240 and Herwig version 6521 to simulate. We choose two event classes, inelastic non-diffractive events (INEL) and high multiplicity events with dN ch /dy ≥ 15, to check the multiplicity dependence of the prediction. In Pythia8 simulations, we further check the prediction with default string fragmentation tune (marked as Pythia8 in Fig. 3) and that with rope hadronization mechanism (marked as Pythia8 rope in Fig. 3). Panels (a)-(c) show the scaled spectra of Ω and φ where coefficient κ is chosen so that two spectra are coincident with each other at small p T . Panel (d) shows the ratio of two scaled spectra. We see that the constituent quark number scaling property in two event generators with current tunes is violated by more than 20% at p T 1.5GeV/c.   In this paper, we apply a specific quark combination model proposed in our recent works [17,21] to systematically study the production of identified hadrons in pp collisions at √ s = 13 TeV. We mainly calculate the p T distributions and yields of identified hadrons and focus on various ratios or correlations for hadronic yields and p T spectra. We compare our results with the available experimental data to systematically test the quark combination hadronization in pp collisions at LHC energies. Predictions are made for the further test in future.
The paper is organized as follows: Sec. II briefly introduces a specific model of quark (re)combination mechanism under equal velocity combination approximation. Sec. III and Sec. IV present our results and relevant discussions in inelastic events and different multiplicity classes, respectively. A summary and discussion is given at last in Sec. V.

II. QUARK COMBINATION MODEL UNDER EQUAL VELOCITY COMBINATION APPROXIMATION
Quark (re-)combination/coalescence mechanism was proposed in 1970s [32] and has many applications both in elementary e − e + , pp and heavy-ion collisions [33][34][35][36][37][38][39]. In particular, ultra-relativistic heavy-ion collisions create the deconfined hot quark matter of large volume whose microscopic hadronization process can be described by QCM naturally [40][41][42][43][44][45]. In this section, we briefly introduce a quark combination model proposed in previous works [17,21] within QCM framework under the equal velocity combination approximation. We take the constituent quarks and antiquarks as the effective degrees of freedom of the soft parton system created in collisions just at hadronization. The combination of these constituent quarks and antiquarks with equal velocity forms the identified baryons and/or mesons.

A. Hadron production at given numbers of quarks and antiquarks
The momentum distributions of identified baryon and meson are denoted as Here p B and p M are momenta of baryon B i and meson M i , respectively. N Bi and N Mi are the momentumintegrated multiplicities of B i and M i , respectively. The superscript (n) denotes the distribution function is normalized to one. Under the equal velocity combination approximation, also called comoving approximation, momentum distributions of baryon and meson can be simply obtained by the product of those of constituent quarks and/or antiquarks. We have, for where f (n) q (p) ≡ dn q /dp is the momentum distribution of quarks normalized to one.
q2 (x 2 p) are normalization coefficients for baryon B i and meson M i , respectively. Momentum fractions x are given by recalling momentum p = mγv ∝ m, where indexes i, j = 1, 2, 3 for baryon and i, j = 1, 2 for meson. Quark masses are taken to be constituent masses m s = 500 MeV and m u = m d = 330 MeV.

Multiplicities of baryon and meson are
Here N q1q2q3 is the number of all possible three quark combinations relating to B i formation and is taken to be for cases of three different flavors, two identical flavor and three identical flavor, respectively. Factors 6 and 3 are numbers of permutation relating to different quark flavors. N q1q2 = N q1 Nq 2 is the number of all possible q 1q2 pairs relating to M i formation.
Considering the flavor independence of strong interaction, we assume the probability of q 1 q 2 q 3 forming a baryon and the probability of q 1q2 forming a meson are flavor independent, the combination probability can be written as Here N B /N qqq denotes the average (or flavor blinding) probability of three quarks combining into a baryon. N B is the average number of total baryons and N qqq = N q (N q − 1)(N q − 2) is the number of all possible three quark combinations with N q = f N f the total quark number. C Bi is the probability of selecting the correct discrete quantum number such as spin relating to B i as q 1 q 2 q 3 is destined to form a baryon. Similarly, N M /N qq approximately denotes the average probability of a quark and an antiquark combining into a meson and C Mi is the branch ratio to M i as q 1q2 is destined to from a meson. N M is total meson number and N qq = N q Nq is the number of all possible quark antiquark pairs for meson formation.
In this paper, we only consider the ground state J P = 0 − , 1 − mesons and J P = (1/2) + , (3/2) + baryons in flavor SU(3) group. For mesons where we introduce a parameter R V /P which represents the relative production weight of the J P = 1 − vector mesons to the J P = 0 − pseudo-scalar mesons of the same flavor composition. For baryons The parameter R D/O stands for the relative production weight of the J P = (3/2) + decuplet to the J P = (1/2) + octet baryons of the same flavor content. Here, R V /P is taken to be 0.45 by fitting the data of K * /K ratios in pp collisions at √ s = 7 TeV and p-Pb collisions at √ s N N = 5.02 TeV [46] and R D/O is taken to be 0.5 by fitting the data of Ξ * /Ξ and Σ * /Λ [47]. The fraction of baryons relative to mesons is N B /N M ≈ 0.085 at vanishing netquarks [18,21,45]. Using the unitarity constraint of hadronization N M + 3N B = N q , N Bi and N Mi can be calculated using above formulas for the given quark numbers at hadronization.
We summarize the main underlying dynamics of the model. Constituent quarks and antiquarks are assumed to be effective degrees of freedom of soft parton system at hadronization. The combination of these constituent quarks and antiquarks with equal velocity forms baryons and mesons. This is similar to constituent quark model, i.e., the summation of masses (and momenta in motion) of constituent quarks properly constructs the mass (and momentum in motion) of hadron. Model parameters R V /P and R D/O contain unclear non-perturbative dynamics and are obtained by fitting the relevant experimental data and are assumed to be relatively stable in/at different collision systems/energies. Also, the normalization of hadronization process is a prerequisite for quark combination. Quark number conservation is not only globally satisfied via N M + 3N B = N q and N M + 3NB = Nq but also satisfied for each quark flavor via h n qi,h N h = N qi . Here h runs over all hadron species and q i = d, u, s,d,ū,s. n qi,h is the number of constituent quark q i in hadron h. Therefore, this model is a statistical model based on constituent quark degrees of freedom and is different from those popular parton recombination/coalescence models [40,41] which adopt the Wigner wave function method under instantaneous hadronization approximation.
B. Quark number fluctuations and some threshold effects of hadron production As quark numbers at hadronization are small, identified hadron production will suffer some threshold effects. For example, baryon production is forbidden for events with N q < 3. For events with N s < 3, Ω − baryon production is forbidden. In pp collisions at LHC energies, the event-averaged number of strange quarks N s 1 in midrapidity region (|y| < 0.5) for 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 strongly influenced by the distribution of strange quark number. Similar case for Ξ which needs two strange quarks, etc. Here we use P ({N qi } , { N qi }) to denote the distribution of quark numbers around the event average, and obtain the averaged multiplicity of identified hadrons by where N hi is given by Eqs. (9) and (10) and is the function of {N qi }.
For simplicity, we assume flavor-independent quark number distribution where f runs over u, d, s flavors. Here we neglect the fluctuation of net-charges and take N f = Nf in each events. The distribution of u and d quarks is based on the Poisson distribution P oi N u(d) , N u(d) . As aforementioned discussion, we particularly tune strange quark distribution. Because in minimum bias events and small multiplicity classes in pp collisions N s 1 and Poisson distribution P oi (N s , N s ) in this case has a long tail as N s ≥ 3 which may over-weight the events with N s ≥ 3, we distort the Poisson distribution by a suppression factor γ s , that is, is the Heaviside step function and N is normalization constant. γ s is taken to be 0.8 in inelastic (INEL>0) events and various multiplicity classes.
There are other possible effects of small quark numbers. For example, in events with N s = Ns = 1, because s ands are most likely created from the same one vacuum excitation and therefore they are not likely to directly constitute color singlet and therefore φ production in these events is suppressed. In addition, quark momentum distributions are more or less dependent on the quark numbers (in other words, system size) and we neglect such dependence for the given multiplicity classes and its potential effects are studied in the future works.

III. RESULTS IN INELASTIC EVENTS
We apply the above quark combination model to describe the transverse production of hadrons at midrapidity in pp collisions.The approximation of equal velocity combination in the model is reduced to that of equal transverse-velocity combination. Here, we only study one dimensional p T distribution of hadrons by further integrating over the azimuthal angle. The p T distribution functions of quarks at hadronization at midrapidity are inputs of the model and are denoted as f qi (p T ) = N qi f (n) qi (p T ) with q i = d, u, s,d,ū,s. N qi is the number of q i in rapidity interval |y| < 0.5 and f (n) qi (p T ) ≡ dn qi /dp T is quark p T spectrum normalized to one. We assume the iso-spin symmetry between up and down quarks and assume the charge conjugation symmetry between quark and antiquark. Finally we have only two inputs f u (p T ) and f s (p T ) which can be fixed by fitting the data of identified hadrons.

A. Quark pT distribution at hadronization
Using the scaling property in Eq. (3) and experimental data shown in Fig. 2, we can directly obtain the normalized p T distribution of strange quarks at hadronization, which can be parameterized in the form  We emphasize that, by taking advantage of the quark number scaling property, this is the first time we can conveniently extract the momentum distributions of soft quarks at hadronization from the experimental data of hadronic p T spectra. The extracted quark p T spectra carry important information of soft parton system created in pp collisions at LHC energies. First, because of parameters b s and b u in quark distribution function in Eq. (17) are obviously smaller than one, the extracted f (n) u (p T ) and f (n) s (p T ) deviate from Boltzmann distribution in the low p T range. This indicates that thermalization may be not reached for the small partonic system created in pp collisions at LHC. Second, we see that the ratio f Fig. 4 (b), increases at small p T and then saturates (only slightly decreases) with p T . This property is similar to the observation in pp collisions at √ s = 7 TeV [21] and in p-Pb collisions at √ s N N = 5.02 TeV [17] and is also similar to the observation in heavyion collisions at RHIC and LHC [48][49][50]. These information of constituent quarks provides important constraint of developing more sophisticated theoretical models of soft parton system created in high energy collisions.

B. pT spectra of identified hadrons
Among hadrons that are often measured by experiments, pion and kaon are most abundant particles. However, because the masses of pion and kaon are significantly smaller than the summed masses of their constituent (anti-)quarks, the momenta of pion and kaon can not be calculated by the simple combination of those of constituent (anti-)quarks at hadronization [21]. Therefore, momentum spectra of pion and kaon are not the most direct probe of the quark combination model and their results are not shown in this paper. On the other hand, proton, Λ, Ξ − , Ω − , φ and K * 0 can be well constructed by the constituent quarks and antiquarks. These hadrons can be used to effectively test the quark combination model.
In Fig. 5, we show the calculation results of p T spectra of proton, Λ, Ξ − , Ω − , φ and K * 0 in inelastic (INEL>0) events in pp collisions at √ s = 13 TeV using the quark spectra in Fig. 4 and quark numbers N u = 2.8 and N s = 0.86. Here, quark numbers are fixed by globally fitting data of p T -integrated yield densities of these hadrons [27]. Solid lines are QCM results which have included the contribution of strong and electromagnetic decay of resonances. Symbols are preliminary data of hadronic p T spectra at midrapidity measured ALICE collaboration [27]. We see that the data can be well fitted in general by QCM. Our result of K * 0 is slightly smaller than the data. If we multiply the result of K * 0 spectrum by a constant factor and we will see the shape is in good agreement with the data.
Besides the scaling property between the p T spectra of Ω − and φ shown in the introduction (Sec. I), (p +p) /φ ratio as the function of p T can also give an intuitive Solid lines are results of QCM and symbols are preliminary data of ALICE collaboration [27].
picture of microscopic mechanism of hadron production. Proton and φ have similar masses but totally different quark contents. In the central (0-10% centrality) Pb-Pb collisions at √ s N N = 2.76 TeV, the data of (p +p) /φ ratio [51], black squares in Fig. 6, are almost flat with respect to p T . This flat ratio is usually related to similar masses of proton and φ and is usually attributed to the strong radial flow and statistical hadronization under chemical/thermal equilibrium in relativistic heavyion collisions. However, data of (p +p) /φ in inelastic events in pp collisions at √ s = 13 TeV [52], solid circles in Fig. 6, show a rapid decrease with the increasing p T . This is an indication of out of thermal equilibrium in pp collisions. In QCM, p T distributions of identified hadrons are determined by p T spectra of (anti-)quarks at hadronization. (p +p) /φ ratio in QCM reflects the ratio or correlation between the third power of u quark spectrum and the square of s quark spectrum. With quark spectra in Fig. 4 which self-consistently describe data of hadronic p T spectra in Fig. 5, the calculated p/φ ratio in QCM, the solid line in Fig. 6, shows a decreasing behavior with p T and is in good agreement with the data of pp collisions [52].
Hyperons Λ, Ξ − and Ω − contain one, two and three s constituent quarks, respectively. Therefore, ratios Ξ − /Λ and Ω − /Ξ − can reflect the difference in momentum distribution between u(d) quark and s quark at hadronization. Fig. 7 shows our prediction of ratios Ξ − /Λ and Ω − /Ξ − as the function of p T in inelastic events in pp collisions at √ s = 13 TeV. We see that two ratios increase with p T and then tend to be saturate at intermediate p T ∼ 6 GeV/c, which is directly due to the difference be-

IV. RESULTS IN DIFFERENT MULTIPLICITY CLASSES
Using the preliminary data of p T spectra of proton, K * 0 and φ in different multiplicity classes [52,53], we can determine the corresponding p T spectra of constituent quarks at hadronization and predict p T spectra of other identified hadrons. Fig. 8 shows the extracted quark p T spectra (using the parameterized form Eq. (17)) at midrapidity in different multiplicity classes. Because the parameter b q of quark spectrum in high multiplicity classes tends to one, the distribution function in Eq. (17) will asymptotically tend to Boltzmann distribution in low p T range and therefore we see a thermal behavior for quark spectrum. This is related to the increasing multiple parton interactions in these event classes. In small multiplicity classes, parameter b q of quark spectrum is relatively small and the quark spectrum deviates from thermal behavior.

A. Hadronic yields and yield ratios
In Fig. 9, we show the p T -integrated yields of identified hadrons (including kaon) [54] in different multiplicity classes and compare them with the preliminary data in pp collisions at √ s = 13 TeV [52,55]. In general, results of QCM, solid lines, are in good agreement with the data (with maximum deviation about 10%).
Yield ratios of different hadrons can significantly cancel the dependence of model parameters and/or model inputs. Therefore, they are more direct test of the basic physics of the model in confronting with the experimental data. In Fig. 10, we show the yield ratios of hyperons Ω − , Ξ − and Λ to pions divided by the values in the inclusive INEL>0 events. Data of pp collisions at √ s = 7 [5] and 13 TeV [52] and those of p-Pb collisions at √ s N N = 5.02 TeV [56,57] are all presented in order to get a visible tendency with respect to multiplicity of charged particles at midrapidity. Solid lines are numerical results of QCM, which are found to be in agreement with the data. We emphasize that such strangeness-related hierarchy behavior is closely related to the strange quark content of these hyperons during their production at hadronization, which can be understood easily via an analytical relation in QCM. Taking yield formulas Eqs. (9), (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 λ s = N s / N u . Because of complex decay contributions, pion yield has a complex expression [58] and here we can write N π = a π N q with coefficient a π being almost a constant. Then the double ratios in Fig. 10 have simple approximate expressions  Solid lines are results of QCM and symbols are preliminary data of ALICE collaboration [52,55].
The dotted lines in Fig. 10 are results of above analytic formulas with a naively tuned strangeness suppression λ s = λ ′ s 1 + 0.165 log dN ch /dη |η|<0.5 /6.0 with λ ′ s = 0.31. They can well fit the experimental data of these double ratios as dN ch /dη |η|<0. 5 10. In small multiplicity classes dN ch /dη |η|<0. 5 6 small quark number effects are not negligible and the analytic approximation is larger than the experimental data to a certain extent. Our numerical results have included small quark number effects and are found to be closer to the data.
Yield ratio Ξ/φ is also influenced by the small quark number effects. If we neglect small quark number effects, we have and using Eq. (19) we get the ratio which slightly decreases with the increase of λ s and therefore will slightly decrease with the increase of multiplicity dN ch /dη because λ s increases with dN ch /dη . This is in contradiction with the experimental data. However, considering small quark number effects in QCM will predict the correct behavior of the ratio Ξ/φ, see the solid line in Fig. 11. The formation of Ξ − needs not only two s quarks but also a d quark, which is different from the formation of φ needing only a s and as. Therefore, in events of small multiplicity or small quark numbers, the formation of Ξ − will be suppressed to a certain extent (or forbidden occasionally) due to the need of one more light quark, in comparing with the formation of φ. We see that the calculated ratio Ξ/φ using QCM increases with sys-tem multiplicity dN ch /dη and the increased magnitude of the ratio is consistent with the experimental data of pp collisions at √ s = 13 TeV and those of p-Pb collisions at √ s N N = 5.02 TeV [52,59]. Proton and φ have similar masses but different quark contents. Yield ratio φ/p can further test the flavordominated feature of hadron production in QCM. Neglecting small quark number effects, proton yield after taking into account decay contribution of ∆ resonances has a simple expression Using Eq. (24), we get yield ratio which shows a significant dependence on the strangeness suppression factor λ s . We get N φ /N p ≈ 0.22 with λ s = 0.32 in low multiplicity classes and N φ /N p ≈ 0.28 with λ s = 0.36 in high multiplicity classes. The short dashed lines in Fig. 12 are above two values under analytical approximation. They are slightly lower than the experimental data [52,59], symbols in the figure. Small quark number effects will increase the ratio to a certain extent in terms of the suppression of proton yield. We further show the numerical results of our model including small quark number effects, the solid line, and we see a good agreement with the data. TeV, solid circles, are taken from Refs. [52,59]. The solid line is the numerical results of QCM. The short dashed lines are estimation under analytical approximation.

B. pT spectra of identified hadrons
In Fig. 13, we show the fit of data of p T spectra of proton, K * * and φ [52,53] using QCM and the prediction of other identified hadrons in different multiplicity classes in pp collisions at √ s = 13 TeV. Note that classes IV and V are combined for the K * 0 data and the same for our results. Beside directly comparing the prediction of single hadron spectrum with the future data, we emphasize that QCM can be more effectively tested by some spectrum ratios and/or scaling. The first is to test whether the constituent quark number scaling property between the p T spectrum data of Ω − and φ holds in different multiplicity classes. The second is to study the ratio Ω − /φ as the function of p T . Ω − /φ ratio in QCM is solely determined by the strange quark p T spectrum at hadronization and the ratio usually exhibits a nontrivial p T dependence, as shown in Fig. 14(a), which is a typical behavior of 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 and the peak position of Ω − /φ in high multiplicity classes is also enlarged in comparing with that in low multiplicity classes. The third is to study p/φ ratio as the function of p T to clarify the p T dependence of the ratio is flavor originated or mass originated? Results of QCM are shown in Fig.  14(b) which decrease with p T and show relatively weak multiplicity dependence.  9 2 × I( ) 8 2 × II( ) 7 2 × III( ) 6 2 × IV( Solid lines are results of QCM and symbols are preliminary data of ALICE collaboration [52,53].

V. SUMMARY AND DISCUSSION
Taking advantage of available experimental data of hadronic p T spectra and yields at midrapidity, we have systematically studied the production of soft hadrons in pp collisions at √ s = 13 TeV within a framework of quark combination mechanism for hadronization. We applied a quark combination model which assumes the constituent quarks and antiquarks being the effective degrees of freedom for the parton system at hadronization and takes equal velocity combination approximation in hadron formation. We used the model to systematically calculate the p T spectra and yields of soft strange hadrons in inelastic events (INEL>0) and in different multiplicity classes.
We found several interesting results which are sensitive to the hadronization mechanism. (1) Data of p T spectra of Ω − and φ in inelastic events (INEL>0) in pp collisions at √ s = 13 TeV exhibit a constituent quark number scaling property. Data in high multiplicity classes in pp collisions at √ s=7 TeV also show this scaling property. It is the first time that such a scaling property is observed in high energy pp collisions. This is an obvious experimental signal of the quark combination mechanism at hadronization in high energy pp collisions. (2) Data of p/φ ratio in inelastic events (INEL>0) in pp collisions at √ s = 13 TeV show an obvious decrease with the increasing p T , which indicates the statistical hadronization model is not responsible for this observation. We demonstrated that data are naturally explained by the quark combination model. (3) Data of yield ratios Λ/π, Ξ − /π and Ω − /π divided by the values in inelastic events as the function of system multiplicity dN ch /dη 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 during their production via the combination of strange quarks and up/down quarks.
By the quark number scaling property, the p T spectrum of strange quarks was directly extracted from data of Ω and φ. The p T spectrum of up/down quarks was extracted from data of other hadrons containing up/down constituent quarks. These extracted quark momentum distribution functions are important results which describe the property of strongly-interacting partonic system at hadronization in the language of constituent quarks.
To confirm such a new feature of hadronization dynamics in high energy pp collisions, we should carefully study all the related experimental data. We have made plenty of predictions on p T spectra and spectrum ratios of strange hadrons in pp collisions at √ s = 13 TeV to further test our model by the future experimental data. On the other hand, compared to our previous study in pp collisions at √ s = 7 TeV [21] which gives the first indication, the current study provides a stronger suggestion of quark combination hadronization in high energy pp collisions. We still need further systematical studies in pp collisions at other available LHC energies so that we can test the universality of this new hadronization feature and study its relation with the possible creation of mini-QGP in small collision systems.

ACKNOWLEDGMENTS
This work is supported by the National Natural Science Foundation of China (11575100), and by Shandong Province Natural Science Foundation(ZR2019YQ06,ZR2019MA053), and by A Project of Shandong Province Higher Educational Science and Technology Program (J18KA228).