Fully-Heavy Tetraquarks in Strongly Interacting Medium

We study the properties of fully-heavy tetraquarks at finite temperature and their production in high-energy nuclear collisions. We obtain the masses and wave functions of the exotic hadron states $cc\bar c\bar c$ and $bb\bar b\bar b$ by solving the four-body Schr\"odinger equation in vacuum and strongly interacting matter. In vacuum, the tetraquarks are above the corresponding meson-meson mass threshold, and the newly observed exotic state $X(6900)$ might be a $cc\bar c\bar c$ state with quantum number $J^{PC}=0^{++}$ or $1^{+-}$. In hot medium, the temperature dependence of the tetraquark masses and the dissociation temperatures are calculated. Taking the wave function at finite temperature, we construct the Wigner function for the tetraquark states and calculate, with coalescence mechanism, the production yield and transverse momentum distribution of $cc\bar c\bar c$ in heavy-ion collisions at LHC energy. In comparison with nucleon-nucleon collisions, the yield per binary collision is significantly enhanced.

The main difficulty of observing fully-heavy tetraquarks in elementary collisions, such as electronpositron and nucleon-nucleon collisions, is the small production cross-section of heavy quarks. The formation of a fully-heavy tetraquark requires at least two pairs of heavy quarks with small relative momenta, which is very * Electronic address: shuzhe.shi@mcgill.ca rare in an elementary event. This difficulty can be overcame in high-energy nuclear collisions. Since the binding energy among the nucleons of a nucleus can be safely neglected at high energies, a nucleus-nucleus collision contains a number of binary nucleon-nucleon collisions. Therefore, the number of heavy quarks and in turn the number of fully-heavy tetraquarks will be significantly enhanced in high-energy nuclear collisions. From the experimental data [31,32], the charm quark number can reach 10 at the Relativistic Heavy-Ion Collider (RHIC) and even 100 at the Large Hadron Collider (LHC). After the creation in the initial stage of the collisions, the heavy quarks will pass through the new state of matter of light quarks and gluons which is called the quarkgluon plasma (QGP). Due to the strong interaction with the matter, heavy quarks are widely considered as a sensitive probe of the QGP [33,34]. The energy loss of heavy quarks during the evolution in the hot medium makes them be partially or even fully thermalized with the matter before the hadronization. Finally, on the hadronization hypersurface of heavy quarks, tetraquark states are formed via coalescence mechanism [35][36][37][38][39].
The key factor in all coalescence models for light and heavy hadrons is the coalescence probability for quarks to form a hadron state. Considering the big problem of confinement, the coalescence probability or the Wigner function is normally taken as a Gaussian distribution with adjustable parameters [35][36][37][38][39].
Taking into account the fact that charm and bottom quarks are very heavy and their moving velocity is small, there exists a hierarchy of scales in the study of heavy quarks: m mv mv 2 [40,41]. Integrating out the degrees of freedom with momenta larger than m and mv successively in the QCD Lagrangian, one can derive its non-relativistic versions NRQCD and pNRQCD [41]. Furthermore, if neglecting the interaction between colorsinglet and color-octet states, the pNRQCD becomes a potential model [41]. In this case, one can employ the Schrödinger equation to study the properties of hadrons consist of only heavy quarks. The potential model has been successfully applied to open and closed heavy flavors arXiv:2009.10319v1 [hep-ph] 22 Sep 2020 in vacuum and at finite temperature [42][43][44]. It is pointed out that, in comparison with nucleon-nucleon collisions in vacuum the Ξ cc and Ω ccc yields per binary nucleonnucleon collision in heavy-ion collisions at RHIC and LHC will be largely enhanced [45,46]. In this work, we employ the four-body Schrödinger equation to study the properties of fully-heavy tetraquark states cccc and bbbb at finite temperature and their production in high-energy nuclear collisions. Different from light hadrons where the coalescence probability is assumed to be a Gaussian distribution, the probability for fully-heavy tetraquarks is derived from the wave function of the system controlled by the Schrödinger equation. This is essential for predicting the properties of unconfirmed particles.
The structure of the paper is as follows. In Section II we present the theoretical framework of solving the fourbody Schrödinger equation. The tetraquark properties, including mass and size, in vacuum and hot medium are investigated in Sections III and IV. In Section V the total yield and transverse momentum distribution of the fully-charmed tetraquark state cccc in heavy-ion collisions are calculated and compared with its production in nucleon-nucleon collisions. After summarizing in Section VI, we provide supplementary informations about the hyper-spherical harmonic functions in Appendix A, and the method to compute the potentials is described in Appendix B.

