X(3872) Production in Relativistic Heavy-Ion Collisions

Heavy-ion collisions provide a unique opportunity to study the nature of X(3872) compared with electron-positron and proton-proton (antiproton) collisions. We investigate the centrality and momentum dependence of X(3872) in heavy-ion collisions via the Langevin equation and instant coalescence model (LICM). When X(3872) is treated as a compact tetraquark state, the tetraquarks are produced via the coalescence of heavy and light quarks near the quantum chromodynamic (QCD) phase transition due to the restoration of the heavy quark potential at $T\rightarrow T_c$. In the molecular scenario, loosely bound X(3872) is produced via the coalescence of $D^0$-$\bar D^{*0}$ mesons in a hadronic medium after kinetic freeze-out. We employ the LICM to explain both $D^0$ and $J/\psi$ production as a benchmark. Then we give predictions regarding X(3872) production and the nuclear modification factor $R_{AA}^{X(3872)}$. We find that the total yield of tetraquark is several times larger than the molecular production in Pb-Pb collisions. Although the geometric size of the molecule is huge, the coalescence probability is small due to strict constraints on the relative momentum between $D^0$ and $\bar D^{*0}$ in the molecular Wigner function, which significantly suppresses the molecular yield.

Heavy ion collisions provide a unique opportunity to study the nature of X(3872) compared with electron-positron and proton-proton (antiproton) collisions. We investigate the centrality and momentum dependence of X(3872) in heavy-ion collisions via the Langevin equation and instant coalescence model (LICM). When X(3872) is treated as a compact tetraquark state, the tetraquarks are produced via the coalescence of heavy and light quarks near the quantum chromodynamic (QCD) phase transition due to the restoration of the heavy quark potential at T → Tc. In the molecular scenario, loosely bound X(3872) is produced via the coalescence of D 0 -D * 0 mesons in a hadronic medium after kinetic freeze-out. We employ the LICM to explain both D 0 and J/ψ production as a benchmark. Then we give predictions regarding X(3872) production and the nuclear modification factor R X(3872) AA . We find that the total yield of tetraquark is several times larger than the molecular production in Pb-Pb collisions. Although the geometric size of the hadronic molecule is huge, the coalescence probability is small due to strict constraints on the relative momentum between D 0 and D * 0 in the molecular Wigner function, which significantly suppresses the molecular yield.

