Analytical coalescence formula for particle production in relativistic heavy-ion collisions

TL;DR

This paper derives an improved analytical coalescence formula for particle yields in relativistic heavy-ion collisions, incorporating relativistic effects and finite size corrections, aiding predictions of cluster production such as light nuclei and exotic hadrons.

Contribution

It introduces a new analytical coalescence formula that accounts for longitudinal momentum, relativistic corrections, and finite size effects, enhancing yield predictions in heavy-ion collisions.

Findings

The formula accurately predicts light nuclei yields in Pb+Pb collisions.

Application to strangeness factor $S_3$ shows consistency with experimental data.

Predictions for exotic hadron production align with other models.

Abstract

Based on a covariant coalescence model with a blast-wave-like parametrization for the phase-space configuration of constituent particles at freeze-out, we derive an approximate analytical formula for the yields of clusters produced in relativistic heavy-ion collisions. Compared to previous existing formulae, the present work additionally considers the contributions from the longitudinal dimension in momentum space, the relativistic corrections and the finite size effects of the produced clusters relative to the spatial distribution of constituent particles at freeze-out. The new analytical coalescence formula provides a useful tool to evaluate the yield of produced clusters, such as light nuclei from nucleon coalescence and hadrons from quark coalescence, in heavy-ion collisions. As a first application of the new analytical formula, we explore the strangeness population factor $S_3 = ^3_{\Lambda}$H/($^3$He$\times \Lambda/$p) based on nucleon/$\Lambda$ coalescence as well as the production of exotic hadrons based on quark coalescence, in central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV. The results are compared with the predictions from other models.

Tables

Table 1. Table 1: Model parameters of COAL-SH, COAL-Ex and the thermal model, i.e., freeze-out temperature T 𝑇 T (GeV), fireball volume V 𝑉 V  (10 3 fm 3 ), harmonic oscillator frequency w 𝑤 w (GeV), w s subscript 𝑤 𝑠 w_{s} (GeV) and w c subscript 𝑤 𝑐 w_{c} (GeV), baryon chemical potential μ B subscript 𝜇 𝐵 \mu_{B} , strangeness chemical potential μ S subscript 𝜇 𝑆 \mu_{S} , and charm fugacity γ c subscript 𝛾 𝑐 \gamma_{c} , for central Pb+Pb collisions at s = 2.76 𝑠 2.76 \sqrt{s}=2.76  TeV.

	$T$	$V$	Nu/3	Ns/3	Nc/3	$w$	$w_s$	$w_c$
COAL-SH	0.154	13.8	462	196	27	0.184	0.078	0.087
COAL-Ex	0.154	6.06	267	156	20	0.55	0.519	0.385
	$T$	$V$	${\mu }_B$	${\mu }_S$	${\gamma }_c$
Thermal	0.153	6.45	0	0	61.6

Table 2. Table 2: d ​ N / d ​ y 𝑑 𝑁 𝑑 𝑦 dN/dy at midrapidity of p 𝑝 p , Λ Λ \Lambda , ϕ italic-ϕ \phi , Ξ − superscript Ξ \Xi^{-} , Ω − superscript Ω \Omega^{-} and D 0 superscript 𝐷 0 D^{0} in central Pb+Pb collisions at s = 2.76 𝑠 2.76 \sqrt{s}=2.76  TeV predicted from COAL-SH, COAL-Ex and the thermal model. The experimental data and the coalescence factor g c subscript 𝑔 𝑐 g_{c} are also included.

Hadron	$p$	$\varphi$	${\Xi }^{-}$	${\Omega }^{-}$	$\Lambda$	$D^0$
$g_c$	$\frac{2}{3^3\times 2^3}$	$\frac{3}{3^2\times 2^2}$	$\frac{2}{3^3\times 2^3}$	$\frac{4}{3^3\times 2^3}$	$\frac{2}{3^3\times 2^3}$	$\frac{1}{3^2\times 2^2}$
Exp.	$34$	$13.8$	$3.34$	$0.58$	$26$	8.4
	$\pm 3$	$\pm 1.77$	$\pm 0.25$	$\pm 0.1$	$\pm 3$	-
COAL-SH	34.4	13.8	3.23	0.64	26.3	8.4
COAL-Ex	36.2	13.8	3.16	0.71	22.5	8.4
Thermal	31.1	16.7	3.48	0.60	19.3	8.4

Table 3. Table 3: d ​ N / d ​ y 𝑑 𝑁 𝑑 𝑦 dN/dy at midrapidity for some exotic hadrons with various quark configurations in central Pb+Pb collisions at s = 2.76 𝑠 2.76 \sqrt{s}=2.76  TeV predicted from COAL-SH, COAL-Ex and the thermal model.

Models	COAL-SH	COAL-Ex	COAL-SH	COAL-Ex	Thermal
Quark configuration	2q/3q/6q	2q/3q/6q	4q/5q/8q	4q/5q/8q	-
Mesons:
$f_0\left(980\right)$	35.5,3.98($s\overline{s}$)	7.24,1.14($s\overline{s}$)	5.90$\times {10}^{-2}$	5.54$\times {10}^{-2}$	6.85
$a_0\left(980\right)$	72.4	19.6	8.54$\times {10}^{-2}$	13.7$\times {10}^{-2}$	20.6
$K\left(1460\right)$	-	-	4.03$\times {10}^{-1}$	3.35$\times {10}^{-1}$	10.0$\times {10}^{-1}$
$D_s\left(2371\right)$	3.38$\times {10}^{-1}$	1.31$\times {10}^{-1}$	3.76$\times {10}^{-3}$	7.52$\times {10}^{-3}$	2.16$\times {10}^{-1}$
$T_{cc}$	-	-	6.28$\times {10}^{-4}$	9.04$\times {10}^{-4}$	5.17$\times {10}^{-3}$
$X\left(3872\right)$	5.37$\times {10}^{-2}$	6.34$\times {10}^{-3}$	6.28$\times {10}^{-4}$	9.04$\times {10}^{-4}$	3.26$\times {10}^{-3}$
$Z^{+}\left(4430\right)$	-	-	5.35$\times {10}^{-4}$	2.68$\times {10}^{-4}$	1.01$\times {10}^{-4}$
Baryons:
$\Lambda \left(1405\right)$	3.06	0.75	4.35$\times {10}^{-2}$	3.46$\times {10}^{-2}$	1.36
${\Theta }^{+}$	-	-	3.76$\times {10}^{-2}$	0.86$\times {10}^{-2}$	6.75$\times {10}^{-1}$
$\overline{K}KN$	-	-	2.14$\times {10}^{-2}$	0.55$\times {10}^{-2}$	1.44$\times {10}^{-1}$
$\overline{D}N$	-	-	2.64$\times {10}^{-3}$	6.28$\times {10}^{-3}$	2.56$\times {10}^{-2}$
${\overline{D}}^{\ast }N$	-	-	3.86$\times {10}^{-3}$	1.32$\times {10}^{-3}$	2.36$\times {10}^{-2}$
${\Theta }_{cs}$	-	-	1.22$\times {10}^{-3}$	1.12$\times {10}^{-3}$	1.63$\times {10}^{-2}$
Dibaryons:
$H$	6.30$\times {10}^{-4}$	5.32$\times {10}^{-4}$	-	-	5.37$\times {10}^{-3}$
$\overline{K}NN$	3.91$\times {10}^{-3}$	8.46$\times {10}^{-3}$	4.49$\times {10}^{-5}$	3.01$\times {10}^{-5}$	5.71$\times {10}^{-3}$
$\Omega \Omega$	3.49$\times {10}^{-6}$	4.56$\times {10}^{-6}$	-	-	1.47$\times {10}^{-5}$
$H_c^{++}$	1.00$\times {10}^{-4}$	3.52$\times {10}^{-4}$	-	-	1.10$\times {10}^{-3}$
$\overline{D}$NN	-	-	2.12$\times {10}^{-6}$	1.18$\times {10}^{-6}$	8.24$\times {10}^{-5}$