II. THEORETIC FRAMEWORK
For a system of four quarks with the same mass m, the wave function Ψ(r 1 , r 2 , r 3 , r 4 ) and the energy E are controlled by the Schrödinger equation where we have assumed that the interaction potential V = i<j V ij is the summation of the two-body interactions, and the direct three-and four-body potentials are neglected. Taking into account one-gluon-exchange interaction, the two-body potential can be effectively expressed as [47,48], where λ a i (a = 1, ..., 8) are the SU(3) Gell-Mann matrices, the factor 1/4 is from the normalization, V c ij is the spin independent interaction, V ss ij is the strength of the spin-spin interaction, and |r ij | = |r i − r j | is the distance between the two quarks labeled by i and j. In order to solve the four-body Schrödinger equation, we first intro-duce the Jacobi coordinates, where µ is a parameter with dimension of mass and its value does not affect the final result [49]. We take µ = M = 4m in numerical calculations. With such coordinates, the kinetic energy becomes Since the potential depends only on the relative coordinates x i , one can factorize the four-body motion into a center-of-mass motion and a relative motion, Ψ(r 1 , r 2 , r 3 , r 4 ) = Θ(X)Φ(x 1 , x 2 , x 3 ). The bound state properties only relate to the relative motion of the system, and we just need to deal with the nine dimensional wave equation. We then express the relative coordinates x 1 , x 2 and x 3 in the hyper-spherical frame: hyper-radius ρ = x 2 1 + x 2 2 + x 2 3 and hyperangles Ω = {α 2 , α 3 , θ 1 , φ 1 , θ 2 , φ 2 , θ 3 , φ 3 }, where the angles α 2 ≡ arcsin(x 2 / x 2 1 + x 2 2 ) and α 3 ≡ arcsin(x 3 /ρ) are defined within the range [0, π/2], and {x i , θ i , φ i } are the spherical coordinates corresponding to x i . With the hyper-spherical coordinates, the Schrödinger equation governing the relative wave function Φ(ρ, Ω) can be written as where K 3 is the hyper-angular momentum operator, and E r the relative energy (binding energy). As shown in (2), the potential V (ρ, Ω) depends on the color and spin degrees of freedom. We start with constructing the color and spin sector of the wave-function, based on the symmetry properties. We follow the analysis in Ref. [50]. The Pauli exclusion principle requires the wave-function to be anti-symmetric when exchanging two identical fermions, i.e. two quarks or two anti-quarks. From the decomposition in color space, there are two color-singlet states obtained from the first and second terms on the right-hand side. We label them as For the exchange between the two quarks or two antiquarks, |φ 1 is anti-symmetric and |φ 2 symmetric. The decomposition in spin-space is There are two s = 0 states, three s = 1 states, and one s = 2 state where the subscripts denote the spin of the subsystems QQ andQQ and the whole system QQQQ.
As we focus on the tetraquark states consist of identical quarks and anti-quarks, the flavor wave function is symmetric by definition. As a first step, we consider the states with vanishing orbital angular momentum, the space wave function is then symmetric. The Pauli exclusion principle only allows the following combination of color and spin wave functions: |φ 1 χ 2 and |φ 2 χ 1 for the states with J P C = 0 ++ , |φ 1 χ 5 for J P C = 1 +− , and |φ 1 χ 6 for J P C = 2 ++ . While these states are orthogonal to each other, the matrix element φ 2 χ 1 |(λ a i · λ a j )(s i · s j )|φ 2 χ 1 = − 3/2 should be particularly noted in the calculation of the potential.
For the tetraquark states with J P C = 0 ++ , the colorspin wave function is a mixture of |φ 1 χ 2 and |φ 2 χ 1 , and the potential contains diagonal and off-diagonal elements in color-spin space, where for the potentials V c ij and V ss ij we have explicitly labeled the two quarks with indices i, j = 1, 2 and the two anti-quarks with i, j = 3, 4.
For the states with J P C = 1 +− or J P C = 2 ++ , the color-spin wave functions are the eigenstates of both (λ a i · λ a j ) and (λ a i ·λ a j )(s i ·s j ), and the corresponding potentials are The potential V (ρ, Ω) depends not only on the hyperradius but also the eight hyper-angles. In this case, the Schrödinger equation (5) cannot be further factorized into a radial part and an angular part. Instead, one expands the wave function in terms of the hyper-spherical harmonic functions Y κ (Ω) which are the eigenstates of the hyper-angular momentum operator K 2 3 , where κ stands for all the quantum numbers related to the angels, and K is the quantum number describing the magnitude of the angular momentum. Some properties of the hyper-spherical harmonic functions Y κ (Ω) which will be used in the following calculation are presented in Appendix A, and the details can be found in Refs. [49,51,52].
With the above preparations, we now write down the relative wave functions for the 0 ++ states, for the 1 +− states and for the 2 ++ states, where R κ (ρ) is the radial wave function corresponding to the hyper-spherical harmonic function Y κ (Ω). Substituting the above expansions into the relative equation (5), one obtains the coupled radial equations, for 1 +− and 2 ++ states, and − 1 2µ for 0 ++ states, where V κκ is the potential matrix element in angular momentum space with the volume element It is worth noting that computing the potential matrix is non-trivial. In the most general form, (21) is an eight dimensional integral, which is computationally expensive. However, taking the assumption that the total interaction potential is the summation of two-body interaction V ij (|r ij |), one can reduce (21) into a one dimensional integral by performing particle permutation. We show the details of such simplification in Appendix B.
In real calculation, one can only include a finite number of hyper-spherical harmonics, a truncation shall be made according to the symmetry properties of the system. Since we focus in this work on the tetraquark states with vanishing orbital angular momentum, the relevant hyper-spherical harmonics are those corresponding to vanishing total orbital angular momentum L and magnetic quantum number M , i.e. L = M = 0. We choose all such hyper-spherical harmonic functions with hyper-angular quantum number K ≤ 3. This leads to seven coupled differential equations which are numerically solved by using the Inverse Power Method [53,54]. The main advantage to take the Inverse Power Method is its high precision for both ground and excited states.

