Exotic Hadrons from Heavy Ion Collisions

TL;DR

Heavy ion collisions at high energies enable the production and study of exotic hadrons, offering insights into their structures and QCD properties through experimental measurements and theoretical models.

Contribution

This review summarizes current understanding and models of exotic hadron production in heavy ion collisions, highlighting yield differences and detection prospects at RHIC and LHC.

Findings

Compact multiquark states have lower yields than molecular states.

Heavy ion collisions can produce exotic hadrons with significant charm and bottom quarks.

Resonance production estimates are influenced by the resonance-involving coalescence model.

Abstract

Heavy ion collisions (HIC) at high energies are excellent ways for producing heavy hadrons and composite particles. With upgraded detectors at RHIC and LHC, it has become possible to measure hadrons beyond their ground states. Therefore, HIC provide a new method for studying exotic hadrons that are either hadronic molecular states or compact multiquark systems. Because their structures are related to the fundamental properties of QCD, studying exotic hadrons is currently one of the most active areas of research in hadron physics. Experiments carried out at various accelerator facilities have indicated that some exotic hadrons may have already been produced. The present review is a summary of the current understanding of a selected set of exotic particle candidates that can be potentially measured in HIC. It also includes discussions on the production of exotic hadrons in HIC based on the coalescence and statistical models. A more detailed discussion leads to the conclusion that the yield of a hadron is typically an order of magnitude smaller when it is a compact multiquark state than that of an excited hadronic state with normal quark numbers and/or a molecular configuration. Attention is also given to some of the proposed heavy exotic hadrons that could be produced with sufficient abundance in HIC because of the significant numbers of charm and bottom quarks produced at RHIC and LHC, making it possible to study them in these experiments. Further included in the discussion are the general formalism for the coalescence model that involves resonance particles and its implication on the present estimated yield for resonance production. Finally, a review is given on recent studies to constrain the hadron-hadron interaction through correlation measurements in HIC and their implications on the interpretation and the possible existence of exotic states in hadronic interactions.

Figures

Tables

Table 1. Table 2.1: Summary of the charged charmonium-like states. The mass and width are taken from PDG  [ 48 ] . Z c ​ ( 3900 ) subscript 𝑍 𝑐 3900 Z_{c}(3900) , Z c ​ ( 4020 ) subscript 𝑍 𝑐 4020 Z_{c}(4020) and Z c ​ ( 4430 ) subscript 𝑍 𝑐 4430 Z_{c}(4430) are listed in the PDG summary Table.

State	Mass [MeV]	Width [MeV]	Decay mode	$J^P$	Reference
$Z_c\left(3900\right)$	$3886.6\pm 2.4$	$28.1\pm 2.6$	$J/\psi {\pi }^{\pm }$, ${\left(D{\overline{D}}^{\ast }\right)}^{\pm }$	$1^{+}$	[153, 154, 155, 157]
$Z_c\left(4020\right)$	$4024.1\pm 1.9$	$13\pm 5$	$h_c{\pi }^{\pm },{\left(D^{\ast }{\overline{D}}^{\ast }\right)}^{\pm }$	-	[156, 161]
$Z_c\left(4050\right)$	$4051\pm {14}_{-41}^{+20}$	${82}_{-17}^{+21}{}_{-22}^{+47}$	${\chi }_{c1}{\pi }^{\pm }$	-	[151]
$Z_c\left(4055\right)$	$4054\pm 3\pm 1$	$45\pm 11\pm 6$	$\psi \left(2S\right){\pi }^{\pm }$	-	[163]
$Z_c\left(4200\right)$	${4196}_{-29}^{+31}{}_{-13}^{+17}$	$370\pm 70{}_{-132}^{+70}$	$J/\psi {\pi }^{\pm }$	$1^{+}$	[162]
$Z_c\left(4240\right)$	$4239\pm 18{}_{-10}^{+45}$	$220\pm 47{}_{-74}^{+108}$	$\psi \left(2S\right){\pi }^{\pm }$	$0^{-}$	[150]
$Z_c\left(4250\right)$	${4248}_{-29}^{+44}{}_{-35}^{+180}$	${177}_{-39}^{+54}{}_{-61}^{+316}$	${\chi }_{c1}{\pi }^{\pm }$	-	[151]
$Z_c\left(4430\right)$	${4478}_{-18}^{+15}$	$181\pm 31$	$\psi \left(2S\right){\pi }^{\pm }$	$1^{+}$	[146, 148, 149, 150]

Table 2. Table 2.2: Summary of the charged bottomonium-like states.

State	Mass [MeV]	Width [MeV]	Decay mode	$J^P$	Reference
$Z_b{\left(10610\right)}^{+}$	$10607.2\pm 2.0$	$18.4\pm 2.4$	${\pi }^{\pm }{\rm Y}\left(ns\right)$, ${\pi }^{\pm }{h}_b$	$1^{+}$	[165]
$Z_b{\left(10610\right)}^0$	$10609\pm 4\pm 4$	-	${\rm Y}(2,3s){\pi }^0$	$1^{+}$	[166]
$Z_b{\left(10650\right)}^{+}$	$10652.2\pm 1.5$	$11.5\pm 2.2$	${\pi }^{\pm }{\rm Y}\left(ns\right)$, ${\pi }^{\pm }{h}_b$	$1^{+}$	[165]

Table 3. Table 2.3: Summary of exotic hadrons with light flavors. Shown are the mass ( m 𝑚 m ), isospin ( I 𝐼 I ), spin and parity ( J P superscript 𝐽 𝑃 J^{P} ), the quark structure ( 2 ​ q / 3 / q / 6 ​ q 2 𝑞 3 𝑞 6 𝑞 2q/3/q/6q and 4 ​ q / 5 ​ q / 8 ​ q 4 𝑞 5 𝑞 8 𝑞 4q/5q/8q ), molecular configuration (Mol.) and corresponding oscillator frequency ( ω Mol . subscript 𝜔 Mol \omega_{\rm Mol.} ). For the ω Mol . subscript 𝜔 Mol \omega_{\mathrm{Mol.}} , it is fixed by the binding energies B of hadrons ( ω ≃ 6 × B similar-to-or-equals 𝜔 6 B \omega\simeq 6\times\mathrm{B} , marked (B)) or their mean square distances ⟨ r 2 ⟩ delimited-⟨⟩ superscript 𝑟 2 \langle r^{2}\rangle ( ω ≃ 3 / 2 ​ μ ​ ⟨ r 2 ⟩ similar-to-or-equals 𝜔 3 2 𝜇 delimited-⟨⟩ superscript 𝑟 2 \omega\simeq 3/2\mu\langle r^{2}\rangle , marked (R)). In the case of three-body molecular configurations for exotic dibaryons, ω Mol . subscript 𝜔 Mol \omega_{\mathrm{Mol.}} is that for the subsystem.

Particle	$m$ [MeV]	$(I,J^P)$	$q\overline{q}$/$qqq$ ($L$)	multiquark	Mol. ($L$)	${\omega }_{\mathrm{Mol}}$ [MeV]
$f_0\left(980\right)$	980	$(0,0^{+})$	$q\overline{q}$ ($P$) ($s\overline{s}$ ($P$))	$qs\overline{q}\overline{s}$	$\overline{K}K$ ($S$)	67.8(B)
$a_0\left(980\right)$	980	$(1,0^{+})$	$q\overline{q}$ ($P$)	$qs\overline{q}\overline{s}$	$\overline{K}K$ ($S$)	67.8(B)
$K\left(1460\right)$	1460	$(1/2,0^{-})$	—	$qq\overline{q}\overline{s}$ ($P$)	$\overline{K}KK$ ($P$)	69.0(R)
$\Lambda \left(1405\right)$	1405	$(0,1/2^{-})$	$uds$ ($P$)	$udsq\overline{q}$	$\overline{K}N$ ($S$)	20.5(R)
$\Delta \Delta$	2380	$(0,3^{+})$	—	$q^6$	—	—
$\Lambda \Lambda$-$N\Xi$ ($H$)	2245	$(0,0^{+})$	—	$uuddss$	$N\Xi$ ($S$)	73.2(B)
$N\Omega$	2592	$(1/2,2^{+})$	—	$uudsss$	—	—

Table 4. Table 2.4: Summary of exotic hadrons with heavy flavors. The notations are the same as those in Table 2.3 .

Particle	$m$ [MeV]	$(I,J^P)$	$q\overline{q}$/$qqq$ $\left(L\right)$	multiquark	Mol. ($L$)	${\omega }_{\mathrm{Mol}}$ [MeV]
$D_s\left(2317\right)$	2317	$(0,0^{+})$	$c\overline{s}$ ($P$)	$c\overline{s}q\overline{q}$	$DK$ ($S$)	273(B)
$X\left(3872\right)$	3872	$(0,1^{+})$	$c\overline{c}$ ($P$)	$c\overline{c}q\overline{q}$	$D{\overline{D}}^{\ast }$ ($S$)	3.6(B)
$Z_c\left(3900\right)$	3900	$(1,1^{+})$	—	$c\overline{c}u\overline{d}$	—	—
$Z_c\left(4430\right)$	4430	$(1,1^{+})$	—	$c\overline{c}u\overline{d}$	$D_1{\overline{D}}^{\ast }$ ($S$)	13.5(B)
$Z_b\left(10610\right)$	10610	$(1,1^{+})$	—	$b\overline{b}u\overline{d}$	—	—
$Z_b\left(10650\right)$	10650	$(1,1^{+})$	—	$b\overline{b}u\overline{d}$	—	—
$X\left(5568\right)$	5568	$(1,0^{+})$	—	$s\overline{b}u\overline{d}$	—	—
$P_c\left(4380\right)$	4380	${(1/2,3/2^{-})}^b$	—	$c\overline{c}uud$ ($S$)	$\overline{D}{\Sigma }_c^{\ast }$ ($S$)	60(B)
$P_c\left(4450\right)$	4450	${(1/2,5/2^{+})}^b$	—	$c\overline{c}uud$ ($P$)	—	—

Table 5. Table 2.5: Summary of other exotic hadrons. The notations are the same as those in Table 2.3 . More information for each particle is available in Refs.  [ 15 , 174 ] .

Particle	$m$ [MeV]	$(I,J^P)$	$q\overline{q}$/$qqq$	multiquark ($L$)	Mol. ($L$)	${\omega }_{\mathrm{Mol}}$ [MeV]	Ref.
$\Theta \left(1530\right)$	1530	$(0,1/2^{+})$	—	$qqqq\overline{s}$ ($P$)	—	—	[175]
$\overline{K}KN$	1920	$(1/2,1/2^{+})$	—	$qqqs\overline{s}$ ($P$)	$\overline{K}KN$	42(R)	[176]
$\overline{K}NN$	2352	$(1/2,0^{-})$	$q^{5}s$ ($P$)	$q^{6}s\overline{q}$ ($S$)	$\overline{K}NN$	20.5(T)	[177]
$\Omega \Omega$	3228	$(0,0^{+})$	—	$s^6$	—	—	[178]
$T_{cc}^1$	3797	$(0,1^{+})$	—	$ud\overline{c}\overline{c}$	—	—	[179]
$\overline{D}N$	2790	$(0,1/2^{-})$	—	$qqqq\overline{c}$	$\overline{D}N$	6.48(R)	[180]
${\overline{D}}^{\ast }N$	2919	$(0,3/2^{-})$	—	$qqqq\overline{c}\left(D\right)$	${\overline{D}}^{\ast }N$	6.48(R)	[181]
${\Theta }_{cs}$	2980	$(1/2,1/2^{+})$	—	$qqqs\overline{c}$ ($P$)	—	—	[182, 183]
$H_c^{++}$	3377	$(1,0^{+})$	—	$qqqqsc$	—	—-	[184]
$\overline{D}NN$	3734	$(1/2,0^{-})$	—	$q^7\overline{c}$	$\overline{D}NN$	6.48(T)	[185]
${\Lambda }_{c}N$	3225	$(1/2,1^{+})$	—	$cuduud$	${\Lambda }_{c}N$	4.24(R)	[174]
${\Lambda }_{c}NN$	4164	$(0,3/2^{+})$	—	$cuduududd$	${\Lambda }_{c}NN$	33.16(R)	[174]
$T_{cb}^0$	7123	$(0,0^{+})$	—	$ud\overline{c}\overline{b}$	—	—	[184]

Table 6. Table 3.1: Statistical and coalescence model parameters for Scenario 1 and 2 at RHIC (200 GeV), LHC (2.76 TeV) and LHC (5.02 TeV), and those given in Refs.  [ 14 , 15 ] . Quark masses are taken to be m q = 350 ​ MeV subscript 𝑚 𝑞 350 MeV m_{q}=350~{}\mathrm{MeV} , m s = 500 ​ MeV subscript 𝑚 𝑠 500 MeV m_{s}=500~{}\mathrm{MeV} , m c = 1500 ​ MeV subscript 𝑚 𝑐 1500 MeV m_{c}=1500~{}\mathrm{MeV} and m b = 4700 ​ MeV subscript 𝑚 𝑏 4700 MeV m_{b}=4700~{}\mathrm{MeV} . In Refs.  [ 14 , 15 ] , light quark masses were taken to be m q = 300 ​ MeV subscript 𝑚 𝑞 300 MeV m_{q}=300~{}\mathrm{MeV} .

	RHIC		LHC (2.76 TeV)		LHC (5.02 TeV)		RHIC	LHC (5 TeV)
	Sc. 1	Sc. 2	Sc. 1	Sc. 2	Sc. 1	Sc. 2	Refs [14, 15]
$T_H\left(\mathrm{MeV}\right)$	162		156				175
$V_H\left({\mathrm{fm}}^3\right)$	2100		5380				1908	5152
${\mu }_B\left(\mathrm{MeV}\right)$	24		0				20	0
${\mu }_s\left(\mathrm{MeV}\right)$	10		0				10	0
${\gamma }_c$	22		39		50		6.40	15.8
${\gamma }_b$	$4.0\times {10}^7$		$8.6\times {10}^8$		$1.4\times {10}^9$		$2.2\times {10}^6$	$3.3\times {10}^7$
$T_C\left(\mathrm{MeV}\right)$	162	166	156	166	156	166	175
$V_C\left({\mathrm{fm}}^3\right)$	2100	1791	5380	3533	5380	3533	1000	2700
$\omega \left(\mathrm{MeV}\right)$	590	608	564	609	564	609	550
${\omega }_s\left(\mathrm{MeV}\right)$	431	462	426	502	426	502	519
${\omega }_c\left(\mathrm{MeV}\right)$	222	244	219	278	220	279	385
${\omega }_b\left(\mathrm{MeV}\right)$	183	202	181	232	182	234	338
$N_u=N_d$	320	302	700	593	700	593	245	662
$N_s=N_{\overline{s}}$	183	176	386	347	386	347	150	405
$N_c=N_{\overline{c}}$	4.1		11		14		3	20
$N_b=N_{\overline{b}}$	0.03		0.44		0.71		0.02	0.8
$T_F\left(\mathrm{MeV}\right)$	119		115				125
$V_F\left({\mathrm{fm}}^3\right)$	20355		50646				11322	30569
$N_K$	67.5		134				142†	363†
$N_{\overline{K}}$	59.6		134				127†	363†
$N_N$	20		32				62†	150†
$N_{\Delta }$	18		28				—	—
$N_{\Lambda }$	3.8		6.5				—	—
$N_{\Xi }$	2.6		4.4				4.7	13
$N_{\Omega }$	0.37		0.62				0.81	2.3
$N_D=N_{\overline{D}}$	1.5		4.0		5.2		1.0	6.9
$N_{D^{\ast }}=N_{{\overline{D}}^{\ast }}$	2.0		5.4		6.9		1.5	10
$N_{D_1}=N_{{\overline{D}}_1}$	0.20		0.49		0.63		0.19	1.3
$N_B=N_{\overline{B}}$	$8.1\times {10}^{-3}$		0.12		0.20		$5.3\times {10}^{-3}$	0.21
$N_{B^{\ast }}=N_{{\overline{B}}^{\ast }}$	$1.9\times {10}^{-2}$		0.27		0.45		$1.2\times {10}^{-2}$	0.49
$N_{{\Lambda }_c}$	0.17		0.36		0.46		—	—
$N_{{\Sigma }_c}$	0.2		0.41		0.52		—	—
$N_{{\Sigma }_c^{\ast }}$	0.28		0.56		0.71		—	—
$N_{{\Xi }_c}$	0.11		0.25		0.32		0.10	0.65

Table 7. Table 3.2: Estimates of heavy quark pairs d ​ N / d ​ y 𝑑 𝑁 𝑑 𝑦 dN/dy at midrapidity in 0-10% central collision at RHIC and LHC.

	RHIC	LHC @2.76 TeV	LHC @5.02 TeV
Without shadowing
$N_c=N_{\overline{c}}$	4.5	17	23
$N_b=N_{\overline{b}}$	0.034	0.68	1.2
With shadowing
$N_c=N_{\overline{c}}$	4.1	11	14
$N_b=N_{\overline{b}}$	0.031	0.44	0.71

Table 8. Table 3.3: Summary of particle yields for light hadrons (cf. Table  2.3 ).

RHIC
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$f_0\left(980\right)$	2.1 (0.7)	$3.9\times {10}^{-2}$	2.1 (0.7)	$4.0\times {10}^{-2}$	1.7	3.5
$a_0\left(980\right)$	6.4	$1.2\times {10}^{-1}$	6.4	$1.2\times {10}^{-1}$	5.2	10
$K\left(1460\right)$	—	$5.8\times {10}^{-2}$	—	$5.7\times {10}^{-2}$	$1.3\times {10}^{-1}$	$6.3\times {10}^{-1}$
$\Lambda \left(1405\right)$	$4.7\times {10}^{-1}$	$2.3\times {10}^{-2}$	$4.5\times {10}^{-1}$	$2.4\times {10}^{-2}$	$7.3\times {10}^{-1}$	$8.6\times {10}^{-1}$
$\Delta \Delta$	—	$4.2\times {10}^{-3}$	—	$5.3\times {10}^{-3}$	—	$1.8\times {10}^{-2}$
$\Lambda \Lambda$-$N\Xi$ ($H$)	—	$4.7\times {10}^{-4}$	—	$5.0\times {10}^{-4}$	$1.6\times {10}^{-3}$	$4.9\times {10}^{-3}$
$N\Omega$	—	$1.7\times {10}^{-3}$	—	$1.9\times {10}^{-3}$	$1.4\times {10}^{-3}$	$6.7\times {10}^{-3}$
LHC (2.76 TeV)
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$f_0\left(980\right)$	4.3 (1.2)	$5.4\times {10}^{-2}$	4.1 (1.2)	$6.0\times {10}^{-2}$	3.2	6.6
$a_0\left(980\right)$	13	$1.6\times {10}^{-1}$	12	$1.8\times {10}^{-1}$	9.5	$20$
$K\left(1460\right)$	—	$8.2\times {10}^{-2}$	—	$8.0\times {10}^{-2}$	$1.9\times {10}^{-1}$	1.0
$\Lambda \left(1405\right)$	$7.5\times {10}^{-1}$	$2.9\times {10}^{-2}$	$7.0\times {10}^{-1}$	$3.2\times {10}^{-2}$	1.1	1.4
$\Delta \Delta$	–	$5.8\times {10}^{-3}$	–	$1.0\times {10}^{-2}$	—	$1.9\times {10}^{-2}$
$\Lambda \Lambda$-$N\Xi$ ($H$)	—	$5.0\times {10}^{-4}$	—	$6.1\times {10}^{-4}$	$1.8\times {10}^{-3}$	$5.9\times {10}^{-3}$
$N\Omega$	—	$1.8\times {10}^{-3}$	—	$2.3\times {10}^{-3}$	$1.6\times {10}^{-3}$	$7.8\times {10}^{-3}$
LHC (5.02 TeV)
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$f_0\left(980\right)$	4.3 (1.2)	$5.4\times {10}^{-2}$	4.1 (1.2)	$6.0\times {10}^{-2}$	3.2	6.6
$a_0\left(980\right)$	13	$1.6\times {10}^{-1}$	12	$1.8\times {10}^{-1}$	9.5	20
$K\left(1460\right)$	—	$8.2\times {10}^{-2}$	—	$8.0\times {10}^{-2}$	$1.9\times {10}^{-1}$	1.0
$\Lambda \left(1405\right)$	$7.5\times {10}^{-1}$	$2.9\times {10}^{-2}$	$7.0\times {10}^{-1}$	$3.2\times {10}^{-2}$	1.1	1.4
$\Delta \Delta$	—	$5.8\times {10}^{-3}$	—	$1.0\times {10}^{-2}$	—	$1.9\times {10}^{-2}$
$\Lambda \Lambda$-$N\Xi$ ($H$)	—	$5.0\times {10}^{-4}$	—	$6.1\times {10}^{-4}$	$1.8\times {10}^{-3}$	$5.9\times {10}^{-3}$
$N\Omega$	—	$1.8\times {10}^{-3}$	—	$2.3\times {10}^{-3}$	$1.6\times {10}^{-3}$	$7.8\times {10}^{-3}$

Table 9. Table 3.4: Summary of particle yields for heavy hadrons (cf. Table  2.4 ).