I. INTRODUCTION
Since the discovery of X(3872) resonance by Belle in 2003 [1], their properties of X(3872) have been extensively studied both experimentally and theoretically [2][3][4][5]. As the mass of X(3872) is just below the D 0D * 0 (orD 0 D * 0 threshold, it might be a meson-meson molecular state with very small binding energy [6][7][8][9][10]. On the other hand, its constituent quark content is generally believed to be of ccqq type. A superposition of two configurations has also been proposed [11,12]. Whether X(3872) is a loosely bound hadronic molecule, a compact tetraquark or just a kinematic effect such as triangle singularity is still under debate [13][14][15]. In protonproton (pp) collisions, the multiplicity dependence of the yield ratio X(3872)/ψ(2S) at LHCb [16] has been measured, and it seems to disfavor the molecular interpretation of X(3872) [17,18]. In Pb-Pb (AA) collisions at the Large Hadron Collider (LHC) [19], the yield ratio X(3872)/ψ(2S) is clearly greater than that of pp collisions. These experimental studies necessitate more detailed theoretical studies about the internal structure of X(3872). Heavy ion collisions have provided a unique opportunity. With the generation of deconfined matter called "quark-gluon plasma" (QGP) in heavy-ion collisions, most primordially produced X(3872) is melted in QGP due to the strong color screening effect and parton inelastic scatterings from the thermal partons. Additionally, the abundant number of charm pairs produced in the early stage can combine with heavy/light quarks to form new hadrons at the quantum chromody- * baoyi.chen@tju.edu.cn † jiangliu 1997@tju.edu.cn ‡ xiaohai.liu@tju.edu.cn § yunpeng.liu@tju.edu.cn ¶ zhao-jx15@tsinghua.org.cn namic (QCD) phase transition. This contribution is significantly enhanced when charm quark densities become large in QGP. Both theoretical and experimental studies have suggested that this coalescence process becomes essential in the charmonium and multi-charm hadron production at the LHC [20][21][22][23][24][25][26][27]. The coalescence model and statistical hadronization model have been extended to study the hadron production in heavy-ion collisions [28][29][30][31][32][33].
In this work, we employ the Langevin equation to simulate the Brownian motion of charm quarks in QGP and D mesons in a hadronic medium. After the diffusion, heavy and light quarks may combine to form D mesons, charmonia, or tetraquarks (ccqq) in QGP. In a hadronic medium, D mesons may also combine to form a hadronic molecule (D 0D * 0 ). The formation processes are described with the instantaneous coalescence model (ICM). The final production mainly depends on two factors: the charm quark spatial and momentum distributions in the medium, and the Wigner functions of the formed particle. The Wigner function can be determined via the Weyl transform of the wave function. Due to the significant difference between the geometric sizes of a compact tetraquark (mean radius r X ∼ 0.3 − 0.5 fm) [17] and a loosely bound hadronic molecule (mean radius r X ∼ 2.5 − 22 fm, see table I) [10,31], their Wigner functions are significantly different from each other in the two scenarios. Their production is expected to become distinguishable in the two scenarios, which allows us to study whether the X(3872) wave function is close to a molecular state or a tetraquark state in heavy-ion collisions.
This paper is organized as follows. In Section (II), we introduce the Langevin equations for the Brownian motions of heavy quarks and D mesons in the medium. Heavy quark distributions in phase space are also presented with different degrees of momentum thermalization. In Section (III), we introduce the hydrodynamic model to simulate bulk medium expansions. In Section (IV), the formation of tetraquark and molecular states are studied with the ICM. We employ the test particle Monte Carlo method to numerically solve the combined model of the Langevin equation and ICM (LICM). In Section (V), we calculate the D 0 spectrum and J/ψ spectrum as a benchmark. Then, we extend the LICM to the predictions of X(3872) as a compact tetraquark and a hadronic molecule, respectively. Their production as a function of the collision centrality and transverse momentum are analyzed in detail. The final conclusion is given in Section (VI).

II. CHARM QUARK EVOLUTION IN PHASE SPACE
A. Langevin equations for charm diffusion In a hot medium, with the assumption of small momentum transfer in each particle scattering, heavy quark trajectories can be treated as Brownian motion. Heavy quark energy loss/gain in QGP is attributed to two processes: elastic scattering with light partons from the medium and medium-induced gluon radiation. From previous studies on D meson spectra from nucleus-nucleus collisions [34][35][36][37][38][39], heavy quark energy loss is dominated by gluon radiation at high transverse momentum p T [40,41] and elastic scattering at low p T [42]. When the local temperatures of QGP drop to the critical temperature T c of the phase transition, the heavy quark potential is partially restored. Charmonia can be generated via the coalescence of c andc quarks in QGP [43][44][45]. The charmonia from the coalescence process are mainly located at low p T and dominate the final charmonium production. We focus on the production of charmonia and X(3872) at small and moderate p T regions, where the gluon radiation effect is negligible. The evolutions of charm quarks in QGP can be described by the Langevin equation [46][47][48], where p is the heavy quark momentum. η and ξ are the drag and random force terms due to the interactions with the bulk medium, respectively. In the meantime, the evolution of D meson in hadronic medium can also be simulated via Langevin equation. The values of the drag force η can be obtained through the fluctuation-dissipation relation, η = κ/(2T E), where E = m 2 + |p| 2 is the heavy quark (or D meson) energy. m is the mass of the particle. The charm quark mass is set to m c = 1.5 GeV, and the D meson mass is set to m D 0 = 1.875 GeV and mD * 0 = 2.007 GeV [49]. T is the local temperature of the medium where a heavy quark (D meson) is located. The momentum diffusion coefficient κ is connected with the spatial diffusion coefficient D s through the relation D s κ = 2T 2 . Lattice QCD and effec-tive model calculations indicate that D s is approximately D s (2πT ) (4 ∼ 10) at a temperature of approximately T c ; see the following review paper: [50][51][52]. Instead of considering the detailed temperature dependence in D s , we estimate mean values D s (2πT ) = 5 for charm quark diffusion in QGP [48] and D s (2πT ) = 8 for D meson diffusion in a hadronic medium [53]. The κ is connected with the random force in Eq. (1) through with i, j = (1, 2, 3) representing three dimensions. t is the time of heavy quark (D meson) evolution in a hot medium.
The Brownian motion of heavy quarks and D mesons is simulated via the test particle Monte Carlo method. At each time step, Eq. (1) is discretized as where n is an integer. x(t) is the position of the heavy quark at time t. Both the heavy quark momentum and position are updated at each time step with Eq. (3)(4). The random noise term in Eq. (5) is set to a Gaussian distribution with the width κ/∆t.

B. Charm quark distributions in Pb-Pb collisions
The initial momentum distribution of charm quarks in pp collisions at √ s N N = 5.02 TeV is obtained with the FONLL model [54]. In Fig. 1, the normalized momentum distribution of charm quarks d 2 N norm pp /dydp T in the central rapidity |y| < 0.9 is plotted with a black solid line. The charm quark production cross-section is set to dσ cc pp /dy = 1.165 mb [55]. In test particle Monte Carlo simulations, the initial momentum of each test particle is randomly generated based on the black solid line plotted in Fig. 1. Then, the charm quarks exhibit Brownian motion with significant energy loss in the medium. When the charm quarks move to the positions where the local temperature of the medium is lower than the hadronization temperature of a certain hadron such as T cc→J/ψ , charm and anticharm quarks may combine to form a new bound state. In Fig. 1, we plot the momentum distributions of charm quarks before and after the evolution in QGP by taking different spatial diffusion coefficient D s (2πT ) = 1, 2, 5. With smaller value of the spatial diffusion coefficient, charm quarks are closer to the limit of complete momentum thermalization. Heavy quarks at high p T are shifted to moderate and low p T regions. From the experimental and theoretical studies of D mesons [35,56], the charm quark diffusion coefficient in QGP is close to D s (2πT ) = 5. This diffusion coefficient value is used in the X(3872) yield predictions. In Pb-Pb collisions, the initial spatial densities of heavy quarks are proportional to the number of nuclear binary collisions. Therefore, the initial positions of the test particles are randomly generated based on the relative distribution: where T A(B) (x T ) = dzρ(x) is the thickness function of two nuclei and x T is transverse position. The nuclear density ρ(x) is set to the Wood-Saxon distribution. b is the impact parameter, defined as the distance between the centers of two nuclei. With the momentum and spatial distributions given in Fig. 1 and Eq. (6), we can randomly generate the initial momentum and initial position for each test particle and then evolve them event-by-event via Eq. (3-4).

III. HYDRODYNAMIC MODEL FOR BULK MEDIUM EVOLUTION
QGP produced in relativistic heavy-ion collisions is a strong coupling medium. Its expansion can be described with hydrodynamic equations. In this paper, we employ the (2+1)-dimensional hydrodynamic model to characterize the time and spatial dependence of the temperatures and velocities of a hot medium via the MUSIC package [57,58]. The viscosities of the medium are set to zero for simplicity. To close the hydrodynamic equations, the equation of state (EoS) of the medium is needed and can be parametrized via the interpolation between lattice EoS for the deconfined phase and the hadron resonance gas EoS for the hadron phase [59]. The two phases are connected with a crossover phase transition. The hot medium is treated as QGP at T ≥ T c and hadronic gas at T < T c , respectively. The critical temperature T c between QGP and a hadronic medium is taken to be T c = 170 MeV. The initial maximum temperature at the center of the hot medium is T 0 (τ 0 , x T = 0) = 510 MeV in the most central collisions (b=0) [25]. The hot medium reaches local equilibrium at τ 0 = 0.6 fm/c, where the hydrodynamic equations start. In Fig. 2, we plot the time evolution of the QGP local temperatures at the center of the hot medium at different collision centralities. The results are from the MUSIC model [57,58].

IV. THE PRODUCTION OF HEAVY FLAVOR HADRONS
As QGP expansion, the local temperature will drops continuously. When the temperature is lower than the dissociation temperature, where the bound state is disappear due to color screening or scattering, the heavy quark may form a heavy-flavor hadron via color recombination with other quarks in QGP. The hadronization process has been studied with the ICM [60][61][62][63] and resonance recombination model (RRM) [46]. In this section, we introduce the extended model LICM to study the production of compact tetraquark and molecular states in heavy-ion collisions. More discussion on the coalescence model and the statistical model for exotic hadrons has been presented in a previous review [64].

A. Charmonium and compact tetraquark states
In a medium with a low temperature, when the relative distance and relative momentum between charm and anticharm quarks become small, they may combine into a bound state. The coalescence probability is determined by the Wigner function of the formed charmonium. We randomly generate c andc quarks based on the distributions in Section (II B); the ensemble averaged coalescence probability between uncorrelated c andc in the reaction c +c → ψ + g is written as where g M = 1/12 for J/ψ is the statistical factor from color-spin degeneracy.
are the distributions of the test particles representing charm and anticharm quarks. p M and x M are the momentum and the coordinate of the position of the formed charmonium ψ. The momentum carried by the emitting gluon in the coalescence reaction is neglected, yielding the relation p M = p 1 + p 2 . The center of the formed charmonium x M = (x 1 + x 2 )/2 is located at the middle point between the charm and anticharm quarks.
P cc→ψ (x M , p M ) events is the ensemble averaged coalescence probability between uncorrelated c andc quarks. The Wigner function f W M (x r , q r ), which serves as the quark coalescence probability, can be constructed via the Weyl-Wigner transform of the charmonium wave function. It satisfies the normalization dx r dqr (2π) 3 f W M (x r , q r ) = 1. For the ground state J/ψ, we take the wave function as a a simple harmonic oscillator, and the corresponding Wigner function f W M (x r , q r ) becomes [60], where the width can be determined via the mean square radius of the formed particle 65]. The root-mean-square radius of J/ψ is taken as r 2 J/ψ = 0.54 fm based on the potential model [52]. x r and q r are the relative coordinate and relative momentum between two constituent quarks in the center of mass frame. Therefore, the coordinates (x 1,2 , p 1,2 ) in Eq. (7) should be boosted into the center of mass frame (x cm 1,2 , p cm 1,2 ) before substitution into the Wigner function: where E cm i = m 2 c + |p cm i | 2 is the energy of the heavy (anti)quark. Note that in the event-by-event simulation, the uncorrelated c andc quarks are unlikely to move to the QGP fluid cells with the coalescence temperature T cc→ψ at the same time. The value of the Wigner function decreases rapidly when the relative distance x r becomes larger than the typical geometry size of the formed hadron, which guarantees that in events where hadrons are formed, c andc quarks are located close to each other, and their local temperatures are almost the same. Heavy quarks are rare particles in QGP. The combination probability of one c and onec is on the order of magnitude of ∼ 1% in relativistic heavy-ion collisions. This makes charmonium production far below the yield of D mesons produced via heavy-light quark coalescence.
With the coalescence probability of random c andc in the hot medium, we can obtain the J/ψ production in Pb-Pb collisions, which is proportional to the square of the charm pair numbers: where p M is the three-component momentum of the formed charmonium respectively. N AA cc is the number of charm pairs produced in the rapidity range ∆y cc in nucleus-nucleus collisions. As the coalescence probability decreases with increasing relative rapidity between charm and anticharm quarks, only those charm and anticharm quarks with similar rapidities may combine to form a bound state. dσ cc pp /dy is the differential crosssection of charm pairs in pp collisions. R S (b, x T ) is the momentum-averaged nuclear shadowing factor of charm pairs in the central rapidity of Pb-Pb collisions. This term is calculated with the EPS09 model [66,67]. The charm pair number is reduced by approximately ∼ 28% after considering the shadowing effect in the centrality 0-20%. This effect becomes weaker in more peripheral collisions.
The above procedure can be extended to four-body coalescence [68]. We separate the judgement of tetraquark coalescence into two parts: the relative momentum and coordinate between heavy and light quarks need to satisfy the Wigner function, and two "diquark"s also satisfy the Wigner function. This coalescence process happens instantaneously at T ccqq→X , which is the hadronization temperature of the tetraquark. The light quark position is chosen to be the same as the heavy quark in the coalescence process, and its momentum p lrf ζ in the local rest frame (LRF) of the QGP is taken as a normalized Fermi-distribution (ζ represents a light quark): where T is taken to be the coalescence temperature of X(3872). N 0 is the normalization factor. In the event-by-event Monte Carlo simulations, the light quark momentum in the LRF of the QGP fluid cell is randomly generated with the normalized distribution from Eq. (13). It is then boosted to to the lab frame p lab ζ . The diquark momentum is defined as p diquark = p c + p lab ζ . Note that the blueshift of the light quark momentum distribution due to the moving fluid cells is consistently included by boosting the light quark momentum from p lrf ζ to p lab ζ . The four velocities of the QGP fluid cells and the QGP temperature profile T (x T , t) are given by the hydrodynamic model. The light quark thermal mass in the hot medium is set to m ζ = 0.3 GeV in Eq. (13).
With the momenta and coordinates of the diquark and antidiquark, the formation of a tetraquark state is similar to the J/ψ situation. With X(3872) spin of J X(3872) = 1, the tetraquark statistical factor from color-spin degeneracy is g X(3872) = 1/432 in the coalescence equation. The form of the X(3872) Wigner function is that of Eq. (8). With the assumption that the tetraquark root-meansquare radius is approximately 0.3 ∼ 0.5 fm which is similar to J/ψ, we take the coalescence temperature of the tetraquark X(3872) to be T ccqq→X T cc→J/ψ 1.2T c , slightly above the critical temperature T c [52,69,70].