Equations

$$E\frac{\text{d}N}{{\text{d}}^{3}p}=\underset{{\Sigma }^{\mu }}{\int }\text{d}{\sigma }_{\mu }{p}^{\mu }f(x,p)=\int {\text{d}}^{4}xS(x,p),$$

$$S(x,p){\text{d}}^{4}x$$

$$\frac{{\text{d}}^{3}N}{p_T\text{d}{p}_T\text{d}y\text{d}{\varphi }_p}=\frac{g\xi \tau m_T}{(2\pi {)}^3}\int 2K_1(\beta )2\pi I_0(\alpha )rdr,$$

$$\frac{{\text{d}}^{3}N}{p_T\text{d}{p}_T\text{d}y\text{d}{\varphi }_p}$$

$$\frac{dN}{dy}$$

$$K_{\nu }(x)\to \sqrt{\frac{\pi }{2x}}{e}^{-x}(1+\frac{4{\nu }^2-1}{8x}+\mathcal{O}(\frac{1}{x^2})),$$

$$\frac{dN}{dy}$$

$$N=\frac{g}{(2\pi {)}^3}{\xi }^{\prime }{V}^{\prime }{V}_p.$$

$$N_c$$

$$\displaystyle\left(\begin{array}[]{c}\mathbf{R}\\ \mathbf{q}_{1}\\ \cdot\\ \cdot\\ \cdot\\ \mathbf{q}_{N-1}\end{array}\right)=\hat{J}\left(\begin{array}[]{c}\mathbf{x}_{1}\\ \mathbf{x}_{2}\\ \cdot\\ \cdot\\ \cdot\\ \mathbf{x}_{N}\end{array}\right),\left(\begin{array}[]{c}\mathbf{P}\\ \mathbf{k}_{1}\\ \cdot\\ \cdot\\ \cdot\\ \mathbf{k}_{N-1}\end{array}\right)=(\hat{J}^{-1})^{T}\left(\begin{array}[]{c}\mathbf{p}_{1}\\ \mathbf{p}_{2}\\ \cdot\\ \cdot\\ \cdot\\ \mathbf{p}_{N}\end{array}\right),$$

$$\mathbf{R}$$

$${\mathbf{q}}_i$$

$$\prod _{i=1}^N{\text{d}}^3{x}_i{\text{d}}^3{p}_i$$

$$N_c$$

$${\mu }_i=\frac{i+1}{i}\frac{m_{i+1}{\sum }_{k=1}^i{m}_k}{{\sum }_{k=1}^{i+1}{m}_k},(1\le i\le N-1)$$

$$\prod _{i=1}^N{V}^{\prime }{V}_{p,i}$$

$$N_c$$

$$\int {\rho }_{c,i}^W{\text{d}}^3{q}_i$$

$$N_c$$

$$F(\sigma ,\mu ,l,T)$$

$$N_c$$

$$\frac{d{N}_c}{dy}$$

$$g_{\text{rel}}\approx \frac{1+\frac{15}{8}\frac{T}{{\mu }_0}}{{\prod }_{i=1}^N(1+\frac{15}{8}\frac{T}{m_i})}.$$

$$g_{\text{size}}(x,N)\approx \frac{1-(2+0.6(N-1{)}^{0.5})x^{1.3}{e}^{-2x^2}}{1+(N-1)x^2},$$

$$\frac{d{N}_c}{dy}$$

$$\displaystyle\mathcal{S}(l,x)=\Big{(}\frac{x^{2}}{1+x^{2}}\Big{)}^{l}G(l,x),$$

$$\mathcal{S}(l,x)<1,$$

Full text

I Introduction

Particle production in relativistic heavy-ion collisions is of fundamental importance for many issues in nuclear physics, particle physics, astrophysics and cosmology. Coalescence model But61 ; Sat81 ; Cse86 ; Dov91 provides an important approach to describe the particle production in heavy-ion collisions at both intermediate and high energies. For instance, the coalescence model has been successfully and widely applied to describe both light nuclei production from nucleon coalescence Gyu83 ; Mat97 ; ChenLW03 ; Oh09 ; Ste12 ; ChenG12 ; ChenG13 ; Zhu15 ; Zha10 ; Xue15 ; KJSUN15 ; KJSUN16 ; KJSUN16-1 and hadron production from quark coalescence in heavy-ion collisions Lin02 ; Vol03 ; Hwa03 ; Gre03 ; Fri03 ; Mol03 ; Pra05 ; Sha05 ; LWChen06 ; Ko14 ; Fri08 . The quark coalescence provides an important mechanism for the hadronization of partons produced in relativistic heavy-ion collisions and thus is very useful for understanding the partonic dynamics as well as the formation signals and properties of quark-gluon plasma possibly formed in these collisions Fri08 ; Adl03 ; Ada04 . Since the coalescence probability is based on the overlap of density matrix of constituent particles in an emission source with Wigner function of the produced cluster, the predicted cluster yield generally depends on the internal structure wave function of the cluster. Therefore, the coalescence model provides a unique tool to identify the constituent quark structure of some exotic hadrons via quark coalescence from studying their yields in relativistic heavy-ion collisions ChenLW04 ; ChenLW07 ; Cho11 ; Cho11prc .

To describe particle production in heavy-ion collisions, the coalescence model is usually implemented with the phase-space configuration of constituent particles at freeze-out in an emission source obtained from either transport model simulations Gyu83 ; Mat97 ; ChenLW03 ; Oh09 ; Ste12 ; ChenG12 ; ChenG13 ; Zhu15 ; LWChen06 ; Ko14 or empirical parametrization Gre03 ; Fri03 ; ChenLW04 ; ChenLW07 ; Zha10 ; Xue15 ; KJSUN15 ; KJSUN16 ; KJSUN16-1 , and the final results are then obtained by multi-dimension numerical integration. Although the numerical calculations can give exact results in the coalescence model, the approximate analytical formula is extremely useful and has its own merits, e.g., it can significantly simplify the computational work and most importantly it can give more transparent and deep insights into the physics of the coalescence scenario.