RHIC
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$D_s\left(2317\right)$	$2.3\times {10}^{-2}$	$2.4\times {10}^{-3}$	$2.3\times {10}^{-2}$	$2.5\times {10}^{-3}$	$6.5\times {10}^{-3}$	$6.6\times {10}^{-2}$
$X\left(3872\right)$	$5.4\times {10}^{-4}$	$5.0\times {10}^{-5}$	$5.6\times {10}^{-4}$	$5.3\times {10}^{-5}$	$9.1\times {10}^{-4}$	$5.7\times {10}^{-4}$
$Z_c\left(3900\right)$	—	$1.5\times {10}^{-4}$	—	$1.6\times {10}^{-4}$	—	$1.5\times {10}^{-3}$
$Z_c\left(4430\right)$	—	$1.5\times {10}^{-4}$	—	$1.6\times {10}^{-5}$	$5.0\times {10}^{-5}$	$6.5\times {10}^{-5}$
$Z_b\left(10610\right)$	—	$2.0\times {10}^{-9}$	—	$2.1\times {10}^{-9}$	—	$2.1\times {10}^{-8}$
$Z_b\left(10650\right)$	—	$2.0\times {10}^{-9}$	—	$2.1\times {10}^{-9}$	—	$1.6\times {10}^{-8}$
$X\left(5568\right)$	—	$5.1\times {10}^{-5}$	—	$5.2\times {10}^{-5}$	—	$2.3\times {10}^{-3}$
$P_c\left(4380\right)$	—	$2.5\times {10}^{-5}$	—	$2.6\times {10}^{-5}$	$2.9\times {10}^{-5}$	$9.2\times {10}^{-5}$
$P_c\left(4450\right)$	—	$1.5\times {10}^{-5}$	—	$1.5\times {10}^{-5}$	—	$9.1\times {10}^{-5}$
LHC (2.76 TeV)
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$D_s\left(2317\right)$	$5.2\times {10}^{-2}$	$4.3\times {10}^{-3}$	$5.0\times {10}^{-2}$	$4.5\times {10}^{-3}$	$1.4\times {10}^{-2}$	$1.5\times {10}^{-1}$
$X\left(3872\right)$	$1.6\times {10}^{-3}$	$1.1\times {10}^{-4}$	$1.7\times {10}^{-3}$	$1.3\times {10}^{-4}$	$2.7\times {10}^{-3}$	$1.7\times {10}^{-3}$
$Z_c\left(3900\right)$	—	$3.4\times {10}^{-4}$	—	$4.0\times {10}^{-4}$	—	$4.3\times {10}^{-3}$
$Z_c\left(4430\right)$	—	$3.4\times {10}^{-4}$	—	$4.0\times {10}^{-4}$	$1.4\times {10}^{-4}$	$1.7\times {10}^{-4}$
$Z_b\left(10610\right)$	—	$1.3\times {10}^{-7}$	—	$1.5\times {10}^{-7}$	—	$1.9\times {10}^{-6}$
$Z_b\left(10650\right)$	—	$1.3\times {10}^{-7}$	—	$1.5\times {10}^{-7}$	—	$1.5\times {10}^{-6}$
$X\left(5568\right)$	—	$5.0\times {10}^{-4}$	—	$5.2\times {10}^{-4}$	—	$3.1\times {10}^{-2}$
$P_c\left(4380\right)$	—	$5.0\times {10}^{-5}$	—	$5.8\times {10}^{-5}$	$6.4\times {10}^{-5}$	$2.1\times {10}^{-4}$
$P_c\left(4450\right)$	—	$2.9\times {10}^{-5}$	—	$3.2\times {10}^{-5}$	—	$2.0\times {10}^{-4}$
LHC (5.02 TeV)
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$D_s\left(2317\right)$	$6.5\times {10}^{-2}$	$5.4\times {10}^{-3}$	$6.4\times {10}^{-2}$	$5.7\times {10}^{-3}$	$1.8\times {10}^{-2}$	$1.9\times {10}^{-1}$
$X\left(3872\right)$	$2.5\times {10}^{-3}$	$1.8\times {10}^{-4}$	$2.7\times {10}^{-3}$	$2.1\times {10}^{-4}$	$4.5\times {10}^{-3}$	$2.8\times {10}^{-3}$
$Z_c\left(3900\right)$	—	$5.4\times {10}^{-4}$	—	$6.4\times {10}^{-4}$	—	$7.1\times {10}^{-3}$
$Z_c\left(4430\right)$	—	$5.4\times {10}^{-4}$	—	$6.4\times {10}^{-4}$	$2.3\times {10}^{-4}$	$2.8\times {10}^{-4}$
$Z_b\left(10610\right)$	—	$3.4\times {10}^{-7}$	—	$3.9\times {10}^{-7}$	—	$5.0\times {10}^{-6}$
$Z_b\left(10650\right)$	—	$3.4\times {10}^{-7}$	—	$3.9\times {10}^{-7}$	—	$3.9\times {10}^{-6}$
$X\left(5568\right)$	—	$7.9\times {10}^{-4}$	—	$8.2\times {10}^{-4}$	—	$5.0\times {10}^{-2}$
$P_c\left(4380\right)$	—	$7.9\times {10}^{-5}$	—	$9.3\times {10}^{-5}$	$1.0\times {10}^{-4}$	$3.4\times {10}^{-4}$
$P_c\left(4450\right)$	—	$4.7\times {10}^{-5}$	—	$5.0\times {10}^{-5}$	—	$3.4\times {10}^{-4}$

Table 10. Table 3.5: Summary of particle yields for other hadrons (I) (cf. Table  2.5 ).

RHIC
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$\Theta \left(1530\right)$	—	$6.7\times {10}^{-3}$	—	$6.7\times {10}^{-3}$	—	$5.0\times {10}^{-1}$
$\overline{K}KN$	—	$5.0\times {10}^{-3}$	—	$5.1\times {10}^{-3}$	$4.2\times {10}^{-2}$	$1.2\times {10}^{-1}$
$\overline{K}NN$	$7.3\times {10}^{-4}$	$2.7\times {10}^{-5}$	$7.4\times {10}^{-4}$	$2.9\times {10}^{-5}$	$3.9\times {10}^{-3}$	$5.8\times {10}^{-3}$
$\Omega \Omega$	—	$8.2\times {10}^{-6}$	—	$9.4\times {10}^{-6}$	—	$1.5\times {10}^{-5}$
LHC (2.76 TeV)
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$\Theta \left(1530\right)$	—	$8.2\times {10}^{-3}$	—	$8.5\times {10}^{-3}$	—	$6.8\times {10}^{-1}$
$\overline{K}KN$	—	$6.0\times {10}^{-3}$	—	$6.6\times {10}^{-3}$	$5.1\times {10}^{-2}$	$1.5\times {10}^{-1}$
$\overline{K}NN$	$7.9\times {10}^{-4}$	$2.3\times {10}^{-5}$	$8.6\times {10}^{-4}$	$3.0\times {10}^{-5}$	$3.9\times {10}^{-3}$	$6.3\times {10}^{-3}$
$\Omega \Omega$	—	$7.6\times {10}^{-6}$	—	$1.2\times {10}^{-5}$	—	$1.8\times {10}^{-5}$
LHC (5.02 TeV)
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$\Theta \left(1530\right)$	—	$8.2\times {10}^{-3}$	—	$8.5\times {10}^{-3}$	—	$6.8\times {10}^{-1}$
$\overline{K}KN$	—	$6.0\times {10}^{-3}$	—	$6.6\times {10}^{-3}$	$5.2\times {10}^{-2}$	$1.5\times {10}^{-1}$
$\overline{K}NN$	$7.9\times {10}^{-4}$	$2.3\times {10}^{-5}$	$8.6\times {10}^{-4}$	$3.0\times {10}^{-5}$	$3.9\times {10}^{-3}$	$6.3\times {10}^{-3}$
$\Omega \Omega$	—	$7.6\times {10}^{-6}$	—	$1.2\times {10}^{-5}$	—	$1.8\times {10}^{-5}$

Table 11. Table 3.6: Summary of particle yields for other hadrons (II) (cf. Table  2.5 ).

RHIC
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$T_{cc}^1$	—	$5.0\times {10}^{-5}$	—	$5.3\times {10}^{-5}$	—	$8.9\times {10}^{-4}$
$\overline{D}N$	—	$2.6\times {10}^{-3}$	—	$2.6\times {10}^{-3}$	$1.3\times {10}^{-2}$	$1.0\times {10}^{-2}$
${\overline{D}}^{\ast }N$	—	$9.8\times {10}^{-4}$	—	$9.3\times {10}^{-4}$	$1.1\times {10}^{-2}$	$9.6\times {10}^{-3}$
${\Theta }_{cs}$	—	$7.4\times {10}^{-4}$	—	$7.4\times {10}^{-4}$	—	$6.4\times {10}^{-3}$
$H_c$	—	$2.7\times {10}^{-4}$	—	$2.8\times {10}^{-4}$	—	$5.7\times {10}^{-4}$
$\overline{D}NN$	—	$1.8\times {10}^{-5}$	—	$1.8\times {10}^{-5}$	$9.4\times {10}^{-5}$	$5.1\times {10}^{-5}$
${\Lambda }_{c}N$	—	$1.5\times {10}^{-3}$	—	$1.5\times {10}^{-3}$	$5.0\times {10}^{-3}$	$2.9\times {10}^{-3}$
${\Lambda }_{c}NN$	—	$6.7\times {10}^{-6}$	—	$6.7\times {10}^{-6}$	$2.9\times {10}^{-6}$	$9.8\times {10}^{-6}$
$T_{cb}^0$	—	$9.3\times {10}^{-8}$	—	$9.9\times {10}^{-8}$	—	$1.6\times {10}^{-6}$
LHC (2.76 TeV)
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$T_{cc}^1$	—	$1.1\times {10}^{-4}$	—	$1.3\times {10}^{-4}$	—	$2.7\times {10}^{-3}$
$\overline{D}N$	—	$4.3\times {10}^{-3}$	—	$4.2\times {10}^{-3}$	$2.3\times {10}^{-2}$	$1.9\times {10}^{-2}$
${\overline{D}}^{\ast }N$	—	$1.6\times {10}^{-3}$	—	$1.3\times {10}^{-3}$	$2.0\times {10}^{-2}$	$1.8\times {10}^{-2}$
${\Theta }_{cs}$	—	$1.2\times {10}^{-3}$	—	$1.2\times {10}^{-3}$	—	$1.2\times {10}^{-2}$
$H_c$	—	$3.8\times {10}^{-4}$	—	$4.0\times {10}^{-4}$	—	$8.6\times {10}^{-4}$
$\overline{D}NN$	—	$2.0\times {10}^{-5}$	—	$2.0\times {10}^{-5}$	$1.1\times {10}^{-4}$	$6.7\times {10}^{-5}$
${\Lambda }_{c}N$	—	$2.2\times {10}^{-3}$	—	$2.2\times {10}^{-3}$	$7.0\times {10}^{-3}$	$4.3\times {10}^{-3}$
${\Lambda }_{c}NN$	—	$6.7\times {10}^{-6}$	—	$6.5\times {10}^{-6}$	$2.7\times {10}^{-6}$	$9.9\times {10}^{-6}$
$T_{cb}^0$	—	$1.1\times {10}^{-6}$	—	$1.3\times {10}^{-6}$	—	$2.7\times {10}^{-5}$
LHC (5.02 TeV)
Particle	$q\overline{q}/qqq$	multiquark	$q\overline{q}/qqq$	multiquark	Mol.	Stat.
	scenario 1		scenario 2
$T_{cc}^1$	—	$1.8\times {10}^{-4}$	—	$2.1\times {10}^{-4}$	—	$4.4\times {10}^{-3}$
$\overline{D}N$	—	$5.3\times {10}^{-3}$	—	$5.3\times {10}^{-3}$	$3.0\times {10}^{-2}$	$2.4\times {10}^{-2}$
${\overline{D}}^{\ast }N$	—	$2.0\times {10}^{-3}$	—	$1.7\times {10}^{-3}$	$2.6\times {10}^{-2}$	$2.3\times {10}^{-2}$
${\Theta }_{cs}$	—	$1.5\times {10}^{-3}$	—	$1.4\times {10}^{-3}$	—	$1.6\times {10}^{-2}$
$H_c$	—	$4.7\times {10}^{-4}$	—	$4.9\times {10}^{-4}$	—	$1.1\times {10}^{-3}$
$\overline{D}NN$	—	$2.5\times {10}^{-5}$	—	$2.5\times {10}^{-5}$	$1.5\times {10}^{-4}$	$8.6\times {10}^{-5}$
${\Lambda }_{c}N$	—	$2.7\times {10}^{-3}$	—	$2.7\times {10}^{-3}$	$9.1\times {10}^{-3}$	$5.5\times {10}^{-3}$
${\Lambda }_{c}NN$	—	$8.2\times {10}^{-6}$	—	$8.0\times {10}^{-6}$	$3.5\times {10}^{-6}$	$1.3\times {10}^{-5}$
$T_{cb}^0$	—	$2.3\times {10}^{-6}$	—	$2.7\times {10}^{-6}$	—	$5.6\times {10}^{-5}$

Table 12. Table 4.1: Interaction parameters

set	$m\left(\mathrm{GeV}\right)$	$\mu \left(\mathrm{GeV}\right)$	$\lambda \left({\mathrm{GeV}}^4\right)$	$E_c\left(\mathrm{GeV}\right)$
A	$0.125$	$0.5$	$3.0\times {10}^{-3}$	$0.12$
B	$0.125$	$0.5$	$6.25\times {10}^{-3}$	$0.15$
C	$0.125$	$0.5$	$6.25\times {10}^{-4}$	$0.10$

Table 13. Table 4.2: Resonance properties

set	${ℰ}_r\left(\mathrm{GeV}\right)$	$p_r\left(\mathrm{GeV}\right)$	${⟨c|r⟩⟩}^2$
A	$0.080-i0.026$	$0.143-i0.023$	$1.157-i0.116$
B	$0.047-i0.065$	$0.126-i0.064$	$1.758-i0.383$
C	$0.092-i0.0052$	$0.152-i0.0043$	$1.024-i0.019$

Table 14. Table 4.3: Density matrix parameters

$\beta \left({\mathrm{GeV}}^{-1}\right)$	$L\left({\mathrm{GeV}}^{-1}\right)$	$\alpha \left({\mathrm{GeV}}^{-2}\right)$	$n_{\overline{K}N}/n_{\pi \Sigma }$	${\rho }_c$
$8.4$	$1055$	$200$	$1.00$	$0.0025$

Table 15. Table 4.4: The integrated non-resonance and resonance parts divided by ρ c subscript 𝜌 𝑐 \rho_{c}

set	A	B	C
${\Pi }_{ab}^{\left(1\right)nr}/{\rho }_c$	$0.00$	$0.18$	$0.00$
${\Pi }_{ab}^{\left(1\right)r}/{\rho }_c$	$0.86$	$1.39$	$0.12$
${\Pi }_c^{\left(1\right)nr}/{\rho }_c$	$0.00$	$0.11$	$0.00$
${\Pi }_c^{\left(1\right)r}/{\rho }_c$	1.00	$0.89$	$1.00$
${\Pi }^{\left(2\right)nr}/{\rho }_c$	$0.11$	$1.18$	$0.00$
${\Pi }^{\left(2\right)r}/{\rho }_c$	$-0.76$	$-2.78$	$-0.12$

Table 16. Table 4.5: Production probabilities with the resonance wave function, Π ~ r superscript ~ Π 𝑟 \tilde{\Pi}^{r} ’s, and the integrated resonance parts, Π r superscript Π 𝑟 \Pi^{r} ’s, divided by ρ c subscript 𝜌 𝑐 \rho_{c}

set	A	B	C
${\tilde{\Pi }}_c^r/{\rho }_c$	$1.15-0.12i$	$1.76-0.39i$	$1.02-0.00i$
${\Pi }_c^{\left(1\right)r}/{\rho }_c$	$1.00$	$0.89$	$1.00$
${\tilde{\Pi }}_{ab}^r/{\rho }_c$	$-0.11+0.12i$	$-0.38+0.05i$	$-0.01+0.03i$
${\Pi }_{ab}^r/{\rho }_c$	$-0.01$	$-1.36$	$-0.00$
${\tilde{\Pi }}^r/{\rho }_c$	$1.04-0.00i$	$1.36-0.34i$	$1.01+0.03i$
${\Pi }^r/{\rho }_c$	$0.99$	$-0.49$	$1.00$

Table 17. Table 5.2: The binding energy ( E B subscript 𝐸 B E_{\rm B} ), scattering length ( a 0 subscript 𝑎 0 a_{0} ) and effective range ( r eff subscript 𝑟 eff r_{\rm eff} ) with and without the Coulomb attraction in the spin-2 p ​ Ω 𝑝 Ω p\Omega state. Physical masses of the proton and Ω Ω \Omega are used.

Spin-2 $p\Omega$ potentials		$V_{\text{I}}$	$V_{\text{II}}$	$V_{\text{III}}$
	$E_B$ [MeV]	$-$	0.05	24.8
without Coulomb	$a_0$ [fm]	$-1.0$	23.1	1.60
	$r_{\text{eff}}$ [fm]	1.15	0.95	0.65
	$E_B$ [MeV]	$-$	6.3	26.9
with Coulomb	$a_0$ [fm]	$-1.12$	5.79	1.29
	$r_{\text{eff}}$ [fm]	1.16	0.96	0.65

Equations

$$M_{X(3872)}$$

$$M_{D^0}+M_{{\bar{D}}^{\ast 0}}$$

$${\rm Y}+{\pi }_{\mathrm{thermal}}\to Z_b\to {\rm Y}+\pi ,B+{\bar{B}}^{\ast },$$

$$N_{h}^{stat}=V_{H}\frac{g_{h}}{2\pi^{2}}\int_{0}^{\infty}\frac{p^{2}dp}{\gamma_{h}^{-1}e^{E_{h}/T_{H}}\pm 1}\approx\frac{\gamma_{h}g_{h}V_{H}}{2\pi^{2}}m_{h}^{2}T_{H}K_{2}\Big{(}\frac{m_{h}}{T_{H}}\Big{)}\approx\gamma_{h}g_{h}V_{H}\left(\frac{m_{h}T_{H}}{2\pi}\right)^{3/2}e^{-m_{i}/T_{H}},$$

$${\gamma }_h={\gamma }_c^{n_c+n_{\bar{c}}}{\gamma }_b^{n_b+n_{\bar{b}}}{e}^{({\mu }_{B}B+{\mu }_{s}S)/T_H},$$

$$N_c$$

$$N_b$$

$$N_h^{\mathrm{coal}}=g_h\int \left[\prod _{i=1}^n\frac{1}{g_i}\frac{p_i\cdot d{\sigma }_i}{(2\pi {)}^3}\frac{{\mathrm{d}}^3{\mathbf{p}}_i}{E_i}f(x_i,p_i)\right]f^W(x_1,\cdots ,x_n:p_1,\cdots ,p_n),$$

$$\int p_i\cdot d{\sigma }_i\frac{{\mathrm{d}}^3{\mathbf{p}}_i}{(2\pi {)}^3{E}_i}f(x_i,p_i)=N_i,$$

$$f^{W}(x_{1},\cdots,x_{n}:p_{1},\cdots,p_{n})=\int\prod_{i=1}^{n}dy_{i}e^{ip_{i}y_{i}}\psi^{*}\Big{(}x_{1}+\frac{y_{1}}{2},\cdots,x_{n}+\frac{y_{n}}{2}\Big{)}\psi\Big{(}x_{1}-\frac{y_{1}}{2},\cdots,x_{n}-\frac{y_{n}}{2}\Big{)},$$

$$N_h^{\mathrm{coal}}=g_h\prod _{j=1}^n\frac{N_j}{g_j}\prod _{i=1}^{n-1}\frac{\int d^3{y}_i{d}^3{k}_i{f}_i(k_i)f^W(y_i,k_i)}{\int d^3{y}_i{d}^3{k}_i{f}_i(k_i)},$$

$$f_j(p_j)\propto \delta (Y_j-{\eta }_j)\exp \left(-\frac{p_{j,T}^2}{2m_{j}T}\right).$$

$$\prod _{j=1}^n\exp \left(-\frac{p_{j,T}^2}{2m_{j}T}\right)=\exp \left(-\frac{P_T^2}{2MT}\right)\prod _{i=1}^{n-1}{\tilde{f}}_i(k_i),$$

$$\frac{1}{{\mu }_i}=\frac{1}{m_{i+1}}+\frac{1}{{\sum }_{j=1}^i{m}_j},$$

$${\mu }_1=\frac{m_1{m}_2}{m_1+m_2},{\mu }_2=\frac{m_3(m_1+m_2)}{m_1+m_2+m_3},{\mu }_3=\frac{m_4(m_1+m_2+m_3)}{m_1+m_2+m_3+m_4},\mathrm{and}\mathrm{so}on.$$

$$f_s^W(y_i,k_i)$$

$$f_p^W(y_i,k_i)$$

$$f_d^W(y_i,k_i)$$

$$F({\sigma }_i,{\mu }_i,l_i,T)\equiv \frac{\int d^3{y}_i{d}^2{k}_i{\tilde{f}}_i(k_i)f^W(y_i,k_i)}{\int d^2{k}_i{\tilde{f}}_i(k_i)}=\frac{(4\pi {\sigma }_i^2{)}^{3/2}}{1+2{\mu }_{i}T{\sigma }_i^2}\frac{(2l_i)!!}{(2l_i+1)!!}{\left[\frac{2{\mu }_{i}T{\sigma }_i^2}{1+2{\mu }_{i}T{\sigma }_i^2}\right]}^{l_i},$$

$$N_h^{\mathrm{coal}}$$

$$N_{\Lambda (1115)}^{\mathrm{stat},\mathrm{total}}$$

$$N_{\Lambda (1115)}^{\mathrm{coal},\mathrm{total}}=1.60+\frac{1}{3}\times 4.81+1.28+\left(0.87+\frac{0.11}{3}\right)\times 9.62+2.57+0.51=\mathrm{16.3.}$$

$$N_{{\Lambda }_c(2286)}^{\mathrm{stat},\mathrm{total}}$$

$$N_{{\Lambda }_c(2286)}^{\mathrm{coal},\mathrm{total}}=0.064+0.193+0.385+0.67\times 0.051=0.676$$

$$N_{{\Lambda }_b(5620)}^{\mathrm{stat},\mathrm{total}}$$

$$N_{{\Lambda }_b(5620)}^{\mathrm{coal},\mathrm{total}}$$

$$\frac{1}{g_i}\frac{N_i}{V}\frac{(4\pi {\sigma }_i^2{)}^{3/2}}{(1+2{\mu }_{i}T{\sigma }_i^2)}\sim$$

$$\frac{1}{g_i}\frac{N_i}{V}\frac{2}{3}\frac{(4\pi {\sigma }_i^2{)}^{3/2}2{\mu }_{i}T{\sigma }_i^2}{(1+2{\mu }_{i}T{\sigma }_i^2{)}^2}\sim$$