B. Molecular states
As the mass of X(3872) is very close to the threshold mass of D 0D * 0 (orD 0 D * 0 ), it has been suggested that X(3872) is a loosely bound molecular state [10]. The interaction potential between different D mesons can be obtained with the effective Lagrangians, including the contributions from the exchanges of π, η, σ, ρ and ω mesons. The total effective potential for D 0D * 0 (or D 0 D * 0 ) is attractive, as shown in [10]. To regularize the potential, we impose a short-distance cutoff Λ to address the singularity of the effective potential. The cutoff Λ affects the range of the interaction potential. Solving the two-body Schrödinger equation with this potential, we obtain the wave function and the binding energy of the D 0D * 0 molecular state, as shown in table I and Fig. 3. The definition of the binding energy is the mass difference between X(3872) andD 0 D * 0 ; i.e.
From the table, we can see that the binding energy and the mean radius are sensitive to the interaction potential and the cutoff parameter Λ. In Fig. 3, the binding energy is set to 100 keV, the corresponding potential and the wave function are plotted with the mean radius of ∼ 7.6 fm.
In heavy-ion collisions with the production of a hot medium, (anti)charm quarks first form D(orD) mesons at the QCD phase boundary. Then, in the hadronic phase, the D mesons continue diffusing. Due to the low binding energy of the molecular state, molecular X(3872) can be formed only via the coalescence of D 0 andD * 0 mesons after the medium reaches kinetic freeze-out. No molecular states can survive in the hadronic medium above the temperature T kin of the kinetic freeze-out due to the random scattering with surrounding light hadrons.
In this section, we extend the coalescence model to the molecular formation. In the low and moderate p T regions, charm and light quarks can form D mesons at the critical temperature T c . The coalescence formula for D meson is written as (taking D 0 as an example), where H c→D 0 is the hadronization ratio of charm quarks turning into a direct D 0 state (similar for D + , D * 0 , Λ c , etc.) We take the values of the hadronization ratios to be H c→D 0 = 9.5% and H c→D * 0 = 20% [56]. The same hadronization ratios are used in the coalescence ofD 0 andD * 0 mesons. p D is the momentum of the formed D meson. P cq→D 0 (p D ) events is the ensemble-averaged probability of charm quarks turning into D 0 mesons. We assume that all D mesons are produced via the coalescence process and neglect the fragmentation contribu- tion. The integration of Eq. (14) over momentum is 1. dN i /dp i (i=1,2) represents the momentum distributions of two test particles at the positions with the coalescence temperature of D mesons T = T c . The position and momentum of charm quarks at the hadronization surface can be obtained via the Langevin equations given in Eqs. (3)(4)(5). The momentum distribution of the light quark in the LRF of the QGP in Eq. (14) is taken to be the normalized Fermi-distribution; see Eq. (13). The width of the D meson Wigner function is determined in the same way as that of charmonium, and the root-mean-square radius of the D mesons is taken as r 2 D = 0.43 fm [52] for both D 0 and D * 0 .
After the formation of D mesons, they continue to diffuse in the hadronic medium, with a different value for the spatial diffusion coefficient, D s (2πT ) = 8, as discussed in Section II. When D mesons move to the regions where the hadronic medium reaches kinetic freezeout, D 0 andD * 0 mesons may combine to form a loosely bound molecular state. We set the coalescence temperature of molecular X(3872) to be the kinetic freezeout temperature, T mole T kin 0.14 GeV (the kinetic freeze-out temperature can be extracted from experimental data [71]). Due to the uncertainty of the molecular geometric size given in table I, we take its value to be r 2 X = 3.0, 5.5, 9.0 fm in the calculations of hadronic molecule production. Actually, the strategy we used here for the molecular state of X(3872) is similar to the light nuclei production in heavy-ion collisions [72]. The proton and neutron are formed in the QCD phase boundary (also many feed-down contributions in the hadronic phase) and evolve in the hadronic phase. When the system undergoes the chemical freeze-out, the coalescence of light nuclei, such as deuteron and triton, happens.

