Charmed and bottomed hadronic cross sections from a statistical model

In this work, we extended our statistical model with charmed and bottomed hadrons, and fit the quark creational probabilities for the heavy quarks, using low energy inclusive charmonium and bottomonium data. With the finalized fit for all the relevant types of quarks (up, down, strange, charm, bottom) at the energy range from a few GeV up to a few tens of GeV's, the model is now considered complete. Some examples are also given for proton-proton, pion-proton, and proton-antiproton collisions with charmonium, bottomonium, and open charm hadrons in the final state.


Introduction
Low and medium energy (up to 20-30 GeV) inclusive and exclusive cross sections are important inputs for transport calculations in heavy ion collisions [1,2,3,4,5]. Most of these cross sections are coming from measurements [6,7], and/or from some effective model calculations [8] e.g. Walecka type models [9]. For many interesting processes, measurements are nonexistent, or very sparse e.g. for strange vector mesons K * (892) [10,11,12] or charmonium cross sections near threshold [13,14]. As some of these hardly accessible processes could be important in the investigation of strongly interacting nuclear matter formed in heavy ion collisions [15,16,17], phenomenological model calculations or extrapolation from fits to available data are necessary. In [18] we proposed a model, based on the statistical Bootstrap approach [19,20,21], which were able to give reliable estimates to low/medium energy cross sections (few GeV) for exclusive hadronic processes. The model is extended in [22] to be able to give estimatations to inclusive processes as well, with an approximated error bound coming from the model uncertainties. Using the extended model, we were able to fit the charm quark creational probability to existing charmonium data. Using the fitted value, the model had been proved to be able to describe inclusive charmonium production in proton-proton, and in pion-proton collisions. In this paper, we further extend the model by introducing more particles into the calculations, and by fitting the bottom quark creational probability. These extensions will slightly change the previously fitted value for the charm quark creational probability as well. Including heavy quarks into statistical models have been done in the past [23], where usually a suppression factor is introduced into the quarks momentum distribution function, which measures the deviation from chemical equilibrium. In our model the quark creational probabilities are taking over this role, which were fitted to measured cross sections, with final states containing heavy quarks. With the modifications the model is now considered complete, as the energy range, we intend to use is far below the top quark mass, so its inclusion in the model is not necessary. The paper is distributed as follows. In Sec.2 a basic formulation is given, which covers the main ingredients of the model. In Sec.3 the fitting procedure for the charm and bottom quark creational probabilities is covered, while in Sec.4 some charmonium, bottomonium and open charm inclusive cross sections are calculated for proton-proton, proton-antiproton and pion-proton collisions. After these examples Sec.5 concludes the paper.

