Statistical method in quark combination model

TL;DR

This paper introduces a simplified generating function method to solve probability distributions of hadrons in the quark combination model, aiding the study of multiplicity fluctuations and phase transitions in heavy-ion collisions.

Contribution

It presents a new, more straightforward method using generating functions for the quark combination rule, enhancing analysis of hadronization and multiplicity distributions.

Findings

Derived probability distributions for baryons, antibaryons, and mesons.

Analyzed multiplicity cumulant ratios of net-protons.

Discussed potential for studying phase transitions in heavy-ion collisions.

Abstract

We present a new method of solving the probability distribution for baryons, antibaryons and mesons at the hadronization of constituent quark and antiquark system. The hadronization is governed by the quark combination rule in the quark combination model developed by the Shandong Group. We use the method of the generating function to derive the outcome of the quark combination rule, which is much simpler and easier to be generalized than the original method. Furthermore, we use the formula of the quark combination rule and its generalization to study the property of multiplicity distribution of net-protons. Taking a naive case of quark number fluctuations and correlations at hadronization, we calculate ratios of multiplicity cumulants of final-state net-protons and discuss the potential applicability of quark combination model in studying hadronic multiplicity fluctuations and the underlying phase transition property in relativistic heavy-ion collisions.

Equations

$$q_1{\overline{q}}_2{q}_3{q}_4{q}_5{\overline{q}}_6{q}_7{q}_8{\overline{q}}_9{\overline{q}}_{10}{q}_{11}{\overline{q}}_{12}{\overline{q}}_{13}{\overline{q}}_{14}{\overline{q}}_{15}{q}_{16}{q}_{17}{q}_{18}{\overline{q}}_{19}{\overline{q}}_{20}$$

$$N_M+3N_B+N_r$$

$$N_M+3N_{\bar{B}}+N_{\bar{r}}$$

$$(a)$$

$$(b)$$

$$(c)$$

$$(d)$$

$$S(N_q,N_{\bar{q}})=$$

$$\times {\delta }_{N_M+3N_B+N_r,N_q}{\delta }_{N_M+3N_{\bar{B}}+N_{\bar{r}},N_{\bar{q}}}$$

$$\displaystyle=\left(\begin{array}[]{c}N_{q}+N_{\bar{q}}\\ N_{q}\end{array}\right)\equiv\frac{(N_{q}+N_{\bar{q}})!}{N_{q}!N_{\bar{q}}!},$$

$$F(N_M,N_B,N_{\bar{B}},N_r,N_{\bar{r}})=\sum _{i=0}^{N_M}F(i,N_B,N_{\bar{B}},0,0).$$

$$F(N_M,N_B,N_{\bar{B}},1,0)$$

$$=F(N_M,N_B,N_{\bar{B}},0,0)+F(N_M-1,N_B,N_{\bar{B}},2,0)$$

$$=F(N_M,N_B,N_{\bar{B}},0,0)+F(N_M-1,N_B,N_{\bar{B}},1,0).$$

$$F(N_M,N_B,N_{\bar{B}},0,0)$$

$$=F(N_M-1,N_B,N_{\bar{B}},1,0)+F(N_M-1,N_B,N_{\bar{B}},0,1)$$

$$+F(N_M,N_B-1,N_{\bar{B}},2,0)+F(N_M,N_B,N_{\bar{B}}-1,0,2).$$

$$F(N_M,N_B,N_{\bar{B}},0,0)$$

$$=3F(N_M-1,N_B,N_{\bar{B}},0,0)$$

$$+F(N_M,N_B-1,N_{\bar{B}},0,0)+F(N_M,N_B,N_{\bar{B}}-1,0,0).$$

$$F(N_M,0,0,0,0)=3F(N_M-1,0,0,0,0),$$

$$F(N_M,0,0,0,0)=2\times 3^{N_M-1},$$

$$A\left(x,y,z\right)$$

$$=\sum _{N_M=0}^{\infty }\sum _{N_B=0}^{\infty }\sum _{N_{\bar{B}}=0}^{\infty }F(N_M,N_B,N_{\bar{B}},0,0)x^{N_M}{y}^{N_B}{z}^{N_{\bar{B}}},$$

$$A(x;N_B,N_{\bar{B}})=\sum _{N_M=0}^{\infty }F(N_M,N_B,N_{\bar{B}},0,0)x^{N_M}.$$

$$\left(1-3x\right)A(x;N_B,N_{\bar{B}})-F(0,N_B,N_{\bar{B}},0,0)$$

$$=A(x;N_B-1,N_{\bar{B}})+A(x;N_B,N_{\bar{B}}-1)$$

$$-F(0,N_B-1,N_{\bar{B}},0,0)-F(0,N_B,N_{\bar{B}}-1,0,0).$$

$$A(x;0,N_{\bar{B}})$$

$$=\frac{2x}{(1-3x{)}^{N_{\bar{B}}+1}}+\frac{1}{(1-3x{)}^{N_{\bar{B}}}},$$

$$A(x;N_B,0)=\frac{2x}{(1-3x{)}^{N_B+1}}+\frac{1}{(1-3x{)}^{N_B}}.$$

$$\sum _{N_B=0}^{\infty }A(x;N_B,0)y^{N_B}$$

$$\sum _{N_{\bar{B}}=0}^{\infty }A(x;0,N_{\bar{B}})z^{N_{\bar{B}}}$$

$$A(x,y,z)=\frac{1-x}{1-3x-y-z}.$$

$$C(z^{N_{\bar{B}}})=\frac{1-x}{(1-3x-y{)}^{N_{\bar{B}}+1}}.$$

$$C(z^{N_{\bar{B}}}y^{N_{B}})=\left(\begin{array}[]{c}N_{B}+N_{\bar{B}}\\ N_{B}\end{array}\right)\frac{1-x}{(1-3x)^{N_{B}+N_{\bar{B}}+1}},$$

$$(1\pm w)^{-n}=\sum_{k=0}^{\infty}\left(\begin{array}[]{c}n-1+k\\ k\end{array}\right)(\mp 1)^{k}w^{k},$$

$$\displaystyle F(N_{M},N_{B},N_{\bar{B}},0,0)=3^{N_{M}}\left(\begin{array}[]{c}N_{B}+N_{\bar{B}}\\ N_{B}\end{array}\right)\times$$

$$\displaystyle\left[\left(\begin{array}[]{c}N_{M}+N_{B}+N_{\bar{B}}\\ N_{M}\end{array}\right)-\frac{1}{3}\left(\begin{array}[]{c}N_{M}-1+N_{B}+N_{\bar{B}}\\ N_{M}-1\end{array}\right)\right],$$

Full text

5.75.-q, 12.38.Mh, 24.60.Ky

1 Introduction