C. Numerical simulations
We employ the test particle Monte Carlo method to numerically solve the LICM. In each event, two test particles are randomly generated with uncorrelated initial positions and initial momenta. Their dynamical evolution is described with two independent Langevin equations. When they move to the regions where local medium temperatures drop to the coalescence temperature, their relative distance and relative momentum are calculated; these parameters are used in the Wigner functions to calculate the probability of coalescence that forms a new bound state. With the coalescence probability between two test particles, we generate a random number between 0 and 1 and compare it with the coalescence probability. If the coalescence probability is larger than this random number, the new bound state can be formed. Otherwise, the test particles continue independently evolving. In event-by-event simulations, the particle distributions in Eqs. (7 and 14) become delta functions. For example, the charm quark distribution before the coalescence process can be written as where (x c , p c ) includes the coordinate and momentum of the charm quark at the moment of coalescence.

V. D MESON, CHARMONIUM AND X(3872) OBSERVABLES
In the above sections, we introduced the LICM to describe the diffusions of charm quarks and D mesons in the hot medium and the coalescence process. Now, we calculate the spectra of prompt D 0 and J/ψ mesons in Pb-Pb collisions as a benchmark of X(3872) production. Due to different binding energies of D and J/ψ, they are decoupled with hot medium by different temperatures.
In the prompt D 0 spectrum in Fig.4, as we do not include radiative energy loss, the theoretical calculations with D s (2πT ) = 5 (solid lines) underestimate the energy loss of D 0 mesons at high p T . As a compensation, D s (2πT ) = 2 (dotted-dashed lines) is also taken. This can significantly change the spectrum of D 0 mesons at high p T but becomes negligible at low p T , as we expected. The ratio of prompt D 0 meson over total charm number is determined with the ratio given in pp collisions N prompt D 0 /N cc = 39% [55]. We focus on the p T -integrated yields of J/ψ and X(3872), which are dominated by the coalescence process at low and moderate p T . For J/ψ experimental data at high p T , the inclusive production is dominated by the primordial production and B-decay contributions, which are absent in the theoretical calculations (color bands) [74]. This explains why our J/ψ calculations are lower than the experimental data at p T 4 GeV/c. At p T 4 GeV/c, our theoretical calculations explain the experimental data well for both prompt D 0 and J/ψ. The lower and upper limits of the color bands in the J/ψ calculations correspond to the situations with The collision centrality is 0-20%. The experimental data are the J/ψ inclusive production from the ALICE Collaboration [73]. The theoretical results are the regeneration production from the coalescence. The two bands correspond to different values of charm quark spatial diffusion coefficient in QGP. The lower and upper limits of the theoretical bands correspond to situations with and without the shadowing effect, respectively. and without the nuclear shadowing effect.
In the formation of a tetraquark, first, a charm quark combines with a light quark to form a diquark, and then the diquark and an antidiquark combine to form a tetraquark at the coalescence temperature T ccqq→X . As light quarks are abundant in QGP, tetraquark production is mainly determined by the density of charm pairs and the Wigner function of the tetraquark state. Different from pp collisions, most primordially produced tetraquarks are melted in QGP due to the strong color screening effect. The final production of tetraquarks mainly comes from the coalescence process. We plot tetraquark production as a function of centrality in Fig.  5, and J/ψ production is plotted as a comparison. The band of J/ψ production represents the situations with and without the nuclear shadowing effect. In tetraquark production, different values of the width in the Wigner function are considered by setting the root-mean-square radius of the tetraquark to r 2 X = 0.3 fm and 0.54 fm (the latter is the same as J/ψ). First, we can see that the J/ψ production is much larger than the tetraquark production. This is mainly induced by the different statistical factors in the coalescence equation. Our predictions for the tetraquark yield are consistent with Ref. [64]. In Fig. 5, when the geometric size of the tetraquark is increased, its production increases by approximately 40% in the central collisions. However, in peripheral collisions, due to the smaller volume and shorter lifetime of the QGP, charm quarks experience less energy loss in the medium, which increases the relative momentum between uncorrelated c andc. Considering the relative momentum part of the Wigner function given in Eq. (8), with a larger mean radius, the tetraquark yield is more reduced for centrality 60-80%, as shown in Fig. 5. If X(3872) is a molecular state, its binding energy is on the order of ∼ keV. X(3872) is then produced via the coalescence of D 0 -D * 0 orD 0 -D * 0 in the hadronic medium at the temperature at which the medium reaches kinetic freeze-out T mole = 0.14 GeV. The molecular geometric size is much larger than that of the compact tetraquark. Its mean radius and the binding energy are calculated based on the potential model in table I. Determining the exact value of the X(3872) geometric size is beyond the scope of this work. Instead, we take different geometric sizes for the hadronic molecule and calculate the X(3872) production. The root-mean-square radius of the molecular state is set to r 2 X = 3.0, 5.5, 9.0 fm. In Fig. 6, the molecular production with r 2 X = 3.0 fm is at the same order as the tetraquark production. When the molecular geometric size increases, σ in the Wigner func-tion also increases. This gives strict momentum conditions in the coalescence of D 0 andD * 0 mesons. Only when D 0 andD * 0 mesons are separated by a large distance but also carry almost the same momentum can they form a molecular state. This constraint significantly suppresses the molecular yield. In the limit of the molecular binding energy approaching zero, the mean radius of the loosely bound hadronic molecule goes to infinity. This means that D 0 andD * 0 mesons must carry almost the same momentum to form a molecular state, which makes the coalescence probability between D mesons very small in a hadronic medium. Both the tetraquark and molecular yields in Fig. 5-6 show clear centrality dependence. They are proportional to the square of the heavy flavor densities in the hot medium. In more central collisions, more charm pairs and X(3872) are produced. This centrality dependence of X(3872) production is qualitatively consistent with the rate equation model [75]. If the relative momentum between D mesons in the center of mass frame of the hadronic molecule is a few pion mass [76], the the molecular root-mean-square radius is taken as ∼ 3 fm. Molecular production will be strongly enhanced and become comparable with the tetraquark production. If the molecular geometric size is larger, the molecular production is several times lower; see Fig. 5-6. One of the main reasons for this is the Wigner functions used for X(3872) production. In this model, the molecular formation conditions in both physical and momentum space are closely connected via one parameter: the width in the Wigner function. With a very large geometric size for molecular X(3872), more D andD mesons satisfy the spatial formation conditions, but this also results in a strict momentum constraint on the momentum part of the Wigner function. The value of the momentum part of the Wigner function exp(−σ 2 q 2 r ) is significantly reduced when the relative momentum q r between D 0 andD * 0 mesons increases. The consistent constraints from both spatial and momentum formation conditions result in molecular production not increasing with geometric size. The freezeout temperature of the molecular state is a parameter in this work and depends on the collision centrality [77]. We assume that the values of the freeze-out temperature change between 0.16 GeV and 0.10 GeV in different collision centralities. The molecular production with different kinetic freeze-out temperature is plotted with bands in Fig.6. With higher freeze-out temperature, molecular production is enhanced due to the larger spatial density of the D mesons in the hot medium.
With the production of tetraquark and molecular states in Fig.5-6, we can obtain the nuclear modification factor R X(3872) AA of X(3872). First, we calculate the nuclear modification factor of J/ψ at 5.02 TeV Pb-Pb collisions. Take the differential cross section of prompt J/ψ to be dσ  [25,45] and the experimental data [78]. For the production cross section of X(3872) in pp collisions, the yield ratio N is 0.24 ∼ 0.46 in the scenario of r 2 X = 5.5 fm. The ratio between X(3872) and ψ(2S) production in Pb-Pb collisions has also been measured by CMS Collaboration [19]. ψ(2S) prompt production can be estimated via a simple thermal weight factor (m ψ(2S) /m J/ψ ) 3/2 exp(−(m ψ(2S) − m J/ψ )/T ) [60]. The temperature in the exponential factor is taken as the J/ψ coalescence temperature. The yield ratio of ψ(2S) to J/ψ is ∼ 7.3%. Then we obtain the value of the ratio to be around N X(3872) AA /N ψ(2S) AA 0.30 (tetraquark scenario with r 2 X = 0.54 fm) and 0.06 (hadronic molecule scenario with r 2 X = 5.5 fm), respectively. If the geometry size of the molecular state becomes smaller by taking r 2 X to be or smaller than 3.0 fm, the yield of the molecular state can become larger than the tetraquark production. The final production of X(3872) depends on its wave function which is characterized by the parameter r 2 X . Note that the yield ratio from above theoretical calculations are in the low p T region where X(3872) and ψ(2S) are mainly from the coalescence process, while the experimental data in Ref. [19] are located in high p T region where X(3872) are produced by the primordial parton hard scatterings. In Fig. 7, the p T spectra of X(3872) as a tetraquark and hadronic molecule are plotted. The uncertainty bands in the theoretical calculations are due to the different choices for the width in the X(3872) Wigner functions. With an increasing value for the width, the tetraquark and hadronic molecule production values show different changes. Tetraquark production is enhanced, but hadronic molecule production is reduced. This is due to the combined effects from the spatial and momentum formation conditions that are consistently given via the Wigner function in Eq. (8). The peak of the molecular p T spectrum is shifted to larger p T compared with that of the tetraquark spectrum. This is because the molecular state is produced in the later stage of the hot medium expansion and the p T of the hadronic molecule can be shifted by the radial flows of the expanding hot medium. With the violent expansion of the hot medium, its radial flows increase with time, which will be picked up by charm quarks and D mesons via random scatterings with the medium. The p T spectra of different particles produeced at different stages of hot medium expansion will be sequentially modified [81]. At very high p T , in-spired by J/ψ studies, the production of exotic heavy flavor hadrons (or hadronic molecules) in the coalescence process is believed to become negligible compared with the primordial production.
It is interesting to compare the production of X(3872) with other hadronic molecules (light nuclei), such as deuteron (d), helium-3 ( 3 He), and hypertriton ( 3 Λ H). In the high multiplicity p-p collisions, the comparison has been made, and the results indicate any loosely bound hadronic molecule interpretation of X(3872) is questionable [33]. Here, we focus on the Pb-Pb collisions. Due to the lack of experimental data in 5.02 TeV Pb-Pb collisions, we add the results of d, 3 He, and 3 Λ H in 2.76TeV Pb-Pb collisions in Fig. 7. We can see the yield of d is about 2 orders of magnitude larger to the X(3872) production. And the yields of 3 He and 3 Λ H are comparable with the molecular-like X(3872) in heavy-ion collisions. Even the production mechanism of light nuclei and molecular-like X(3872) are similar in relativistic heavyion collision at low p T region. But the abound of protons and neutrons in the hadronic phase enhance the yield of two-body molecular state, d. For the three-body molecular state 3 He and 3 Λ H, the coalescence probability constrains the phase-space distribution of protons, neutrons, and hyperons, which in turn reduce their production. Due to the coalescence production, this behavior is much different from the case in pp collisions, especially in the high p T region [33]. Our results show, in relativistic heavy-ion collision, the hadronic molecule interpretation of X(3872) is not excluded in the low p T region so far.
We also check the sensitivity of X(3872) production with the different choices of parameters. When the coalescence temperature of the tetraquark state is shifted to the critical temperature T c , heavy quarks diffuse to a larger volume in QGP before forming a tetraquark state. The tetraquark yield is fractionally suppressed due to the smaller spatial density of heavy quarks in the medium. This effect is similar in the molecular scenario. Different degrees of heavy quark kinetic thermalization can also affect the final production values of tetraquarks and hadronic molecules. In the limit of charm quark kinetic thermalization, both tetraquark and molecular production can be enhanced by approximately ∼ 2 times compared with the situations in Fig. 5-6. Different from D mesons, the production of X(3872) depends on the square of the charm pair number in heavy-ion collisions. The uncertainty in the charm pair production cross-section dσ cc pp /dy is amplified in X(3872) production. The scope of the work is to distinguish the nature of X(3872) via the geometric size of its wave function, which is one of the most important differences between the compact tetraquark and the loosely bound hadronic molecule.