Basic formulation
The model is based on the assumption that during a collision process a so called fireball is formed, which after a short time, will hadronize into a specific final state. The collision process therefore can be factorized into an initial dynamical part forming the fireball, and a mixed dynamical and statistical part, which describes the hadronization, as it is shown in Eq.1 where σ n→k is a generalized transitional probability for an n → k process, p i is the 3-momenta of the incoming particles, E is the center of momentum energy of the nbody collision, and q i is the 3-momenta of the outgoing particles. The function R(E, p 1 , ...p n ) describes the initial stage of the collision, while w(E, q 1 , ..., q k ) describes the hadronization of the fireball. Throughout this paper, we only consider two-body collisions, so σ 2→k becomes the cross section of the collision. It is assumed that the dynamical part is described by the inclusive cross section of the colliding particles, so the integral in the first bracket simply becomes the inelastic cross section of the two-body reaction. This assumption is widely used in many statistical models in the past [24,25]. Although it is worth to mention that it is possible to include more than two-body reactions as well, where only the initial dynamical part R(E, p 1 , ..., p n ) has to be changed accordingly to some theory. The integral in the second bracket describes the fireball hadronization, and will be denoted by W (E). The hadronization probability depends on the following factors: (1) Fireball decay scheme probability (2) Fireball density of states (3) Phase space of the final state (4) Quark combinatorial factor (5) Quark creational probabilities. In the followings a short description is given to each of the ingredients. After the fireball is formed, there is a possibility that the fireball with invariant mass M hadronizes into some final state or decays into smaller fireballs with invariant masses M 1 , M 2 . Each of these smaller fireballs are then again choose to hadronize or decay into smaller fireballs. At the end of the chain there will be a specific number of fireballs, with invariant masses M 1 , M 2 , . . . M k , and each of them ultimately decay into a specific hadronic final state. The probability of the number of fireballs is given by kinematical considerations and is described in detail in [18]. The hadronization procedure is based on the statistical bootstrap approach, which sets the probability for the number of hadrons coming from one fireball to P f b 2 = 0.69 and P f b 3 = 0.24 [26]. More particles have negligible probabilities compared to the two-, and three body cases, therefore in our model the possible number of hadrons coming from one fireball is two, or three. However it is straightforward that due to the many number of smaller fireballs the number of hadrons in the full final state could be much more than three.
Following the bootstrap approach the hadronization into final state particles should depend on the density of states (DOS), and also the two-, and three body phase spaces with a slight difference on how we handle non-stable and stable particles. The density of states is given by: where T 0 is the hadronization temperature, E 0 = 500 MeV is a cut off parameter, and a is an irrelevant normalization factor, which disappears in the final results. These parameters are fitted using experimentally measured particle multiplicities and momentum spectra for protons, pions and kaons in [27], and in [19]. It is also mentioned that T 0 could be essentially anything from 130 MeV to 170 MeV (with some corresponding changes in E 0 and a), giving almost the same mass spectrum, therefore, as our model is not really sensitive to E 0 and a, we did our own fit to T 0 , using exclusive cross section ratios in [18], giving T 0 = 160 MeV. In the model, stable and non-stable hadrons are distinguished, which is manifested in the phase spaces considered. A hadron is considered a resonance, if it has a non-negligible width and is able to decay into a lower lying hadronic state via the strong interaction. Taking into consideration the two types of particles, the phase space is written as: where q j is the three momenta of the j'th particle, E is the energy, k = 2, 3 the number of hadrons in the final state of a fireball with invariant mass M , V is the interaction volume, set to a corresponding interaction radius of r = 0.5 fm, and F BR r is a Breit-Wigner factor given by: where M r is the invariant mass of the resonance, with a corresponding decay width of Γ r , and a pole position m r . The product r∈R goes through all of the resonant particles, while for stable hadrons only the momentum integral remains, with F BR r = 1. The last ingredient to the hadronization probability is the quark combinatorial factor and the corresponding quark creational probabilities. This is, in some way, similiar to the parton model [28,29], where the partonic cross sections are integrated out with the parton density functions giving a hadronic cross section. In this model the number of quark-antiquark pairs are estimated from simple phase space considerations and is given by [30]: where T 0 is the hadronization temperature. For each quark there is a creational probability P i , i = u, d, s, c, b and the quark distribution function is then given by the multinomial distribution: where n i is the number of quarks of type i = u, d, s, c, b.
The distribution function gives the probability of a specific quark number configuration. In our model, we take the configuration, where the probability is at maximum, which corresponds to a quark number n i = P i N . With the number of quarks, the number of possibilities of a final state is then calculated with simple combinatorics, and multiplied by a color factor given by the number of colorless combinations, then normalized by the possible two-, and three body final state factors, giving a final quark combinatorial probability. Finally this probability is multiplied by the maximum value of the distribution function, giving the final quark combinatorial factor C Qi . When the expected number of quarks are very small, which is the case for the heavy quarks with small quark creational probabilities, a suppression factor is introduced. This is described in Sec.3, where we fit the charm and bottom quark creational probabilities.
If the calculated process respects the conservation laws and the quarks/antiquarks from the colliding particles are also included into the combinatorial factors, meaning n i = n i +n initial will be the number of a specific quark/antiquark, there will be no quarks or antiquarks left unpaired at the end, and every quark, which is not used in the final state, could annihilate with a corresponding antiquark.
Putting together the previously described factors the normalized hadronization probability for k fireball is given by Eq.7.
where P f b k (E) is the fireball decay scheme probability for k-fireballs, N i1,...i k is a symmetry factor counting the fireballs with the same hadrons after hadronization, x min and x max limits are given by kinematical considerations, and Z k (E) is the energy-dependent k-fireball normalization factor given by: where the subscript < l 1 ..l k >∈ S means that the summation goes for all of the possible processes having quantum numbers S, where S includes baryon number, electric charge, strangeness, charmeness, and bottomness. Also, for simplicity the T i (E) = C Qi (E)P H,i ni (E) notation is introduced, where C Qi (E) is the quark combinatorial factor, and P H,i ni (E) is given by: where n = 2, 3 the number of hadrons coming from one fireball, s l is the spin of the hadron, P d n is the probability of two-, or three hadrons, and N I is a symmetry factor given by the number of indistinguishable hadrons in the final state of one fireball. In the normalization, only processes with the quantum numbers of the initial state is considered and for many fireballs only the full final state, with all the produced hadrons has to respect the conservation laws, and not each fireball separately.