Most hadronization models such as the Lund string model [1, 2, 3, 4] or the coalescence/recombination models [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] in electron-positron and hadron-hadron collisions assume that a hadron is formed by quarks and antiquarks in the neighborhood of phase space. Normally the phase space can be decomposed into the longitudinal and transverse direction. The longitudinal phase space plays a special role. For example in the Lund model a string is formed between a quark and an antiquark moving back to back, it is an object with only one spatial dimension, i.e. with only the longitudinal phase space. A quark-antiquark pair is produced as the result of the string break into two pieces. The process of string breaks continues until it ends at a lowest energy scale. The transverse momentum space of the excited quark-antiquark pair is rather limited and can be described by the exponentially suppressed function of the transverse momentum. In relativistic heavy-ion collisions, the longitudinal phase space is still dominant although the transverse expansion is significant.

The quark combination model developed by the Shandong group (SDQCM) [22, 9, 23, 10, 11, 24, 25, 28, 26, 27, 29, 30, 31, 32, 33, 34, 35] is a kind of the exclusive or statistical hadronization model which is different from the Lund string model or the coalescence model. The model firstly takes the constituent quark degrees of freedom as an effective description for the strongly-interacted quark gluon system at hadronization. Then, the model adopts a quark combination rule (QCR) to combine the quarks and antiquarks in the neighborhood of the longitudinal phase space into baryons and mesons. Since the longitudinal phase space is easily described by the momentum rapidity, the correlation in rapidity is the basis of the QCR. Such kind of QCR in SDQCM has successfully explained many data of hadronic production in $e^{+}e^{-}$ and $pp$ collisions [36, 37, 24, 25] and also rapidity distributions of hadrons in heavy-ion collisions [26, 38, 27, 31].

The probability distributions of the particle numbers of baryons, antibaryons and mesons is an important observable in high energy collisions and are closely related to the hadronization dynamics. Especially the particle number distributions are essential to search for the critical end point of the QCD phase diagram in relativistic heavy-ion collisions [39, 40, 41, 42, 43]. In this paper, we propose a new method to solve the baryon and meson number distribution in the context of the QCR in the SDQCM. The new method is much simpler and easier to be generalized to more sophisticated cases than the original one [44]. Then we use these formula to calculate the multiplicity distribution and the ratios of cumulants for net protons in relativistic heavy-ion collisions.

The paper is organized as follows. In Sec. 2, we briefly present the original QCR in the SDQCM. Then we use the generating function method to solve the particle number probability of baryons, antibaryons and mesons for given numbers of quarks and antiquarks. In Sec. 3 we generalize the original QCR and derive the corresponding particle number probability of baryons, antibaryons and mesons using the generating function method. In Sec. 4 we compare the numerical difference between the original QCR and the generalized QCR in terms of moments of the antibaryon number. In Sec. 5, we formulate the fluctuation and correlation of identified hadrons. In Sec. 6 and 7, we study the effects of the quark number fluctuation and resonance decays. In Sec. 8, we give an illustrative example of applying the QCR and generalized QCR to calculate the ratios of cumulants for net protons in heavy-ion collisions and compare with data. Finally we give a summary and make discussions in Sec. 9.

2 Baryon and meson formation in original quark combination rule

Due to the non-perturbative QCD feature, the transition from quarks and/or gluons to hadrons is not yet solved from first principle and is only described by the phenomenological models at present. Inspired by the simple and effective description of constituent quark model in explaining the static property of hadrons, SDQCM assumes the constituent quarks and antiquarks as the effective degrees of freedom for the strongly-interacted quarks and gluons at hadronization and therefore builds a simple hadronization phenomenology by the combination of these constituent quarks and antiquarks into hadrons. We emphasize that there exist explicit experimental signals for such constituent quark degrees of freedom in high energy collisions, such as the quark number scaling property of elliptic flows and transverse momentum spectra for hadrons observed in recent years [45, 46, 47, 48, 49, 50]. A phenomenological quark combination rule was proposed in Ref. [22] to describe the combination of these constituent quarks and antiquarks in one-dimensional phase space and can successfully describe the data of yields and momentum distributions of hadrons in high energy reactions [36, 37, 26, 38, 27, 31, 24, 25]. In this section, we use the generating function method to solve the probability distribution for the number of baryons, antibaryons and mesons formed by QCR.

2.1 Original quark combination rule

In the QCM $N_{q}$ quarks and $N_{\bar{q}}$ antiquarks produced in an event are put into a queue and then allow them to combine into hadrons one by one following a rule called the quark combination rule (QCR)[22]. QCR is based on the basic property of QCD.

A $q\overline{q}$ pair pair may be in a color octet with a repulsive interaction or a singlet with a attractive interaction. If $q\overline{q}$ is adjacent in phase space, $q\overline{q}$ will have enough time or opportunity to be in a color state and hadronizes into a meson. For a $qq$ pair, it may be in a sextet or an antitriplet. If its nearest neighbor is a $q$ in phase space, they can hadronizes into a baryon. If the neighbor of $qq$ is a $\overline{q}$, $q\overline{q}$ will win the competition to form a meson and leave a $q$ alone to combine with other quarks and/or antiquarks. This is because the attraction strength of the singlet for $q\bar{q}$ is two times that of the antitriplet for $qq$ (via counting color factor in one-gluon exchange case).

The original QCR proposed in Ref. [22] reads:

Check if there are partons in the queue. If there are no partons, the process ends. Otherwise, start from the first parton ($q$ or $\overline{q}$) in the queue and go to the next step. 2. 2.

Look at the second parton. If there is no second parton, the process ends. Otherwise, if the baryon number of the second parton in the queue is different from the first one, i.e. the first two partons are either $q\overline{q}$ or $\overline{q}q$, they combine into a meson and are removed from the queue, go to step 1; Otherwise they are either $qq$ or $\overline{q}\overline{q}$, go to the next step. 3. 3.

Look at the third parton. If there is no third parton, the process ends. Otherwise, if the third parton is different in baryon number from the first one, the first and third parton form a meson and are removed from the queue, go to step 1; Otherwise the first three partons combine into a baryon or an antibaryon and are removed from the queue, go to step 1.

The following example show how the above QCR works

[TABLE]

In relativistic heavy-ion collisions, the longitudinal rapidity space is predominant and the rapidity density of quarks and antiquarks is quite large. Therefore, it is suitable and straightforward to apply such a QCR in one-dimensional longitudinal rapidity space. However, it is quite complicated in three dimensional phase space because one can not easily define a particular hadronization order to apply QCR [51].

2.2 Recursive relation for $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$

We consider the system consisting of $N_{q}$ quarks and $N_{\bar{q}}$ antiquarks stochastically populated in one-dimensional phase space. After combination by QCR, there are $N_{B}$ baryons, $N_{\bar{B}}$ antibaryons, and $N_{M}$ mesons formed and $N_{r}$ quarks and $N_{\bar{r}}$ antiquarks left. There are only five different configurations with $(N_{r},N_{\bar{r}})=(0,0),(0,1),$ $(1,0),(2,0),(0,2)$. The quark number conservation gives