$$\frac{1}{g_i}\frac{N_i}{V}\frac{8}{15}\frac{(4\pi {\sigma }_i^2{)}^{3/2}(2{\mu }_{i}T{\sigma }_i^2{)}^2}{(1+2{\mu }_{i}T{\sigma }_i^2{)}^3}\sim$$

$${\tau }_{scatt}^i=\frac{1}{{\sum }_j⟨{\sigma }_{ij}{v}_{ij}⟩n_j},$$

$${\tau }_{exp}=\frac{1}{\partial \cdot u},$$

$${\tau }_{exp}={\tau }_{scatt}\to \frac{R}{3dR/dt}=\frac{1}{n\sigma ⟨v⟩}.$$

$$\frac{N}{R_{fo}^2}=\frac{4\pi }{{\sigma }_{fo}},$$

Full text

††journal: Progress of Particle and Nuclear Physics

1. 1 Introduction

2. 2 Current status of exotic hadrons

3. 2.1 Light hadrons

4. 2.1.1 Scalar mesons

5. 2.1.2 $\Lambda(1405)$

6. 2.1.3 Dibaryons

7. 2.2 Heavy hadrons

8. 2.2.1 $D_{s0}^{\ast}(2317)$

9. 2.2.2 Charmonium-like states

10. 2.2.3 Charged bottomonium-like states

11. 2.2.4 New states: $P_{c}$ and $X(5568)$

12. 3 Yields of particles

13. 3.1 Statistical model

14. 3.2 Coalescence model

15. 3.2.1 Quark coalescence

16. 3.2.2 Hadron coalescence

17. 3.2.3 Heavy quark pair production

18. 3.3 Freeze-out conditions for molecular states

19. 3.4 Yields of hadrons

20. 4 Coalescence model for resonances

21. 4.1 Coalescence model

22. 4.2 Model for S-wave resonance

23. 4.3 Coalescence model for scattering states

24. 4.4 Numerical examples

25. 4.5 Summary and discussions

26. 5 Hadron-hadron interactions from two particle momentum correlations

27. 5.1 General property of the two-particle momentum correlation function

28. 5.1.1 Formalism

29. 5.1.2 Correlations from strong interactions and quantum statistics

30. 5.1.3 Lednický and Lyuboshits Model

31. 5.1.4 Effect of collectivity

32. 5.1.5 Feed-down contribution

33. 5.2 Non-exotic channels

34. 5.2.1 $pp$ and $\bar{p}\bar{p}$ correlation

35. 5.2.2 $p\Lambda$ and $p\bar{\Lambda}$ correlations

36. 5.3 Exotic channels

37. 5.3.1 $\Lambda\Lambda$ correlation

38. 5.3.2 $p\Omega$ correlation

39. 5.3.3 $K^{-}p$ Correlation

40. 6 Summary

1 Introduction

High-energy heavy-ion collisions provide a unique opportunity to study the properties of the high energy density QCD matter formed at the initial stage of the collisions. It is now well established from the study of jet quenching and anisotropic flow [1] that the produced high energy density matter is strongly interacting with a very small specific viscosity. These collisions at the same time provide the suitable conditions for producing weakly bound hadronic states such as the light nuclei, hypernuclei [2] and anti-nuclei [3]. An interesting experimental finding is that the yields of these nuclei follow the statistical model predictions with temperature and chemical potentials that are fitted to the yields of the ground state particles [4]. These results thus suggest that the final abundance of ground state hadrons and light nuclei are already determined near the energy density at which quarks and gluons hadronize. However, when analyzing particles beyond the ground states, one finds that their yields sometimes deviate from the statistical model predictions, which, on the other hand, may reveal valuable information about the evolution of the hadronic system after the hadronization of quarks and gluons and/or the structure of these particles. For example, resonances with large width that are reconstructed from daughter particles are found to be less produced than the statistical model prediction [5], suggesting that hadronic interactions and their freeze-out conditions play important roles in determining the final yields. Moreover, depending on their quantum numbers, the yields of excited states can be either suppressed [6] or enhanced [7] relative to those predicted by the statistical model. These results suggest the importance of the structure of a hadron on its yield in heavy ion collisions.

During the past decade, there has also been a revival of hadron physics research due to the observation of many exotic state candidates. These findings started from the observation of $D_{sJ}$(2317) [8] and X(3872) [9] and continued on to the recently observed pentaquark states at the LHC [10]. Hadrons of such exotic structures have been proposed since the early days of bag models as the color confinement does not rule out the existence of a multiquark configuration in a color singlet state [11, 12]. A multiquark configuration can be either a flavour exotic state or a cryptoexotic state. For the latter, one cannot distinguish its structure among an internally excited state, compact multiquark state or a molecular configuration, as shown in Fig. 1.1 based just on its quantum numbers. Even for flavor exotic states, a multiquark configuration can always be decomposed into meson and baryon states, making it impossible to discriminate a compact multiquark state from a hadronic molecular state just from its quantum numbers. Moreover, some observed structures from an elementary process that have been interpreted as exotic states could be caused by kinematical effects. It is worth to note that the existence of stable tetraquark configurations can also influence the properties of QCD at finite temperature and density as the tetraquark condensation may lead to a second chiral phase transition [13].

In Refs. [14, 15] published by the authors of present review, it is proposed that heavy-ion collisions at ultrarelativistic energies provide a unique opportunity to study exotic particle candidates. Due to the abundant number of heavy quarks and antiquarks produced in these collisions, various exotic hadrons could be formed, and their yields have been estimated by using the quark coalescence model. When the parameters in the model are fit to the yields of ground state hadrons as predicted from the statistical model, it has been found that the yield of a hadron is typically an order of magnitude smaller when it is a compact multiquark state than that of an excited hadronic state with normal quark numbers and/or a molecular configuration. Combined with the fact that only resonances with large natural width are affected by the subsequent hadronic evolution, it is thus possible to determine if an exotic hadron produced in relativistic heavy ion collisions is a compact multiquark state or has a hadronic molecular configuration. Also, owing to the significant numbers of charm and bottom quarks produced at RHIC and even larger numbers at LHC, some of the proposed and recently measured heavy exotic hadrons could be produced with sufficient abundance for experimental detection, making it possible to study these new exotic hadrons in heavy ion collisions. Moreover, the structures in the invariant mass spectrum of a specific decay channel of, for example, a B meson, which are generated by kinematic effects, cannot be produced statistically and thus will not appear in heavy ion collisions. Therefore, heavy ion collisions also make it possible to discriminate such effects from real resonances.

Another important recent development is the measurement of two-particle momentum correlations in relativistic heavy ion collisions. In recent experiments, correlations have been measured for particle pairs such as $p\bar{p}$ [16], $\bar{p}\bar{p}$ [16], $p\Lambda$ [17, 18, 19], $\bar{p}\Lambda(p\bar{\Lambda})$ [18], and $\Lambda\Lambda$ [20, 21, 22, 23, 24] in addition to $\pi\pi$, $KK$ and $pp$. The two-particle correlation is generated by quantum statistics and final state interactions, and it also depends on the size and lifetime of the emission source [25, 26, 27, 28, 29, 30]. Therefore, one can use, on the one hand, the two-particle correlation to determine the source size if the interaction between the two particles is known. On the other hand, it is possible to extract information on the interactions between two particles by using the experimental data on their correlation functions [31, 32, 33, 34, 35, 36, 37], if the property of their emission source is known. The latter provides, in particular, a unique opportunity to explore interactions between short-lived hadrons, such as the meson-meson, hyperon-hyperon and antibaryon-antibaryon interactions, which can serve as crucial inputs for understanding possible exotic hadronic states, such as the hadronic molecules and dibaryons. While the $\Lambda\Lambda$ interaction is accessible via the double $\Lambda$ hypernuclei [38], there is no other ways to access other short-lived hadron-hadron interactions experimentally. In heavy ion collisions, the freeze-out conditions are well studied and the chaotic source assumption is known to work reasonably. It is thus possible to constrain various hadron-hadron interactions from the experimental data on two-particle correlations.

This review is organized as follows. It begins with a survey of the current status on our understanding of exotic hadrons. This is followed by a discussion on the general issues related to particle production in heavy-ion collisions. An updated account is then given on the yields of some potentially exotic hadrons that could be measured in experiments on heavy-ion collisions. Next, we will present the new development in the coalescence model for resonances. Also, reviews are given on the basic theoretical framework and the current status on the study of hadron-hadron interactions using two-particle correlations.

2 Current status of exotic hadrons

2.1 Light hadrons

2.1.1 Scalar mesons

The scalar mesons $f_{0}(980)$ and $a_{0}(980)$, together with $f_{0}(500)~{}(\sigma)$ and $K_{0}^{\ast}(800)~{}(\kappa)$, have been thought to have some exotic structures, since they exhibit an inverted mass spectrum compared to what is expected if they have simple $q\bar{q}$ configurations. Model studies have suggested that $f_{0}(980)$ and $a_{0}(980)$ could be compact $qq\bar{q}\bar{q}$ systems [11, 12], $K\bar{K}$ molecules in $s$ wave [39, 40], and dynamically generated states in the $K\bar{K}$ and $\pi\pi$/$\pi\eta$ coupled-channel chiral dynamics [41, 42, 43]. However, the RHIC data seem to rule out a dominant tetraquark configuration for $f_{0}(980)$ as shown in previous studies by the authors [14, 15]. Therefore, the structure of these scalar mesons is still highly controversial.

Recently the so-called compositeness has been introduced to hadron physics for investigating the internal structure of hadrons in terms of the hadronic molecular configuration [44, 45]. The compositeness is defined as the two-body composite part of the normalization of the total wave function, and with this quantity it has been found that the $f_{0}(980)$ resonance in the coupled-channel chiral dynamics is indeed dominated by the $K\bar{K}$ composite state [45]. The $K\bar{K}$ compositeness of $f_{0}(980)$ and $a_{0}(980)$ has been evaluated in Ref. [46] from experimental observations, and is found to have a large absolute value for $f_{0}(980)$ and a small but nonnegligible absolute value for $a_{0}(980)$. Moreover, the spatial structure, which reflects the hadronic spatial expanse inside the system, has been theoretically studied for $f_{0}(980)$ by using the finite volume method [47], and the distance between $K$ and $\bar{K}$ inside $f_{0}(980)$ has been estimated as $2.6$–$3.0\text{ fm}$.

2.1.2 $\Lambda(1405)$

The $\Lambda(1405)$ resonance is a negative parity excited state of $\Lambda$ baryon in the energy region between the $\pi\Sigma$ and $\bar{K}N$ thresholds [48]. While it should be described as a $p$-wave excited baryon in constituent quark models, the observed mass is too small in comparison with other negative parity baryons [49]. Rather, the $\Lambda(1405)$ is considered to be a meson-baryon molecule, as initiated in a study using the phenomenological model with vector meson exchange [50, 51]. In recent studies, the meson-baryon molecule picture is understood on the basis of chiral SU(3) dynamics [52, 53, 54, 55, 56, 57]. The $\bar{K}N$ molecular nature of the $\Lambda(1405)$ is also supported by the lattice QCD simulation [58], as well as the analysis with the weak-binding relation [59]. A remarkable feature of this resonance is the two-pole nature [60], which stems from the attraction in the singlet and octet channels in the SU(3) basis [61] and the $\bar{K}N$ and $\pi\Sigma$ channels in the isospin basis [62]. A recent striking experimental achievement is the precise measurement of the kaonic hydrogen by SIDDHARTA [63, 64], which significantly reduces the uncertainty of the subthreshold extrapolation of the $\bar{K}N$ amplitude in the $\Lambda(1405)$ region [65, 66].

The spatial structure of the $\Lambda(1405)$ has been studied by evaluating its form factor in the chiral unitary approach [67, 68]. Here the form factor is evaluated at the higher energy pole, which gives a dominant contribution to the $\Lambda(1405)$. By switching off the decay into the $\pi\Sigma$ channel, the mean distance between $\bar{K}$ and $N$ in the $\Lambda(1405)$ is obtained as $\langle r^{2}\rangle\sim 2.8$ fm2, which is in fair agreement with the estimation by the finite volume method [47]. By using the effective single-channel $\bar{K}N$ potential constrained by the SIDDHARTA data, the mean distance is evaluated as $\langle r^{2}\rangle\sim 2.1$ fm2 [69].

From the experimental viewpoint, the $\Lambda(1405)$ has been observed in the low energy exclusive reactions. Traditionally, the kaon and pion beams have been used [70, 71], and several high-statistics data are recently available in photoproductions by the LEPS collaboration [72] and the CLAS collaboration [73, 74], and in the proton-proton collisions by the HADES collaboration [75]. It is remarkable that the spin-parity $1/2^{-}$ is experimentally determined by the CLAS collaboration [76].

So far, the $\Lambda(1405)$ is not observed in high-energy inclusive processes, such as heavy ion collisions. This is mainly because its decay is dominated by the $\pi\Sigma=\{\pi^{+}\Sigma^{-},\pi^{0}\Sigma^{0},\pi^{-}\Sigma^{+}\}$ mode, except for a tiny fraction of the radiative decays into $\gamma\Lambda$ and $\gamma\Sigma^{0}$. In order to measure the $\pi\Sigma$ modes, it is necessary to detect at least one neutral particle ($\pi^{0}$ or neutron) in the weak decay of $\Sigma$, which is in general difficult. If the detection of neutrons is possible in heavy ion collisions, the $\Lambda(1405)$ can be reconstructed from the $\pi^{+}\Sigma^{-}$ mode where $\Sigma^{-}$ decays into $n\pi^{-}$ by almost 100%. Alternatively, the $\pi^{-}\Sigma^{+}$ mode is also feasible, where the $\Sigma^{+}$ has 50% probability of decaying into $n\pi^{+}$.

2.1.3 Dibaryons

The $S=-2$ H dibaryon was first predicted in Ref. [77] as the color spin interaction is the most attractive and the Pauli principle between quarks does not operate in the $J=0,I=0$ channel in flavour SU(3). Extensive experimental searches ruled out the possibility of a deeply bound H dibaryon [78, 38]. The two-$\Lambda$ binding energy in the double $\Lambda$ hypernucleus observed in the Nagara event is determined to be $B_{\Lambda\Lambda}({}_{\Lambda\Lambda}^{6}\mathrm{He})=6.91\pm 0.16$ MeV [38, 79], which sets the H mass range as $M_{H}>2M_{\Lambda}-B_{\Lambda\Lambda}({}_{\Lambda\Lambda}^{6}\mathrm{He})$. There is still a possibility that the H particle exists as a shallow bound state or as a resonance. Experimentally, a bump structure was observed between the $\Lambda\Lambda$ and the $N\Xi$ threshold ($10-15$ MeV above the $\Lambda\Lambda$ threshold) in the $\Lambda\Lambda$ invariant mass spectrum from the $(K^{-},K^{+})$ reaction [21], while no clear signal was found so far in the $\Lambda p\pi^{-}$ invariant mass analysis from heavy ion collisions by ALICE [23] as well as in the $\Lambda\Lambda$ invariant mass spectrum from $\Upsilon(1S)$ and $\Upsilon(2S)$ decays by Belle [24]. Lattice calculations show that the H dibaryon becomes bound in the massive pion cases [80, 81], and it may evolve to a resonance near the $\Xi N$ threshold at a smaller pion mass and with the SU(3)f breaking effects, by combining physical hadron masses with the SU(3) symmetric potentials at unphysical quark masses [82]. The latter result is consistent with the old prediction of the quark cluster model calculation [83] that takes account of the instanton-induced interaction effects [84, 85]. Chiral extrapolations of the lattice QCD results to the physical point indicate an unbound H dibaryon with respect to the $\Lambda\Lambda$ threshold [86, 87, 88]. A constituent quark model calculation also shows that a compact multiquark configuration would be highly unlikely [89]. The $\Lambda\Lambda$ correlation in heavy-ion collisions and its implication to the $\Lambda\Lambda$ interaction are discussed in Sec. 5.3.1.

As in the H particle case, stable or resonance dibaryon states may appear in those channels in which the color spin interaction is attractive [90, 91] and the Pauli principle does not operate, since the repulsive core of the baryon-baryon interaction is mainly due to the Pauli principle between quarks in the presence of the color spin interaction in the quark-quark force [92]. The SU(3)f breaking and channel coupling effects are also important, since a bound or resonance state is sensitive to the threshold. The channels with the most attractive color spin interactions with $S=0,-1,-2,-3$ are $\Delta\Delta(I=0,J^{\pi}=3^{+})$, $N\Sigma^{*}-\Delta\Sigma(I=1/2,J^{\pi}=2^{+})$, $\Lambda\Lambda-N\Xi-\Sigma\Sigma(I=0,J^{\pi}=0^{+})$, and $N\Omega-\Lambda\Xi^{*}-\Sigma^{*}\Xi-\Sigma\Xi^{*}(I=1/2,J^{\pi}=2^{+})$, respectively [91]. Among these channels, the third ($S=-2$) is the H particle channel, and recent lattice and experimental dibaryon studies suggest that $N\Omega$ and $\Delta\Delta$ may have bound states.

The $N\Omega$ channel was also found to be attractive in the quark model mainly due to the changes in the spatial wave function as the two baryons merge and the absence of any repulsion from the color spin interaction [93]. In fact, recent lattice calculations find an attractive potential in the S-wave spin 2 channel that allows for a bound state with binding energy in the order of 18.9 MeV [94], although the pion mass in the simulation is still larger than the physical one. The $N\Omega$ correlation function can tell if the $N\Omega$ bound state exists [37] as discussed in Sec. 5.3.2.

The nature of the low-mass enhancement in the $\pi\pi$-invariant mass spectrum from the double-pion fusion reaction [95, 96], called the Abashian-Booth-Crowe (ABC) effect, is recently found to be related to a possible resonance structure called the $d^{*}(2380)$ [97, 98, 99, 100, 101, 102]. The resonance structure has a width of only 70 MeV and about 80 MeV below the $2\Delta$ threshold [103]. There are works that claim the state is a dibaryon state [104] but from a constituent quark model analysis, it is problematic if a compact multiquark configuration can be stable in that channel [105].

2.2 Heavy hadrons

2.2.1 $D_{s0}^{\ast}(2317)$

Recent progress in experimental facilities, particularly the B factories [106], has provided an opportunity to investigate the hadron spectroscopy involving heavy charm and bottom quarks. One of the most important tasks in the heavy-hadron spectroscopy is to understand how the properties of the exotic hadron candidates change with the heavy quark masses because one expects the existence of heavy exotic states that are analogous to the corresponding light exotic hadrons, such as the light scalar mesons discussed in Sec. 2.1.1.

In this respect, the charmed and strange scalar meson $D_{s0}^{\ast}(2317)$, with quantum numbers $(I,J^{P})=(0,0^{+})$, is of special interest. The $D_{s0}^{\ast}(2317)$ was first observed by BaBar through its isospin violating $\pi^{0}D_{s}^{+}$ decay mode [8], which was followed by its confirmation by CLEO [107], Belle [108], and FOCUS [109]. An important feature of $D_{s0}^{\ast}(2317)$ is that its mass $M_{D_{s0}^{\ast}(2317)}=2317.7\pm 0.6\text{ MeV}$ [48] is about 160 MeV below the prediction of a quark model for the charmed meson [110, 111], implying that it has some exotic configuration besides an ordinary $q\bar{q}$ configuration. There have thus been extensive discussions on the $D_{s0}^{\ast}(2317)$, such as a $c\bar{s}$ state [112, 113, 114, 115, 116], two-meson molecular state [117, 118, 119, 120, 121, 122], $D$-$K$ mixing [123], tetraquark state [124, 125, 126, 127], or a mixture of two-meson and tetraquark states [128]. Also, in terms of the compositeness, the $KD$ molecular component of the $D_{s0}^{\ast}(2317)$ has been studied experimentally via some observables [129], and a lattice simulation [130] has indicated the existence of a dominant $KD$ molecular component of about $70\%$ in the $D_{s0}^{\ast}(2317)$ wave function.

2.2.2 Charmonium-like states

In recent years, there has been a remarkable progress in the heavy meson spectroscopy [131, 132, 133, 134, 48]. Plenty of new states, called $XYZ$, are observed above the open charm/bottom thresholds. Because the properties of these states are not well described in the conventional constituent quark model, the $XYZ$ states are expected to have an exotic structure. Among many interesting states, here we summarize the current status of $X(3872)$ and the charged charmonium-like states, $Z_{c}^{\pm}$.

One of the most intensively studied states is the $X(3872)$. It is firstly observed by the Belle collaboration in the $B$ decay [9]. Subsequently, the $X(3872)$ is confirmed by the CDF collaboration [135], the D0 collaboration [136], the BaBar collaboration [137], the LHCb collaboration [138], and the CMS collaboration [139]. The mass $M_{X(3872)}$ and width $\Gamma_{X(3872)}$ are given by [48]

[TABLE]

It is worth noting that the mass is very close to the threshold of the $D^{0}\bar{D}^{*0}$ state

[TABLE]

The $C$-parity $C=-1$ is determined by measuring the $J/\psi\gamma$ decay [140, 141]. The spin-parity can be studied from the angular distribution of the final states. Eventually, the LHCb collaboration determines $J^{PC}=1^{++}$ based on the angular correlations [142]. No charged partners are observed in the $J/\psi\pi^{\pm}\pi^{0}$ mode.

While the quantum numbers of the $X(3872)$ can in principle be given by the $\bar{c}c$ configuration, the proximity of the $D^{0}\bar{D}^{*0}$ state urges us to consider the hadronic molecule interpretation. In fact, the phenomenological Lagrangian analysis of the radiative decay [143] indicates the dominance of the $D^{0}\bar{D}^{*0}$ molecular component. On the other hand, an analysis of the prompt production cross section [144] shows that the $\bar{c}c$ component is required in addition to the molecule component. In addition, lattice study in Ref. [145] indicates that a candidate for the $X(3872)$ is found only when both the $\bar{c}c$ and $D\bar{D}^{*}$ operators are included in the analysis.