Quark creational probabilities for bottom and charm quarks
At high energies it is expected that the quark creational probabilities will be asymptotically the same, so P u = P d = P s = P c = P b , with P u + P d + P s + P c + P b = 1, however in the energy range we are interested in (2-20 GeV) the equality is only valid for the up and down quarks. For the strange, charm and bottom quarks P s , P c , P b has to be fitted from experiments or has to be estimated from some underlying theory. We choose to fit the parameters from experiments, where for the strange quark, the exclusive pp → K + K − process was used, assuming a constant P s [18]. For the up, down and strange quarks the values are P u = 0.425, P d = 0.425, and P s = 0.15. For the charm and bottom quarks the fit is not that straightforward due to their considerable masses and therefore their huge suppression. Due to the very small charmed and bottomed probabilities it is assumed that only one charm/bottom quark-antiquark pair is created, but the following arguments can be made to more charm and bottom quarks as well. Due to this assumption the constrained probability (one charm/bottom quark is certainly created) obtained from the quark number distribution Eq.6, will be much smaller than the non-constrained maximal probability. The suppression for both heavy quarks is coming from the ratio of the maximal and the constrained probability, obtained from the quark number distribution mass function, and can be expressed as: where it is assumed that only one heavy quark is created, replacing one strange quark, meaning n ′ c,b = 1, and n ′ s = n s −1, while n ′ u = n u and n ′ d = n d . It is evident that the suppression is energy dependent even when the quark creational probability is constant and the ratio of the global and the constrained maximum tends to unity as the energy goes higher. After this point the global maximum value of the quark number distribution mass function is used and the expected number of heavy quarks, which corresponds to the global maximum are n c,b = P c,b N , so the transition between the constrained and the global maximum case is continous.
A constant quark creational probability is however not too realistic, because at high energies these parameters are essentially have to be equal, therefore we assume a simple linear relationship between the energy and the quark creational probability P c,b = a c,b E, for the charm and bottom quarks. For the strange quarks, we keep the constant P s value, as it was sufficient to describe exclusive K, Λ, and inclusive strange vector meson production cross sections in the energy range of a few tens of GeV's, where we intend to use the model. It is possible that at higher energies these values will change, but at this relatively low energy range the linear relationship for the heavy quarks and the constant values for the up, down and strange quarks seem to be good approximations. It is also possible that the rise of the strange quark creational probability is so small that P s can be approximated by a constant value at the energy range considered. Using Eq.10 and the energy dependent P c,b , the heavy quark suppression can be expressed as: where a c,b slope parameters have to be determined from experiments. The suppression is energy dependent through the number of quarks, and P c,b . This is of course only valid until γ c,b reaches unity, as after that point the global maximum will naturally give one heavy quark, and it is not necessary to constraint the quark number distribution anymore. The actual energy when this happens, depends on the slope parameters a c,b of the charm and bottom quarks and can be calculated setting γ c,b = 1 and solving the following equation: where E L is the limit energy. On Fig.1 the limit energy is shown as the function of the slope parameter a c,b . As the quark creational probability depends only on a c,b , this is the only parameter that has to be fitted from experiments. In [22] a very simple fitting procedure is used where, a c was fitted only to one data point from the π − p → J/Ψ X reaction and checked if it gives the correct values to different measurement points, in different collisions as well. In this paper a least squares fitting method is used, where an optimum point is sought in the error function, defined as: where meas i is a measurement point, model i is the corresponding model calculation, and N is the number of measurement points used for fitting. For the charm quarks the inclusive π − N → J/Ψ X, and pp → J/Ψ X reactions were used, where the measured values were collected from [31,32,33,34,35,36,37,38,39,40,41,42]. To calculate the inclusive cross section from the model, apart from the prompt J/Ψ production, the radiative and hadronic decays from the excited charmonium states Ψ (3686), χ c1 , and χ c2 are included as well, using a two fireball decay scheme. In Fig.2 the calculated error for the J/Ψ production cross section is shown, where a clear minimum can be seen at a c = 8.5 · 10 −4 GeV −1 , which is really close to our previous fit in [22]. The cause of the difference is that now more particles are introduced in the model e.g. new charmed mesons/baryons, bottomed mesons/baryons, which in exchange slightly changes the normalization values as well.
For the bottom quarks the fitting procedure is the same as it was for the charm quarks. The process π − p → Υ (1S)X is used for fitting, where the measured points are taken from [43,44]. In the model calculations apart from the direct Υ (1S) meson, the excited onium states Υ (2S), Υ (3S), χ b1 (1P ), χ b2 (1P ), χ b1 (2P ), and χ b2 (2P ) are used as well. The results for the error function can be seen in Fig.3, where a minimum is obtained, now at a b = 1.05 · 10 −5 GeV −1 .
The energy dependent quark creational probabilities and the suppressions can be seen in Fig.4, where the limiting energy -the energy, where the suppression from the constrained quark number distribution breaks down -is also shown. After the limiting energy, P c,b still increases until a full equilibrium is reached with the other type of quarks, however, we do not intend to use the model until that energy. The probabilities are energy dependent, therefore the previously fitted P u , P d , and P s values have to be changed accordingly, to  bottom probabilities, the corrections are very small at the energy range we are interested in, so in practical calculations P ′ i ≈ P i can be used. In the next section the results for different onium state production, as well as some open charm calculations are shown for proton-proton, pion-proton, and proton-antiproton collisions.