In the literature, there indeed exist some analytical coalescence formulae (see, e.g., ChenLW07 ; Cho11 ; Cho11prc ) for cluster yields in relativistic heavy-ion collisions, obtained under some approximations. In particular, the analytical formula reported recently by the ExHIC collaboration in Refs. Cho11 ; Cho11prc (denoted as COAL-Ex in the present work) is quite general for $N$-body cluster production with the $N$ constituent particles at various orbit angular momentum states. In the derivation of the COAL-Ex formula, for simplicity, the integration of Wigner function in longitudinal ($z$-) direction is neglected, namely, the integration in momentum space is only treated in two-dimensional transverse plane although the integration in coordinate space is treated three-dimensionally, which leads to a problematic feature that the yield of the produced cluster in the COAL-Ex formula increases as a power function of its root-mean-square (RMS) radius $r_{\text{rms}}$ when $r_{\text{rms}}$ becomes large. Even though the single-particle momentum distribution for a particle in the emission source may only depend on its transverse momentum due to boost invariance Bjo83 in the longitudinal direction, the many-body coalescence process cannot be treated merely in two-dimensional transverse plane. It is thus of great interest to include the contribution of longitudinal momentum dimension and check its importance.

In the present work, based on a blast-wave-like parametrization, which is inherently longitudinal boost invariant by design, for the emission source of constituent particles at freeze-out, we derive a new analytical coalescence formula (denoted as COAL-SH in the following) by treating the integration of Wigner function in full phase-space for the many-body coalescence process. We find that the COAL-SH formula possesses a nice saturation property that when $r_{\text{rms}}$ of the produced cluster is large, its yield converges to a constant. In addition, we consider the relativistic corrections to leading order in the COAL-SH formula. Furthermore, an empirical expression is proposed for the corrections of finite size effects of the produced clusters relative to the size of the emission source for constituent particles. As a first application of the new analytical coalescence formula, we explore the strangeness population factor $S_{3}=^{3}_{\Lambda}$H/(3He$\times\Lambda/$p) as well as the production of exotic hadrons in central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, and compare the results with the predictions from other models.

II Theoretical formalism

II.1 Blast-wave-like parametrization for emission source of constituent particles at freeze-out

We assume that particles are emitted from a hypersurface $\Sigma^{\mu}$, and then the Lorentz invariant one-particle momentum distribution is given by

[TABLE]

where $\sigma_{\mu}$ denotes the normal vector of the hypersurface $\Sigma^{\mu}$ and $p^{\mu}$ is the four-momentum of the emitted constituent particle. For the cluster production at midrapidity in relativistic heavy-ion collisions that we are considering in this work, we adopt the longitudinal boost invariance assumption Bjo83 and assume the constituent particles are emitted at a fixed proper time $\tau_{0}$, and the emission function can then be expressed as Ret04

[TABLE]

where we use the longitudinal proper time $\tau=\sqrt{t^{2}-z^{2}}$, space-time rapidity $\eta=\frac{1}{2}\text{ln}\frac{t+z}{t-z}$, polar coordinates ($r$, $\phi_{s}$), rapidity $y=\frac{1}{2}\ln(\frac{E+p_{z}}{E-p_{z}})$, transverse momentum ($p_{T},\phi_{p}$) and transverse mass $m_{T}=\sqrt{m^{2}+p_{T}^{2}}$. The statistical distribution function $f(x,p)$ is given by Coo74 $f(x,p)=g(2\pi)^{-3}[\exp(p^{\mu}u_{\mu}/kT)/\xi\pm 1]^{-1}$ where $g$ is spin degeneracy factor, $\xi$ is the fugacity, $u_{\mu}$ is the four-velocity of a fluid element in the fireball of the emission source and $T$ is the corresponding local temperature. In addition, we note $p^{\mu}u_{\mu}=m_{T}\cosh\rho\cosh(\eta-y)-p_{T}\sinh\rho\cos(\phi_{p}-\phi_{s})$ gives the energy in local rest frame of the fluid, and $p^{\mu}\text{d}^{3}\sigma_{\mu}=\tau m_{T}\cosh(\eta-y)\text{d}\eta r\text{d}r\text{d}\phi_{s}$. The symbol $\rho$ represents the transverse flow rapidity distribution of the fluid element in the fireball with a transverse radius $R$.

Assuming the temperature of fireball is much smaller than the mass of constituent particles, we can then use the Boltzmann approximation for $f(x,p)$, and Eq. (1) can then be analytically obtained as Ret04

[TABLE]

where we have $\beta=m_{T}\cosh[\rho(r)]/T$ and $\alpha=p_{T}\sinh[\rho(r)]/T$; $I_{0}(x)$ and $K_{1}(x)$ are the first and second kind of modified Bessel functions, respectively. The radial expansion of the fireball would lead to a blue-shift on the transverse spectrum of the emitted particles, and thus effectively increase the temperature. As a result, we can approximately take $T$ as an effective temperature and set radial flow rapidity $\rho$ to be zero, i.e., $\alpha=0$ and $I_{0}(0)=1$, and in this case the above formula can be further simplified to be

[TABLE]

where we denote $V=\pi R^{2}\tau$ as an effective volume. Furthermore, the multiplicity of the constituent particles can be integrated out to be

[TABLE]

In the case of $m\gg T$, according to the asymptotic behavior of Bessel function $K_{\nu}(x)$ in large $x$ limit, i.e.,

[TABLE]

we can keep only the first term to make a non-relativistic approximation, and then Eq. (5) can be written as

[TABLE]

where $V_{p}=(2\pi Tm)^{\frac{3}{2}}$ reflects the effective volume of the emission source in momentum space. If we denote $y_{L}=y_{\text{max}}-y_{\text{min}}$, $V^{\prime}=y_{L}V$ and $\xi^{\prime}=\xi e^{-\frac{m}{T}}$, then the number of the emitted constituent particle in rapidity region [$y_{\text{min}},~{}y_{\text{max}}$] is given by

[TABLE]

II.2 Analytical coalescence formula

We consider the case that $N$ constituent particles are coalesced into a cluster, and the total multiplicity of the cluster can be expressed as

[TABLE]

where $\rho_{c}^{W}(x_{1},...,x_{N};p_{1},...,p_{N})$ is the Wigner density function which gives the coalescence probability, and $g_{c}$ is the coalescence factor Sat81 . Eq. (8) can be used to normalize the coalescence probability and so Eq. (9) can be rewritten as

[TABLE]

where $N_{i}=\frac{g_{i}}{(2\pi)^{3}}\xi_{i}^{\prime}V^{\prime}V_{p,i}$ is the multiplicity of the $i$-th constituent particle. The Wigner function in the above formula cannot be analytically integrated out and the result is usually obtained by a multi-dimension numerical integration method (see, e.g., Ref. KJSUN15 ).

Following Refs. Cho11 ; Cho11prc , we consider the produced clusters are in midrapidity region where one has $\tau\simeq t$, and the volume element on the hypersurface $p_{i}^{\mu}d\sigma_{i\mu}\text{d}^{3}p_{i}/E_{i}$ can then be approximated to be $\text{d}^{3}x_{i}\text{d}^{3}p_{i}$. Furthermore, the non-relativistic Boltzmann approximation for the distribution function $f(x,p)$ is adopted in the following derivation, and thus one has $V_{p,i}=(2\pi Tm_{i})^{\frac{3}{2}}$. We will discuss the corrections of relativistic effects later. Therefore, Eq. (II.2) can be approximated as