III. TETRAQUARKS IN VACUUM
We start with computing the cccc and bbbb bound states in vacuum. We employ the Cornell potential to describe the spin independent central interaction V c ij between two quarks and the lattice result [48] for the spinspin coupling, The parameters α, σ, β, γ and the quark mass m in the potential model are fixed by fitting the experimental data of charmonium and bottomonium masses. We calculate the quarkonium states QQ via the two-body Schrödinger equation with the potential where the factor 4/3 is the color factor for color-singlet states. With the model parameters presented in Table I, we obtain the quarkonium masses shown in Table II. One can see that, the potential model is effective in describing the heavy quarkonia. With the known parameters, we then solve the radial equations (19) and (20) for fullyheavy tetraquarks. It should be worth noting that similar calculations have been done in literatures [20][21][22][23][24], by taking different potentials and employing different numerical method, e.g. variational method based on Gaussian expansion [55]. We emphasize that the variational method is more appropriate for ground states. The numerical algorithm we use in this work, namely the Inverse Power Method, can be applied to any bound state, including excited states. The tetraquark mass comes from the summation of the constituent masses M = 4m and the binding energy E r which is determined by the radial equations, The root-mean-squared radius of the tetraquark state can be expressed as [34] The root-mean-squared radius is µ independent. Taking µ = 4m in our numerical calculations makes the prefactor be equal to unity, and the hyper-radius ρ can be considered as the radius of the tetraquark state. The calculated mass and mean radius for the ground and radial-excited states 1S, 2S and 3S with quantum numbers J P C = 0 ++ , 1 +− and 2 ++ are shown in Table III, where we have used the orthogonal and normalized condition for the hyper-spherical harmonic functions Y κ (Ω) and the normalization κ |R κ (ρ)| 2 ρ 8 dρ = 1 for the   radial functions R κ (ρ). The mass spectrum is also plotted in Figure 1. All the fully-heavy tetraquark states lay above the meson-meson mass threshold, 2m J/ψ or 2m Υ , shown as dotted lines in Figure 1. The J P C dependence of the mass is weak, and the small difference comes mainly from the spin-spin interaction.
For the 0 ++ tetraquarks, there are two possible colorspin states |φ 1 χ 2 and |φ 2 χ 1 for any ground and radialexcited state. The mixture between the two color-spin states, see the coupling between the two radial functions R (1) and R (2) in (20), will modify the tetraquark mass. The two modified masses are listed in Table III and shown in Figure 1. The left and right ones in Table 3 and lower and higher ones in Figure 1 correspond to the modified results based on the states |φ 2 χ 1 in the representation 6 c ⊗6 c and |φ 1 χ 2 in3 c ⊗ 3 c . To see the modification from the coupling between the two colorspin states, we show also the radial probability fractions in Figure 2 for the ground state 1S of the fully-charmed tetraquark cccc. The thick and thin lines represent the fractions |R (1) κ | 2 ρ 8 and |R (2) κ | 2 ρ 8 . The small difference between the two shows a strong coupling, and the very The calculated tetraquark mass MT and the root-mean-squared radius rrms for the ground and radial-excited states, 1S, 2S and 3S of cccc and bbbb with quantum numbers J P C = 0 ++ , 1 +− , and 2 ++ .  small contributions from larger κ indicate a very fast convergence in the numerical calculation. It is easy to understand that, for excited states the contributions from larger κ should increase. To guarantee a good convergence for both ground and excited states, κ runs from 1 to 7 in our calculation. A big problem in the study of multi-quark states is to distinguish the multi-quark states from molecular states. In the case of tetraquarks, while |φ 1 and |φ 2 form a complete and orthonormal set of color-singlet states, there can be other representations, especially the mesonmeson states. For instance, a quark pair forms a meson state, and the other pair the other meson state, and then the two combine into a meson-meson state. Considering all the possible combinations, there are two meson-meson states, These two states are not orthogonal to each other but orthogonal to the possible color-octet states The above four molecular states can be expressed as a linear combination of the color-singlet states |φ 1 and |φ 2 , and |φ 4 = − 1/3|φ 1 + 2/3|φ 2 , Therefore, a tetraquark state with quantum number J P C = 0 ++ can be expanded in color space in terms of either |φ 1 and |φ 2 or |φ 3 and |φ 5 or |φ 4 and |φ 6 . The projection probabilities of each 0 ++ state are shown in Table IV. Finally, we look at the exotic hadron X(6900) recently observed by the LHCb Collaboration [30]. The current experiment measures only the mass and the width, and is not able to determine the spin and parities J P C . Our theoretical result indicates that, X(6900) maybe the first radial excited state 2S of cccc with J P C = 0 ++ (6908 MeV) or 1 +− (6896 MeV). No matter what the quantum number J P C is, X(6900) is likely to be a tetraquark state, instead of a meson-meson state.