Results
After the quark creational probabilities for heavy quarks are obtained, it is possible to calculate inclusive cross sections for charmed and bottomed hadrons as well. The first reaction shown is the energy dependence of the pp → J/Ψ X process, which is used to fit the charm quark creational probability. The decays from the excited Ψ (3686),  [31,32,33,34,35,36,37,38,39,40,41,42].
χ c1 , χ c2 states are also included in the calculations. The results can be seen in Fig.5, where as expected, the agreement with the data is really good. The second example is the π − p → J/Ψ X process, which is expected to give a few times larger cross section compared to the protonproton case. The results can be seen in Fig.6, where again the agreement is still satisfactory. The third example is the π − p → Υ (1S)X process, where apart from the direct bottomonium production, the decays from the excited onium states Υ (2S), Υ (3S), χ b1 (1P ), χ b2 (1P ), χ b1 (2P ), and χ b2 (2P ) are included. The results can be seen in Fig.7.
The agreement between measurement and model in the previous three processes are, however expected as these were the ones used to fit the quark creational probabilities. The fourth process is the pp → Υ X, where Υ = Υ (1S)+ Υ (2S) + Υ (3S). The measured data points are at √ s = 29.1 GeV, with proton beam and five different target nuclei (Be, Al, Cu, Ag, W), taken from [45]. The results can be seen in Fig.9, where on the upper side all the five measured data points are shown with the model calculation at √ s = 29.1 GeV, while at the bottom the energy dependence of the model calculation is shown.
For the fifth and sixth processes, we checked the protonproton collision suppression to proton-antiproton collisions in the inclusive charmonium and bottomonium channels. In these calculations, we excluded the one body final states in the proton-antiproton collisions e.g. pp → J/Ψ , which only gives a significant contribution near threshold. The results can be seen in Fig.9. The proton-antiproton collision is therefore favorable in constrast to the protonproton case, however its production rate is still lower than in pion-proton collisions. This is not just a curiosity on its own, because in heavy ion collisions an interesting problem is the propagation of onium states in dense nuclear matter, as they could have a significant mass shift, which  Fig. 7. Model calculation for the π − p → Υ X inclusive cross section. The black and red lines are the model result and its error respectively, while the circles are measurement points from [43,44].
is proportional to the gluon condensate [47,48], so if one could measure the mass shift of e.g. Ψ (3686) the nuclear expectation value of the gluon condensate at a specific nuclear density can be deduced. Due to the really small cross sections it is necessary to create the most onium states one can get, which is one reason to use antiprotons or pions in these investigations.
Some of our transport calculations [49,50] requires the determination of the direct charmonium production, without the decays from higher states as well. As there are only just a few measured values for e.g. direct Ψ (3686) production, we checked the ratios defined in Eq.14 and in Eq.15  for the χ c1 , χ c2 and Ψ (3686) particles: where A and B represents the colliding particles, in this case π − and p, and Br is the branching fraction of the particles to decay into a J/Ψ meson. The measured ratios  Table 2. Total J/Ψ production ratio to the direct χc1, χc2 production cross sections. The measured values are collected from [51,52,53,54,55,56].
are collected from [51,52,53,54,55,56] and the comparison with the model calculations at different collision energies are shown in Table.1 and in Table.2, where a really good match with the data is achieved. This allows us to calculate the full energy dependence of the direct production of these states and put it into our transport simulations.
The following few processes are somewhat more complicated than the previous ones. These are the inclusive open charm productions in π − p and in pp collisions. The main problem is that our model is heavily based on the known resonances, excited states, which could decay into hadrons we are interested in, and if the branching fractions or the quantum numbers of such particles are not well known, then it could not be included in the calculations. , Λ c (2860) + , Λ c (2880) + , Λ c (2940) + are included, whose quantum numbers are at least fairly well known in the PDG. However, the branching fractions are not known at all, and in the excited Λ c baryon resonances even less is known than the D mesons. To be able to give an estimation, we introduced an asymmetry parameter r defined as: with the assumption that every other decay is negligible, so Br(X → D + /D − ) + Br(X → D 0 /D 0 ) ≈ 1 .
The ratio thus gives the fraction that a resonance 'X' e.g. D * (2300) 0 giving a D + to the expense of D 0 , and varied it between [0. 11,9], which for the branching ratios means Br i ∈ [0.1, 0.9]. This applies to the differently charged and neutral resonances with different possible final states as well. This is reasonable as for example in the D * (2007) + decays, where the branching fractions are known, the asymmetry in the decay products is r ≈ 0.46. For the Λ resonances, where the final state products are also not well known, we assumed it could decay as Λ + c → pD 0 , Λ + c → nD + , and for the antiparticles Λ c → pD 0 , Λ c → nD − according to the quark flow diagrams in Fig.10, with the same asymmetry variations as in the meson resonances. In [57] an effective model calculation estimated the asymmetry to r ≈ 1, and it is highly possible that each resonance has different branching fractions, so the best we can do at this point is to stay at the varied asymetries and give an order of magnitude estimation. It could be an interesting further step to include all the possible resonances and branching ratios from e.g. effective model calculations. Using the mentioned assumptions, the results are shown in Fig.11-Fig.14, where only cross section intervals are shown due to the varied asymmetries. It can be seen that an order of magnitude estimation is very much achievable, even with the high uncertainty in the parameters.