[TABLE]

Since the Wigner function does not depend on the center-of-mass coordinate of the cluster, a Jacobi transformation Mat97 ; ChenLW03 ; ChenLW04 ; LWChen06 ; ChenLW07 can be performed to separate the center-of-mass coordinate and the relative coordinates of constituent particles of the cluster, i.e.,

[TABLE]

where $\mathbf{x}_{j}$ ($\mathbf{p}_{j}$) is the spatial (momentum) coordinate of the $j$-th constituent particle and $\hat{J}$ is Jacobi matrix whose property can be found in Appendix D. In particular, the center-of-mass position vector of the cluster $\mathbf{R}$ and the relative spatial coordinate vectors $\mathbf{q}_{i}$ can be expressed as

[TABLE]

Correspondingly, in the momentum space, $\mathbf{P}$ is the total momentum of the cluster and $\mathbf{k}_{i}$ are the relative momentum vectors. It should be noted that the Jacobi matrix $\hat{J}$ we use here is different from that in Ref. Cho11prc , but the physical results should be independent of the choice of coordinate transformation. The modulus of determinant of the Jacobi matrix is $|\hat{J}|=N^{-\frac{1}{2}}$, and one then has the following identity

[TABLE]

From the above identity, Eq. (II.2) can be simplified as

[TABLE]

where

[TABLE]

is the reduced mass related to the relative coordinates (see Appendix D for more details); $\mu_{0}=\sum_{i=1}^{N}m_{i}$ is the total mass of constituent particles inside the cluster, which is equal to the rest mass of the cluster if the binding energy of the cluster is neglected; $V^{\prime}$ and $(2\pi\mu_{0}T)^{\frac{3}{2}}$ are the spatial and momentum effective volumes, respectively, for the center-of-mass motion of the cluster. To obtain Eq. (II.2), following Refs. ChenLW07 ; Cho11 ; Cho11prc , we have assumed $\eta_{i}\approx y_{i}$, and this can be justified from the fact that the coalescence probability is highly suppressed when the relative position or the relative momentum of the constituent particles inside the cluster is large. One sees that the integration with respect to momentum in Eq. (II.2) is three-dimensional, which is different from Refs. ChenLW07 ; Cho11 ; Cho11prc where the momentum integration is assumed to be two-dimensional in transverse plane.

The denominator in Eq. (II.2) can be re-expressed as

[TABLE]

where the factor $N^{-\frac{3}{2}}$ comes from the coordinate transformation and will be canceled out at last. It should be noted that to obtain Eq. (43), Eq. (79) has been used.

We note that the contributions of center-of-mass coordinate in the numerator and denominator (i.e., Eq. (43)) of Eq. (II.2) cancel out, and Eq. (II.2) can be written as

[TABLE]

Furthermore, for simplicity, we assume that the Wigner function of the cluster can be factorized into the Wigner functions of each relative coordinate, i.e., $\rho_{c}^{W}=\prod_{i=1}^{N-1}\rho_{c,i}^{W}$, then Eq. (44) can be recast into

[TABLE]

If we adopt the harmonic oscillator assumption for the potential of constituent particles (see Appendix A), and consider a mixed ensemble with a definite orbital angular momentum state ($l$) in the lowest energy state with $n=l$, then the spatial integration of Wigner function can be obtained as

[TABLE]

where $(2\pi)^{3}$ comes from the normalization of Wigner function $\int\rho_{c,i}^{W}\text{d}^{3}q_{i}\text{d}^{3}k_{i}=(2\pi\hbar)^{3}$ with the reduced Plank constant $\hbar$ taken to be unity, and the ensemble averaged probability distribution is represented by $P(k_{i})=\frac{(4\pi\sigma_{i}^{2})^{\frac{3}{2}}}{(2\pi)^{3}}\frac{(2\sigma_{i}^{2}k_{i}^{2})^{l}}{(2l+1)!!}e^{-\sigma_{i}^{2}k_{i}^{2}}$ with $\sigma_{i}^{2}=\frac{1}{\mu_{i}w}$ where $w$ is the frequency of the harmonic oscillator (see Appendix A) and it can be determined by the RMS radius ($r_{\text{rms}}$) of the cluster (see Eq. (87) in Appendix D). After the momentum integration, the multiplicity of the cluster can be obtained as

[TABLE]

where $F(\sigma,\mu,l,T)$ is expressed as (see Appendix B)

[TABLE]

with the function $G(l,x)$ defined as a ratio of two hypergeometric function (see Appendix C).

Finally. we obtain the following analytical coalescence formula, i.e.,

[TABLE]

Since $V^{\prime}=Vy_{L}$, we can express the multiplicity per unit rapidity as

[TABLE]

II.3 Relativistic corrections

In the derivation of Eqs. (II.2) and (II.2), a non-relativistic approximation of the momentum distribution has been adopted. However, it should be noted that the relativistic effect could decrease the coalescence probability as shown in the study for the exotic $\Theta^{+}$ production in relativistic heavy-ion collisions ChenLW04 . It is thus interesting and important to include this effect in the analytical coalescence formula. We use $g_{\text{rel}}$ to denote the relativistic correction factor.

For a cluster, it is reasonable to assume that the center-of-mass motion should obey relativistically covariant one-particle momentum distribution, while the relative motions of constituent particles inside the cluster should be approximately non-relativistic because the relative momentum and spatial coordinates have to be small, otherwise the coalescence probability of forming a cluster is highly suppressed. This means that the relativistic correction of the Wigner function integration should be small, and thus the leading order relativistic correction should mainly come from the one-particle momentum distribution. From Eqs. (5) and (6) and keeping the leading order correction in the latter, one obtains the effective volume in momentum space as $V_{p}=(2\pi Tm)^{\frac{3}{2}}(1+\frac{15T}{8m})$. Therefore, the relativistic correction factor $g_{\text{rel}}$ can be expressed as

[TABLE]

On one hand, it is easy to prove that $g_{\text{rel}}$ is always less than one and this leads to a general conclusion that the relativistic effect tends to decrease the cluster yield. On the other hand, one can see that when $\frac{m_{i}}{T}$ is much larger than one, $g_{\text{rel}}$ approaches to unity and the relativistic correction can then be neglected.

II.4 Cluster finite size effects

In evaluating the spatial part of the integration in Eq. (45), we have assumed that the radius $R$ of the fireball is much larger than the size of the produced cluster. However, for some loosely bound clusters such as ${}^{3}_{\Lambda}$H, its $r_{\text{rms}}$ can be as large as $4.9$ fm Nem00 . It has been shown in Ref. Sch99 that the large value of $\frac{r_{\text{rms}}}{R}$ will decrease the value of integration. In this work, we denote the correction factor of this cluster finite size effect by $g_{\text{size}}$. When $\frac{r_{\text{rms}}}{R}$ approaches to [math], $g_{\text{size}}$ should approach to $1$.