IV. TETRAQUARKS IN HOT MEDIUM
It is widely accepted that there exists a deconfinement phase transition from hadron gas to quark matter at high temperature and baryon density. From the lattice QCD simulations [57,58] and many effective model studies [59,60], the critical temperature of the transition is about T c = 165 MeV at zero baryon density. Considering the fact that heavy-quark mass is much larger than the temperature scale, the tightly bound states of heavy quarks, such as J/ψ and Υ, can survive in the quark matter and be used to probe the properties of the new matter. In this Section, we study the temperature dependence of the tetraquark properties and their dissociation temperatures in the QGP.
In the color-deconfined QCD medium, the heavy-quark potential is screened, and the long-range interaction is strongly weakened when the temperature is high enough. The lattice QCD simulations indicate that, the finitetemperature potential between a pair of heavy quarks can be approximated by the free energy F (r, T ) [34,61,62]. For the heavy-quark bound-states in the hot QCD medium, we take the free energy F (r, T ) as the central potential V c ij (|r ij |, T ), and neglect the finite-temperature corrections to the spin-spin interaction, where Γ and K 1/4 are the Gamma functions and modified Bessel function of the second kind. The temperature dependent screening mass µ(T ) is extracted from fitting the lattice data. We solve again the coupled radial equations (19) and (20) with the central interaction (31) and obtain the binding energy E r (T ) and relative wave function Φ(ρ, Ω, T ) as functions of temperature. The radial probabilities for the ground and excited states of the fullycharmed tetraquark cccc with J P C = 1 +− at the critical temperature and the comparison with the vacuum result are shown in Figure 3. At finite temperature, the long range confinement force is suppressed and the interaction strength is weakened due to the color screening. As a result, the wave function expands outward, and the averaged size of the tetraquark becomes larger in comparison with vacuum, especially for the excited states. The temperature effect changes also the radial symmetry of the system. The random thermal motion of the heavy quarks will smear the angle dependence of the wave function, and the asymmetric components with larger values of the hyper-angular quantum number κ are suppressed. These features for tetraquarks are consistent with the properties of quarkonia and heavy-flavor baryons at finite temperature [44,63].
From Figure 3 the wave function for the second radial excited state 3s expands with temperature very fast, and the peaks in vacuum almost disappear at the critical temperature T c . This means that the bound state is close to the disappearance. Similar to the definition for quarkonium dissociation, the tetraquark dissociation temperature T d is defined as the divergence of the size and the vanish of the binding energy, The dissociation temperatures for different tetraquark states are shown in Table V. Considering the very weak J P C dependence of the tetraquark mass M T shown in Figure 1, the dissociation temperature is almost independent of the quantum numbers J P C , and we have neglected this small difference in Table V. Different from the cccc states which are already dissociated a little bit above the critical temperature, the bbbb states can survive in the QGP phase at very high temperature, due to the extremely large mass of b quark. It is clear that, the excited states will disappear first. The deconfinement phase transition can be realized in the early stage of relativistic heavy-ion collisions at RHIC and LHC when the temperature of the system is above the critical temperature T c . The appearance of QGP significantly changes the production mechanism of hadrons. In particular, the production of low-momentum hadrons are dominantly contributed by the coalescence of partons when the QGP cools down due to the expansion of the colliding system and the temperature reaches T c . The coalescence model [35] has successfully described the light hadron production in heavy-ion collisions, especially the quark number scaling of the elliptic flow [36,37] and the enhancement of the baryon to meson ratio [38,39]. Since heavy quarks are rare particles in the QGP, their hadronization is more in line with the spirit of the coalescence mechanism. The production of quarkonia and multi-charmed baryons in heavy-ion collisions are well studied in different coalescence models [45,46,64]. It shows that their production in heavy-ion collisions is largely enhanced due to the combination of uncorrelated charm quarks in the QGP [43]. This provides a most probable way to discover those multi-charmed baryons like Ξ cc and Ω ccc in heavy-ion collisions at the RHIC and LHC energies. On the other hand, the previous studies on exotic hadron production in heavy-ion collisions show that, the yields of exotic hadrons are expected to be strongly affected by their structures [65,66]. Con-sidering the fact that, bottom quarks are rarely created even in heavy-ion collisions at LHC energy, we discuss only the production of fully-charmed tetraquark cccc in this section. Taking into account the cccc dissociation temperature which is almost the same as the critical temperature, the initially produced tetraquarks via nucleonnucleon collisions will be dissociated in the QGP phase and all the tetraquarks measured in the final state are from the coalescence at the freeze-out of the QGP.
In the coalescence model, the differential production cross-section of a tetraquark state is given by where R µ = i (t i , r i )/4 is the four dimensional centerof-mass coordinate of the tetraquark, and P µ represents its four-momentum with P 0 = M 2 T + P 2 being the energy, P = i q i the total three-momentum, P T the transverse-momentum, and y the rapidity. The ninedimensional coordinate and momentum x and p are the shorthands of three relative coordinates and momenta x i and p i (i = 1, 2, 3). They are defined in the rest-frame of the tetraquark. The factor C = 1/1296, 1/432 and 5/1296 for J P C = 0 ++ , 1 +− and 2 ++ states are from the statistics determined by the intrinsic symmetry, i.e. color, spin, and isospin, and σ inel N N is the inelastic crosssection of the corresponding nucleon-nucleon collisions.
In heavy-ion collisions at RHIC and LHC energies, the heavy quarks in the QGP phase are almost all created through the initial nucleon-nucleon collisions. For the four heavy quarks Q, Q,Q andQ to form a tetraquark state, they can be from two, or three or four binary collisions. In a heavy-ion collision (AA) with fixed number N coll of binary nucleon-nucleon (N N ) collisions, the averaged number of combinations to have four quarks reads where n N N QQ represents the averaged pair number of heavy quarks created in a nucleon-nucleon collision. At the colliding energy √ s = 2.76 TeV, the inelastic crosssection is σ inel N N = 65 mb [67], and the production crosssections of heavy quarks are dσ cc /dy = 0.7 mb [32] and dσ bb /dy = 15 µb [68]. Therefore we have the averaged number n N N cc = (dσ cc /dy)/σ inel N N = 1.1 × 10 −2 and n N N bb = (dσ bb /dy)/σ inel N N = 2.3 × 10 −4 . In the calculation we have neglected the probability of producing two pairs of heavy quarks in a nucleon-nucleon collision.
The integration region Σ in the coalescence model (33) is the isothermal hadronization hyper-surface controlled by the critical temperature T (R µ ∈ Σ) = T c , and the integration element dσ µ is the normal four-vector to Σ. Such an isothermal hyper-surface can be extracted from hydrodynamic calculations. The space-time evolution of the QGP phase can be successfully described by relativistic hydrodynamics [69]. The theory is based on the conservation laws of the matter. For ideal hydrodynamics without considering the dissipation of the fluid, the evolution of the QGP is governed by the energy-momentum conservation and baryon number conservation, where T µν = ( + P )u µ u ν − P g µν is the energymomentum tensor with being the energy density, P the pressure and u µ the fluid velocity, and n µ = nu µ is the baryon current with n being the baryon number density. , P and n are functions of temperature T and baryon number n, given by the equation of state of the hot medium. To compute the equation of state, we treat the deconfined phase at high temperature as an ideal gas of gluons, massless u and d quarks, as well as s quarks with mass m s = 150 MeV, and the hadron phase at low temperature as an ideal gas of all known hadrons and resonances with mass up to 2 GeV [70]. The phase transition temperature is chosen as T c = 165 MeV. The initial condition of the equation (35) at proper time τ 0 = 0.6 fm/c is determined by the colliding energy and nuclear geometry, which lead to a maximum initial temperature T 0 = 484 MeV for central Pb-Pb collisions at √ s N N = 2.76 TeV.
For such extremely high-energy nuclear collisions, the baryon number density approaches to zero. By solving the hydrodynamic equations, we obtain the space-time profiles of temperature T (t, r) and fluid velocity u µ (t, r). Based on the hydrodynamic profiles, we can determine the isothermal hadronization hyper-surface Σ and its normal four-vector dσ µ at the hadronization temperature T c .
There are two key ingredients in the coalescence model (33). One is the phase space distribution of heavy quarks f Q (t, r, q) and fQ(t, r, q), and the other is the coalescence probability W (x, p) (Wigner function) for four heavy quarks to combine into a tetraquark. We first consider f Q and fQ. Since fully-bottomed tetraquarks are extremely rarely produced in heavy-ion collisions at RHIC and LHC energies, we will only calculate here charmed tetraquarks. In Pb-Pb collisions at energy √ s N N = 2.76 TeV, the experimental measurement on D-meson elliptic flow [71] indicates that, charm quarks reaches kinetic equilibrium with the QGP. Therefore, we can take the normalized Fermi-Dirac distribution f F D (q) = A/(e uµq µ /T + 1) as the charm quark and anti-quark momentum distribution, where u µ (t, r) is the local fluid velocity of the matter determined by the hydrodynamics (35), and A(t, r) is the normalization factor. The number density n c (t, r) in coordinate space is controlled by the charm conservation law for the charm current n µ c = n c u µ , with the initial condition where T A (r T +b/2) and T B (r T −b/2) are thickness functions of the two colliding nuclei, r T is the transverse coordinate, and b the impact parameter of the collision. Combining the momentum and spatial distributions, we obtain the phase-space distribution function for charm quarks f c (t, r, q) = n c (t, r)f F D (t, r, q).
We now come to the Wigner function which reflects the dynamics of the hadronization of heavy quarks in hot medium. For light hadrons and light-heavy systems, the non-perturbative (confinement) properties of hadronization makes it difficult to theoretically calculate the Wigner function. It is usually to take a double Gaussian distribution [35][36][37][38][39] in the phase space with adjustable parameters which can be fixed by fitting the experimental data. For fully-heavy tetraquark systems, however, one can non-perturbatively solve the corresponding Schrödinger equation with confinement potential and obtain the wave function Φ(x, T ) of the system and in turn the Wigner function via a Fourier transformation, Taking all the ingredients discussed above for the coalescence model (33), we calculated numerically the fullycharmed tetraquark cccc yield and transverse momentum distribution in Pb-Pb collisions at LHC energy. The It is easy to understand the strong tetraquark enhancement in heavy-ion collisions in comparison with nucleon-nucleon collisions, because a nuclear collision consists of N coll nucleon-nucleon collisions. In a most central Pb-Pb collision, N coll reaches 1937. Let us consider the tetraquark yield in a binary nucleon-nucleon collision. From the previous study [19,72,73], the tetraquark production cross-section in p-p collisions is dσ cccc pp /dy = 78 pb at √ s = 7 TeV. In a most central Pb-Pb collision the effective cross-section per binary collision is dσ cccc AA /dy/N coll = 0.77 nb, which is almost 10 times larger than that in corresponding p-p collisions! The reason for this nontrivial enhancement is from the many combinations for having four quarks to form a tetraquark state. It is highly nonlinear in N coll , see (34).
The difference between nucleus-nucleus and nucleonnucleon collisions is not only the yield but also the momentum distribution. In nucleon-nucleon collisions, the initially created heavy quarks via hard processes carry high momentum, and the produced tetraquarks will inherit the high momentum. In nucleus-nucleus collisions, the heavy quarks lose energy when they pass through the medium and get thermalized before the hadronization. Therefore, the formed tetraquarks via coalescence mechanism are mainly distributed in low momentum region, see Figure 5.