The charged charmonium-like states $Z_{c}^{\pm}$ are of particular interest, because a light quark-antiquark pair is required in addition to $\bar{c}c$ as the valence component. At present, altogether eight charged charmonium-like states have been reported, although not all the states are firmly established. Below we briefly overview the experimental status of these charged charmonium-like states. The basic properties are summarized in Table 2.1.

The first candidate is called $Z_{c}(4430)$, observed in the $\pi^{\pm}\psi(2S)$ spectrum of the $B$ decay into $K\pi^{\pm}\psi(2S)$ by the Belle collaboration [146]. Although the BaBar collaboration did not confirm the $Z_{c}(4430)$ [147], the Dalitz plot analysis of the $B$ decay by Belle reconfirmed the original findings [148]. The spin-parity $J^{P}=1^{+}$ is favored, based on the full amplitude analysis by Belle [149]. The $Z_{c}(4430)$ is finally confirmed by LHCb [150] with more than ten times higher statistics than the original observation. The analysis of LHCb [150] also suggests the possibility of additional state $Z_{c}(4240)$ in the $\pi^{-}\psi(2S)$ spectrum. Two other charged states $Z_{c}(4050)$ and $Z_{c}(4250)$ were found in the $\pi^{+}\chi_{c1}$ mass distribution in $\bar{B}^{0}\to K^{-}\pi^{+}\chi_{c1}$ by Belle [151]. These states were not confirmed in the analysis by BaBar [152].

The $Z_{c}(3900)$ state is observed in the $\pi^{\pm}J/\psi$ spectrum of the process $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ by the BESIII collaboration [153]. The state was confirmed in the same decay mode by Belle [154] and by the analysis of CLEO-c data [155]. BESIII also looked for other decay modes. The $Z_{c}(3900)$ signal was not found in $h_{c}\pi^{\pm}$ spectrum [156], while the structure is observed in the $(D\bar{D}^{*})^{\pm}$ mode [157]. Angular analysis by BESIII favors the spin-parity assignment of $J^{P}=1^{+}$ [157]. A recent lattice QCD study finds no convincing signal for the $Z_{c}(3900)$ candidate [158, 159]. The result by the HAL QCD collaboration indicates the interpretation of $Z_{c}(3900)$ as a threshold cusp effect [160].

A narrow state was observed by Belle in the $h_{c}\pi^{\pm}$ spectrum, which is called $Z_{c}(4020)$ [156]. The state is also found in the $(D^{*}\bar{D}^{*})^{\pm}$ mode [161]. In the study of the $J/\psi\pi^{+}$ spectrum of the $\bar{B}^{0}\to K^{-}J/\psi\pi^{+}$ decay by Belle, the $Z_{c}(4200)$ was observed [162]. Favored spin-parity assignment turns out to be $J^{P}=1^{+}$. Belle also found a signal of $Z_{c}(4055)$ in the $J/\psi\pi^{+}$ channel [163].

Among many charged states, we will concentrate on the structure of the most established states, $Z_{c}(3900)$ and $Z_{c}(4430)$. In the molecular interpretation, these are considered to be $\bar{D}D^{*}$ and $D_{1}\bar{D}$ molecules, respectively. A novel idea based on a strong diquark-antidiquark correlation is also considered for the structure of X and Z states [164].

2.2.3 Charged bottomonium-like states

As bottom analogue of charged charmonium-like state, $Z_{c}^{+}$, charged bottomonium-like state called $Z_{b}$ was reported in Belle [165]. They are $Z_{b}(10610)^{+}$ and $Z_{b}(10650)^{+}$ with spin-parity $J^{P}=1^{+}$. The neutral state, $Z_{b}(10610)^{0}$, was also reported [166]. Because the masses of the two $Z_{b}$’s are very close to $B\bar{B}^{\ast}$ ($B^{\ast}\bar{B}$) and $B^{\ast}\bar{B}^{\ast}$ thresholds, respectively, it might be natural to regard those states as a $B\bar{B}^{\ast}$ ($B^{\ast}\bar{B}$) hadronic molecule and a $B^{\ast}\bar{B}^{\ast}$ hadronic molecule, respectively [167, 168, 169]. However, $Z_{b}(10610)^{+}$ and $Z_{b}(10650)^{+}$ decay to $B\bar{B}^{\ast}$ ($B^{\ast}\bar{B}$) and $B^{\ast}\bar{B}^{\ast}$, and hence they cannot be simple $B\bar{B}^{\ast}$ ($B^{\ast}\bar{B}$) and $B^{\ast}\bar{B}^{\ast}$ bound states. The hadronic molecule structure can be studied in the decay properties of $Z_{b}^{+}\rightarrow\Upsilon\pi^{+}$ in view of the heavy quark symmetry [167, 168, 170, 171] and within a phenomenological Lagrangian approach [172]. However, the situation may be much different, because the observed experimental peaks of $Z_{b}$’s could be explained by the cusp effect [173].

In relativistic heavy ion collisions, enhancement of the scattering amplitude of $\Upsilon+\pi$ can be seen via

[TABLE]

with the intermediate $Z_{b}$ states, where $\pi_{\mathrm{thermal}}$ is the thermal pion in the hadron phase.

2.2.4 New states: $P_{c}$ and $X(5568)$

The existence of hidden-charm pentaquark states has been predicted in many theoretical works [186, 187, 188, 189, 190]. Therefore, it was with great excitement that the theoretical community heard about the LHCb observation of two hidden-charm pentaquark-like structures $P^{+}_{c}(4380)$ and $P^{+}_{c}(4450)$ in the $J/\psi p$ invariant mass distribution in the decay $\Lambda^{0}_{b}\to J/\psi pK^{-}$ [10]. They used an amplitude analysis of the three-body final-state, and the masses and widths obtained are $M_{P_{c}(4380)}=(4380\pm 8\pm 29)$ MeV, $\Gamma_{P_{c}(4380)}=(205\pm 18\pm 86)$ MeV, $M_{P_{c}(4450)}=(4449.8\pm 1.7\pm 2.5)$ MeV, and $\Gamma_{P_{c}(4450)}=(39\pm 5\pm 19)$ MeV. The significance of the lower mass and higher mass states is 9 $\sigma$ and 12 $\sigma$, respectively. The preferred spin-parity assignments are $J^{P}=3/2^{\pm}$ or $5/2^{\mp}$.

Before the LHCb observation, the predicted masses for genuine pentaquark states, with both negative and positive parities [187], or meson-baryon bound states with $J^{P}=3/2^{-}$ [186, 188, 189], cover the observed masses of the two $P^{+}_{c}$ structures. After the observation of the two $P^{+}_{c}$ structures, many theoretical works appeared proposing various explanations for these structures  [191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209]. It was even suggested that the observed structures could be due to kinematical triangle singularities [210, 211, 212]. Of course this possibility needs to be examined by future experiments.

Since the masses of these two $P^{+}_{c}$ structures are very close to the mass thresholds of the $\bar{D}\Sigma_{c}^{*}$ and $\bar{D}^{*}\Sigma_{c}$, a very natural explanation is that the observed structures could be meson-baryon molecules [191, 192, 193, 194, 195, 196, 197, 198]. Other possible explanations are: diquark-diquark-antiquark pentaquarks [199, 200, 201, 202, 203], compact diquark-triquark pentaquarks [204, 205, 206, 207], topological soliton model [208], genuine multiquark states other than molecules [209], etc. For more comprehensive discussions, see Refs. [213, 214, 215].