The cluster size correction factor $g_{\text{size}}$ cannot be analytically obtained, but it can be calculated by a multi-dimensional numerical integration KJSUN15 . Based on the multi-dimensional numerical integration, we propose the following empirical expression for $g_{\text{size}}$ for $N$-particle ($2\leq~{}N\leq~{}6$) coalescence, i.e.,

[TABLE]

where $x=\frac{r_{\text{rms}}}{R}$ is assumed to be in range [0, 0.5]. For the same $x$, we note $g_{\text{size}}$ is larger for 3-particle coalescence than for 2-particle coalescence, but the difference is relatively small compared with their own values. This empirical expression can be applied to central heavy-ion collisions at both RHIC and LHC energies. For coalescence of constituent particles with different masses, we notice that the mass difference tends to decrease the value of $g_{\text{size}}$, but the correction is not significant.

III Discussions

Combining the relativistic correction factor $g_{\text{rel}}$ and the cluster size correction factor $g_{\text{size}}$ with Eq. (II.2), we obtain the following full analytical coalescence formula for the cluster multiplicity per unity rapidity, i.e.,

[TABLE]

We denote the above formula as COAL-SH. For comparison, the COAL-Ex formula (i.e., Eq. (18) of Ref. Cho11prc ) reads

[TABLE]

One can see that compared with the formula COAL-Ex, besides the additional factors $g_{\text{rel}}$ and $g_{\text{size}}$, the new formula COAL-SH displays very similar structure but different dependence on the $w$ parameter and the orbital angular momentum quantum numbers, which will be discussed in the following.

III.1 Suppression from orbital angular momentum

In the COAL-SH formula (Eq. (III)), the suppression factor from non-zero orbital angular momentum quantum number $l$ is

[TABLE]

with $x=(\frac{2T}{w})^{\frac{1}{2}}$. From the property of $G(l,x)$ (see Appendix C), one can easily show

[TABLE]

indicating the suppression feature due to the finite orbital angular momentum quantum number $l$.

For $x\gg 1$, one can show that $\mathcal{S}(l,x)$ is very close to $1$ for all $l$. Taking $l=1$ as an illustration, one has

[TABLE]

for $x\gg 1$. This feature is consistent with the results in Ref. Mul06 where the orbital angular momentum suppression factor is shown to approach to unity in high temperature limit (large $T$). This is a little bit different from the result of the COAL-Ex formula (Eq. (54)) from which one can see the value of the suppression factor is $\frac{2}{3}$ for $l=1$ at the limit of $x\gg 1$.

III.2 Saturation with the cluster size

According to Eq. (87) in Appendix D, the value of harmonic oscillator frequency $w$ is inversely proportional to the $r^{2}_{\text{rms}}$ value of the produced cluster. Under the limit of $\frac{2T}{w}\gg 1$ for large cluster size (but the $r_{\text{rms}}$ of the cluster is still assumed to be smaller than the fireball radius $R$ ), the COAL-SH formula Eq. (III) can be further simplified to be

[TABLE]

The above expression indicates that the $\frac{dN_{c}}{dy}$ will be saturated to a constant at the limit of large (small) $r_{\text{rms}}$ ($w$). In particular, for $R\gg r_{\text{rms}}$, the cluster size correction factor $g_{\text{size}}$ approaches to unity and Eq. (58) becomes

[TABLE]

which means that the size of the produced cluster has no influence on its yields in the limits $2T\gg w$ and $R\gg r_{\text{rms}}$.

For the COAL-Ex formula, in the limit of $\frac{2T}{w}\gg 1$ for large cluster size, the $\frac{dN_{c}}{dy}$ does not saturate with the increment (decrease) of $r_{\text{rms}}$ ($w$). For example, in the case of two-particle coalescence, the $\frac{dN_{c}}{dy}$ is proportional to $r_{\text{rms}}$ Cho11 ; Cho11prc . As we will see in the following, this difference can lead to quite different predictions for the yield ratio ${}^{3}_{\Lambda}\text{H}/^{3}\text{He}$ between the COAL-Ex and COAL-SH formulae since the hypertriton ${}^{3}_{\Lambda}\text{H}$ is a very loosely bound state with a huge radius of about $4.9$ fm Nem00 .

IV Applications

In this section, we first apply COAL-SH and COAL-Ex to investigate the strangeness population factor $S_{3}=^{3}_{\Lambda}$H/(3He$\times\Lambda/$p), and then to study the production of some exotic hadrons, in central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV. The results from thermal (statistical) model are also included for comparison.

IV.1 Strangeness population factor $S_{3}$

The strangeness population factor $S_{3}$ was first suggested in Ref. Arm04 and it is expected to be a good representation of the local correlation between baryon number and strangeness Zha10 ; Koc05 . Therefore, it may provide valuable information of deconfinement in relativistic heavy-ion collisions. The $S_{3}$ was measured to be $0.60\pm 0.13(\text{stat.})\pm 0.21(\text{syst.})$ for central ([math]-$10\%$ centrality) Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV ALICE15H . It should be noted that while there is negligible feed-down from heavier states into ${}^{3}_{\Lambda}$H and 3He, the $\Lambda$ and p are significantly influenced by feed-down from decays of excited baryonic states. When we calculate $S_{3}$ within the coalescence model, the contributions to the yields of $\Lambda$ and $p$ from electromagnetic and weak decays are excluded.

Within the formula COAL-SH, $S_{3}$ can be expressed as

[TABLE]

where $w_{{}^{3}\text{He}}=0.013$ GeV and $w_{{}^{3}_{\Lambda}\text{H}}=1.6\times 10^{-3}$ GeV are the corresponding harmonic oscillator frequencies of ${}^{3}\text{He}$ and ${}^{3}_{\Lambda}\text{H}$, respectively, and their values are obtained from the $r_{\text{rms}}$ of ${}^{3}\text{He}$ (1.76 fm Rop09 ) and ${}^{3}_{\Lambda}\text{H}$ (4.9 fm Nem00 ) through Eq. (87). The relativistic correction of $S_{3}$ due to $g_{\text{rel}}$ (Eq. (51)) can be simply neglected since one has $g_{\text{rel}}(^{3}\text{He})\approx g_{\text{rel}}(^{3}_{\Lambda}\text{H})$. Because ${}^{3}_{\Lambda}$H has a much larger $r_{\text{rms}}$ than 3He, one has to take account of the cluster size correction factor $g_{\text{size}}$. By neglecting the relativistic correction, $S_{3}$ can be expressed as

[TABLE]

Now, the difference for the prediction of $S_{3}$ between the thermal model and the coalescence model becomes clear: the factor $0.845$ in Eq. (61) corresponds to the prediction of thermal model at LHC energy And11 , and the remaining parts in Eq. (61) thus correspond to the corrections from coalescence model.