VI. SUMMARY
In this paper, we solved the four-body Schrödinger equation and investigated the properties of fully-heavy tetraquark states cccc and bbbb in vacuum and at finite temperature. To increase the precision of solving the equation, we expanded the wave functions in series of hyper-spherical harmonics and obtained the eigenstates and eigenvalues by using an iteration algorithm based on the Inverse Power Method. This algorithm allows us to study not only the ground but also excited tetraquark states.
In vacuum, we found that the masses of all the tetraquark states 1S, 2S and 3S with J P C = 0 ++ , 1 +− and 2 ++ are above the 2m J/ψ or 2m Υ threshold. The experimentally observed exotic state X(6900) is likely to be a tetraquark state of cccc, and the possible quantum number is J P C = 0 ++ or 1 +− .
At finite temperature, we determined the tetraquark dissociation temperatures due to the color screening effect on the heavy-quark potential. bbbb can survive in almost all the QGP phase, while cccc is already melted at the critical temperature T c . Taking the wave function at finite temperature, we constructed, without introducing any adjustable parameter, the Wigner function in phase space which is the key ingredient of the coalescence mechanism. In the framework of coalescence model, we calculated the production cross-section and transverse momentum distribution for cccc in heavy-ion collisions. Compared to p-p collisions, the production yield, not only for A-A but also for a binary collision, is extremely enhanced in heavy-ion collisions, and the tetraquarks are mainly distributed at low momentum.
Due to the complicated background in nuclear collisionss, it is a challenge to search for rare particles with low/median p T in heavy-ion collisions. However, for fully-heavy tetraquarks, the four-lepton decay channel X → l + 1 l − 2 l + 3 l − 4 can be well separated from the bulk background and makes it possible to find such exotic states in low p T region [19]. In central collisions, the production cross-section of fully-charmed tetraquarks is around three or four orders of magnitude larger than that in p-p collisions, and the leptons produced in the decay channel are energetic but do not interact with the hot medium. Consequently, we expect that the fully-charmed tetraquark shall be able to be measured by lepton detectors at LHC. Acknowledgement: We thank Guojun Huang and Lu Meng for helpful discussions during the work. This work is supported by the NSFC Grant Number 11890712, and SS is grateful to the Natural Sciences and Engineering Research Council of Canada.

Appendix A: Hyper-spherical harmonic functions
For a four-body system with central two-body interaction, the conserved quantities include the orbital angular momental 1 ,l 2 andl 3 corresponding to the relative coordinates x 1 , x 2 and x 3 and L 1 =l 1 , L 2 =l 1 +l 2 and L = L 3 =l 1 +l 2 +l 3 for the 1-2 sub-system, 1-2-3 subsystem and whole four-body system, and the projections L 1z , L 2z and L 3z . Any two of these operators are commutative and their eigenvalues, L 2 , L, M 2 , M, l 1 , l 2 , l 3 , n 2 and n 3 , form a complete set of quantum numbers.
The hyper-spherical harmonic functions Y κ (Ω) are defined as the eigenstates of the hyper-angular momentum K 2 3 of the system, with the solution K = 2(n 2 + n 3 ) + l 1 + l 2 + l 3 , (A2) where N i is the normalization coefficient P niliKi ≡ P li+1/2,Ki−1+(3j−5)/2 ni is the Jacobi polynomial, Y l k ,m k (θ k , φ k ) are the ordinary spherical harmonic functions, and κ stands for all the quantum numbers.