[TABLE]

The outcome of implementing the QCR to the queue of $N_{q}$ quarks and $N_{\bar{q}}$ antiquarks gives a group of numbers $(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$. The queue $(N_{M},N_{B},N_{\bar{B}},0,0)$ can be reached by one of the following ways of adding one $q$ or $\bar{q}$ to the end of other four queues with smaller number of quarks and antiquarks

[TABLE]

We use $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ to denote the number of different queues for a given group $(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$. We define $F(0,0,0,0,0)=1$ and $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})=0$ for the case that any of $N_{M}$, $N_{B}$, $N_{\bar{B}}$, $N_{r}$ and $N_{\bar{r}}$ are negative. Under constraint Eq. (2), the sum of $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ over all different groups of $(N_{M},$ $N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ should be

[TABLE]

where $\delta_{i,j}=1$ if $i=j$ and $\delta_{i,j}=0$ if $i\neq j$.

For non-zero $N_{r}$ and $N_{\bar{r}}$, $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ has a property

[TABLE]

Proof of Property. We firstly take $(N_{r},N_{\bar{r}})=(1,0)$ for example. We note that the queue giving $(N_{M},N_{B},N_{\bar{B}},1,0)$ can be reached by two ways: (a) add $q$ to the end of the queue with $(N_{M},N_{B},N_{\bar{B}},0,0)$; (b) add $\bar{q}$ to the end of the queue with $(N_{M}-1,N_{B},N_{\bar{B}},2,0)$ which can further be obtained from the queue with $(N_{M}-1,N_{B},N_{\bar{B}},1,0)$ by adding a $q$ to the end. Thus we obtain the recursive relation

[TABLE]

Solving the above equation recursively we get Eq. (7), where we have used $F(0,N_{B},N_{\bar{B}},1,0)=F(0,N_{B},N_{\bar{B}},0,0)$. The proof of the case $(N_{r},N_{\bar{r}})=(2,0)$ is straightforward due to the fact that the queue giving $(N_{M},N_{B},N_{\bar{B}},2,0)$ can be only obtained from the queue with $(N_{M},N_{B},N_{\bar{B}},1,0)$ by adding a $q$ to the end. The proof for the cases $(N_{r},N_{\bar{r}})=(0,1)$ and $(N_{r},N_{\bar{r}})=(0,2)$ is similar.

Using properties Eq. (7) and Eq. (3), we obtain

[TABLE]

We make replacement $N_{M}\rightarrow N_{M}-1$ in the above equation and take the difference between the two, and finally we use Eq. (7) to get the recursive equation for $F(N_{M},N_{B},N_{\bar{B}},0,0)$

[TABLE]

where $N_{M}\geqslant 0$, $N_{B}\geqslant 0$, $N_{\bar{B}}\geqslant 0$ excluding two cases $(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})=(0,0,0,0,0),(1,0,0,0,0)$ which we have $F(0,0,0,0,0)=1$ by definition and $F(1,0,0,0,0)=2$ by simple counting.

2.3 Solution to recursive equation by generating function method

First we consider a special simple case: $N_{M}>1$, $N_{B}=N_{\bar{B}}=0$, $N_{r}=N_{\bar{r}}=0$, Eq. (10) is simplified to

[TABLE]

which immediately leads to the solution

[TABLE]

for $N_{M}\geqslant 1$ with $F(1,0,0,0,0)=2$.

Now we consider the general case for $F(N_{M},N_{B},N_{\bar{B}},0,0)$ with $N_{M}\geqslant 0$, $N_{B}\geqslant 0$, $N_{\bar{B}}\geqslant 0$. We define the generating function

[TABLE]

and $F(N_{M},N_{B},N_{\bar{B}},0,0)$ is the coefficient of $x^{N_{M}}y^{N_{B}}z^{N_{\bar{B}}}$ in the polynomial expansion of $A(x,y,z)$ once solved.

To solve $A(x,y,z)$, we define a partial generating function

[TABLE]

Inserting Eq. (10) for $N_{M}\geqslant 1$, we get

[TABLE]

In a special case with $N_{B}=0$ we obtain

[TABLE]

where we have used Eq. (12) and assume $3|x|<1$. In the same way, we obtain

[TABLE]

We also obtain two sum rules

[TABLE]

where the convergent region is $|y|,|z|<|1-3x|$.

We now multiply Eq. (15) by $y^{N_{B}}z^{N_{\bar{B}}}$ and take a sum over $N_{B}$ from 1 to infinity and $N_{\bar{B}}$ from 1 to infinity. By noticing $A(x,y,z)=\sum_{N_{B}=0}^{\infty}\sum_{N_{\bar{B}}=0}^{\infty}A(x;N_{B},N_{\bar{B}})y^{N_{B}}z^{N_{\bar{B}}}$ and using Eqs. (18) and (19), we solve the generating function as

[TABLE]

To obtain $F(N_{M},N_{B},N_{\bar{B}},0,0)$, we firstly extract the coefficient of $z^{N_{\bar{B}}}$

[TABLE]

Then we extract the coefficient of $y^{N_{B}}$ in $C(z^{N_{\bar{B}}})$ as

[TABLE]

where we have used

[TABLE]

for $n>0$ and $|w|<1$. Finally we extract the coefficient of $x^{N_{M}}$ in $C(z^{N_{\bar{B}}}y^{N_{B}})$, i.e.,

[TABLE]

which give the final result for

[TABLE]

Eq. (32) is the result of QCR in Sect. 2.1 and was first given in Ref. [44] by** **mathematical induction and traditional combination method. The current derivation is a simplified one based on the recursive equation and the method of generating functions. The probability of forming $N_{M}$ mesons, $N_{B}$ baryons and $N_{\bar{B}}$ antibaryons in a quark system with $N_{q}$ quarks and $N_{\bar{q}}$ antiquarks is given by

[TABLE]

Here we have removed the label ’0,0’ from $P(N_{M},N_{B},N_{\bar{B}})$ since this is the real case in hadronization.

3 Baryon and meson formation in generalized quark combination rule

3.1 Generalized quark combination rule

The ratio of baryons to mesons (B/M) given by QCR in Sect. 2.1 is larger than observation of heavy-ion collisions [52, 53]. In order to suppress the B/M ratio, we can generalize the QCR in Sect. 2.1 to decrease the formation probability of baryons relative to mesons. The generalized QCR (gQCR) reads:

Check if there are partons in the queue. If there are no partons, the process ends. Otherwise, start from the first parton and go to the next step. 2. 2.

Look at the second parton. If there is no second parton, the process ends. Otherwise, if the baryon number of the second parton is different from the first one, i.e. the first two partons are either $q\overline{q}$ or $\overline{q}q$, they combine into a meson and are removed from the queue, go to step 1; Otherwise they are either $qq$ or $\overline{q}\overline{q}$, go to the next step. 3. 3.