In general, the freeze-out temperature $T$ lies between $0.1$ to $0.2$ GeV in Pb+Pb collisions that we are considering here, then $0.845(1+\frac{0.013}{2T})^{2}$ is in the range of $0.9\sim 0.95$, and one can see the temperature dependence is quite weak for this part. Compared with the thermal model, the largest correction in the coalescence model comes from the size effect due to different sizes of ${}^{3}_{\Lambda}$H and 3He. The radius $R$ of the fireball at freeze-out in central Pb+Pb collisions at $\sqrt{s_{NN}}$=2.76 TeV is about $19.7$ fm KJSUN16 , and according to Eq. (52), one has $\frac{g_{\text{size}}(^{3}_{\Lambda}\text{H})}{g_{\text{size}}(^{3}\text{He})}=0.60$, which leads to $S_{3}=0.54\sim 0.57$, consistent with the measured value $0.60\pm 0.13(\text{stat.})\pm 0.21(\text{syst.})$. To fit the measured central value $0.60$, one can introduce a multi-freeze-out between nucleon and $\Lambda$ with an earlier $\Lambda$ freeze-out KJSUN16 .

For the formula COAL-Ex, $S_{3}$ can be obtained as

[TABLE]

One can see that compared with the expression (IV.1) based on the COAL-SH formula, the expression (62) from the COAL-Ex formula does not have correction factors $g_{\text{rel}}$ and $g_{\text{size}}$ but includes an additional factor $w_{{}^{3}\text{He}}/w_{{}^{3}_{\Lambda}\text{H}}=8.1$. As a result, $S_{3}$ based on COAL-Ex is larger than about $7$ and thus the COAL-Ex formula significantly overestimates the observed $S_{3}$ factor value around $0.6$ for central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV.

IV.2 Production of exotic hadrons from quark coalescence

In the following, we focus on the production of exotic hadrons from quark coalescence in central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV. Three kinds of exotic hadrons with four quark flavors (i.e., $u$, $d$, $s$, and $c$) are considered, namely, the exotic mesons that could be in four-quark states, the exotic baryons that could be in five-quark states, and the exotic dibaryons that could be in six-quark or eight states. The detailed properties, including mass, isospin, spin, parity, orbital angular momentum quantum number, quark configuration and decay modes, of these exotic hadrons can be found in Table IV of Ref. Cho11prc . For the calculations with the COAL-SH formula, the cluster size correction factor $g_{\text{size}}$ is neglected since the RMS radii of the exotic hadrons we are considering here are small (less than about $2$ fm) based on the harmonic oscillator frequencies determined in the following.

For the quark coalescence model calculations, following Refs. Cho11 ; Cho11prc , the masses of $u$($d$), $s$ and $c$ constituent quarks are taken to be $m_{u,d}=0.3$ GeV, $m_{s}=0.5$ GeV and $m_{c}=1.5$ GeV. The temperature is set to be $T=0.154$ GeV Baz12 ; Baz14 , and we note that a small (e.g., $20\%$) variance of the temperature value does not change our conclusion. In order to apply the formulae COAL-SH and COAL-Ex, one also needs the information on the harmonic oscillator frequencies ($w$ for hadrons with only $u$($d$) quarks, $w_{s}$ for hadrons with strange hadrons, and $w_{c}$ for charmed hadrons), the fireball volume ($V$), the number of constituent quarks ($N_{u(d)}$ for $u(d)$ quarks, $N_{s}$ for $s$ quarks and $N_{c}$ for $c$ quarks).

For COAL-SH, the values of frequencies $w$, $w_{s}$ and $w_{c}$ can be deduced from the hadron size through Eq. (87). In particular, the values of $w$ and $w_{s}$ have been determined to be $w=0.184$ GeV and $w_{s}=0.078$ GeV in Ref. KJSUN16-1 . For the value of $w_{c}$, we obtain $w_{c}=0.087$ GeV by using the value $0.546$ fm for the RMS radius of $\Omega_{ccc}$ Zhu14 . With $w=0.184$ GeV, $w_{s}=0.078$ GeV and $w_{c}=0.087$ GeV, we obtain reasonable RMS radius values for various hadrons, namely, $0.84$ fm for protons ($p$), $0.87$ fm for $\phi$ mesons, $1.1$ fm for $\Xi^{-}$, $1.0$ fm for $\Omega^{-}$, $1.2$ fm for $\Lambda$, and $0.98$ fm for $D^{0}$. Finally, the parameters $N_{u(d)}$, $N_{s}$ and $N_{c}$ can be obtained from fitting the measured yields ($dN/dy$ at midrapidity) of $p$ Abe12 , $\Lambda$ Abe13 , $\phi$ Abe15 , $\Xi^{-}$ Abe14 , $\Omega^{-}$ Abe14 , and $D^{0}$ WangRQ16 .

For the yields of $p$, $\Lambda$, $\phi$, $\Xi^{-}$, $\Omega^{-}$ and $D^{0}$ in central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, it should be noted that the weak decays have already been corrected in the experimental data, but not for the strong decays and electromagnetic decays. In order to compare with experiment results, we have to include the contributions from strong and electromagnetic decays in quark coalescence model. Following Refs. Cho11 ; Cho11prc ; KJSUN16-1 , we assume the relations $N_{\Lambda(1115)}^{\text{measured}}=N_{\Lambda(1115)}+\frac{1}{3}N_{\Sigma(1192)}+(0.87+\frac{0.11}{3})N_{\Sigma(1385)}=7.44N_{\Lambda(1115)}$, $N_{p}^{\text{measured}}=N_{p}+N_{\Delta^{++}(1232)}+\frac{1}{2}N_{\Delta^{+}(1232)}+\frac{1}{2}N_{\Delta^{0}(1232)}=5N_{p}$, $N_{\Xi^{-}}^{\text{measured}}=N_{\Xi^{-}}+\frac{1}{2}N_{\Xi(1530)}=3N_{\Xi^{-}}$, and $N_{D^{0}}^{\text{measured}}=N_{D^{0}}+N_{D^{*0}}+0.677N_{D^{+*}}=6.0N_{D^{0}}$. For $\phi$ and $\Omega^{-}$, we assume no strong and electromagnetic decay corrections, and thus $N_{\phi}^{\text{measured}}=N_{\phi}$ and $N_{\Omega^{-}}^{\text{measured}}=N_{\Omega^{-}}$. Here $N_{\Lambda(1115)}$, $N_{p}$, $N_{\Xi^{-}}$, $N_{\phi}$, $N_{\Omega^{-}}$ and $N_{D^{0}}$ represent the corresponding hadron multiplicity obtained directly from the quark coalescence model.

For COAL-SH, by fitting the measured yields of $p$, $\Lambda$, $\phi$, $\Xi^{-}$, $\Omega^{-}$ and $D^{0}$, we obtain the fireball volume $V=1.38\times 10^{4}$ fm3, the constituent quark numbers $N_{u(d)}=3\times 462$, $N_{s}=3\times 196$ and $N_{c}=3\times 27$.