VI. SUMMARY
In this work, we develop the Langevin equation and instant coalescence model (LICM) to study the produc-tion of open and hidden charm flavors including prompt D 0 , J/ψ and X(3872) in heavy-ion collisions. Calculations regarding J/ψ and D 0 mesons are the benchmark of our predictions regarding X(3872) as a tetraquark state and a hadronic molecule, respectively. The realistic diffusions of charm quarks in quark-gluon plasma (QGP) and D mesons in the hadronic medium are described with the Langevin equation. The spatial and momentum formation conditions of X(3872) are consistently given in the Wigner function, which encodes the internal structure of the formed particle. The compact tetraquark and loosely bound hadronic molecule are produced at different medium temperatures: a tetraquark is formed in QGP above the critical temperature, while a hadronic molecule is formed only in the hadronic medium after the kinetic freeze-out. With the constraints of color-spin degeneracy, X(3872) production as a tetraquark state becomes much smaller than J/ψ production. The geometric size of molecular state is very large, and its binding energy is almost zero. This requires the relative momentum between D 0 andD * 0 mesons to be small to form a loosely bound hadronic molecule. Strict constraints on the relative momentum in the Wigner function significantly suppress the molecular yields. Nuclear modification factor R X(3872) AA of X(3872) as a tetraquark and molecular states are also calculated. Its value becomes R X(3872) AA > 1 and < 1 respectively in the scenarios of tightly bound state and weakly bound state, which is characterized by the parameter of the root-mean-square r 2 X . The ratio N X(3872) AA /N ψ(2S) AA between X(3872) and ψ(2S) production in Pb-Pb collisions can be enhanced and become larger than the value in pp collisions N X(3872) pp /N ψ(2S) pp 0.1 when treating X(3872) as a tightly bound state. Otherwise, the yield ratio is suppressed if X(3872) is a weakly bound state. Different degrees of charm quark kinetic thermalization are studied. It is nonnegligible in X(3872) production, which demonstrates the necessity of realistic heavy quark evolution in the study of X(3872) in heavyion collisions. The coherent treatment of charm quark and D meson evolution in a hot medium and the coalescence process are necessary and meaningful for studies on exotic candidates in heavy-ion collisions.