Look at the third parton. If there is no third parton, the process ends. Otherwise, if the baryon number of the third parton is different from the first one, the first and third parton form a meson and are removed from the queue, go to step 1; Otherwise the first three partons are either $qqq$ or $\overline{q}\overline{q}\overline{q}$, go to next step. 4. 4.

Look at the fourth parton. If there is no fourth parton, the first three partons form a baryon or an antibaryon and the process ends. Otherwise, if the baryon number of the fourth parton is different from the first one, the first and fourth parton form a meson and are removed from the queue, go to step 1; Otherwise the first three partons combine into a baryon or an antibaryon and are removed from the queue, go to step 1.

The outcome of implementing the gQCR to the queue of stochastically populated $N_{q}$ quarks and $N_{\bar{q}}$ antiquarks gives a group of numbers $(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$. Same as in the QCR case, the special queue with $(N_{M},N_{B},N_{\bar{B}},0,0)$ can be reached by one of four ways in (3). Other queues with $(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ where $N_{r}\neq 0$ or $N_{\bar{r}}\neq 0$ can be reached by adding $q$ or $\bar{q}$ to the end of queues with smaller $N_{M}$, $N_{B}$ and $N_{\bar{B}}$. Queues marked as $(N_{M},N_{B},N_{\bar{B}},1,0)$ can be reached by

[TABLE]

In manner $(a)$, the special queue $(N_{M},N_{B},N_{\bar{B}}-1,0,3)$ (which is included in $(N_{M},N_{B},N_{\bar{B}},0,0)$) is excluded because of the step 4 in gQCR. Queues marked as $(N_{M},N_{B},N_{\bar{B}},0,1)$ can be reached by

[TABLE]

Queues marked as $(N_{M},N_{B},N_{\bar{B}},2,0)$ can be reached by

[TABLE]

Queues marked as $(N_{M},N_{B},N_{\bar{B}},0,2)$ can be reached by

[TABLE]

The special queues $(N_{M},N_{B},N_{\bar{B}},3,0)$ and $(N_{M},N_{B},$ $N_{\bar{B}},0,3)$ can be build by normal ones

[TABLE]

Properties (34-39) lead to following recursive equations,

[TABLE]

For non-zero $N_{r}$ and $N_{\bar{r}}$, $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ can be obtained with the help of two properties in Appendix A. Using them, we can derive the recursive equation for $F(N_{M},N_{B},$ $N_{\bar{B}},0,0)$

[TABLE]

The detailed derivation of Eq. (41) is given in Appendix B.

3.2 Solution to recursive equation by generating functions for gQCR

We start with the most simple case $F(N_{M},0,0,0,0)$ with $N_{M}\geqslant 1$, Eq. (41) is simplified as

[TABLE]

With the initial values $F(0,0,0,0,0)=1$, $F(1,0,0,0,0)=2$, $F(2,0,0,0,0)=6$, and $F(3,0,0,0,0)=20$, we obtain the solution

[TABLE]

Now we multiply Eq. (41) by $x^{N_{M}}$ and sum over $N_{M}\geqslant 3$, we obtain the recursive equation Eq. (90) for the partial generating function $A(x;N_{B},N_{\bar{B}})$ defined in Eq. (14). In the special case with $N_{B}=0$ in Eq. (90) we obtain

[TABLE]

where we have used

[TABLE]

In the same way, we derive $A(x;N_{B},0)$ whose result is given by Eq. (43) by the replacement $N_{\bar{B}}\rightarrow N_{B}$. We also obtain two summation properties

[TABLE]

We multiply Eq. (90) by $y^{N_{B}}z^{N_{\bar{B}}}$ and take sums over $N_{B}\geqslant 1$ and $N_{\bar{B}}\geqslant 1$, then we can solve $A(x,y,z)$ as

[TABLE]

The detailed derivation of (47) is shown in (88-93).

We can extract $F(N_{M},N_{B},N_{\bar{B}},0,0)$ as the coefficient of $x^{N_{M}}y^{N_{B}}z^{N_{\bar{B}}}$ in polynomial expansion of $A(x,y,z)$. The result is a sum of four terms,

[TABLE]

where these terms are given by Eq. (173) in Appendix D.

4 Particle number distribution for baryons and mesons

In a system of $N_{q}$ quark and $N_{\bar{q}}$ antiquark, we obtain in Eq. (33) the probability $P(N_{M},N_{B},N_{\bar{B}})$ to form $N_{M}$ mesons, $N_{B}$ baryons and $N_{\bar{B}}$ antibaryons. A general form of the raw moments of meson and baryon numbers is

[TABLE]

where we use an overline to denote the average at fixed quark and antiquark numbers. The central moments for baryons and mesons are related by the quark number conservation

[TABLE]

where $\delta N_{i}\equiv N_{i}-\overline{N_{i}}$ with $i=M,B,\bar{B}$.

In Fig. 4, we show the ratios of the cumulants for antibaryons as functions of quark-antiquark asymmetry

[TABLE]

with the original QCR and gQCR at two values of total quark number $x=N_{q}+N_{\bar{q}}$. The cumulants for antibaryons are given by

[TABLE]

Note that the first order cumulant for baryons is related to that for antibaryons by $C_{1}^{B}=C_{1}+zx/3$. As $x$ is not small ($x\gtrsim 200$ ), $C_{1}/x$, $C_{2}/C_{1}$, $C_{3}/C_{2}$ and $C_{4}/C_{2}$ are almost independent of $x$ (i.e., the system size) and are only mainly dependent on the quark-antiquark asymmetry $z$ (proportional to the net-baryon density). Therefore, we plot their $z$ dependence under the original and generalized QCR respectively in Fig. 4. Panel (a) shows antibaryons are less produced with gQCR than that the QCR. The antibaryon production with the gQCR is more suppressed than the QCR for larger $z$. We can parameterize the antibaryon number [54] as $\overline{N_{\bar{B}}}/x=\left(z/3\right)\left(1-z\right)^{a}/\left[\left(1+z\right)^{a}-\left(1-z\right)^{a}\right]$ with $a=3$ for the QCR and $a=5$ for the gQCR. We see that the cumulant ratios $C_{2}/C_{1}$, $C_{3}/C_{2}$ and $C_{4}/C_{2}$ tend to unity at large $z$. This is because the antiquark number is small at large $z$ and the aggregation of quarks and antiquarks in phase space to form baryons and antibaryons is more stochastic and follows the Poisson distribution. Because antibaryons are less produced in the gQCR and thus the feature of the Poisson distribution is more obvious, which can be seen by that the cumulant ratios in the gQCR approaches unity more quickly at large $z$ than in the QCR. At small $z$, $C_{3}/C_{2}$ and $C_{4}/C_{2}$ are small and approach the Gaussian distribution.

\figcaption

The ratios of cumulants for baryon or antibaryon number as functions of quark-antiquark asymmetry $z$ with the original QCR and gQCR at two different values of total quark number $x$: $C_{1}/x$, $C_{2}/C_{1}$, $C_{3}/C_{2}$ and $C_{4}/C_{2}$.

5 Multiplicity property of identified hadrons

Following the method of Ref. [35, 55], we can obtain some multiplicity properties of identified hadrons by taking advantage of the stochastic combination rule. In this paper, we only consider the production of octet baryons with $J^{P}=(1/2)^{+}$ and decuplet baryons with $J^{P}=(3/2)^{+}$, pseudo-scalar mesons with $J^{P}=0^{-}$ and vector mesons $J^{P}=1^{-}$ in the flavor SU(3) ground state. The mean values of the multiplicities for identified baryons and mesons are given by

[TABLE]

with

[TABLE]

where $f$ denotes the quark or antiquark flavors in the baryon $B_{i}$ and the meson $M_{i}$, $n_{f,B_{i}}$ and $n_{f,M_{i}}$ are the number of flavor $f$ quarks in $B_{i}$ and $M_{i}$ respectively, $S_{B_{i}}$ counts the number of different permutations in quarks and antiquarks, $g_{B_{i}}$ and $g_{M_{i}}$ are spin selection factors for $B_{i}$ and $M_{i}$ respectively. We introduce a parameter $R_{V/P}$ to denote the relative weight of vector mesons to pseudo-scalar mesons with the same quark content, and introduce parameter $R_{D/O}$ to denote the relative weight of decuplet baryons to octet baryons with the same quark content. We take $R_{V/P}=0.45$ and $R_{D/O}=0.5$ from studying the hadronic yields in relativistic high energy collisions [56, 58]. Because the values of the two parameters are extracted from the experimental data of hadronic yields, we emphasize that the two parameters can absorb the effects/contribution of various excited states and higher-mass resonances to a certain extent.

Taking the proton as an example, we have $N_{\mathrm{p}}^{\left(q\right)}=3N_{u}\left(N_{u}-1\right)N_{d}$ for all possible combinations of $uud$. Then the ratio $N_{\mathrm{p}}^{\left(q\right)}/N_{q}\left(N_{q}-1\right)\left(N_{q}-2\right)$ stands for the formation probability of the proton having the quark flavor $uud$. The spin selection factor for the proton is $g_{\mathrm{p}}=1/(1+R_{D/O})$.

We introduce the configuration of multiple hadrons as $\left\{k_{\mathrm{h}_{1}}\mathrm{h}_{1},\,k_{\mathrm{h}_{2}}\mathrm{h}_{2},\,...\right\}\equiv\left\{k_{\mathrm{h}}\mathrm{h}\right\}$, where $\mathrm{h}_{i}$ denotes one hadron and $k_{\mathrm{h}_{i}}$ its number in the configuration. We define the joint moment of the multiplicity for such a configuration as $N_{\left\{k_{h}h\right\}}=\prod_{\mathrm{h}_{i}}N_{\mathrm{h}_{i}}^{\underline{k_{\mathrm{h}_{i}}}}$, where we have used the falling factorial of the hadron multiplicity $N_{\mathrm{h}}^{\underline{k_{\mathrm{h}}}}=\prod_{n=1}^{k_{\mathrm{h}}}\left(N_{\mathrm{h}}-n+1\right)$. The average value of the joint moment of the multiplicity reads

[TABLE]

where $k_{B}=\sum_{\mathrm{h}_{i}}k_{\mathrm{h}_{i}}Q_{B,\mathrm{h}_{i}}$ counts the number of baryons in the multi-hadron configuration with $Q_{B,\mathrm{h}_{i}}=1$ for $\mathrm{h}_{i}$ being a baryon and $Q_{B,\mathrm{h}_{i}}=0$ for $\mathrm{h}_{i}$ being a meson and an antibaryon, similarly $k_{\bar{B}}$ counts the number of antibaryons and $k_{M}$ the number of mesons in the multi-hadron configuration, $k_{q}=\sum_{\mathrm{h}_{i}}k_{\mathrm{h}_{i}}Q_{q,\mathrm{h}_{i}}$ counts the number of constituent quarks in the multi-hadron configuration with $Q_{q,\mathrm{h}_{i}}=3$ for $\mathrm{h}_{i}$ being a baryon, $Q_{q,\mathrm{h}_{i}}=0$ for $\mathrm{h}_{i}$ being an antibaryon and $Q_{q,\mathrm{h}_{i}}=1$ for $\mathrm{h}_{i}$ being a meson, similarly $k_{\bar{q}}$ counts the number of antiquarks. The numerator in (55) is defined as

[TABLE]

denotes the number of all possible combinations out of all quarks with specific flavors in the hadrons in the multi-hadron configuration, where $f$ runs over $u,d,s,\bar{u},\bar{d},\bar{s}$ and $n_{f}=\sum_{\mathrm{h}_{i}}k_{\mathrm{h}_{i}}n_{f,\mathrm{h}_{i}}$ counts the number of $f$ flavor quarks or antiquarks in the multi-hadron configuration.

We take a few examples of multi-hadron configurations to illustrate Eq. (55). The first example is the configuration with two protons, we have $k_{\mathrm{p}}=2$ and $N_{\mathrm{p}}^{\underline{2}}\equiv N_{\mathrm{p}}\left(N_{\mathrm{p}}-1\right)$, so we have

[TABLE]

where $\overline{N_{B}\left(N_{B}-1\right)}$ denotes the average for all possible baryon pairs, $N_{\mathrm{pp}}^{\left(q\right)}=3^{2}N_{u}^{\underline{4}}N_{d}^{\underline{2}}$ is the number of all possible combinations out of six quarks with specific flavors in two protons. The ratio $N_{\mathrm{pp}}^{\left(q\right)}/N_{q}^{\underline{6}}$ gives the probability of the six-quark combination having the particular flavor structure $(uud)(uud)$. The second example is the configuration $\mathrm{pp\bar{n}}$, we have

[TABLE]

where the first parentheses in the second line gives the probability of two baryons having the flavor and spin structure of two protons and the second parentheses gives the probability of an antibaryon having the flavor and spin structure of the antineutron.

The all orders of moments and correlation functions of hadron multiplicity can be built from Eq. (55). The two-body correlation function reads

[TABLE]

where $\alpha$, $\beta$ denote two hadrons, and we have used $\overline{N}_{\alpha\beta}\equiv\overline{N_{\alpha}N_{\beta}}$ for $\alpha\neq\beta$ and $\overline{N}_{\alpha\beta}\equiv\overline{N_{\alpha}^{\underline{2}}}$ for $\alpha=\beta$. The three-body correlation function reads

[TABLE]

where $\overline{N_{\alpha}N_{\beta}N_{\gamma}}$ can be expressed by falling factorials,

[TABLE]

Here we have used $\overline{N}_{\alpha\beta\gamma}\equiv\overline{N_{\alpha}N_{\beta}N_{\gamma}}$ for $\alpha\neq\beta\neq\gamma$, $\overline{N}_{\alpha\beta\gamma}\equiv\overline{N_{\alpha}^{\underline{2}}N_{\gamma}}$ for $\alpha=\beta\neq\gamma$ and $\overline{N}_{\alpha\beta\gamma}\equiv\overline{N_{\alpha}^{\underline{3}}}$ for $\alpha=\beta=\gamma$. Note that $\overline{N}_{\alpha\beta\gamma}$ is symmetric for any permutation of $\alpha$, $\beta$ and $\gamma$. The four-body correlation function can be written as

[TABLE]

where $\overline{N_{\alpha}N_{\beta}N_{\gamma}N_{\epsilon}}$ can be expressed by falling factorials

[TABLE]

Here we have used $\overline{N}_{\alpha\beta\gamma\epsilon}\equiv\overline{N_{\alpha}N_{\beta}N_{\gamma}N_{\epsilon}}$ for $\alpha\neq\beta\neq\gamma\neq\epsilon$, $\overline{N}_{\alpha\beta\gamma\epsilon}\equiv\overline{N_{\alpha}^{\underline{2}}N_{\gamma}N_{\epsilon}}$ for $\alpha=\beta\neq\gamma\neq\epsilon$, $\overline{N}_{\alpha\beta\gamma\epsilon}\equiv\overline{N_{\alpha}^{\underline{3}}N_{\epsilon}}$ for $\alpha=\beta=\gamma\neq\epsilon$ and $\overline{N}_{\alpha\beta\gamma\epsilon}\equiv\overline{N_{\alpha}^{\underline{4}}}$ for $\alpha=\beta=\gamma=\epsilon$. Note that $\overline{N}_{\alpha\beta\gamma\epsilon}$ is symmetric for any permutation of $\alpha$, $\beta$, $\gamma$ and $\epsilon$.

The cumulants of the net proton number $N_{\mathrm{p}}-N_{\bar{\mathrm{p}}}$ can be calculated by combinations of above multi-body correlation functions (except $C_{1}$ which is just the mean value of the net proton number)

[TABLE]

In Fig. 5, we show the cumulant ratios for the net proton number with the QCR and the gQCR at given total quark number $x=2000$ as functions of the quark-antiquark asymmetry $z$. We have checked that the cumulant ratios of the net proton number are independent of $x$ for large $x$. Here we assume the number of strange quarks is $N_{s}=N_{\bar{s}}=0.45N_{\bar{u}}$, where the strangeness conservation is satisfied and strangeness suppression factor 0.45 is consistent with the observation in relativistic heavy-ion collisions [27]. Because the net-baryon number is fixed at given numbers of quarks and antiquarks, Fig. 5 only show the fluctuations of the net proton number brought by the quark combination process. $C_{1}/C_{2}$ and $C_{3}/C_{2}$ of net proton number increase with $z$ and approach to 1.5 and $1/3$ as $z\to 1$, respectively. $C_{4}/C_{2}$ decreases with $z$ and approaches to $-1/3$ as $z\to 1$. These results are different from those of the statistical model for hadron resonance gas with thermal equilibrium [57], in which $C_{1}/C_{2}$ and $C_{3}/C_{2}$ increase to one and $C_{4}/C_{2}$ almost keeps a constant of one at large baryon number chemical potential.

\figcaption

The cumulant ratios of the net proton number directly-produced through the QCR and the gQCR as functions of quark-antiquark asymmetry $z$ at $x=2000$.

6 Effects of quark number fluctuations at hadronization

The numbers of quarks and antiquarks may fluctuate in heavy energy collisions in event by event, so hadronic observables should be influenced by the fluctuations at the quark level. We denote the distribution of quark numbers at hadronization as

[TABLE]

and the event average of a hadronic observable $A_{\mathrm{h}}$ is given by

[TABLE]

where $\overline{A}_{\mathrm{h}}\left(\{N_{f}\}\right)$ denote the average at fixed quark and antiquark numbers which is obtained in the previous section.

In practical evaluation, it is more convenient to take the expansion of $\overline{A}_{\mathrm{h}}\left(\{N_{f}\}\right)$ in $N_{f}$ around the event-average numbers of quarks and antiquarks $\{\langle N_{f}\rangle\}$ as

[TABLE]

where $\delta N_{f}=N_{f}-\langle N_{f}\rangle$ and the subscript ’0’ denotes the values at $\langle N_{f}\rangle$. Substituting it into Eq. (66), we obtain

[TABLE]

which involves multi-body correlations for the quark and antiquark number with the quark number distribution $P\left(\{N_{f}\}\right)$

[TABLE]

where we have used the same symbols $C$ as in Sec. 5 to denote the quark number correlation functions but with quark flavor indices.

Using the above moment expansion method, we can study, for the selected phase-space window such as midrapidity region $|y|<0.5$, the influence of different quark and antiquark distributions in the window on the production of hadrons at hadronization. This extends the applicability of QCR described in Sec. 2 and 3 where only the stochastic quark-antiquark population is considered.

7 Influence of resonance decays

The long life (stable) hadrons measured in experiments include contributions from resonance decays. In this section we consider the influence of resonance decays to the multiplicity correlations of long life hadrons. For a resonance $\mathrm{h}_{i}$, its stable daughters are denoted as $\alpha$, $\beta$, $\gamma$, etc., and the corresponding decay branching ratios are $D_{i\alpha}$, $D_{i\beta}$, $D_{i\gamma}$, etc., respectively. The distribution function of $\left\{N_{\alpha},N_{\beta},N_{\gamma},\cdots\right\}$ from decays of the resonance $\mathrm{h}_{i}$ with the particle number $N_{\mathrm{h}_{i}}$ is denoted as $f\left(N_{\mathrm{h}_{i}},\left\{N_{\alpha}^{i}\right\}\right)$ where $N_{\alpha}^{i}$ denotes number of the stable hadron $\alpha$ from the resonance $\mathrm{h}_{i}$. We take the multi-nominal distribution for selection of decay channels. Convoluting with the joint distribution of directly produced hadrons $P\left(\left\{N_{\mathrm{h}_{i}}\right\},\left\{\left\langle N_{f}\right\rangle\right\}\right)$ at hadronization, we get the joint distributions of stable hadrons $F\left(\left\{N_{\alpha},N_{\beta},N_{\gamma},\cdots\right\}\right)$ [35].

From the joint distribution functions of stable hadrons, we get the average yield of a stable hadron $\alpha$

[TABLE]

where the index $i$ runs over all directly produced hadron species including stable hadrons (we define $D_{\alpha\alpha}=1$), and we have used shorthand notation $N_{i}\equiv N_{\mathrm{h}_{i}}$. The two-body multiplicity correlation function is

[TABLE]

The three-body multiplicity correlation function of stable hadrons is

[TABLE]

The four-body multiplicity correlation function of stable hadrons is

[TABLE]

The higher order multiplicity correlation functions can similarly be derived from the joint distribution functions of stable hadrons.

8 Application: cumulants for net protons in heavy ion collisions

In this section, we take a simple example of a quark system and calculate the cumulant ratios for net protons in the final state with gQCR. We will discuss our results in the context of the experimental observation in Au+Au collisions at RHIC.

We consider a quark system having the property $S\sigma\equiv C_{3}/C_{2}=z$ and $\kappa\sigma^{2}\equiv C_{4}/C_{2}=1$ for the baryon quantum number, which is similar to the base line in the grand canonical ensemble in statistical model [59]. A simple case for the quark number correlation functions satisfying the above property is

[TABLE]

for diagonal elements and non-vanishing off-diagonal elements ($f_{1}\neq f_{2}$)

[TABLE]

We also assume following properties for correlation functions of the strangeness as a result of strangeness conservation

[TABLE]

Because total quark number $x$ in midrapidity region $|y|<0.5$ in central heavy-ion collisions at RHIC energy is usually large ($x\gtrsim$500), the cumulant ratios for net protons in the final state in our model is not sensitive to $x$ (the system size). Here, we just take $x=5000$ in the calculation. We use Eq. (70) to fit the data of the $\bar{\mathrm{p}}/\mathrm{p}$ yield ratio in Au+Au collisions and determine $z$ in the collision energy range $\sqrt{s_{NN}}\in[7,200]$ GeV.

In Fig. 8, we show the results for the cumulant ratios for net protons in the final state at different collisional energies. Our results only incorporate the contributions of quark number fluctuations and correlations up to fourth order as in Eq. (68). The auxiliary horizontal axis shows the corresponding quark-antiquark asymmetry parameter $z$. Results including the contributions from strong and electromagnetic (S&EM) decays of resonances are in dashed lines, while results with full decay contributions including weak decays are in solid lines. The STAR data are shown in solid circles with error bars.

\figcaption

The cumulant ratios for net proton at different collision energies. Solid circles with error bars are experimental data [43, 60]. Solid and dashed lines are theoretical results of SDQCM.

We see in the figure that the cumulant ratios for net protons as functions of collisional energies in our model can describe the experimental data: $M/\sigma^{2}$ and $S\sigma$ increase with $z$ or equivalently decrease with collisional energies, while $\kappa\sigma^{2}$ decreases with collisional energies until it reaches a minimum at $\sqrt{s_{NN}}=$20 GeV and then increases toward unity at high energies. Our results for $\kappa\sigma^{2}$ are the consequence of the competition between two effects: (1) The cumulant ratio $C_{4}/C_{2}$ from directly produced net protons by quark combination, as shown in Fig. 5, always decreases with increasing $z$; (2) The cumulant ratio $C_{4}/C_{2}$ for baryons or antibaryons, as shown in Fig. 4, rapidly increases as $z\gtrsim 0.2$ (corresponding to $\sqrt{s_{NN}}\lesssim 20$ GeV).

We now make some remarks about the results for $\kappa\sigma^{2}$ of net protons. Although our results for net protons can reproduce the nontrivial dependence on collisional energies, we always have $\kappa\sigma^{2}=1$ for net baryons at all collisional energies due to the flavor/charge conservation in quark combination. Therefore, our results indicate that cumulant ratios of net proton number in this simple case do not exactly follow those of the net baryon number. This is different from that of the statistical model [57] for hadron resonance gas with thermal equilibrium which predicts the quite similar cumulant ratios for net-protons and net-baryons.

We emphasize that the current results are preliminary and are mainly used to show the potential application of our model in hadronic fluctuations and the related phase transition in relativistic heavy-ion collisions. There are some limitations in the current calculations that should be clarified. In this paper, we only consider the production of ground-state hadrons. Effects of higher-mass resonances are only partially absorbed by parameters $R_{V/P}$ and $R_{D/O}$. According to our previous study [35], the production of hadrons with small yields tends to follow Poisson distribution. Including those higher-mass resonances with small yields may enlarge the fourth moments of protons to a certain extent. Our model is a static one and does not consider the diffusion of hadrons/quarks during the finite hadronization time, which will cause some influence for the calculation of net-proton fluctuations in a specific window. In addition, in order to make final comparison with data of net protons, we should consider more realistic quark number fluctuations and correlations in the studied rapidity window, which may be obtained by grand-canonical statistics of thermal quark system or by canonical statistics with Bernoulli trial selection of the specific window. We should also consider other effects related to the finite acceptance window such as the diffusion/blur of charges during the hadronic scatterings stage as well as that caused by resonance decay [61, 62, 63, 64, 57, 65, 66, 67, 68, 69, 70, 55, 71, 72, 73]. We will study these effects in this framework in the future.

9 Summary

We developed a new statistical method to solve the probability distribution for the numbers of baryons, antibaryons and mesons formed in hadronization of a constituent quark and antiquark system governed by the quark combination rule (QCR) in the quark combination model developed by Shandong Group. We use a set of numbers $(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ to classify the outcome of implementing the QCR for a queue of $N_{q}$ quarks and $N_{\bar{q}}$ antiquarks, where there are $N_{M}$ mesons, $N_{B}$ baryons and $N_{\bar{B}}$ antibaryons formed by combination but with $N_{r}$ quarks and $N_{\bar{r}}$ antiquarks left without forming hadrons. The number of different ways of a given configuration $(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ is denoted as $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$. We build the recursive relation for $F(N_{M},N_{B},N_{\bar{B}},0,0)$. We adopt the generating function method to solve the recursive relation and give the analytical expression of $F(N_{M},N_{B},N_{\bar{B}},0,0)$. This method is much simpler than the previous one and is easy to generalize. To accommodate the baryon yield in experimental data in relativistic heavy-ion collisions, we consider a generalized combination rule (gQCR). We give the solution of baryon and meson probability distribution function under the new rule. Because the (anti)baryon production is more suppressed in the gQCR, we find that the cumulant ratios of the antibaryon number reach the Poisson statistics more rapidly in the gQCR than in the QCR with increasing baryon number density.

We studied the multiplicity fluctuation and correlation functions for identified baryons directly produced in collisions. We also studied correlation functions of final state hadrons including contributions from resonance decays. As an illustrative example we consider a quark-antiquark system having the property $S\sigma=z$ and $\kappa\sigma^{2}=1$ for the baryon quantum number, which is similar to the base line in the grand-canonical ensemble in the statistical model. We calculate the cumulant ratios in the quark combination model and find that $S\sigma$ for net protons in the model increases with decreasing collisional energies, which is consistent with the experimental observation. More interesting is that $\kappa\sigma^{2}$ for net protons has a nontrivial energy behavior consistent with the data.

Appendix A Derivation of two properties for gQCR

We present two properties for $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ in this appendix.

Property 1. For $(N_{r},N_{\bar{r}})=(1,0),(0,1)$, $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ can be expressed in terms of $F(N_{M}^{\prime},N_{B}^{\prime},N_{\bar{B}}^{\prime},0,0)$ with $N_{M}^{\prime}\leqslant N_{M}$, $N_{B}^{\prime}\leqslant N_{B}$ and $N_{\bar{B}}^{\prime}\leqslant N_{\bar{B}}$,

[TABLE]

where the coefficients $C_{n_{B},n_{\bar{B}}}^{10}(n_{M})$ are given by $C_{0,0}^{10}(n_{M})=1$ and

[TABLE]

for $a+b>0$. The coefficients $C_{n_{B},n_{\bar{B}}}^{01}(n_{M})$ are given by $C_{0,0}^{01}(n_{M})=1$ and

[TABLE]

for $a+b>0$. Note that the sums over $n_{B}$, $n_{\bar{B}}$ and $n_{M}$ in Eq. (77) continues until any of the baryon, antibaryon and meson number become negative since $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})=0$ for the case that any of $N_{M}$, $N_{B}$, $N_{\bar{B}}$, $N_{r}$ or $N_{\bar{r}}$ are negative.

Property 2. For $(N_{r},N_{\bar{r}})=(2,0),(0,2)$, $F(N_{M},N_{B},N_{\bar{B}},N_{r},N_{\bar{r}})$ can be expressed in terms of $F(N_{M}^{\prime},N_{B}^{\prime},N_{\bar{B}}^{\prime},0,0)$ with $N_{M}^{\prime}\leqslant N_{M}$, $N_{B}^{\prime}\leqslant N_{B}$ and $N_{\bar{B}}^{\prime}\leqslant N_{\bar{B}}$,

[TABLE]

where the coefficients are identical to Eqs. (78,79): $C_{n_{B},n_{\bar{B}}}^{20}(n_{M})=C_{n_{B},n_{\bar{B}}}^{10}(n_{M})$ and $C_{n_{B},n_{\bar{B}}}^{02}(n_{M})=C_{n_{B},n_{\bar{B}}}^{01}(n_{M})$.

In order to prove the two properties in Eq. (77) and Eq. (80), we can solve $F(N_{M},N_{B},N_{\bar{B}},2,0)$ by substituting the first equation into the third one in Eq. (40),

[TABLE]

In the same way, we can solve $F(N_{M},N_{B},N_{\bar{B}},0,2)$ by substituting the second equation into the fourth one in Eq. (40)

[TABLE]

From Eqs. (81,82) we obtain Eq. (80). Then from Eq. (80) and the last two equations of (40) we obtain Eq. (77).

Appendix B Derivation of recursive equation Eq. (41)

We derive the recursive equation for $F(N_{M},N_{B},N_{\bar{B}},0,0)$ with the help of Eq. (40). We start from Eq. (9) which also holds for gQCR. Using Eq. (40), Eq. (9) can be rewritten as

[TABLE]

where $G$ is defined by

[TABLE]

We need to define another auxiliary function

[TABLE]

where we have used the first two equalities of Eq. (40) to obtain second line and used Eq. (83) to obtain the last line. On the other hand, we can also derive another form of $H(N_{M},N_{B},N_{\bar{B}})$ by applying the last two equalities and then the first two equalities of Eq. (85) and finally applying Eq. (40), the result is

[TABLE]

From Eq. (85), we also obtain

[TABLE]

By equating Eq. (86) and Eq. (87) we derive the recursive equation for $G$ and then for $F(N_{M},N_{B},N_{\bar{B}},0,0)$ through Eq. (83) as in Eq. (41).

Appendix C Derivation of generating functions for gQCR

We multiply Eq. (41) by $x^{N_{M}}$ and sum over $N_{M}\geqslant 3$ (the equation is well defined for $N_{M}\geqslant 3$), we obtain

[TABLE]

One can verify that a complete cancellation occurs for terms proportional to $x^{2}$ and those proportional to $x$ after applying Eq. (41) for $N_{M}=2$ and $N_{M}=1$ respectively. The constant terms in $F$ reads

[TABLE]

where we have used the initial value $F(0,N_{B},N_{\bar{B}},0,0)=\delta_{N_{B}N_{\bar{B}},0}$. This means that if there are no mesons, baryons and antibaryons cannot coexist, and if there are only quarks or antiquarks in the system the number of different queues is 1. Then Eq. (88) is simplified as

[TABLE]

We now multiply Eq. (90) by $y^{N_{B}}z^{N_{\bar{B}}}$ and take a sum over $N_{B}\geqslant 1$ and $N_{\bar{B}}\geqslant 1$ to obtain

[TABLE]

To simplify the above equation, we use

[TABLE]

and

[TABLE]

as well as Eqs. (45) and (46), we finally obtain Eq. (47).

Appendix D Derivation of Eq. (48) in gQCR

From Eq. (47) we can extract $F(N_{M},N_{B},N_{\bar{B}},0,0)$, the coefficient of $x^{N_{M}}y^{N_{B}}z^{N_{\bar{B}}}$ in the polynomial expansion of $A(x,y,z)$. To this end we rewrite $A(x,y,z)$ as

[TABLE]

where we have used the multinomial theorem

[TABLE]

with multinomial coefficients

[TABLE]

Then we can extract the coefficient of $y^{N_{B}}z^{N_{\bar{B}}}$ as

[TABLE]

where $C_{1,2,3,4}$ are given by

[TABLE]

Note that we have expressed the summation indices $j$ and $k$ in terms of $N_{B}$ and $N_{\bar{B}}$ for a given summation index $i$. The sum over $i$ involves a finite number of terms. For example, in $C_{1}$, we have $N_{B}+N_{\bar{B}}\geqslant i\geqslant\mathrm{Max}(N_{B},N_{\bar{B}})$, since any terms in the summation with the factorial of a negative integer in the denominator are vanishing.

We can factorize the following two polynomials as

[TABLE]

and expand $C_{1}$, $C_{2}$, $C_{3}$, $C_{4}$ into a power series of $x$ with the help of the binomial theorem. Then we can extract the coefficient of $x^{N_{M}}$ in $C(y^{N_{B}}z^{N_{\bar{B}}})$. This gives the coefficient of $x^{N_{M}}y^{N_{B}}z^{N_{\bar{B}}}$ in $A(x,y,z)$ or $F(N_{M},N_{B},N_{\bar{B}},0,0)$ as the sum of the following four terms

[TABLE]

which give the final result of gQCR in Sec. 3.