For the formula COAL-Ex, we find it cannot reasonably fit the measured yields of $p$, $\Lambda$, $\phi$, $\Xi^{-}$, $\Omega^{-}$ and $D^{0}$ if the frequency parameters $w$, $w_{s}$ and $w_{c}$ are taken to have the same values as those in the COAL-SH formula. In Refs. Cho11 ; Cho11prc , for central Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV and central Pb+Pb collisions at $\sqrt{s_{NN}}=5$ TeV, these frequencies were assumed to be free parameters and determined to be $w=0.55$ GeV, $w_{s}=0.519$ GeV and $w_{c}=0.385$ GeV by reproducing the yields of $w(782)$, $\rho(770)$, $\Lambda(1115)$ and $\Lambda_{c}(2286)$ predicted from the thermal model. If we take these frequency values for the formula COAL-Ex, we then obtain the fireball volume $V=6.06\times 10^{3}$ fm3, the constituent quark numbers $N_{u(d)}=3\times 267$, $N_{s}=3\times 156$ and $N_{c}=3\times 20$.

For the thermal model (see Refs. Cho11 ; Cho11prc for details), the baryon and strange chemical potentials are set to be zero. And by fitting the measured yields of $p$, $\Lambda$, $\phi$, $\Xi^{-}$, $\Omega^{-}$ and $D^{0}$, we obtain the temperature $T=0.153$ GeV, the fireball volume $V=6.45\times 10^{3}$ fm3, and the charm fugacity $\gamma_{c}=61.6$.

Table 1 summarizes the model parameters for COAL-SH, COAL-Ex and the thermal model. Table 2 displays the predictions from COAL-SH, COAL-Ex and the thermal model, for the yields of $p$, $\Lambda$, $\phi$, $\Xi^{-}$, $\Omega^{-}$ and $D^{0}$ in central Pb+Pb collisions at $\sqrt{s}=2.76$ TeV. Also included in Table 2 are the corresponding experimental data and the coalescence factors $g_{c}$ due to spin and color degrees of freedom. From Table 2, one can see COAL-SH gives a best fit, COAL-Ex also gives a nice fit, while the thermal model does not fit the $\Lambda$ yield well. We note that COAL-SH can give a similar good description for the yields of $p$, $\Lambda$, $\phi$, $\Xi^{-}$ and $\Omega^{-}$ as the more realistic multi-dimensional numerical integration method, and the latter has been shown to describe very well the spectra and yields of these hadrons KJSUN16-1 . This implies that the COAL-SH formula can be served as a very useful tool to evaluate the yield of clusters produced in relativistic heavy-ion collisions.

Shown in Tab. 3 are the predicted exotic hadron yields from COAL-SH, COAL-Ex and the thermal model. The properties of these exotic hadrons can be found in Table IV of Refs. Cho11prc . In Tab. 3, most of the predicted yields from COAL-SH and COAL-Ex are very close and the differences are within about factor two. However, there are also some cases that the difference between the COAL-SH and COAL-Ex predictions is larger than factor two. This is mainly due to different suppressions of non-zero orbital angular momentum states of the exotic hadrons within these two formulae. Compared with the predictions from thermal model, most of the results of COAL-SH are much smaller by about one order, and even smaller by about two orders for some cases. However, when the masses of the exotic hadrons, such as $a_{0}(980)$, $X(3872)$, $Z^{+}(4430)$ and $\Lambda(1450)$, are much larger than the total mass of their constituent quarks, their yields from COAL-SH are larger than those from the thermal model.

In particular, it is suggested that the mesons $f_{0}(980)$, $a_{0}(980)$, $D_{s}(2371)$ and $X(3872)$ could be either in normal two-quark states or in exotic four-quark states. Similarly, the baryon $\Lambda(1405)$ could be either in normal three-quark state (i.e., $uds(L=1)$) or in exotic five-quark state (i.e., $qqqs{\bar{q}}$). It is interesting to see that the predicted yields of these hadrons in exotic four(five)-quark states are smaller by about two orders than those in normal two(three)-quark states, indicating that measuring the yields of these hadrons in central Pb+Pb collisions at $\sqrt{s}=2.76$ TeV will be potentially useful to identify the internal quark configuration of these hadrons.

Furthermore, the exotic dibaryon $\bar{K}NN$ could be either in six-quark state (i.e., $qqqqqs(L=1)$) or in eight-quark state (i.e., $qqqqqqs{\bar{q}}$) with the latter (eight-quark state) having a yield smaller than that of the former (six-quark state) by two orders. In particular, based on the COAL-SH calculations, the yield ($dN/dy$ at midrapidity) of the exotic dibaryon $\bar{K}NN$ is $3.91\times 10^{-3}$ for six-quark state and $4.49\times 10^{-5}$ for eight-quark state. These yields are measurable in central Pb+Pb collisions at $\sqrt{s}=2.76$ TeV at LHC via the strong decay $\bar{K}NN\rightarrow\Lambda N$.

The above results based on COAL-SH and COAL-Ex demonstrate that the yields of various exotic hadrons sensitively depend on the quark configuration and therefore the yield measurement of these hadrons is potentially useful to determine the quark configuration of these exotic hadrons. It should be noted that these exotic hadrons could be in molecular states and the formulae COAL-SH and COAL-Ex can also be used to evaluate the molecular state yield if the freeze-out information is known for nucleons, $\Xi$, $\Omega$, $\Xi_{c}$, $K$, $\bar{K}$, $D$, $\bar{D}$, $D^{*}$, $\bar{D}^{*}$ and $D_{1}$, and this is beyond the scope of the present work and may be pursued in future.

V conclusion

Based on a blast-wave-like parametrization for phase-space configuration of constituent particles at freeze-out in relativistic heavy-ion collisions, we have derived a new approximate analytical formula COAL-SH for cluster yield within the covariant coalescence model. The new analytical coalescence formula COAL-SH improves some aspects of the existing formulae by treating the integration of Wigner function in the many-body coalescence process in full phase-space, and thus has a good saturation property when the $r_{\text{rms}}$ of the produced clusters increases. The corrections of relativistic effects and the finite size effects of the produced cluster relative to the emission source are also considered in the COAL-SH formula.

We have applied COAL-SH to investigate the strangeness population factor $S_{3}=^{3}_{\Lambda}$H/(3He$\times\Lambda/$p) from nucleon/$\Lambda$ coalescence and the exotic hadron production from quark coalescence in central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV. The results have been compared with the prediction from the analytical coalescence formula COAL-Ex derived by the ExHIC collaboration and the thermal model. Our results indicate that the COAL-SH formula can reasonably describe the measured $S_{3}$ factor while the COAL-Ex formula predicts a too large value of the $S_{3}$ factor. It is also indicated that the cluster finite size effect is important for the $S_{3}$ evaluation due to the large size of the hypernucleus ${}^{3}_{\Lambda}$H.

For the exotic hadron production in central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, we have determined the model parameters by fitting the yields of $p$, $\Lambda$, $\phi$, $\Xi^{-}$, $\Omega^{-}$ and $D^{0}$. We have found that COAL-SH can nicely fit the experiment data of these yields but COAL-Ex cannot reasonably fit the data if the harmonic oscillator frequencies are determined from the RMS radii of normal hadrons. By adjusting the values of the harmonic oscillator frequencies, the COAL-Ex can reasonably fit the measured yields. The thermal model can also fit the measured yields except that the $\Lambda$ yield cannot be well reproduced.