Conclusions
We extended our model with the inclusion of charmed and bottomed hadrons and estimated the charm and bottom quark creational probabilities, assuming a simple linear relationship between the probability and the energy P c,b = a c,b E. A fit is made for the two slope parameters Fig. 11. Estimated cross section interval for the π − p → D + /D − process, with varied branching fraction asymetries. The data is taken from [58,59,60]. a c , and a b , using inclusive charmonium and bottomonium data in proton-proton and in pion-proton collisions. The suppression of the higher charmonium states Ψ (3686), χ c1 , and χ c2 to the direct J/Ψ production is also calculated and compared to the measured ratios in π − p and pp collisions, giving a really good match with the data. The model is further validated through the processes pp → J/Ψ X, π − p → J/Ψ X, pp → Υ X, and π − p → Υ X. Estimations were made to proton-antiproton collisions and the ratios σ pp→J/Ψ X /σ pp→J/Ψ X , σ pp→Υ X /σ pp→Υ X are calculated. It can be concluded that at low energies proton-antiproton and pion-proton collisions are very much favorable if one wants to produce charmonium particles. For further validation open charm production cross sections were also estimated, namely the π − p → D + /D − X, pp → D + /D − X, π − p → D 0 /D 0 X, and the pp → D 0 /D 0 X processes.  The main problem with these calculations are the missing branching fractions in the PDG, which are needed to make reliable estimates, however introducing an asymmetry parameter, an order of magnitude estimation was made, which is in agreement with the measurements. With these extensions the model is now considered complete regarding non-exotic states. It is however a further question how to include e.g. tetra-quarks into the model, as there are no measurements for these states at low energies, where the model works.