The D0 Collaboration has recently announced the observation of a narrow enhancement of the experimental data in the $B_{s}^{0}\pi^{\pm}$ mass spectrum in the energy around 5.6 GeV. The enhancement was interpreted as a new state: $X^{\pm}(5568)$ [216]. The mass and width for this state have been found to be $m=5567.8\pm 2.9(\mbox{sta})^{+0.9}_{-1.9}(\mbox{syst}$ MeV/c2 and $\Gamma=21.9\pm 6.4(\mbox{sta})^{+5.0}_{-2.5}(\mbox{syst})$ MeV/c2, respectively [216]. The isospin of $X(5568)$ is clearly one. Its spin-parity is not yet known although a scalar four quark interpretation has been suggested in Ref. [216]. The $X(5568)$ would be a very interesting addition to the list of undoubtedly exotic mesons, since its wave function consists of four different flavors: $u$, $b$, $d$ and $s$ quarks. However, the LHCb Collaboration has not confirmed the observation of the $X(5568)$. In their analysis [217] no structure is found in the $B_{s}^{0}\pi^{\pm}$ mass spectrum from the $B_{s}^{0}\pi^{+}$ threshold up to $M_{B_{s}^{0}\pi^{+}}\leq 5700$ GeV.

The announcement of the exotic state $X(5568)$ stimulated the theoretical interest and several theoretical works have been done to investigate the properties of such state. There are studies based on QCD sum rules [218, 219, 220, 221, 222, 223, 224, 225, 226, 227], quark models [228, 229, 230, 231, 232, 233, 234, 235], rescattering effects [236], coupled-channel analysis [237] and more general arguments [238, 239]. In Ref. [238] various interpretations for the $X(5568)$ signal were considered and the authors concluded that threshold, cusp, tetraquark and molecular models were all unfavored. In Ref. [239] additional arguments, based on general properties of QCD, were provided to question the existence of the $X(5568)$. In Refs. [219, 228, 233, 234] although it was possible to find a tetraquark state with $J^{P}=0^{+}$ and the same quark content as the $X(5568)$, the obtained masses were around 200 MeV higher than the announced $X(5568)$ mass [228, 234], or too high to be candidate of $X(5568)$ [219, 232, 233]. However, the width of the state can be of the same order as the width reported by D0 [223]. In Ref. [233] no molecular structure could be obtained to explain the $X(5568)$ state, while in Ref. [227] molecular and tetraquark states could be obtained but with masses around 5200 MeV.

In all other calculations it was possible to explain the properties of the $X(5568)$. In particular, in Refs. [218, 220, 221, 222, 224, 225, 226, 235] the results of the calculations for the properties of $X(5568)$ are in excellent agreement with the experimental value. In Refs. [218, 220, 222, 224, 225, 226, 235] $J^{P}=0^{+}$ was assumed while in Ref. [221] scalar as well as axial tetraquark currents were considered. In Ref. [228] a model using multiquark interactions has been used and a 150 MeV higher mass is found for $X(5568)$, although the systematic errors still allow their state to be related to $X(5568)$. Another multiquark model calculations using the color-magnetic interaction has been presented in Refs. [229, 231] with very good agreement with the experimental value. The possibility of explaining the enhancement in the data as near threshold rescattering effects has been studied in Ref. [236]. The $B\bar{K}$ and $B^{*}\bar{K}$ molecular interpretations have been suggested in Ref. [230]. In Ref. [237], using a coupled-channel analysis, it is possible to find a pole that can be associated with the $X(5568)$ state, although the cutoff used is much larger than the normal one.

Clearly, more analysis are required to clarify this situation from the experimental side as well as from the theoretical side.

3 Yields of particles

In relativistic heavy ion collisions, a variety of hadrons and their resonances are produced. To describe the yields of these particles, a number of approaches have been used. These include the microscopic transport model [240, 241] and the macroscopic hydrodynamic model [242, 243]. In the transport model, hadron production from the produced quark-gluon plasma is treated either kinetically via parton to hadron reactions [240] or via the quark coalescence [241]. In the hydrodynamic approach, hadrons are produced through the statistical hadronization model [244].

On the other hand, the schematic model, which is based on isentropic boost invariant longitudinal and accelerated transverse expansions, has been adopted in Ref. [245] to understand the production of exotic hadrons in heavy ion collisions and to quantitatively evaluate their yields from the dynamically expanding quark-gluon plasma. In this approach, both the initial QGP and final hadronic matter are treated as noninteracting free gases. For the crossover transition between these two phases of matter, it starts at the critical temperature $T_{C}$ and ends at the hadronization temperature $T_{H}$. During this phase transition, the system expands from the critical volume $V_{C}$ at $T_{C}$ to the hadronization volume $V_{H}$ at $T_{H}$ while maintains a constant entropy.

In the present review on the yields of hadrons in relativistic heavy ion collisions, we update the results reported in Refs. [14, 15] by using the schematic model of Ref. [245]. Following this reference, we take the same hadronization temperature $T_{H}=162~{}(156)$ MeV and volume $V_{H}=2100~{}(5380)$ fm3 for RHIC (LHC) as used in the statistical hadronization model analysis of the experimental data [246, 247]. The critical temperature $T_{C}$ and volume $V_{C}$ are then determined by tracing back in time and using the entropy conservation condition $s_{H}V_{H}=s_{C}V_{C}$, where $s_{C}$ and $s_{H}$ are the entropy density of the system at $T_{C}$ and $T_{H}$, respectively. The information on the entropy density of the system $s=(\epsilon+p)/T$ at different temperatures is available from the lattice results $p(T)/T^{4}$ and $(\epsilon-3p)(T)/T^{4}$ for the pion mass 135 MeV [248]. Figure 3.1 shows the temperature dependence of the volume of the matter produced at RHIC (LHC) during its isentropic expansion.

To determine the critical temperature $T_{C}$ and volume $V_{C}$, we consider here two possibilities. The first scenario is to take the critical temperature to be the same as the hadronization temperature, $T_{C}=T_{H}$, and also for the critical volume $V_{C}=V_{H}$. Since hadrons are continuously produced during the crossover phase transition, we require that the number of hadrons produced via the statistical hadronization is equal to that produced via coalescence of constituent quarks at the end of the phase transition.

In the second scenario, we assume that the critical temperature $T_{c}$ and the oscillator frequency, used in the harmonic oscillator model for the wave functions of composite systems, are the same at RHIC and LHC, and determine the respective volumes by requiring that the numbers of various hadrons produced via the quark coalescence at the beginning of phase transition ($T_{C}$) are equal to those produced at the end of transition ($T_{H}$) via the statistical hadronization.

Although hadron yields in the statistical hadronization model only depends on the hadronization temperature $T_{H}$ and volume $V_{H}$, those in the quark coalescence model also depend on the the number of quarks in the system and the oscillator frequency in the wave function of the produced hadron. The critical temperature $T_{C}$ and volume $V_{C}$ can thus be determined by choosing an appropriate value for the oscillator frequency. Using the rho and omega mesons as examples, we obtain the critical temperature $T_{C}=166$ MeV at RHIC and LHC. Hence, in the second scenario, we find that the critical temperature at LHC drops from 166 MeV to 156 MeV, while that at RHIC changes slightly from 166 MeV to 162 MeV during the crossover transition.

In Table 3.1, we tabulate the critical temperature $T_{C}$ and volume $V_{C}$ at the beginning of the crossover transition between the quark-gluon plasma and hadronic matter, and the hadronization temperature $T_{H}$ and volume $V_{H}$ at the end of the mixed phase or hadronization for the two scenarios. Also given in Table 3.1 are the baryon and strange chemical potentials evaluated in the statistical hadronization model [246]. The small strange chemical potential at RHIC, where the produced QGP has small net baryons or baryon chemical potential $\mu_{B}$, indicates that there is an approximate chemical equilibrium of strangeness at RHIC, and this is due to the short equilibration time and net zero strangeness in the QGP. Because of the higher energy, longer lifetime, and almost zero net baryons of QGP at LHC, a complete chemical equilibrium is reached for strangeness, resulting in a zero strange chemical potential.

3.1 Statistical model

Both the statistical hadronization model and the coalescence model have been used to evaluate the yields of exotic hadrons produced from heavy ion collisions. The statistical model has been very successful in explaining the relative yields of normal hadrons in relativistic heavy ion collisions. In the statistical hadronization model, the number of produced hadrons of a given type $h$ is given by [244]

[TABLE]

with $g_{h}$ being the degeneracy of the hadron, $\gamma_{h}$ the fugacity, $K_{2}$ the modified Bessel function of the second kind, and $V_{H}$ and $T_{H}$, respectively, the volume and temperature of the source for the statistical production of hadrons. The fugacity $\gamma_{h}$ of hadron species $h$, is generally expressed as

[TABLE]

where $B$, $S$, $n_{c}(n_{\bar{c}})$ and $n_{b}(n_{\bar{b}})$ are the baryon number, strangeness, (anti-)charm quark number, and (anti-)bottom quark number of the hadron, respectively. Values of the strangeness chemical potentials for heavy ion collision at RHIC and LHC are listed in Table 3.1. Since charm and bottom quarks are mostly produced from initial hard scattering, their numbers are much larger than those expected from a chemically equilibrated QGP. Therefore, we consider the fugacity $\gamma_{h}>1$ for both charmed and bottom hadrons. The fugacity of $n$-multiple charm quark hadrons is the product of $n$ charm quark fugacity $\gamma_{c}$, or $\gamma_{c}^{n}$, and also same for bottomed hadrons.

As shown in Ref. [15], the charm and bottom fugacities $\gamma_{c}$ and $\gamma_{b}$ can be determined by requiring that the total yield of charm or bottom hadrons in the statistical hadronization model to be the same as the total charm $N_{c}$ or bottom $N_{b}$ quark number from initial hard nucleon-nucleon scattering. With the values $N_{c}=4.1$ and $N_{b}=0.03$ for heavy ion collisions at RHIC at 200 GeV, which will be explained in detail in the next session, we obtain $\gamma_{c}=22$ and $\gamma_{b}=4.0\times 10^{6}$. For examples, when $N_{c}=4.1$ and $N_{b}=0.03$,

[TABLE]

In the above evaluation, the average yield of heavy anti-strange and strange mesons as well as that of heavy anti-baryons and baryons have been considered in order to average out the effect of strangeness and baryon chemical potentials. A similar analysis for LHC at 2.76 (5.02) TeV based on the charm and bottom quark numbers $N_{c}=11~{}(14)$ and $N_{b}=0.44~{}(0.71)$ gives the charm and bottom fugacities $\gamma_{c}=39~{}(50)$ and $\gamma_{b}=8.6\times 10^{8}~{}(1.4\times 10^{9})$, respectively. In Table 3.1, we give the fugacities needed for the evaluation of exotic hadron yields at RHIC and LHC in heavy ion collisions.

Also shown in Table 3.1 are the yields of various kinds of hadrons obtained in the statistical hadronization model, Eq. (3.1). We have evaluated here only the yield of directly produced hadrons, and have not taken into account the feed-down contributions. For examples, the yield of nucleons $N_{N}$ includes only directly produced $N$ at chemical freeze-out, and those from the strong decay of $\Delta$ have not been included.

3.2 Coalescence model

The coalescence model describes the production of hadrons through the recombination or coalescence of constituents. It is based on the sudden approximation by considering the overlap between the Wigner function of the produced particle and the density matrix of the constituents. The model has been successful in explaining the enhanced production of baryons compared to that of mesons at mid-rapidity in the intermediate transverse momentum [249] and the quark number scaling of elliptic flows [250] in heavy ion collisions. It has also been extensively adopted to investigate hadron production from the quark-gluon plasma in relativistic heavy ion collisions [251, 252, 249, 253, 254]. In this model, the number of hadron of type $h$ produced from the coalescence of $n$ constituents is given by [252]

[TABLE]

where $g_{h}$ is the degeneracy of the hadron, $g_{i}$ is that of its $i$th constituent, and d$\sigma_{i}$ is an element of a space-like hypersurface. The covariant phase-space distribution function of the constituents, $f(x_{i},p_{i})$ are normalized to their numbers, i.e.,

[TABLE]

and the Wigner function of the produced hadron, $f^{W}(x_{1}...x_{n}:p_{1}...p_{n})$ is defined by

[TABLE]

in terms of its wave function $\psi(x_{1},\cdots,x_{n})$.

In the non-relativistic limit and using a spherically symmetric harmonic oscillator function for the hadron wave function, the above equation becomes [255, 256],

[TABLE]

where $f^{W}(y_{i},k_{i})$ is the Wigner function associated with the internal (relative) wave function with its internal (relative) spatial and momentum coordinates $y_{i}$ and $k_{i}$. For heavy ion collisions at the RHIC and LHC energies, one can assume that the momentum distributions $f_{j}(p_{j})$ are Boltzmann at temperature $T$ for the transverse momentum $p_{j,T}$ and the strong Bjorken correlation of equal momentum-energy rapidities, $Y_{j}=\log[(E_{j}+p_{j,z})/(E_{j}-p_{j,z})]/2$ and space-time rapidities, $\eta_{j}=\log[(t_{j}+z_{j})/(t_{j}-z_{j})]/2$ for the longitudinal momentum, i.e.,

[TABLE]

Using the relation

[TABLE]

with $P_{T}$ and $M$ being the total transverse momentum and mass, respectively, the 2-dimensional momentum distribution functions of the constituents become $\tilde{f}_{i}(k_{i})\propto e^{-k_{i}^{2}/(2\mu_{i}T)}$ in the Jacobi coordinates. In the above, the reduced constituent masses $\mu_{i}$ are defined by

[TABLE]

or, explicitly

[TABLE]

The rapidity variables can be simplified at midrapidities, or $Y=\eta\sim 0$ as $\eta_{j}\simeq z_{j}/t_{j}$ and $Y_{j}\simeq p_{j,z}/m_{j}$ in the non-relativistic limit. Therefore, as long as the time $t_{j}$ when the coalescence occurs after the collision is large compared with the internal time scale of the hadron, $1/\omega$, or $t_{j}\gg 1/\omega$, one can omit in the Wigner function $f^{W}$ the contribution from the longitudinal momentum. In this case, the 3-dimensional momentum integrations in Eq. (3.7) reduces to 2-dimensional ones over the Wigner functions $f^{W}(y_{i},k_{i})$ and $\tilde{f}_{i}(k_{i})$ in transverse momentum $k_{i}$.

For a uniform distribution of particles in the emission source, the Wigner functions for the $s$-wave, $p$-wave, and $d$-wave are given explicitly as

[TABLE]

with the parameters $\sigma_{i}=1/\sqrt{\mu_{i}\omega}$ related to the reduced constituent masses $\mu_{i}$ and the oscillator frequency $\omega$.

After carrying out the phase-space integrals in Eq. (3.7), one gets the coalescence factor for each relative coordinate,

[TABLE]

with $l_{i}$ being the angular momentum of the wave function associated with the relative coordinate $y_{i}$. The final expression for the yield of hadrons in the coalescence model is then

[TABLE]

where $l_{i}$ is [math], $1$, and $2$ for an $s$-wave, a $p$-wave and a $d$-wave constituent, respectively, and $L=\sum_{i=1}^{n-1}l_{i}$, and $M=\sum_{i=1}^{n}m_{i}$. In the above, the relation $\mu_{i}\sigma_{i}^{2}=1/\omega$ has been used to convert the main dependence on $l_{i}$ into the form of the orbital angular momentum sum $L$. When $L\geq 2$, the factor in Eq. (3.14) depends on the way $L$ is decomposed into $l_{i}$, e.g., when $L=2$ and $n=3$ a factor $4/9$ has to be considered for the combination $(l_{1},l_{2})=(1,1)$, while a factor $8/15$ for $(l_{1},l_{2})=(2,0)$.

3.2.1 Quark coalescence

To evaluate the yields of hadrons produced from the QGP at the critical temperature $T_{c}$ when the volume is $V_{C}$ in the coalescence model, one needs to find the appropriate oscillator frequency $\omega$. One can choose the oscillator frequencies for light, strange, charmed, and bottom hadrons, denoted by $\omega_{q}$, $\omega_{s}$, $\omega_{c}$ and $\omega_{b}$, respectively, in the quark coalescence to reproduce the yields of certain normal hadrons in the statistical model. These oscillator frequencies can then be used to predict the yields of exotic hadrons.

For hadrons composed of only up and down quarks, one fixes the oscillator frequency $\omega_{q}$ to obtain in the coalescence model a similar yield of omega mesons as in the statistical model. This is shown in Table 3.1 together with the numbers of light quarks, $N_{u}=N_{d}$, and that of strange quarks, $N_{s}$, in the quark-gluon plasma, which are obtained from the statistical model based on the hadronization temperature $T_{H}$ and volume $V_{H}$ given in the table.

The parameter $\omega_{s}$ needed to evaluate the yield of hadrons composed of light and strange quarks in the coalescence model is determined by fitting the statistical model prediction for $\Lambda(1115)$ including the contribution from resonance decays [15]. Considering states that decay dominantly to $\Lambda(1115)$, one obtains the following result for heavy ion collisions for scenario 1 at RHIC:

[TABLE]

where 0.87 and $0.11/3$ in the parentheses represent the branching ratios of $\Sigma(1385)\rightarrow\Lambda+\pi$ and $\Sigma(1385)\rightarrow\Sigma^{0}+\pi$ in the $\Sigma(1385)$ decay, respectively. All numbers are calculated based on values in Table 3.1 for $T_{H}$ and $V_{H}$ with $\mu_{s}=10$ MeV and $\mu_{B}=24$ MeV at RHIC. With the constituent quark masses $m_{u,d}=350$ MeV and $m_{s}=500$ MeV for light and strange quarks, the above result can be reproduced by the coalescence model with $\omega_{s}=431$ MeV after taking into account the contributions from the decay of same resonances as in Eq. (3.15). Specifically, one has in the coalescence model

[TABLE]

The same method is used to determine the oscillator frequencies $\omega_{s}$ for the two scenarios at both RHIC and LHC, and the values are given in Table 3.1.

The oscillator frequency for charmed hadrons, $\omega_{c}$, is determined from reproducing the yield of $\Lambda_{c}(2286)$ including the feed-down contributions in the statistical model [257] but without considering the effect of diquarks [258]. Resonances included in the feed-down contributions include $\Sigma_{c}(2455)$, $\Sigma_{c}(2520)$, and $\Lambda_{c}(2625)$ as states of higher masses are negligible. For scenario 1 at RHIC, one then obtains

[TABLE]

from the statistical model and

[TABLE]

from the coalescence model, if the charm quark mass $m_{c}=1500$ MeV and the value $\omega_{c}=222$ MeV are used. The oscillator frequencies $\omega_{c}$ for other scenarios at both RHIC and LHC are summarized in Table 3.1, where it is seen that the different total charm quark number 11 at LHC 2.76 TeV energy and 14 at 5.02 TeV energy yield similar oscillator frequencies in both scenarios.

The oscillator frequency for bottom hadrons is fixed by considering the yield of $\Lambda_{b}(5620)$ and contributions from $\Sigma_{b}(5810)$ and $\Sigma_{b}^{*}(5830)$ decays in the statistical model and fitting the results to that in the coalescence model. That is,

[TABLE]

yielding $\omega_{b}=183$ MeV for RHIC in scenario 1 using the bottom quark mass $m_{b}=4700$ MeV. The oscillator frequencies for bottom hadrons for the other scenario at both RHIC and LHC are obtained using the same method, and they are summarized in Table 3.1. Again, the different total bottom quark number 0.44 at LHC 2.76 TeV energy and 0.71 at 5.02 TeV energy lead to almost same oscillator frequencies in both scenarios.

Using the values of $\omega_{q}$ for normal hadrons, it has been confirmed from Eq. (3.14) that the addition of an $s$-wave, a $p$-wave, or $d$-wave quark yields, respectively, a coalescence or suppression factor for the yields [15]

[TABLE]

for scenario 1 at RHIC. The coalescence factors for the other scenario are similar to those in Eq. (3.21). They show that the $d$-wave coalescence is more suppressed than the $p$-wave coalescence, which is further suppressed relative to the $s-$wave coalescence [6]. Since the production of multiquark hadrons involves more $s$-, $p$-, and $d$-wave coalescence factors, their yields are therefore generally suppressed.

3.2.2 Hadron coalescence

Since weakly bound hadronic molecules are expected to be continuously produced from the constituent hadrons and dissociated by interactions with other hadrons during the hadronic stage in heavy ion collisions, their yields in the coalescence model are determined at the end of the hadronic evolution at the kinetic freeze-out temperature $T_{F}$ and volume $V_{F}$. For the oscillator frequency in the wave function of a hadronic molecule, it can be determined from its relation to the mean square distance $\langle r^{2}\rangle$ between the two constituent hadrons. For a hadronic molecule in the relative $s$-wave state, the oscillator frequency is given by $\omega=3/(2\mu_{R}\langle{r^{2}}\rangle)$ with the reduced mass $\mu_{R}=m_{1}m_{2}/(m_{1}+m_{2})$. The mean square distance of the hadronic molecule is also related to its binding energy $B$ via the scattering length $a_{0}$ of the two interacting constituent hadrons, that is, $B\simeq\hbar^{2}/(2\mu_{R}a_{0}^{2})$ and $\langle{r^{2}}\rangle\simeq a_{0}^{2}/2$. These relations are valid when the binding energy is small and the scattering length is large compared to the range of the hadronic interaction. For a weakly bound two-body state, this leads to the simple relation $\omega=6B$. One notes that $\langle{r^{2}}\rangle$ is not the mean squared radius from the center of mass but is rather the mean square distance in the relative coordinate between the two hadrons. The oscillator frequency for $f_{0}(980)$, for instance, obtained from its binding energy of $B^{f_{0}}=M_{K^{+}}+M_{\bar{K}_{0}}-M_{f_{0}(980)}=493.7+497.6-980=11.3$ MeV, is $\omega_{f_{0}(980)}=6\times B^{f_{0}}=67.8$ MeV.

The oscillator frequencies for all hadronic molecules evaluated in the above described method are summarized in the tables for exotic hadrons. The number of constituent hadrons used in calculating the yield of hadronic molecules in the hadron coalescence are determined from the statistical model at the hadronization temperature and volume, and are given in Table 3.1.

Also shown in Table 3.1 are the temperature $T_{F}$ and volume $V_{F}$ of the hadronic matter at kinetic freeze-out, which are determined from requiring the yield of well-known hadronic molecules, such as the deuteron, from the hadron coalescence at kinetic freeze-out temperature and volume to be equal to that from the statistical hadronization model at hadronization temperature $T_{H}$ and volume $V_{H}$. Since the statistical hadronization model explains very well the yield of deuterons [247], it is necessary for the hadron coalescence model to explain also the deuteron production in order for it to be able to predict the production of exotic hadron molecules.

Because the temperature and volume are related in the isentropic expansion of the system in heavy ion collisions as already shown in Fig. 3.1, the kinetic freeze-out temperature and volume can thus be simultaneously determined by fitting the yield of deuteron from the hadron coalescence model to that in the statistical hadronization model. This is achieved by taking the oscillator frequency in the deuteron wave function to be $\omega_{d}=6\times B^{d}=6\times 2.2=13.2$ MeV and using the number of nucleons in Table 3.1. As shown in the table, the kinetic freeze-out at RHIC is found to take place at a higher temperature but a smaller system volume than those at LHC.

3.2.3 Heavy quark pair production

Charm quark pairs are produced through nucleon-nucleon binary collisions in relativistic heavy-ion collisions. Their numbers can be determined by using the Monte Carlo Glauber model to simulate nucleon-nucleon binary collisions [259, 260]. In a given binary collision, the probability to produce a charm quark pair is given by the ratio of the charm production cross section to its total inelastic scattering cross section. For the charm production cross section, it can be obtained from fitting to the experimental data at various collision energies [259]. The energy-momentum of each charm quark pair produced from a nucleon-nucleon collision can be obtained using the PYTHIA event generator [261]. Although heavy quark production in PYTHIA is based on the leading-order calculations in pQCD, albeit taking into account the effects due to initial and final parton showers, it can be tuned to reproduce the transverse momentum spectrum and rapidity distribution of charm quarks from Fixed-Order Next-to-Leading Logarithm (FONLL) calculations [262, 263]. As for the bottom quark, it can be simply determined from the ratio of the cross section for bottom production to that for charm production in the FONLL calculations.

An effect that needed to be considered for heavy quark production in heavy ion collisions is the modification of the parton distribution functions in heavy nucleus. For example, the parton distribution functions are known to decrease at small momentum fraction $x$ compared to those in a single nucleon. This so-called shadowing effect suppresses the production of heavy quark pair. Since the small-$x$ region in the parton distribution function can also contribute to heavy quark production as the collision energy increases, it is important to take into account the shadowing effect on heavy quark production in heavy ion collisions, and this can be included by using the EPS09 package [264]. It has been shown that recent experimental data on charm production at the LHC are much better described after including the shadowing effect [260, 265].

Using the above described method, the number of charm pairs produced at midrapidity in 0-10 % central collision at RHIC and LHC are given as follows. Without the shadowing effect, there are 4.5, 17, 23 pairs, respectively at 200 GeV, 2.76 TeV, 5.02 TeV. With shadowing, the numbers are 4.1, 11, 14 pairs, respectively, which are the ones used in Table 3.2. For the bottom quarks, they are estimated by using the ratio of the bottom production cross section to that of charm given by FONLL calculations, which are 0.75 %, 4 %, 5.1% at collision energies of 200 GeV, 2.76 TeV, 5.02 TeV, respectively.

3.3 Freeze-out conditions for molecular states

In this subsection, we will closely follow the discussion given in Ref. [245] on the freeze-out condition of a hadron from an expanding system without further elastic collisions.

Kinetic freeze-out of a particle of species $i$ from a matter occurs when its scattering time $\tau^{i}_{scatt}$ becomes larger than the expansion time of the system $\tau_{exp}$ [266]. The scattering time scale depends on the elastic scattering cross section with other particles as follows:

[TABLE]

with $\langle\sigma_{ij}v_{ij}\rangle$ being the thermal average of the product of the cross section times the relative velocity between particle species $i$ and $j$, and $n_{j}$ the density of particle $j$. The expansion time is defined as

[TABLE]

with $u$ being the expansion velocity of the system and can be approximated by the ratio of the fireball volume $V$ to its change in time, $V/(dV/dt)$.

For a spherically symmetric expanding fireball with radius $R$, the expanding time scale is reduced to $R/(3dR/dt)$ [267]. If the radius expands with a constant velocity, $R=vt$, one finds that $\tau_{exp}=\frac{1}{3}t$. Assuming for simplicity that the system is composed of only one species and that the cross section is independent of velocity, the freeze-out condition then becomes

[TABLE]

Although there is no general relation between $dR/dt$ and $\langle v\rangle$, particularly in the presence of a collective flow, the condition for the kinetic freeze out becomes simple if the rate of change in the radius is close to the average velocity of the particles in the system. In this case, it is simply given by [267]

[TABLE]

where the subscript ”$fo$” stands for physical quantities at kinetic freeze-out and $N$ is the total number of particles. It is seen that the two dimensional density determines the condition for freeze-out, and this is because the transverse total cross section determines whether a particle still interacts when it escapes from the medium.

On the other hand, the three dimensional density at the freeze-out goes as

[TABLE]

and it decreases with the square root of the total number of particles. This suggests that for higher collision energies and/or larger initial temperature and/or number of particles, the three dimensional density at which freeze-out takes places becomes smaller [245].

The above result is a general one not restricted to spherically symmetric expansion. For a system that follows the boost invariant Bjorken picture $R_{L}=c\tau$ with transverse expansion $R_{T}=v\tau$ of constant velocity $v$, one again finds that at large time, $\tau_{exp}=\frac{1}{3}\tau$ with $\tau$ being the invariant time, which has been explicitly confirmed in a hydrodynamical calculation [245]. In this case, using $V=\pi R_{T}^{2}R_{L}$ leads to the freeze-out condition

[TABLE]

This then leads to the following density at kinetic freeze-out:

[TABLE]

where $\langle v\rangle=(v^{2}c)^{1/3}$, with $c$ being the velocity of light, has been assumed for simplicity. As can be seen from Eq. (3.26) and Eq. (3.28), the relation between the freeze-out density and the cross section and/or the total number of particles seems to have a universal behaviour in three dimensions. The freeze-out conditions of the constituents in the hadronic matter for light nuclei or hadronic molecules, which are bound, plays an important role in determining their yields in the coalescence model.

For resonances with large decay width compared to the inverse lifetime of the hadronic phase, the freeze-out condition of its daughter particles will determine the decrease of its yield from the statistical model prediction. This can be understood by considering, for example, a simple rate equation for the $K^{*}$ meson during the hadronic stage,

[TABLE]

with $1/\tau^{K^{*}}_{scatt}=\sum_{i}\langle\sigma_{K^{*}i}v_{K^{*}i}\rangle n_{i}+\langle\Gamma_{K^{*}}\rangle$, and $1/\tau^{K}_{scatt}=\sum_{j}\langle\sigma_{Kj}v_{Kj}\rangle n_{j}$. Here $i$ and $j$ stand for mostly the light mesons such as the pion and $\rho$ meson, i.e., $1/\tau^{K^{*}}_{scatt}=\langle\sigma_{K^{*}\rho\to K\pi}v_{K^{*}\rho}\rangle n_{\rho}+\langle\sigma_{K^{*}\pi\to K\rho}v_{K^{*}\pi}\rangle n_{\pi}+\langle\Gamma_{K^{*}}\rangle$ and $1/\tau^{K}_{scatt}=\langle\sigma_{K\pi\to K^{*}\rho}v_{K\pi}\rangle n_{\pi}+\langle\sigma_{K\rho\to K^{*}\pi}v_{K\rho}\rangle n_{\rho}+\langle\sigma_{K\pi\to K^{*}}v_{K\pi}\rangle n_{\pi}$ with $\langle\Gamma_{K^{*}}\rangle$ being the thermally averaged decay width of the $K^{*}$ meson [268]. In the above, non-linear terms originated from the interaction between $K^{*}$ mesons or kaons, like $K\bar{K}\to\rho\pi$, are not considered.

Consider a simple picture where the total number of light mesons and $K$ mesons are fixed as the system expands. The equilibrium number of $K^{*}$ mesons is given by the asymptotic value obtained by taking the right hand side of Eq. (3.29) to be zero, given as

[TABLE]

At chemical freeze-out, this value should correspond to that given by the statistical model. As the system expands, while the total number of light hadrons and $K$ meson remain fixed, the $K^{*}$ number decreases due to decay as the freeze-out volume $V(\tau)$ increases, leading to a suppression factor that depends on the freeze-out volume, a result borne out in the measured $K^{*}$ number in heavy ion collision [5]. This mechanism becomes relevant only for particles that have natural decay width, which leads to terms in the rate equation that are proportional to their numbers and thus scale with the volume of the system. For bound states composed of hadrons, they do not have natural decay widths and are thus not affected by this suppression mechanism. Although the kinetic freeze-out condition for hadrons depend on their elastic scattering cross sections, a universal kinetic freeze-out temperature is used in Ref. [245], and it is determined by requiring the deuteron yield from the coalescence model at this temperature to reproduce the experimental value, which has been found to follow the statistical model prediction at the chemical freeze-out point.

3.4 Yields of hadrons

This section summarizes the expected yields of exotic hadrons in central Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV at RHIC, central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV at LHC, and central Pb+Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV at LHC. We show the results for normal quark ($\bar{q}q$, $3q$), multiquark, and hadronic molecule configurations, calculated from the coalescence model in addition to those estimated from the statistical model. These results are shown in Tables  3.3,  3.4,  3.5 and  3.6. We also give some discussions on the obtained results.

In Refs. [14, 15], it was found that for most of the hadronic states, the yield from the coalescence model for the compact multiquark state is smaller than that for the usual quark configuration as a result of the suppression due to the coalescence of additional quarks. For the same state, the yield from the coalescence model for a molecular configuration is larger than that from the statistical model prediction. This is in contrast to high energy pp collisions, where molecular configurations with small binding energies are difficult to produce at high transverse momentum $p_{T}$ [144]. The upper panel of Fig. 3.2, shows the coalescence-statistical yield ratio, $R_{cs}=N_{\mathrm{coal}}/N_{\mathrm{stat}}$, given in Refs. [14, 15] using parameters given in Table 3.1 and assuming that the hadron coalescence takes place at $T_{F}=125~{}\mathrm{MeV}$ as well as including the resonance decay contributions to the $K(\bar{K})$ and $N$ yields. In this treatment, however, the coalescence model overestimates the deuteron yield, which is known to follow the statistical model prediction at the chemical freeze-out temperature.

In Fig. 3.3, we show the freeze-out temperature dependence of the coalescence-statistical yield ratio $R_{h}^{CS}$ for deuteron and $\Lambda(1405)$ at RHIC. Requiring $R_{h}^{CS}=1$ for the deuteron leads to a freeze-out temperature of deuteron $T_{F}=119~{}\mathrm{MeV}$. According to Eq. (3.28), the density at which a particle freezes out is inversely proportional to its scattering cross section with other particles in the medium. Since the elastic cross section of a particle is related to its size, one would expect that the freeze-out density for a particle decreases as its size increases. This result suggests that the freeze-out temperature will be smaller for particles of larger size. From the relation between the radius and the oscillator frequency $\omega$ in the wave function of a hadron, the root mean square radii for $\Lambda(1405),d,X(3872)$ are found to be 1.71, 1.77 and 2.36 fm, respectively. Because of their different radii and thus sizes, these particles are expected to freeze-out at different temperatures. Such subtleties are, however, neglected in previous studies, and using a common freeze-out temperature seems to still give results for hadronic molecules that are roughly consistent with the statistical model results, i.e., $0.2<R_{h}^{CS}<2$.

In the lower panel of Fig. 3.2, we show the updated results in the present treatment described in this review. Similar to the old results, the yield from the coalescence model for a compact multiquark state is generally suppressed due to the coalescence of additional quarks, as discussed above and in Refs. [14, 15]. The yield from the coalescence model for a molecular configuration strongly depends on its size. Loosely bound hadronic molecules are more easily produced, and tightly bound molecules have smaller size and their production is suppressed. Thus, the yield of a hadronic molecule in heavy ion collisions can be used as a measure of its spatial size.

4 Coalescence model for resonances

4.1 Coalescence model

One of the commonly used tools for calculating production probabilities of composite particles in high energy collisions is the coalescence model [14, 15, 269]. It is based on a sudden approximation and gives the production probability of a composite particle by the probability of finding it in the particle source formed by the collision just before the freeze-out stage. For a composite particle, $\gamma$, which is a bound state, $|\gamma\rangle\!\rangle$, of two particles, $a$ and $b$, the probability is given by

[TABLE]

where $\hat{\rho}$ is the density matrix for the two particles in the source.

Most of the composite particles produced in high energy collisions are, however, not stable bound states but resonances with non-negligible widths. When the particle $\gamma$ is a resonance due to the interaction between $a$ and $b$, a straightforward extension of Eq. (4.1) would be to replace the bound state wave function by the resonance wave function with the appropriate modification of the complex conjugation [45]. The calculated probability, $P_{\gamma}$, would be complex, however, and its meaning becomes unclear in the case of broad resonances as we will see later in the numerical examples. Experimentally, the resonance is observed as a peak in the invariant mass spectrum for the $a$ and $b$ scattering system.

In this section, the formulation for the resonance particle production in the way it is experimentally observed is discussed.

4.2 Model for S-wave resonance

To make the discussion clear, we consider here a simple (Lee type) model of a S-wave resonance [270, 44], where a particle $c$ is coupled to two particles, $a$ and $b$, giving rise to a resonance in the two particle scattering system. In the center-of-mass system, the Hamiltonian is given by

[TABLE]

where $H_{0}$ describes the system without the coupling which is given by $V$. The free two particle state with relative momentum ${\bf{k}}$ is denoted by $|{\bf{k}}\rangle$, while the one particle state is denoted by $|c\rangle$, and they are eigenstates of $H_{0}$ with eigenvalues $E_{k}$ and $E_{c}$, respectively, i.e.

[TABLE]

The non-zero matrix elements of $V$ are expressed as

[TABLE]

where $g$ and $v(k)$ are assumed to be real. The one-particle state, $|c\rangle$, together with the two-particle states, $|{\bf{k}}\rangle$, constitutes a complete set of the model space.

[TABLE]

The scattering state with the asymptotic relative momentum, ${\bf{p}}$, can be expressed as

[TABLE]

where $E_{p}^{\pm}=E_{p}\pm i\eta$ specifies the asymptotic boundary condition and the T-matrix, $T(E)$, satisfies the Lippmann Schwinger equation,

[TABLE]

for generally complex $E$. If the system has no bound state, the scattering states $|{\bf{p}}\rangle\!\rangle^{\pm}$ give the complete set and thus

[TABLE]

The Lippmann-Schwinger equation for $T$, Eq. (4.7), can be easily solved and the relevant matrix elements are given by

[TABLE]

where $\Sigma(E)$ is the self-energy of the particle $c$,

[TABLE]

Eq. (4.6) then gives the scattering state, $|{\bf{p}}\rangle\!\rangle^{\pm}$, as

[TABLE]

A resonance state, $|r\rangle\!\rangle$, is an eigenstate of $H$ with a complex eigenvalue, $\mathcal{E}_{r}$, and satisfies

[TABLE]

leading to the eigenvalue equation

[TABLE]

which confirms that $\mathcal{E}_{r}$ is a pole of T-matrix in the complex $E$ plane. In this model, the conjugate state, $|\tilde{r}\rangle\!\rangle$, which is an eigenstate of $H$ with the eigenvalue $\mathcal{E}_{r}^{*}$, is related to $|r\rangle\!\rangle$ by complex conjugation, i.e.

[TABLE]

The normalization condition becomes

[TABLE]

and determines $\langle c|r\rangle\!\rangle^{2}$ as

[TABLE]

The T-matrix can now be decomposed into the resonance and non-resonance terms as

[TABLE]

Depending on the boundary condition specified by $E_{p}^{\pm}$, one can choose the resonance wave function $|r\rangle\!\rangle$ satisfying the corresponding boundary condition with ${\rm{Im}}\mathcal{E}_{r}=\mp\Gamma_{r}/2$ (half width of the resonance). The resonance term of the T-matrix introduced in Eq. (4.20) has the matrix elements

[TABLE]

4.3 Coalescence model for scattering states

A straightforward extension of Eq. (4.1) to scattering states gives the probability, $P({\bf{p}})$, of finding the two particles, $a$ and $b$, with the relative momentum ${\bf{p}}$ as

[TABLE]

Using Eq. (4.6), one can decompose it into the free (background) term, $P^{(0)}$, the interaction term, $P^{(1)}$, and the interference term, $P^{(2)}$, as

[TABLE]

Since $P({\bf{p}})$ gives the invariant mass spectrum for the $(a,b)$ pairs, one expects the resonance structures to appear in $P^{(1)}({\bf{p}})$ and $P^{(2)}({\bf{p}})$. In order to get information on the resonance state from the production processes, one decomposes these quantities further into resonance and non-resonance parts, using the decomposition of the T-matrix, Eq. (4.20).

[TABLE]

where Eq. (4.14) is used to get Eqs. (4.29) and (4.30). It is seen that the resonance parts, $P^{(1),r}({\bf{p}})$ and $P^{(2),r}({\bf{p}})$, carry the information on the properties of the resonance through its energy, $\mathcal{E}_{r}$, and wave functions, $|r\rangle\!\rangle$ and $|\tilde{r}\rangle\!\rangle$. The main problem is how to extract these quantities from the experimentally observed invariant mass spectrum given by $P({\bf{p}})$. As to be discussed later with numerical examples, it might be possible to obtain $P^{(0)}({\bf{p}})$, $P^{(1)}({\bf{p}})$ and $P^{(2)}({\bf{p}})$ separately through their angular distributions depending on the nature of the density matrix $\hat{\rho}$ describing the source. Once the interaction term, $P^{(1)}({\bf{p}})$, and the interference term, $P^{(2)}({\bf{p}})$, are separately obtained, the further decomposition of them into the resonance and non-resonance parts can in principle be done through their $p=|{\bf{p}}|$ dependences, since the non-resonance parts have smooth energy dependence and are also small in magnitude compared with the resonance parts in the relevant region.

As for the density matrix $\hat{\rho}$ describing the source, one can assume that it has no matrix elements between $|c\rangle$ and $|{\bf{k}}\rangle$ and denote its non-zero matrix elements as

[TABLE]

$\rho_{c}$ then contributes only to $P^{(1)}$ in the decomposition (4.24) of $P$ and one has

[TABLE]

and

[TABLE]

Applying Eq. (4.29) to calculate the resonance part, one notes that the contribution of $\rho_{c}$ contains divergences at $E_{p}=E_{c}$, which are of course cancelled by the corresponding divergences in the non-resonance part. In order to avoid this problem and get a more reasonable resonance part, one modifies Eq. (4.21) by replacing $\mathcal{E}_{r}$ with $E$ in the numerator, i.e.

[TABLE]

so as to eliminate the divergences without changing the residue at the pole. One then obtains for the resonance parts

[TABLE]

where $F^{(1)}$ and $F^{(2)}$ are defined by Eqs. (4.37) and (4.40), respectively.

It is noticed here that the completeness of $|{\bf{p}}\rangle\!\rangle^{-}$, Eq. (4.5), leads to the following sum rules for the integrated probabilities denoted by $\Pi$s.

[TABLE]

where

[TABLE]

obtained from Eqs. (4.23) and (4.25) have been used. The sum rules imply that the number of $(a,b)$ pairs is not affected by their mutual interactions but is increased by the decay of $c$.

4.4 Numerical examples

To see how the resonance appears in the probability, $P({\bf{p}})$, which gives the invariant mass spectrum of the $(a,b)$ pairs, calculated by the above formalism, some numerical examples are presented below.

With the non-relativistic kinetic energy and a monopole form factor, i.e.

[TABLE]

$\Sigma(E)$ is then given by

[TABLE]

where $\lambda=\frac{m^{2}g^{2}}{2\pi\mu}$ and $p_{E}=\pm\sqrt{2mE}$ with the appropriately chosen sign. For example, $p_{E_{p}^{\pm}}=\pm\sqrt{2mE_{p}}=\pm p$ for real positive $E_{p}$ and $p$. The full off-shell T-matrix becomes

[TABLE]

where $p_{c}^{2}=2mE_{c}$. The poles of $T(E)$ are given by the solutions of a 4th order equation for the variable $p_{E}$, and one chooses one of them with positive real part and negative imaginary part for the resonance pole which will be denoted as $p_{r}$, i.e. $p_{r}^{2}=2m\mathcal{E}_{r}$. In this notation, one has

[TABLE]

and $\langle c|r\rangle\!\rangle^{2}$ given by Eq. (4.19) becomes

[TABLE]

The density matrix given by Eq. (3.8) and used in the ExHIC papers [14, 15] is in the present notation

[TABLE]

where $k_{z}$ and ${\bf{k}}_{T}$ are the longitudinal and transverse components of ${\bf{k}}$, respectively, and the normalization $N$ is determined by the condition,

[TABLE]

where $(2\pi)^{3}\delta({\bf{0}})$ has been replaced by the volume of the emission source $V$ and $n_{ab}$ is the number of pairs of the particles $a$ and $b$. Inserting Eq. (4.52) into Eqs. (4.37) and (4.40), however, one sees that $O(1/\eta)$ terms appear in both equations. The first delta function in Eq. (4.52), which gives rise to this problem, is based on the assumption that the size of the emission source is much larger than the sizes of hadrons. Although the assumption is justified for the bound state formation, it is not the case for scattering states. One therefore chooses a Gaussian distribution in the transverse direction, i.e.

[TABLE]

and the normalization condition (4.53) determines $\tilde{N}$ as

[TABLE]

where $2\pi\delta(0)$ is replaced by the longitudinal length $L$ of the source. Since for large $\alpha$, $\exp(-\alpha({\bf{k}}_{T}-{\bf{k}}^{\prime}_{T})^{2})\rightarrow\frac{1}{4\pi\alpha}(2\pi)^{2}\delta({\bf{k}}_{T}-{\bf{k}}^{\prime}_{T})$, $4\pi\alpha L$ corresponds to the volume $V$ of the source in Eq. (4.53). Substituting Eq. (4.54) into Eqs. (4.37) and (4.40), and carrying out the angular parts of the integrals, one obtains

[TABLE]

where $\alpha_{\pm}=\alpha\pm\frac{\beta}{8m}$ and the variable $k_{T}$ has been changed to $k$ since they are the same for $k_{z}=0$. It is seen that the two terms have very different angular dependences. $P^{(1)}_{ab}$ is isotropic reflecting the S-wave nature of the interaction, while the $P^{(2)}$ contributes only on the $p_{z}=0$ plane due to the strong anisotropy of the density matrix (Eq. (4.54)). $P^{(1)}_{c}({\bf{p}})$ given by Eq. (4.35) is

[TABLE]

and has the isotropic angular dependence as $P_{ab}^{(1)}$. Thus, in the direction of $p_{z}=0$, the interference term, $P^{(2)}$ is dominant while in the other direction, the probability is given by the interaction term, $P^{(1)}=P_{ab}^{(1)}+P^{(1)}_{c}$. The resonance parts can be similarly calculated by Eqs. (4.42), (4.43) and (4.44).

For numerical examples, three sets of interaction parameters are considered with the set A generating a resonance similar to $\Lambda(1405)$ while the sets B and C giving typical examples of broad and narrow resonances, respectively. Thus, $(a,b)$ is $(\pi,\Sigma)$ and $c$ is the $\bar{K}N$ bound state with its coupling to $\pi\Sigma$ switched off. The interaction parameters are given in Table 4.1 and the resulting properties of the resonances (pole positions in the complex $E$ and $p$ planes, $\mathcal{E}_{r}$ and $p_{r}$, and the square of the overlap between the resonance state and the one-particle state, $\langle c|r\rangle\!\rangle^{2}$ (Eq. (4.51))) are given in Table 4.2. The square of the overlap for a bound state would represents the probability of the one particle state to remain in the bound state and therefore would be a real positive number less than or equal to 1. In the case of a resonance state, however, the probability interpretation is not applicable and the magnitude of deviation from 1 gives a measure of the contribution from $(a,b)$ scattering states in forming the resonance.

The density matrix parameters are chosen to be similar to those used in 3.4 with the parameters given in Table 3.1 for RHIC [14, 15] and are given in Table 4.3. Thus the temperature $1/\beta$ and the volume $V=4\pi\alpha L$ are the freeze-out temperature $T_{F}(0.119{\rm{GeV}})$ and volume $V_{F}(20355{\rm{fm}}^{3})$, respectively, and $n_{\pi\Sigma}$ is estimated by the statistical model at the hadronization temperature, $T_{H}(0.162{\rm{GeV}})$, and volume, $V_{H}(2100{\rm{fm}}^{3})$. $\rho_{c}$ is calculated by the coalescence model with the $\bar{K}N$ density matrix chosen analogously to that for $\pi\Sigma$ and the bound state wave function taken to be that of the harmonic oscillator ground state with the oscillator frequency $\omega=20.5{\rm{MeV}}$. To see the relative importance of the three terms, one takes here $n_{\pi\Sigma}=1$ and show only the ratio $n_{\bar{K}N}/n_{\pi\Sigma}$ and the resulting $\rho_{c}$.

Shown in the left panels of Fig.4.1 are the two terms, $P_{ab}^{(1)}$ and $P_{c}^{(1)}$, which have the same isotropic angular dependence, and the sum $P^{(1)}=P_{ab}^{(1)}+P_{c}^{(1)}$, multiplied by the phase space factor $p^{2}/2\pi^{2}$ as functions of $p$, for the interaction parameter sets A , B and C.

One sees that $P_{ab}^{(1)}$ and $P_{c}^{(1)}$ give similar but slightly different spectra, the peak position of the former being visibly shifted downward in the cases of the sets A ($\Lambda(1405)$ ) and B (broad resonance). The integrated values of the former, $\Pi_{ab}^{(1)}$, are $0.0017$, $0.0041$ and $0.00031$ for the sets A, B and C, respectively, to be compared with that of the latter, $\Pi_{c}^{(1)}$, which is $\rho_{c}=0.0025$ (cf. Eq. (4.46)). There is a strong correlation between these integrated probabilities and the widths of the resonances ($-2{\rm{Im}}\mathcal{E}_{r}$). They are also correlated with the contributions of the $(a,b)$ scattering states to the resonance wave function expressed by $|\langle c|r\rangle\!\rangle^{2}-1|$ (see Table 4.2).

On the other hand, the interference term, $P^{(2)}$, has the same angular dependence as the background term, $P^{(0)}$, which is larger by nearly two orders of magnitude, and thus can be observed only as a tiny structure of the spectrum in the direction of $p_{z}=0$. The right panels of Fig.4.1 show this term multiplied by the phase space factor, $p/2\pi$, for the sets A, B and C, respectively, together with the background term divided by 50 for the illustrative purpose. The interference term is seen to have an oscillatory behavior near the resonance giving rise to a tiny but peculiar structure in the observed spectrum. Since the background term is assumed to be known, its subtraction is in principle possible. According to Eq. (4.45), the net (integrated) contribution of the interference term, $\Pi^{(2)}$, is negative and cancels that of the interaction term, $\Pi^{(1)}_{ab}$.

As for the resonance and non-resonance parts, the latter for the interaction terms, $P_{ab}^{(1)nr}$ and $P_{c}^{(1)}$, are generally small in this model and are visible in the invariant mass spectra only in the case of the set B. The structure of the resonance can thus be learnt directly from the observed spectra. On the other hand, the non-resonance part for the interference term, $P^{(2)nr}$, is not small. It is therefore more complicated to analyze the resonance properties through this term, though, as discussed before, its extraction from the observed spectra is difficult anyway. The left panel of Fig. 4.2 shows the interaction terms, $P_{ab}^{(1),r}$ (red), $P_{ab}^{(1),nr}$ (orange). $P_{c}^{(1)r}$ (blue) and $P_{c}^{(1)nr}$ (green), while the right panel of Fig. 4.2 shows the interference terms, $P^{(2),r}$ (red) and $P^{(2),nr}$ (orange) for the set B.

The contributions of the non-resonance parts can be seen more quantitatively from their contributions to the integrated probabilities, $\Pi^{nr}$’s, together with $\Pi^{r}$’s in Table 4.4. They are divided by $\rho_{c}$, which is the total number of $(a,b)$ pairs produced by the interaction in the present model. One sees that for the interference terms, the resonance parts, $\Pi^{(2)r}$ is larger in magnitude than $\Pi^{(2)}$ and thus the non-resonance parts have opposite signs. Note that the resonance parts need not satisfy the sum rules Eqs. (4.45) and (4.46). In the case of set B, however, $\Pi^{(2)r}$ is negative and too large in magnitude and makes the summed resonance part, $\Pi^{r}=\Pi_{ab}^{(1)r}+\Pi_{c}^{(1)r}+\Pi^{(2)r}$ negative. Thus the concept of “integrated (or total) resonance production probability” becomes unclear for broad resonances.

Once the resonance parts are obtained from the observed spectra, they are expressed by the resonance wave function, $|r\rangle\!\rangle$, the non-resonance T-matrix, $T^{nr}$, and the density matrix, $\hat{\rho}$ as given in Eqs. (4.42), (4.43) and (4.44). The structure of the resonance in the present model is determined by $|{\bf{k}}\rangle$ and $|c\rangle$ components of its wave function, $\langle{\bf{k}}|r\rangle\!\rangle$ and $\langle c|r\rangle\!\rangle$, and $\rho_{c}$ which reflects the structure of the particle $c$ representing the $\bar{K}N$ bound state with the coupling to the $\pi\Sigma$ channel switched off, and they can be extracted by the analysis described above.

With the resonance wave function, $|r\rangle\!\rangle$, and the density matrix given by Eq. (4.33), the straightforward extension of Eq. (4.1) gives the production probability in the coalescence calculation as

[TABLE]

where Eq. (4.15) has been used to get the last expression. Note that these quantities are all complex and thus their physical meaning is not clear unless the real parts are dominant. As discussed in the previous subsections, what can be extracted from the observed spectra to study the structure of the resonance are the resonance parts, $P^{(1)r}$ and $P^{(2)r}$, defined by Eqs. (4.29) and (4.30). In the present model, they can be calculated by Eqs. (4.42), (4.43) and (4.44), with the decomposition, $P^{(1)r}=P_{c}^{(1)r}+P_{ab}^{(1)r}$.

The production probabilities calculated by the above expressions using the resonance wave function are compared with the integrated resonance parts, $\Pi_{c}^{(1)r}$, $\Pi_{ab}^{r}=\Pi_{ab}^{(1)r}+\Pi^{(2)r}$ and $\Pi^{r}$ in Table 4.5. It is seen that the correlation between $\tilde{\Pi^{r}}$’s and $\Pi^{r}$’s is very good for the narrow resonance, tolerable for $\Lambda(1405)$ but very poor for the broad resonance.

4.5 Summary and discussions

The coalescence model for resonance particle production in high energy heavy ion collisions is formulated in the way it is observed in experiments. A simple S-wave (Lee type) model where two particles, $a$ and $b$, couple to a particle, $c$, forming a resonance in the $(a,b)$ scattering states, is used to clarify the discussion. The probability, $P({\bf{p}})$, of finding the two particles with a relative momentum ${\bf{p}}$, which is simply related to the invariant mass spectrum of the two particle system, is calculated in the coalescence model for scattering states. It consists of the background term, $P^{(0)}$, the interaction term, $P^{(1)}$, and the interference term, $P^{(2)}$ (see Eq. (4.24)). In this model, $P^{(1)}$ can be decomposed into the contribution of the two particle states, $P_{ab}^{(1)}$, and that of the $c$ particle, $P_{c}^{(1)}$, and the former together with the interference term satisfy the sum rule, Eq. (4.45), while the latter satisfies Eq. (4.46). Thus the integrated contribution of the interaction and interference terms gives $\rho_{c}$, the probability of finding the particle $c$ in the source. In the numerical examples, one has in mind $\Lambda(1405)$ as the resonance and $a,b$ is $\pi\Sigma$. $\rho_{c}$ is calculated assuming $c$ to be a $\bar{K}N$ bound state with its coupling to $\pi\Sigma$ states switched off. This is the procedure used in 3.4 and the previous ExHIC papers to estimate the production yield of $\Lambda(1405)$. In the case of the interaction parameters corresponding to $\Lambda(1405)$ (set A) in this example, the non-resonance parts are small and the integrated sum of $P^{(1)}$ and $P^{(2)}$ can be regarded as the production probability of the resonance $\Lambda(1405)$. Thus the ExHIC procedure is justified in this example. $P^{(1)}$ and $P^{(2)}$ can be themselves considered as the resonance terms and tell us further details of the resonance such as its $\pi\Sigma$ components.

In a more general case where the non-resonance parts are not small, another step of separating the resonance parts becomes necessary. This is the case in the example of a broad resonance (set B), where the non-resonance part is about the same in magnitude as the resonance part in the interference term (see Fig. 4.1 and Table 4.4) and the integrated resonance part, $\Pi^{r}$, becomes negative. Even in such a case, however, the non-resonance part is still small (about $10\%$) in the interaction term and the information on the structure of the resonance can be extracted through this term. It is then crucial to get the interaction and interference terms separately from the observed invariant mass spectrum. This is in principle possible for anisotropic sources as seen in the discussed examples.

Although the coalescence calculation for bound states can be straightforwardly extended to resonances, it leads however, to complex values for the production probabilities (see Tables  4.5). In this case, the calculated probabilities have quantitative meaning for narrow resonances, but this becomes doubtful for broad ones.

It is straightforward to extend the above formulation to more general cases of multichannels with non-separable interactions, though the expressions for the production probabilities become more involved and the numerical calculations based on them would be time-consuming. Such extensions are necessary to use realistic models of resonances and study how their structures are reflected in the production probabilities. For quantitative analyses of the experimental data, the applicability of the coalescence model itself should also be examined in detail.

5 Hadron-hadron interactions from two particle momentum correlations

Pairwise hadronic interactions as well as quantum statistics produce a correlation at low relative momenta in multiparticle production from elementary to heavy ion collisions [25, 26, 27, 28, 29]. The momentum correlation of identical particles from quantum statistics, known as the Hanbury-Brown and Twiss (HBT) [271] or Goldhaber-Goldhaber-Lee-Pais (GGLP) effect [272], can give information on the size of the emission source through the (anti-)symmetrization of the two-boson (fermion) wave function. The quantum statistical effect of stable hadrons, particularly pions, has been used to estimate the source sizes created in relativistic nucleus-nucleus collisions [273, 29]. By comparison, one expects substantial interaction effects on the correlation function for particle pairs whose interaction is sufficiently strong in the range comparable to the effective source size [26, 27]. In particular, the correlation function of non-identical pairs is directly related to the pairwise interaction due to the absence of the quantum statistical effect [27]. Thus, high statistics measurement of the correlation function might provide information on the pairwise interaction of any measurable channel, including those difficult to perform the scattering experiments. In this section, a brief review is given on recent activities on constraining hadron-hadron interactions through momentum correlations in heavy-ion collisions and their implications for the interpretation and the possible existence of exotic states.

5.1 General property of the two-particle momentum correlation function

5.1.1 Formalism

The two-particle momentum correlation function is defined as the ratio of the two-particle spectrum to the product of single particle inclusive momentum spectra [25, 26, 27, 29],

[TABLE]

where $P$ and $q$ are the center-of-mass and the relative momentum of the pair, respectively, and $E_{i}$ ($i=1,2$) is the energy of the hadron $i$. In the last equality in Eq. (5.2), the relative momentum is expressed in the center-of-mass frame of the pair (the pair rest frame), where $E^{\prime}_{i}(i=1,2)$ is the energy of the hadron $i$ in this frame and $M_{\mathrm{inv}}=E^{\prime}_{1}+E^{\prime}_{2}$ is the invariant mass. In the non-relativistic limit, $E^{\prime}_{1}\to M_{1}$ and $E^{\prime}_{2}\to M_{2}$, the definition of the relative momentum reads $q=(M_{2}p_{1}-M_{1}p_{2})/(M_{1}+M_{2})$.

Assuming independent (chaotic) emission from the source, i.e., particles are produced with random phases, and the correlation function can be expressed in terms of the single particle source function $S(x_{i},\bm{p_{i}})$, which describes the emission probability from a space-time point $x_{i}$ with momentum $\bm{p_{i}}$, and the weight factor $|\varphi^{(-)}(\bm{r},\bm{q})|^{2}$, which depends on the relative coordinate $\bm{r}$ and momentum $\bm{q}$,

[TABLE]

The weight factor $|\varphi^{(-)}(\bm{r},\bm{q})|^{2}$ can be identified as the relative wave function of the pairs in the outgoing state, provided that the difference between the emission times of the two particles is small [26]. In general, the emission time difference modifies the relative coordinate $\bm{r}$ in the relative wave function from the position difference $\bm{x}_{1}-\bm{x}_{2}$. By using $P$ and $q$, particle momenta are given as $p_{1}=E^{\prime}_{1}P/M_{\mathrm{inv}}+q$ and $p_{2}=E^{\prime}_{2}P/M_{\mathrm{inv}}-q$. Then the free two-particle wave function is given as

[TABLE]

When the interaction between the two particles is switched on, the relative wave function is modified into a superposition of $\exp(i\bm{q}\cdot\bm{r})$, which is a function of $\bm{r}$, while the center-of-mass wave function $\exp(-iP\cdot X)$ is kept unchanged. The modification of the relative coordinate in Eq. (5.5) serves as a generalization of the formula derived in Refs. [274, 27].

The above formula can be reduced to a convenient form in the pair rest frame where $P$ and $q$ becomes temporal and spatial four-vectors, $P=(M_{\mathrm{inv}},\bm{0})$ and $q=(0,\bm{q})$, respectively. In this case, the center-of-mass coordinate $X$ and relative time $t$ can be integrated out to obtain the relative source function

[TABLE]

and the correlation function is then given by the Koonin-Pratt (KP) formula

[TABLE]

In frames different from the pair rest frame, one needs to put $x=(t,\bm{r}+\bm{v}t)$ in Eq. (5.6). In the KP formula, the relative source function $S_{12}(\bm{r})$ can be interpreted as the relative source distribution integrated over time in the pair rest frame. In particular, for the free (anti-)symmetrized wave function, this formula reduces to the three-dimensional Fourier transformation of the source function. It should be noted that the relative source function $S_{12}(\bm{r})$ depend on $P$ when the emission point and momentum are correlated. In the later discussions, the KP formula, Eq. (5.7), is mainly used. While there are several conditions under which the KP formula works [275, 29], they seem to be satisfied in high-energy heavy-ion collisions.

5.1.2 Correlations from strong interactions and quantum statistics

Both the pairwise interaction and the quantum statistics modify the relative wave function from the simple plain wave $e^{i\bm{q\cdot r}}$. Since the strong interaction is of short range, the modification of the relative wave function appears mainly in the $s$-wave. In this case, one can write the relative wave function in the two-body outgoing state with an asymptotic relative momentum $\bm{q}$ as

[TABLE]

where $q=|\bm{q}|$, $j_{0}$ is the spherical Bessel function and $\psi(r)$ is the relative wave function in the $s$-wave, which is regular at $r\to 0$ and has an asymptotic form,

[TABLE]

with $\delta$ being the phase shift. It should be noted that the above wave function $\psi$ is different from that appearing in the two-particle scattering by a factor $e^{2i\delta}$. In the two-body outgoing state, the coefficient of the outgoing wave is unity and the incoming spherical wave is modified in contrast to the scattering of two particles where the coefficient of the incoming wave is unity and the outgoing spherical wave is modified.

For illustration, let us consider a spherical and static Gaussian source, $S_{i}(x_{i},\bm{p}_{i})\propto\delta(t_{i}-t_{0})\,\exp(-\bm{x}_{i}^{2}/2R_{i}^{2})$ and simultaneous emission of the pairs. Then the correlation function from the non-symmetrized wave function $\varphi^{(-)}$ is obtained as

[TABLE]

where

[TABLE]

Since the sum of first two terms in Eq. (5.8) does not contain the $s$-wave components, the cross term involving the third term disappears for a spherical source. As a result, the effect of interaction on $C(\bm{q})$ appear as the deviation from unity by the difference of squared wave functions between the free and $s-$wave, as seen from Eq. (5.12).

For identical pairs, the relative wave function is symmetric or antisymmetric with respect to the exchange of the two-particle spatial coordinates ($\bm{r}\to\,-\bm{r}$). For spin-$1/2$ pairs, the spatial part of the wave function is symmetric for the spin-singlet (${}^{1}S_{0}$) state and antisymmetric for the spin-triplet (${}^{3}S_{1}$) state, which does not have $s-$wave interaction. Then the wave function is given as

[TABLE]

The wave functions $\varphi_{E}$ and $\varphi_{O}$ have even and odd parities, respectively.

For (anti)symmetrized wave functions, the KP equation is reduced to the Fourier transform of the source function, then the correlation function is given by,

[TABLE]

where $\widetilde{S}$ denotes the Fourier transform, $\widetilde{S}(2\bm{q})=\int d\bm{r}S(\bm{r})\exp(-2i\bm{q}\cdot\bm{r})$. The second equality in Eqs (5.17) and (5.19) is obtained for the Gaussian source (5.13). The Gaussian term in Eqs. (5.17) and (5.19) represents the effects from the quantum statistics. For the symmetric (asymmetric) wave function, the correlation function exhibits enhancement (reduction) from unity and its width in $q$ is inversely proportional to the size of the source. In the case of identical interacting particles, Eq. (5.16), the effect of the interaction appears as deviation not from unity but from the free correlation function. This fact provides an intuitive understanding of the correlation function as follows [36, 37];

For large $q$, the wave function rapidly oscillates to give $\Delta C(\bm{q})\simeq 0$. Thus one needs to look at small $q$ to get information on the interaction.

• 2.

Weakly attractive interaction gives $|\psi(r)|>j_{0}(qr)$ in the range of the interaction and thus leads to small enhancement of the correlation function.

• 3.

Strongly attractive interaction having an bound state gives a node to $\psi(r)$. Since the contribution from the integrand with $r\simeq 0$ is suppressed by $r^{2}$ in $d\bm{r}$, the correlation function is also suppressed due to $|\psi(r)|<|j_{0}(qr)|$. Repulsive interaction also leads to the similar behavior.

In reality, the hadron-hadron correlation function is expressed as combinations of the above correlation functions. Considering the spherical Gaussian source and neglecting the Coulomb potential and channel coupling effects, one may have the following classification for non-identical spinless meson pairs, identical spinless meson pairs, pairs of a spinless meson and a spin-half baryon, non-identical spin-half baryon pairs, and identical spin-half baryon pairs;

[TABLE]

Interaction generally depends on the spin of the pair, and so does the interaction dependent part of the correlation function, $\Delta C(\bm{q})$, as found in the $BB^{\prime}$ pair, Eq. (5.23). The correlation function of the spin-half baryon pairs is obtained as the spin-average over the spin-singlet and the triplet states. As a result, the correlation at $\bm{q}=0$ is not zero like Eq. (5.19) but $1/2$ for the non-interacting case ($\Delta C(\bm{q})=0$).

5.1.3 Lednický and Lyuboshits Model

In order to examine the interaction dependence of the correlation function, an analytic model developed by Lednický and Lyuboshits (LL) [26] is useful. In the LL model, the correlation function is obtained by using the asymptotic wave function together with the shape-independent approximation in the scattering phase shift. Then the correlation function is given in terms of the scattering amplitude and the effective range.

The asymptotic wave function Eq. (5.10) can be rewritten in the following form,

[TABLE]

where $f(q)=(\mathcal{S}-1)/2iq$ is the scattering amplitude and $S=e^{2i\delta}$ is the S-matrix. With the Gaussian source (5.13), the integral in the KP formula for $\psi_{\text{asy}}$ is reduced to

[TABLE]

where $x=2qR$, $F_{1}(x)=\int_{0}^{x}dte^{t^{2}-x^{2}}/x$ and $F_{2}(x)=(1-e^{-x^{2}})/x$. The use of the asymptotic wave function is well justified when the source size is sufficiently large compared to the range of the interaction [276]. In the single channel case, the deviation from the asymptotic wave function at small $q$ can be accounted for by using the effective range formula [277],

[TABLE]

The integral in the left hand side of Eq. (5.27) gives the correction to Eq. (5.26), when the integrand is multiplied by the factor $e^{-r^{2}/4R^{2}}$. By using Eqs. (5.26) and (5.27), one arrives at the interaction dependent part of the correlation function in the LL model [26],

[TABLE]

where $x=2qR$ and the effective range correction appears in $F_{3}(r_{\mathrm{eff}}/R)=1-r_{\mathrm{eff}}/2\sqrt{\pi}R$. In the formula given in Ref. [26], one assumes $\psi^{(-)}_{\bm{q}}=(\psi^{(+)}_{-\bm{q}})^{*}$ and $|\mathcal{S}|=1$, then the last term in Eq. (5.28) does not exist.

Figure 5.1 displays the interaction dependence of the correlation function in the LL model, $C(q)=1+\Delta C^{\mathrm{LL}}(q)$. The correlation function is given in terms of the scattering amplitude $f(q)$, which is known to be well described by the scattering length $a_{0}$ and the effective range $r_{\mathrm{eff}}$ at low energy.

[TABLE]

Here the effects of Coulomb potential, channel coupling and the imaginary part of the potential (absorption) have been ignored. It should be noted that the above discussion is based on the “nuclear physics” convention for the scattering length, which leads to $\delta\simeq-a_{0}q$ at low energy.

The behavior of the correlation function can be understood on the same footing as discussed above; now the properties of the wave function is represented by the corresponding scattering length. At negative scattering length $a_{0}<0$, the correlation function is always enhanced by the attractive interaction. At positive but small scattering length $R/a_{0}\gtrsim 0.453$, the correlation function is always suppressed. The positive $a_{0}$ means that there is a bound state or the interaction is repulsive, and the scattering wave function has a node at $r\simeq a_{0}$ at low energy. Then the wave function squared is generally suppressed compared with the free wave function. At around the unitary limit $1/a_{0}\simeq 0$, the correlation function is strongly enhanced at low energy, while it becomes close to unity at $qR\simeq 1$. If this qualitative behavior survives other effects, the correlation function measurement can provide useful information on the scattering length.

5.1.4 Effect of collectivity

In high-energy heavy ion collisions, the hot matter created in the collisions undergoes hydrodynamic expansion, which eventually affects final particle spectra. Assuming a static and spherical Gaussian source, Eq. (5.13), is thus a crude approximation, this ignores the dynamical property of the particle emission sources.

The effect of the collective expansion can be taken into account by properly modeling the source function $S(x,\bm{p})=EdN/d\bm{p}d^{4}x$. Denoting the four-velocity of the collective expansion at space-time point $x$ as $u^{\mu}=\gamma(1,\bm{v})$, the particle energy at the local rest frame reads $u\cdot p$. Then, the thermal distribution is modified into $\exp(-u\cdot p/T)$ with the chemical potential ignored for simplicity. This factor causes a correlation between the freeze-out point $x$ and particle momentum $\bm{p}$. As naturally expected, a fast moving source can produce particle with high momentum easier than a static source. Also in microscopic approaches for the collective phenomena such as transport models, the position-momentum correlation is produced through the multiple scattering of particles [278].

So far the effect of the collective expansion on the correlation function has been mainly discussed in $\pi^{\pm}\pi^{\pm}$ correlation where the final state interaction is negligible except for the repulsive Coulomb force. The apparent dependence of the correlation function on the momentum of the pion pairs can be understood as a consequence of the expansion and gives stringent constraints on the property of the hot QCD matter, such as equation of state and transport coefficients, as well as the dynamics of the collisions [30]. In short, one may regard the effect of the expansion as a modification of the source size into momentum dependent effective source size, often called “length of homogeneity” [279, 273]. For instance, taking longitudinally expanding boost-invariant source in which $u^{\mu}=(\cosh\eta_{s},0,0,\sinh\eta_{s})$ with $\eta_{s}=\ln\sqrt{(t+z)/(t-z)}$, one may have an effective source size in the longitudinal direction [279], $R_{L}\simeq\tau\sqrt{\frac{T}{m_{t}}}$, where $m_{t}=\sqrt{p_{t}^{2}+m^{2}}$ is the transverse mass. Thus, the source size of heavy or high momentum particles becomes effectively small. Because one may obtain information on the detailed source shape by measuring $C(\bm{q})$ as three-dimensional function of $\bm{q}$, the corresponding source sizes are differently affected by the profile of the expansion. In the present review, however, we concentrate on the one-dimensional correlation function of $q=|\bm{q}|$ since the current statistics of experiments are not sufficient for such analyses in the specific channels we are interested in. Then, the effect of the collective expansion might be regarded as an effective reduction of the source size $R$. Although rescatterings and resonance decays produce a non-Gaussian tail in the source function [280], this effect can be safely ignored in the analyses of the final-state interaction since the dominant part of the source function in the pair rest frame can be well approximated by a Gaussian [281] and only small distance pairs are important.

5.1.5 Feed-down contribution

In the discussions above, two particles are assumed to be directly emitted from the hot matter. This assumption would be valid if one could remove contributions from the decay of parent particles. In reality, substantial fraction of observed particles come from decay of resonances. Strongly decaying short-lived resonances with lifetimes of $\mathcal{O}(\text{several fm})$ will give the source function an effective long lifetime and a tail of the spatial distribution, and thus might influence low $q$ behavior of $C(q)$ through the change of the source geometry.

However, it has been known that long-lived parents give a sharp correlation near $q\simeq 0$, which cannot be resolved, and thus cause an apparent reduction of the intercept $C(q=0)$ [282, 283, 284]. With $N_{\text{tot}}^{A}$ being the total number of measured particle $A$ of interest and $N_{\text{res}}=\sum_{i}N_{i\rightarrow A}$ being the long-lived parent contribution decaying into $A$, the effective intercept $\lambda$ is given by

[TABLE]

In this case the observed correlation function is expected to take the following form;

[TABLE]

Since $\lambda\leq 1$, the long-lived resonance decay dilutes the strength of the correlation. Practically the same correction should also be applied even without long-lived resonances because the experiments cannot perfectly identify the particles. Thus the $\lambda$ parameter is often called the “purity” parameter, indicating the purity of particle identification in data samples. When the percentage of the misidentification is known and estimation of the long-lived resonance decay contribution is feasible, one may construct a purity-corrected correlation function by inverting (5.31) as

[TABLE]

This correction serves as a crucial input when one tries to extract the pairwise interaction, since the overall magnitude is sensitive to the scattering length as demonstrated in Fig. 5.1 via the inverse of its ratio to the source size.

Decay parents may induce residual correlations to the observed ones. For example, the observed $pp$ correlation may have been affected by $p\Lambda$ correlation before the $\Lambda$ decays into the proton. Introducing a pair fraction (or pair purity) $x_{ij}$, which is defined as a fraction of $(i,j)$ pairs to total number of the pairs of interest, one may include such residual correlations as [285]

[TABLE]

The fraction parameters $x_{i,j}$ as well as the effective intercept $\lambda$ can be estimated from experimental data and production models, and the relative momentum of the parent $Q_{i,j}$ can be obtained from decay kinematics [285]. It should be noted that, however, this is also affected by the purity of the particle identification.

5.2 Non-exotic channels

To begin with, it is instructive to discuss some examples for the correlation functions of non-exotic channels whose interaction is not expected to produce exotic states.

5.2.1 $pp$ and $\bar{p}\bar{p}$ correlation

The correlation method to extract the pairwise interaction, albeit limited to $s$-wave, is applicable to any measurable particle species. In particular, high-energy heavy-ion collisions at the top RHIC energy and the LHC energies produce as many antimatters as matters. The STAR experiments reported measurements of $\bar{p}\bar{p}$ correlation as well as $pp$ correlation in Ref. [16]. The measured $\bar{p}\bar{p}$ correlation is consistent with $pp$ correlation within errors, so are the extracted scattering parameters.

The data were analyzed within the LL model (Sec. 5.1.3), but extended to include the appropriate quantum statistics effect (Eq. (5.24)), residual correlation from $p\Lambda$ and $\Lambda\Lambda$ ($\bar{p}\bar{\Lambda}$ and $\bar{\Lambda}\bar{\Lambda}$ for $\bar{p}\bar{p}$ correlation), and the Coulomb repulsion. The Coulomb interaction can be taken into account by replacing the plane-waves in Eq. (5.8) with the corresponding Coulomb wave functions and by applying the effective range formula with the Coulomb interaction. The pair fractions $x_{pp}$, $x_{p\Lambda}$, and $x_{\Lambda\Lambda}$, are adopted from the THERMINATORS2 model [286] which is an extended version of one of the implementations of statistical models. The extracted low energy scattering parameters of $\bar{p}\bar{p}$ interactions are $a_{0}=-7.41\pm 0.19(\text{stat.})\pm 0.36(\text{sys,})$ fm and $r_{\text{eff}}=2.14\pm 0.27(\text{stat.})\pm 1.34(\text{sys,})$ fm, which are consistent with the known values for protons, $a_{0}^{pp}=-7.82$ fm and $r_{\text{eff}}^{pp}=2.78$ fm. The Gaussian radius was also obtained as $R_{\bar{p}\bar{p}}=2.75$ fm and $R_{pp}=2.8$ fm. These small radii indicate the influence of the collective expansion. Indeed, similar measurements at LHC [287] show that the Gaussian radii scale with $m_{t}$ as indicated in Sec. 5.1.4.

5.2.2 $p\Lambda$ and $p\bar{\Lambda}$ correlations

$p\Lambda$ correlation has been measured in several experiments [17, 18, 19]. Since the $p\Lambda$ interaction is rather known from scattering experiments and hypernuclear data, the $p\Lambda$ correlation measurements in $pA$ and $AA$ have been used to constrain dynamics of the collisions [288]. Nevertheless, high-energy collisions allows for measuring $p\bar{\Lambda}$ and $\bar{p}\Lambda$ correlations [18] which are not known and serve as inputs for transport model calculations, and precise analysis to extract the scattering parameters with modern facilities [18, 19] provides the cross-check with the scattering experiments.

In Ref. [18], the purity-corrected $p\Lambda$ and $p\bar{\Lambda}$ correlations and their anti-particle pairs were reported. The LL model (5.28) was used with the known scattering lengths and effective ranges of the $p\Lambda$ interaction in the spin-triplet $(t)$ and spin-singlet $(s)$ channels ($a_{0}^{t}=-1.66$fm, $a_{0}^{s}=-2.88$fm, $r_{\text{eff}}^{t}=3.78$ fm, and $r_{\text{eff}}^{s}=2.92$ fm [288]) to extract the source size. The measured $p\Lambda$ correlation function is found to fairly reflect the weakly attracting nature, and the extracted source size is found to follow the same trend with that in the $pp$ correlation. On the other hand, $p\bar{\Lambda}$ and $\bar{p}\Lambda$ correlations exhibit small suppression below unity in intermediate range of $q$, $0<q<0.15$ GeV. The LL model was again used with complex scattering length which accounts for the annihilation effects, under the assumption that the effective range is zero and spin dependence is neglected. While the extracted $\text{Im}a_{0}$ is found to be comparable with that of $p\bar{p}$ channel, the source size was found to be significantly smaller than the result from $p\Lambda$ correlation. Later refined analyses in Refs. [34, 281], which also apply a sophisticated dynamical model incorporating hydrodynamic expansion and hadronic cascades, pointed out that including residual correlations of decay parents (see (5.33)) would cure this problem, although current statistics in the data do not allow for a precise determination of the scattering lengths.

5.3 Exotic channels

5.3.1 $\Lambda\Lambda$ correlation

The measurements of the $\Lambda\Lambda$ correlation function in heavy ion collisions provide another constraint on the $\Lambda\Lambda$ interaction. The correlation function is also complementary to direct search of the exotic $H-$dibaryon from the invariant mass of the decay products discussed in Sec. 2.1.3. Indeed STAR collaboration reported the first measurement of the $\Lambda\Lambda$ correlation in Au+Au collisions at the top RHIC energy [22]. The data were corrected only for pair purity via Eq. (5.32) for the identification (92%) after rejecting most of $\Lambda$s from weak decay of higher mass hyperons by using the distance of closest approach to the primary vertex. The long-lived resonance contribution from $\Sigma^{0}$ and a part of $\Xi$ is still supposed to reduce the correlation strength via Eq. (5.31).

In Ref. [22], the data were analyzed within the LL model Eq. (5.28) with an intercept parameter $\lambda$. Furthermore, a Gaussian term with two parameters taking account of the residual correlation at large $q$ is included, although its origin has not been understood. Therefore, a six-parameter fit to the data is made with

[TABLE]

where $\Delta C^{\text{LL}}(q)$ is given by Eq. (5.28). Optimized parameters given in Ref. [22] are summarized in Table 5.1.

Although the quality of the fit is quite well ($\chi^{2}/N_{\text{dof}}\simeq 0.56$), the obtained scattering length222The opposite sign convention of the scattering length is adopted in Ref. [22], $a_{0}=1.10\pm 0.37^{+0.68}_{-0.08}$ fm, seems to conflict with the results from the observed double hypernucleus. Indeed, the $\Lambda\Lambda$ bond energy in ${}^{~{}~{}6}_{\Lambda\Lambda}\mathrm{He}$ is found to be $\Delta B_{\Lambda\Lambda}=B_{\Lambda\Lambda}(^{~{}~{}6}_{\Lambda\Lambda}\mathrm{He})-2B_{\Lambda}(^{5}_{\Lambda}\mathrm{He})\simeq 1.01~{}\mathrm{MeV}$ [38]. From $\Delta B_{\Lambda\Lambda}(^{~{}~{}6}_{\Lambda\Lambda}\mathrm{He})$, the scattering length and the effective range in the $\Lambda\Lambda$ ${}^{1}\mathrm{S}_{0}$ channel are suggested as $(a_{0},r_{\mathrm{eff}})=(-0.77~{}\mathrm{fm},6.59~{}\mathrm{fm})$ [289] or $(a_{0},r_{\mathrm{eff}})=(-0.575~{}\mathrm{fm},6.45~{}\mathrm{fm})$ [290]. Recent update of the bond energy due to the update of the $\Xi^{-}$ mass [291] gives $\Delta B_{\Lambda\Lambda}(^{~{}~{}6}_{\Lambda\Lambda}\mathrm{He})=0.67\pm 0.17~{}\mathrm{MeV}$ [79], which suggests $(a_{0},r_{\mathrm{eff}})=(-0.44~{}\mathrm{fm},10.1~{}\mathrm{fm})$ [292].

A detailed investigation of the $\Lambda\Lambda$ correlation function by making use of the KP formula Eq. (5.7) with various $\Lambda\Lambda$ interaction potentials and source functions including collective expansion in both longitudinal and transverse directions has been carried out in Ref. [36], It was found that after taking into account the correction of electromagnetic decays from $\Sigma^{0}$, the scattering length is found to be consistent with the double hypernuclei. The detailed comparison of the methods is discussed in Ref. [35], which concludes that it is crucial to determine the value of $\lambda$. Here we briefly outline the above points.

First, we clarify the difference between the $C(q)$ obtained from the LL formula Eq. (5.28) and the KP formula Eq. (5.7). In the left panel of Fig. 5.2, $C(q)$ with the fss2 $\Lambda\Lambda$ interaction is displayed. The corresponding values $a_{0}=-0.81$ fm and $r_{\text{eff}}=3.99$ fm are used as inputs for the LL formula. The difference between the two is small, thus confirming previous studies [276] that indicate insensitivity of the correlation to the detailed shape of the wave function within the interaction range. The difference of $C(q)$ between the static spherical source (thin red, circles) and the expanding source (thin green, triangles) indicates the effect of the collective expansion. The existence of the fast boost-invariant longitudinal expansion deforms the source function such that the correlation function takes a different shape in the best fit to the data [36]. Note that such a difference does not take place in the case of non-identical pairs; as seen in Eqs. (5.16) and (5.18), the quantum statistics effect makes $C(q)$ more sensitive to the source shape through the Fourier transformation.

Second, we estimate the the contribution to $N_{\text{tot}}^{\Lambda}$ with the help of the statistical model and experimental data, to correct the data for the long-lived resonance decay via Eqs. (5.30) and (5.31). Here $\Sigma^{0}$ and $\Xi$ are treated as long-lived resonances, since other decay parents have much shorter lifetime thus only change the effective source size or have a negligible contribution. Adopting data from $p+$Be collisions at $p_{\text{lab}}=28.5$GeV [295], we take $N_{\Sigma^{0}}/N_{\Lambda}=0.278$, which is also consistent with thermal model calculations. Taking into account the fact that the $\Xi$ yield in Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV has been shown to be 15% of total $\Lambda$ [296] and the STAR candidate selection with the distance of closest approach less than 0.4 cm may exclude a part of $\Xi$ decay contributions to $\Lambda$, we estimate $\lambda=(0.67)^{2}$. If we take account of the $\Xi$ contribution into the total yields, $\lambda=(0.572)^{2}$. It has been confirmed that the lower value $\lambda=(0.572)^{2}$ only leads to small quantitative changes in the following analyses.

The right panel of Fig. 5.2 compares the results in the fixed $\lambda=(0.67)^{2}$ (solid line) and free $\lambda$ (dashed line) cases in the LL model formula with the residual correlation term Eq. (5.34). Table 5.1 summarizes the results of the fit to the $\Lambda\Lambda$ correlation data. The free $\lambda$ case confirms the result obtained by the STAR in Ref. [22], while the fixed $\lambda$ case shows the opposite sign of the scattering length. In the free $\lambda$ case where the optimal value is found to be $\lambda\simeq 0.18$, quantum statistics and the pair purity give $C(q\to 0)=1-\lambda/2\sim 0.91$ while the data show $C(q\to 0)\simeq 0.82$. Thus $C(q)$ needs to be reduced at small $q$ by the $\Lambda\Lambda$ interaction and positive $a_{0}$ is favored. By contrast, for a fixed $\lambda=(0.67)^{2}$, the corresponding quantum statistical correlation $C_{\scriptscriptstyle\Lambda\Lambda}(q\to 0)=1-\lambda/2\simeq 0.78$ is slightly smaller than the observed correlation. With the residual source contribution, $a_{\text{res}}\sim-0.06~{}\mathrm{fm}$, the difference from the data becomes more evident. The $\Lambda\Lambda$ interaction needs to enhance the correlation, and the optimal $a_{0}$ value is found in the negative region, as concluded in Ref. [36].

One should note that the best fit result of the LL formula in the fixed $\lambda$ case differs from the KP formula result from the cylindrical source including flow effects. This result may indicate the importance of fixing not only the purity $\lambda$ but also the source geometry including the flow effects.

Figure 5.3 summarizes the constraints from the $\Lambda\Lambda$ correlation data at the present stage and its dependence on the assumptions made. Also shown is the boundary of the favored region, given by $\chi^{2}/\text{DOF}=0.65$ ($0.56$), in the fixed (free) $\lambda$ case in the LL model. The region in the free $\lambda$ case is consistent with that by the STAR collaboration [22]. As shown in the previous subsection, negative and positive scattering lengths are favored in the fixed and free $\lambda$ cases, respectively. It is found that negative scattering lengths are more favored in the pair purity probability range of $\lambda>0.35$. Namely, the $\chi^{2}/\mathrm{DOF}$ at the negative $a_{0}$ local minima is smaller than that at the positive $a_{0}$ local minima when $\lambda$ is fixed at a value $\lambda>0.35$. The low energy scattering parameters $(1/a_{0},r_{\text{eff}})$ of several $\Lambda\Lambda$ interactions are also shown; Boson exchange potentials (ND, NF, NSC89, NSC97, ESC08, Ehime) [297, 298, 299, 300, 301, 303, 302] and Nijmegen-based potentials fitted to the Nagara data (FG,HKMYY) [289, 290, 292], in addition to the quark model potential (fss2) [293, 294]. It should be noted that the fixed $\lambda$ region covers recently proposed $\Lambda\Lambda$ potentials, fss2 and ESC08 [294, 303].

The shaded areas in Fig. 5.3 show the favored region in the analysis using the KP formula [36]. The dark (grey) shaded area shows the region with $\chi^{2}/\mathrm{DOF}<5$ from the cylindrical source including flow effects but without feed-down and residual correlation effects. The light (yellow) shaded area shows the region with $\chi^{2}/\mathrm{DOF}\lesssim 1$ under the condition $R>r_{\mathrm{res}}$ with flow, feed-down and residual correlation effects. The light shaded area includes the favored region in the fixed $\lambda$ case in the LL model analysis.

On the basis of the scattering length and the effective range of the $\Lambda\Lambda$ interaction obtained in the present analyses, the existence of $H$ particle as a bound state of $\Lambda\Lambda$ is not preferred. This can be understood from the enhanced $\Lambda\Lambda$ correlation function observed in the data compared with the free case. If there was a bound state in $\Lambda\Lambda$ with $H$ being the dominant component, the correlation function would be suppressed from the free case, as seen in the scattering length dependence of the correlation function (Fig. 5.1).

The existence of $H$ as a resonance pole above the $\Lambda\Lambda$ threshold is another interesting possibility, as discussed in Sec. 2.1.3. The strength of the resonance $H$ signal in the correlation function is shown in Fig. 5.4 [36]. Here the $H$ and $\Lambda$ yields per event per unit rapidity are evaluated by using the statistical model, $N_{H}\simeq 1.3\times 10^{-2}$ and $N_{\Lambda}\simeq 30$ [15], and the mass distributions for the $H$ and $\Lambda\Lambda$ pair are assumed to be the Breit-Wigner and thermal distributions. Then the resonance $H$ contribution to the $\Lambda\Lambda$ relative momentum spectrum is evaluated to be,

[TABLE]

where $f_{\mathrm{BW}}(E)=\Gamma_{H}/[(E-E_{H})^{2}+\Gamma_{H}^{2}/4]/2\pi$ is the Breit-Wigner function, $N_{\Lambda\Lambda}=N_{\Lambda}^{2}$ and $\mu=M_{\Lambda}/2$. The bump structures in the data are roughly explained by adding $\Delta C_{H}$ multiplied by 10 and 2, to a simple smooth function fitting the STAR data for $(E_{H},\Gamma_{H})=(14~{}\mathrm{MeV},4.5~{}\mathrm{MeV})$ and $(E_{H},\Gamma_{H})=(1.8~{}\mathrm{MeV},1.5~{}\mathrm{MeV})$, respectively. These bumps may come from the statistical fluctuations 333We thank N. Shah for this information., but it would be possible to confirm or rule out the existence of resonance $H$ with higher statistics.

5.3.2 $p\Omega$ correlation

As discussed in Sec. 2.1.3, the spin-2 nucleon-Omega ($N\Omega$) state with $S=-3$ [93] is the most promising candidates for bound or resonant dibaryons besides the $H$. The measurement of the $p\Omega$ correlation in order to determine the $N\Omega$ interaction has been recently proposed by Morita et al. [37].

In Ref. [37], the $p\Omega$ correlation function is calculated through the KP formula. Since the $p\Omega$ state has either spin-1 or 2, the wave function can be expressed by the statistical average,

[TABLE]

The interaction in the ${}^{3}S_{1}$ channel is assumed to be a complete absorption at short distance $r<r_{0}$, because there would be a strong coupling to the octet-octet channels when the spatial distance between $N$ and $\Omega$ becomes small. This can be modeled by an imaginary potential $V(r;\ ^{3}S_{1})=-iV_{0}\theta(r_{0}-r)$ with $V_{0}\rightarrow+\infty$ for the strong interaction part. We choose $r_{0}=2$ fm, because the Coulomb potential dominates over the strong interaction for $r>2$ fm.

The interaction in the ${}^{5}S_{2}$ channel was described by the following potential with an attractive Gaussian core + an attractive (Yukawa)2 tail with a form factor; $V_{N\Omega}(r)=b_{1}e^{-b_{2}r^{2}}+b_{3}(1-e^{-b_{4}r^{2}})(e^{-b_{5}r}/r)^{2}$, which well fits the lattice QCD data with heavy quarks ($m_{\pi}$=875 MeV and $m_{K}$=916 MeV) [94] with $b_{1,3}<0$ and $b_{2,4,5}>0$. Assuming that the qualitative form of this attractive potential does not change even for physical quark masses, a series of potentials can be generated by varying the range-parameter at long distance, $b_{5}$. Shown in the following are the results for three typical examples: $V_{\rm I}$ with weaker attraction, $V_{\rm II}$ with a shallow bound state, and $V_{\rm III}$ with stronger attraction. The binding energies, scattering lengths and effective ranges in the ${}^{5}{\rm S}_{2}$ $p\Omega$ channel with and without the Coulomb potential are summarized in Table 5.2.

Figure 5.5 shows the $p\Omega$ correlation function for the various cases. When the Coulomb interaction is turned off (left panel), the behavior of the $C(q)$ only with ${}^{5}S_{2}$ interaction follows the description given in Sec. 5.1.2 and 5.1.3 (see Fig. 5.1). The effect of absorption in the ${}^{3}{\rm S}_{1}$ channel tends to suppress the particle correlation as expected from the vanishing wave function inside the interaction range and resultant scattering length $a_{0}=r_{0}$. The absorption effect is not negligibly small, but is not significantly large enough to change the qualitative behavior of $C(q)$ obtained by the ${}^{5}{\rm S}_{2}$ scattering alone. One notes that when the imaginary part of the potential is finite, the suppression of the correlation becomes the strongest at non-zero $q$, as seen in the analyses of $p\bar{\Lambda}$ correlation (see Sec. 5.2.2).

When the Coulomb interaction is switched on, the long-range attraction gives a strong enhancement of $C(q)$ at small $q$. This is in contrast to the $pp$ and $\bar{p}\bar{p}$ collisions where the attractive strong interaction is separated from the Coulomb repulsion, as represented by the solid lines in the right panel of Fig. 5.5. The different ordering of the three curves ($R_{p}=R_{\Omega}=5$ fm) is due to the large reduction of the scattering length for $V_{\rm II}$ by the Coulomb effect (Table 5.2). For larger source size, $R_{p,\Omega}=5$ fm, the correlation function is more sensitive to the long-range part of the interaction as found for the proton-proton correlation [26, 27]. As a result, the ordering of the correlation function is further changed such that $C(q)$ for $V_{\text{II}}$ becomes the lowest.

In principle, one may try a full Coulomb correction with source-size dependence to isolate the effect of strong interaction. Instead, it has been proposed in Ref. [37] to take an “SL (small-to-large) ratio” of the correlation functions for systems with different source sizes,

[TABLE]

as an alternative and model-independent way to handle the Coulomb effect.

As shown in Fig. 5.6, an advantage of this ratio is that the effect of the Coulomb interaction for small $q$ is largely canceled, so that it has a good sensitivity to the strong interaction without much contamination from the Coulomb interaction. Moreover, taking the ratio of $C(q)$ reduces the apparent reduction of its sensitivity to the strong interaction due to the purity factor. There are in principle two ways to extract $C_{\rm SL}(q)$ experimentally in ultrarelativistic heavy ion collisions at RHIC and LHC: (i) Comparison of the peripheral and central collisions for the same nuclear system, and (ii) comparison of the central collisions with different system sizes, e.g. central Cu+Cu collisions and central Au+Au collisions at RHIC.

By using the SL ratio data, it would be possible to guess the existence or nonexistence of a $N\Omega$ bound state. Suppressed (slightly enhanced) $C_{\mathrm{SL}}$ from unity suggests an attractive $N\Omega$ interaction with (without) a bound state, and strongly enhanced $C_{\mathrm{SL}}$ implies a large scattering length $|a_{0}|$ of the $N\Omega$ interaction.

5.3.3 $K^{-}p$ Correlation

For the elucidation of the structure of $\Lambda(1405)$ and the properties of nuclei consisting of $\bar{K}$, the study of the $\bar{K}N$ interaction is essential. Historically, it has been analyzed from the $\bar{K}N$ scattering amplitude constrained by the $K^{-}p$ scattering experiments. While the recent kaonic hydrogen measurement [63, 64] reduces the uncertainty of the subthreshold extrapolation of the $\bar{K}N(I=0)$ amplitude [65, 66], the accuracy in the $I=1$ component is still to be improved. Here we discuss the $K^{-}p$ correlation in heavy ion collisions using the potential developed in Ref. [69], and demonstrate that the $K^{-}p$ correlation can be used as a complementary observable to the existing $\bar{K}N$ data.

The correlation function is calculated with the KP formula, as discussed in the previous sections. It should be noted that the $K^{-}p$ channel couples with the $\bar{K}^{0}n$ channel, in contrast to the single-channel problems studied so far. The $K^{-}p$ wave function is therefore obtained by solving the coupled-channel Schrödinger equation,

[TABLE]

where $\mu$, $V^{\rm strong}$, and $V^{\rm Coulomb}$ respectively represent the reduced mass, the strong interaction, and the Coulomb interaction. The mass difference between $K^{-}p$ and $\bar{K}^{0}n$ is neglected. It should be noted that the Coulomb interaction acts only in the charged $K^{-}p$ channel, and the strong interaction induces the off-diagonal channel coupling between $K^{-}p$ and $\bar{K}^{0}n$.

To extract the physical meaning of the $K^{-}p$ correlation function, it is instructive to start from the case without the Coulomb interaction. In this case, the isospin basis, $\bar{K}N(I=0)$ and $\bar{K}N(I=1)$, is adequate as well as the physical basis, $K^{-}p$ and $\bar{K}^{0}n$. Considering the relation between these two basis444The phase convention is chosen to be $|K^{-}\rangle=-|I=1/2,I_{3}=-1/2\rangle$.,

[TABLE]

the $\bar{K}N$ interactions in the physical basis can be represented by those in the isospin basis $V^{I=0,1}$ as

[TABLE]

Using $V^{I=0,1}$ from Ref. [69], one can construct the $\bar{K}N$ interaction with the physical basis.

A general form of the $\bar{K}N$ wave function $\Psi^{(-)}_{{\scriptscriptstyle\bar{K}N},\ell=0}$ can be written as the superposition of the isospin wave function $\psi_{I}(r)$, which has the asymptotic form $e^{-i\delta_{I}}\sin(qr+\delta_{I})/(qr)$,

[TABLE]

where $\chi(K^{-}p)$ and $\chi(\bar{K}^{0}n)$ represent the isospin wave function of the physical state. For the wave function used in the correlation function, the $K^{-}p$ channel should satisfy the outgoing boundary condition as in Eq. (5.10). On the other hand, the outgoing wave in the $\bar{K}^{0}n$ channel should disappear. From these conditions, the coefficients $C_{0}$ and $C_{1}$ are determined as $C_{0}=-C_{1}=1/\sqrt{2}$. Thus, the asymptotic $K^{-}p$ wave function is found to be

[TABLE]

Because of the characteristic boundary condition for the coupled-channel correlation function, the obtained $\tilde{\mathcal{S}}_{\scriptscriptstyle K^{-}p}$ is different from the $S$-matrix in the $K^{-}p$ channel $\mathcal{S}_{\scriptscriptstyle K^{-}p}=(\mathcal{S}_{0}+\mathcal{S}_{1})/2$ for usual scattering experiments.

The left panel of 5.7 (left) shows the $K^{-}p$ correlation function without the Coulomb interaction. The source size of nonidentical particle pairs can be estimated as $R=\sqrt{(R_{K}^{2}+R_{p}^{2})/2}$. Considering that the kaon source size in Au+Au collisions at $\sqrt{s_{{\scriptscriptstyle NN}}}=200~{}\mathrm{GeV}$ is estimated as $R_{K}=2$-$5~{}\mathrm{fm}$ [304, 305] and the proton source size is expected to be similar, $R=3.0$ fm is used in this study. Because of the small interaction range of the $\bar{K}N$ potential (0.4 fm [69]) owing to the absence of the $\pi$ exchange, the short range details of the $\bar{K}N$ interaction does not affect the correlation function for the source size $R=3.0$ fm. Actually, the correlation function is well reproduced by the LL model explained in Sec. 5.1.3, as shown by the dashed line in Fig. 5.7 (left), which assumes a zero range interaction and uses the asymptotic behavior for the wave function.

There is another interesting feature, i.e., the existence of the bump and dip structures, around $q\sim 0.05$-$0.15$ GeV/$c$, which does not appear in the $K^{-}p\to K^{-}p$ scattering amplitude. Its origin seems to be the characteristic isospin combination of $\tilde{\mathcal{S}}_{K^{-}p}$ in Eq. (5.66). Especially, the dip structure around $q\sim 0.15$ GeV/$c$ is a good example, because the $K^{-}p$ correlation function is smaller than unity, though both of the $\bar{K}N(I=0)$ and $\bar{K}N(I=1)$ correlation functions are larger than unity in the corresponding energy region [see dashed ($I=0$) and dash-dotted ($I=1$) lines in Fig. 5.7 (left)], reflecting the attractive $\bar{K}N(I=0,1)$ interaction. Thus, the coupled-channel correlation function gives us information complementary to that from the $K^{-}p$ scattering.

For the direct comparison with future experiments, the $K^{-}p$ correlation with the Coulomb interaction is shown by the solid line in Fig. 5.7 (right). Similar to the $p\Omega$ correlation in Sec. 5.3.2, the $K^{-}p$ correlation is largely enhanced by the Coulomb interaction in the small $q$ region ($q\lesssim 0.1$ GeV/$c$). On the other hand, in the relatively higher energy region, the correlation function is determined by the strong interaction. As a result, the interesting dip structure in Fig. 5.7 (left) is kept in the case with the Coulomb interaction in Fig. 5.7 (right).

It should be noted that the $\Lambda(1520)$ effect, which appears in the $d$-wave $\bar{K}N(I=0)$ scattering, is not included in the above results. Because the $\Lambda(1520)$ energy region corresponds to $q\sim 0.24$ GeV/$c$ and the width of $\Lambda(1520)$ is not very large ($\sim 15$ MeV), the inclusion of the $\Lambda(1520)$ would not affect very much the dip structure around $q\sim 0.15$ GeV/$c$. Thus, the interesting feature of the isospin interference is expected to be seen in actual measurements.

6 Summary

High-energy heavy-ion collisions provide a unique opportunity to study the properties of the high energy-density QCD matter formed at the early stage of the collisions. At the same time, these experiments can also be used to investigate the hadronic interactions at low energies because the final stages of heavy ion collision is an excellent environment where heavy hadrons and composite particles including the light (anti)nuclei can be produced.

In this article, we have summarized the present status of recently observed exotic hadrons that can potentially also be measured in a heavy ion collision. We have also reviewed the current understanding on the production mechanisms of hadronic resonances and bound states as well as hadrons of multiquark configurations in heavy ion collision. We have further reported the yields of these particles in relativistic heavy ion collisions and how they can be used to discriminate between different configurations for their structures that otherwise would not be possible by simply considering their quantum numbers. We have finally discussed the current status of two-particle correlation measurements in relativistic heavy ion collisions and how these studies can shed light on the interactions between the particles involved in the measurements. Further theoretical and experimental studies along these directions will open up a new window for understanding the properties of QCD at low energy from high energy heavy ion collisions.

Acknowledgments

This work was supported in part by the Grants-in-Aid for Scientific Research on Innovative Areas from MEXT (Grants No. 24105008 and No. 24105001), by JSPS KAKENHI (the Grant-in-Aid for Scientific Research from Japan Society for the Promotion of Science (JSPS)) with Grants No. 24740152 and No. 16K17694 (Tetsuo Hyodo), No. 16K05349 (Kenji Morita), No. 25247036 and No. 15K17641 (Shigehiro Yasui), No. 15K05079, No. 15H03663, No. 16K05350 (Akira Ohnishi), by the Yukawa International Program for Quark-Hadron Sciences (YIPQS), by the Korea National Research Foundation under the grant number KRF-2011-0030621 and the Korean ministry of education under the grant number 2016R1D1A1B03930089 (Su Houng Lee), by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No. 2016R1C1B1016270) and 2015 Research Grant from Kangwon National University (Sungtae Cho), and by FAPESP and CNPq-Brazil, by US Department of Energy under Contract No. $\mathrm{DE}$-$\mathrm{SC0015266}$ and the Welch Foundation under Grant No. A-1358.