For the yields of exotic hadrons, most of the predictions from COAL-SH and COAL-Ex (with adjusted harmonic oscillator frequencies) are very close to each other within about factor two. COAL-SH and COAL-Ex may give significantly different predictions for the exotic hadrons with non-zero angular momentum states as the two formulae have different suppressions of non-zero angular momentum states. Compared with the predictions from thermal model, most of the results of COAL-SH are about one order smaller or even more except for $a_{0}(980)$, $X(3872)$, $Z^{+}(4430)$ and $\Lambda(1450)$, for which the hadron masses are much larger than the total mass of their constituent quarks.

The new formula COAL-SH can be served as a useful tool to study particle production in relativistic heavy-ion collisions. It should be pointed out that for the COAL-SH formula, we have assumed the size of the emission source for constituent particles at freeze-out should be larger than that of the produced cluster, and thus it cannot be simply applied to some large-size cluster production in small collision system like $pp$ collision. This interesting issue will be explored in the near future.

Acknowledgements.

We are grateful to Che Ming Ko for helpful discussions. This work was supported in part by the Major State Basic Research Development Program (973 Program) in China under Contract Nos. 2015CB856904 and 2013CB834405, the National Natural Science Foundation of China under Grant Nos. 11625521, 11275125 and 11135011, the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning, Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education, China, and the Science and Technology Commission of Shanghai Municipality (11DZ2260700).

Appendix A Quantum 3-dimension isotropic harmonic oscillators

We consider a 3-dimension isotropic quantum harmonic oscillator with the potential of $V(r)=\frac{1}{2}\mu\omega^{2}r^{2}$, where $\mu$ ($\omega$) is the mass (frequency) of the oscillator. The wave function can be obtained by solving Sch$\ddot{o}$dinger equation and the solution reads

[TABLE]

where $n=2s+l$ is the principal quantum number related to the energy $E=\hbar\omega(n+\frac{3}{2})$, $l$ is the orbital angular momentum quantum number, $m=-l,...,l$ is the magnetic quantum number, and $s$ is a non-negative integer representing radial excitation. $N_{sl}=\sqrt{\sqrt{\frac{1}{\pi}}\frac{2^{s+l+2}s!\sigma^{-2l-3}}{(2s+2l+1)!!}}$ is a normalized constant with $\sigma^{2}=\frac{1}{\mu w}$, and $L_{s}^{l+\frac{1}{2}}(x)$ is the generalized Laguerre polynomials with order of $s$. $Y_{lm}(\theta,\phi)$ is the spherical harmonic function which is normalized as $\int|Y_{lm}|^{2}d\Omega=1$. The RMS radius $r_{\text{rms}}$ can then be obtained via the following expression

[TABLE]

The wave function in momentum ($k$) representation reads

[TABLE]

For an ensemble having the lowest energy with a given $l$, the averaged probability distribution in momentum space can be obtained as

[TABLE]

and $P(k)$ is normalized as $\int P(k)\text{d}^{3}k=1$. To obtain Eq. (66), the following Uns$\ddot{o}$ld’s theorem has been used

[TABLE]

Appendix B Function $F(\sigma,\mu,l,T)$

The function $F(\sigma,\mu,l,T)$ is defined as

[TABLE]

where $k_{T}=k\sin\theta$ is the transverse momentum. With Eq. (66), the $k$-dependent part can be integrated out as follows

[TABLE]

where the $\Gamma$-function identity $\Gamma(n+\frac{1}{2})=\frac{\sqrt{\pi}(2n)!}{4^{n}n!}$ has been used. With the function $G(l,x)$ (see Appendix C), $F$ can finally be obtained as

[TABLE]

Appendix C Orbital angular momentum factor

The orbital angular momentum factor $G(l,x)$ is an integration which can be done analytically as follows

[TABLE]

where ${}_{2}F_{1}$ is the hypergeometric function. $G(l,x)$ can also be expanded as

[TABLE]

From the series expansion, the following properties of $G(l,x)$ can be easily proved, i.e.,

[TABLE]

For $x^{2}\gg 1$, $G(l,x)$ is thus very close to unity.

Appendix D Jacobi transformation and the RMS radius of clusters

The Jacobi matrix $\hat{J}$ for the coordinate transformation, defined in Eqs. (36), (38) and (39), has some special properties. The Jacobi matrix $\hat{J}$ is introduced to separate the center-of-mass and the relative coordinates of a many-body system since the Wigner function does not depend on the center-of-mass coordinate and it can be expressed in terms of the relative coordinates.

Considering $N$ independent constituent particles with mass $m_{i}(i=1,2,3,...,N)$ in ($3N$-dimension) harmonic oscillators with same frequency $w$, the total potential energy of the system can be expressed as

[TABLE]

where $\hat{M}$ is the mass matrix which is diagonal and $\mathbf{r}^{\prime}$ is the new coordinate vector. The elements of transformation matrix $\hat{J}^{-1}$ is

[TABLE]

The new mass matrix $\hat{M}^{\prime}$ is also diagonal and the reduced mass can be expressed as

[TABLE]

from which one can see that $\mu_{0}$ is the total mass and $\mu_{i}$ is the reduced mass related to the constituent particles in the corresponding relative coordinates, i.e.,

[TABLE]

From the above equations, the following identity can then be obtained

[TABLE]

After Jacobi transformation, the $N$ independent constituent particles with mass $m_{i}$ $(i=1,2,3,...,N)$ in harmonic oscillators with frequency $w$ are transferred into $N$ independent particles with mass $\mu_{0}$ (total mass of the cluster) and $\mu_{j}$ $(j=1,2,3,...,N-1)$ (the reduced mass related to the $N-1$ relative coordinates) in harmonic oscillators with frequency $w$. At the same time, the coordinates $\mathbf{r}_{i}$ $(i=1,2,3,...,N)$ are correspondingly transferred into the center-of-mass coordinate $\mathbf{R}$ and the $N-1$ relative coordinates $\mathbf{q}_{j}$ $(j=0,1,2,3,...,N-1)$ as in Eq. (36). The total potential energy of the system can be re-expressed as

[TABLE]

The above expression is the basis of deriving the Wigner function in Section II.2. The wave functions for harmonic oscillators with mass $\mu_{j}$ $(j=1,2,3,...,N-1)$ and frequency $w$ can be found in Appendix A.

Furthermore, the mean-square radius $r_{\text{rms}}^{2}$ can be obtain as

[TABLE]

where we define $\hat{S}=\hat{L}^{T}\hat{L}$ with the matrix $\hat{L}$ defined via the following relations, i.e.,

[TABLE]

With Eq. (76) and Eq. (D), one can obtain the elements of matrix $\hat{S}$, i.e.,

[TABLE]

and then one has

[TABLE]

Furthermore, from Eq. (64), one has

[TABLE]

where $\sigma_{i}^{2}=1/(\omega\mu_{i})$, and $l_{i}$ and $s_{i}$ are, respectively, the quantum numbers of orbital angular momentum and radial excitation for the $i$-th relative coordinate.

In particular, for $l=s=0$, we have

[TABLE]

From Eqs. (78), (83), (84) and (86), one can then obtain the analytical formula for mean-square radius of the cluster as

[TABLE]