Production of exotic tetraquarks $QQ\bar{q}\bar{q}$ in heavy-ion collisions at the LHC

We investigate the production of exotic tetraquarks, $QQ\bar{q}\bar{q} \equiv T_{QQ}$ ($Q=c$ or $b$ and $q=u$ or $d$), in relativistic heavy-ion collisions using the quark coalescence model. The $T_{QQ}$ yield is given by the overlap of the density matrix of the constituents in the emission source with the Wigner function of the produced tetraquark. The tetraquark wave function is obtained from exact solutions of the four-body problem using realistic constituent models. The production yields are typically one order of magnitude smaller than previous estimations based on simplified wave functions for the tetraquarks. We also evaluate the consequences of the partial restoration of chiral symmetry at the hadronization temperature on the coalescence probability. Such effects, in addition to increasing the stability of the tetraquarks, lead to an enhancement of the production yields, pointing towards an excellent discovery potential in forthcoming experiments. We discuss further consequences of our findings for the search of exotic tetraquarks in central Pb+Pb collisions at the LHC.


I. INTRODUCTION
There is a long-standing prediction that flavor-exotic four-quark states with two units of heavy flavor, QQqq ≡ T QQ (Q = c or b and q = u or d), are stable against decay into two Qq mesons, the binding energy increasing with the heavy-to-light quark-mass ratio M Q /m q [1][2][3]. The critical value of M Q /m q for binding is somewhat model dependent, but there is nowadays a broad theoretical consensus in the literature -see Ref. [4] for a recent compendium-about the existence of a deeply-bound doubly-bottom tetraquark, T bb , with quantum numbers (I)J P = (0)1 + , strong-and electromagnetic-interaction stable with a binding energy that might be as large as 100 MeV or more [4][5][6][7][8][9][10][11][12][13]. This exciting perspective is further reinforced by recent calculations predicting the stability of tetraquarks with distinguishable heavy quarks, QQ ′qq ≡ T QQ ′ [6,14,15].
Let us review the different recent theoretical studies leading to the stability of the T bb tetraquark. A novel lattice QCD calculation [5] employing a nonrelativistic formulation to simulate the bottom quark finds unambiguous signals for a strong-interaction-stable (0)1 + tetraquark, 189(10) MeV below the corresponding two-meson threshold,BB * . With such binding, the tetraquark will be stable also with respect to electromagnetic decays. Ref. [6] uses the mass of the doubly-charmed baryon Ξ ++ cc recently discovered by the LHCb Collaboration [16] to calibrate the binding energy of a QQ diquark. Assuming that the bb diquark binding energy in a T bb is the same as that of the cc diquark in the Ξ ++ cc , the mass of the (0)1 + doubly-bottom tetraquark is estimated to be 215 MeV below the strong decay thresholdBB * . Combining heavy-quark-symmetry (HQS) mass relations of heavy-light and doubly-heavy-light mesons and baryons with leading-order corrections for finite heavy-quark mass, corresponding to hyperfine spin-dependent terms and kinetic energy shift that depends only on the light degrees of freedom, Ref. [7] predicts that the T bb state is stable against strong decays. More specifically, using as input the masses of the doubly-bottom baryons (not yet experimentally measured) obtained by the model calculations of Ref. [8], Ref. [7] finds an axial-vector tetraquark bound by 121 MeV. Ref. [9] solves the Schrödinger equation with a potential extracted from a lattice QCD calculation for static heavy quarks. Using pion masses of m π ∼ 340 MeV they find evidence for an isoscalar doubly-bottom axial-vector stable tetraquark.
When extrapolated to physical pion masses it has a binding energy of 90 +43 −36 MeV. The robustness of these predictions comes reinforced by detailed few-body calculations using phenomenological constituent models based on quark-quark Cornell-like interactions [4,10], which predict that the isoscalar axial-vector doubly-bottom tetraquark is strong-and electromagnetic-interaction stable with a binding energy ranging between 144-214 MeV for different realistic quark-quark potentials. Recent studies using a simple color-magnetic model come to similar conclusions [11]. The QCD sum rule analysis of Ref. [12] also points to the possibility of a stable doubly-bottom isoscalar axial-vector tetraquark. Finally, the recent phenomenological analysis of Ref. [13] also presents evidence in favor of the existence of a stable T bb state. In summary, the theoretical evidence seems to be very compelling and entices one to claim that the T bb tetraquark is an unavoidable hadron ! On the other hand, the theoretical evidence for the existence of a T cc tetraquark is not so convincing [4,17], the results depending on the dynamical model.
Tetraquarks have the simplest multiquark configuration among the exotic states reported by experiments up to now [18]. The tetraquark picture was first introduced in the light-quark sector [19] as an attempt to explain the inverted mass spectrum (inverted in comparison to the simple quark-antiquark structure favored by the naive quark model) exhibited by the low-lying scalar mesons: a 0 (980), f 0 (980), f 0 (500) and K * 0 (800) [20]. For the heavy tetraquarks, all observed candidates fit to the substructure QQqq (see Refs. [21] for a recent compendium). However, there are not yet heavy tetraquarks with a QQqq configuration reported by experiment. If such exotic states do exist, producing and identifying them is an extraordinary experimental challenge. Most of the discoveries of exotic hadrons in the last decade were made in e + e − collisions, starting with the charmonium-like state X(3872) observed in B → K π ± ψ ′ decays by the Belle Collaboration [22]. In recent years, protonproton collisions at the LHC have shown an enormous potential by confirming some earlier discoveries and also revealing new states. Note that a major difficulty in the production of a T QQ state in e + e − collisions is that two heavy-quark pairs, QQ, produced in hard scatterings must rearrange into QQ andQQ diquarks, which makes it a much rarer event than the production of hadrons with QQ content [7]. Despite these difficulties, the recent estimates in Ref. [23] for the production cross sections of T bb and T bc tetraquarks based on Monte Carlo event generators point towards an excellent discovery potential in ongoing and forthcoming proton-proton collisions at the LHC.
An alternative that circumvents those rare rearrangement processes is the production of quarks by coalescence in the environment of the matter produced in heavy-ion collisions at ultra-relativistic energies, the quark-gluon plasma (QGP), since the number of heavy quarks available for producing such structures is appreciable [24][25][26]. Along with the large array of applications offered by relativistic heavy-ion collisions [27], search of exotic hadrons in the QGP is an exciting new direction in our quest to understanding their structure. In generic terms, the coalescence model is based on an adiabatic approximation, in which the probability for the production of, for example, a tetraquark from deconfined quarks is given by the overlap of the density matrix of the quark distribution with the Wigner function of the tetraquark. Ref. [28] gives a review on applications of the model to hadron formation from a QGP, in which the underlying assumptions of the model and its successes in reproducing hadron yields in relativistic heavy-ion collisions are thoroughly discussed.
Within this perspective, in this work we intend to study the production by coalescence of T QQ tetraquarks in central Pb + Pb collisions at the LHC at √ s N N = 2.76 TeV and √ s N N = 5.02 TeV -√ s N N is the total collision energy per nucleon-nucleon pair in the center of mass (c.m.) frame. We employ the dynamical coalescence model extensively used in exotic-hadron production reviewed recently in Ref. [29]. Previous studies invariably make use of a single Gaussian for the hadron wave functions which, although serving to obtain simple expressions for the yields, are far from being realistic and might be an important source of uncertainty. We avoid such approximation and calculate the Wigner function of the tetraquark employing the four-body wave function obtained from constituent models that correctly reproduce the low-lying meson and baryon spectra. In particular, we use the chiral constituent quark model, χCQM, of Ref. [30] and the Cornell-like interaction, AL1, of Ref. [31]. Both, the χCQM and the AL1 models, predict a binding energy for the T bb tetraquark [4,10] comparable to the recent HQS and Lattice QCD estimates [5][6][7]9].
We also address the important question of how a partial restoration of chiral symmetry affects the coalescence process [28]. Since the coalescence happens at nonzero temperature, at which the coalescing (light) constituent quarks have properties different from those in vacuum, the tetraquark wave function is expected to be modified. The importance of such effects has already been investigated in transport [32] and molecular-dynamics [33] descriptions of hadron production. In the present context, this issue becomes particularly relevant for the stability of the produced tetraquark against two-meson decays, as not only the tetraquark mass is changed from its vacuum value, but the threshold energy, which is given by the sum of the masses of two mesons, is also modified. The in-medium stability of a T QQ state is of central importance for assessing the effects of the interactions of the tetraquark with other particles during the expansion of the system before kinetic freeze-out [34].
The paper is organized as follows. In Sect. II we briefly review the basic features of the coalescence model. In Sect. III we discuss the structure of the tetraquark wave function employed in this work and its use in the computation of the Wigner function. In Sect. IV we present and discuss our results in comparison with previous estimates in the literature.
We will concentrate our discussion on the most promising exotic tetraquark candidate, the isoscalar doubly-bottom axial-vector tetraquark T bb . Finally, in Sect. V we summarize the most important findings of our work.

II. COALESCENCE OF TETRAQUARKS
According to the coalescence model, the probability of producing tetraquark hadrons from quarks in the medium formed in a QGP is given by the overlap of the Wigner function of the produced hadron with the phase-space distribution of the constituents in the medium. Here we follow the developments in Refs. [25,26,29], in which the coalescence model was employed to study the production of exotic hadrons in heavy-ion collisions. The implementation of the coalescence model in those references is particularly suitable for the present investigation since the Wigner function of the produced exotic hadrons is motivated by a nonrelativistic constituent quark model. Explicitly, the number of T QQ hadrons is given by where N j is the total number of quarks of flavor j produced in the collision and g j its degeneracy, p i⊥ is the transverse momentum of a quark with flavor i, T is the hadronization temperature, and ρ W P (x 1 , · · · , x 4 ; p 1 , · · · , p 4 ) is the Wigner function of the tetraquark. P = p 1 + p 2 + p 3 + p 4 is the c.m. momentum. Finally g T QQ is the degeneracy factor of the tetraquark given by (2J T QQ + 1)(2I T QQ + 1) 1 . 1 We will be interested in an isoscalar axial-vector tetraquark, thus J TQQ = 1 and I TQQ = 0. The degeneracy factor g j of quarks of flavor j would correspond to its color-spin degeneracy, i.e., (3 × 2) for each constituent [24,35].
Ref. [29] reviews the details and discusses the different hypotheses made in arriving to the expression of N T QQ . The most important ones are: neglect of transverse flow of the produced matter, consideration of only the central unit rapidity assuming uniform rapidity quark distributions, use of nonrelativistic approximations, and use of a Boltzmann distribution for the transverse quark momenta for the phase-space distribution of the quarks.
In addition, it is further assumed that the time in which the coalescence occurs after the collision is large compared with the internal time scale of the hadron, what allows to omit the contribution from the longitudinal relative momenta. It is worth to note that Ref. [36] has derived an alternative implementation of the coalescence model for the study of exotic hadrons overcoming some of the approximations mentioned above, as it could be to consider relativistic effects or finite-size effects of the produced cluster relative to the emission source.
However, it is explicitly stated in Ref. [36] that the alternative implementation gives significantly different predictions for exotic hadrons with nonzero orbital angular momentum, which is not the case of the T bb state of our interest. In all other cases it gives results very close to the original derivation of Refs. [25,26,29] that we follow in the present work.
The integration over P in Eq. (1) can be done by expressing the tetraquark wave function in terms of the following Jacobi coordinates [37,38]: where M = 4 i m i , with m i being the masses of the constituent quarks. In terms of these coordinates, the Wigner function is given by where ρ W int (r 1 , r 2 , r 3 ; k 1 , k 2 , k 3 ) is given in terms of the tetraquark wave function ψ(r 1 , · · · , r 3 ) as with k i being the conjugate momenta relative to r i . Integrating over P and performing the phase-space integrals in the denominator of Eq. (1), one obtains where V is the volume of the source and we have defined the temperature-dependent overlap function F T QQ (T ): with the reduced masses µ i given by Note that in addition to the explicit T dependence due to the presence of the Boltzmann distribution in Eq. (6), F T QQ (T ) might also acquire an implicit T dependence through the parameters of the constituent model when chiral symmetry restoration effects are taken into account, as it is discussed in the forthcoming sections.

III. EVALUATION OF THE TETRAQUARK WIGNER FUNCTION
We evaluate the tetraquark Wigner function with a four-body wave function obtained from realistic constituent quark models by means of a generalized Gaussian variational method. As mentioned in the introduction, we will present results for two different constituent quark models, χCQM and AL1, to check the robustness of our predictions. Chiral symmetry restoration effects will be addressed by means of the χCQM model, that has already been used for various studies of hadron masses and hadron-hadron interactions inmedium [39,40].
Let us first of all briefly summarize the most important features of the constituent quark models. The χCQM takes into account short-distance perturbative QCD effects through a one-gluon-exchange potential. In addition to the masses for the constituent quarks, dynamical chiral symmetry breaking generates (pseudo) Goldstone bosons, introduced as explicit degrees of freedom via π and σ fields. This aspect makes the model ideally suited to study the effects of partial restoration of chiral symmetry. Quark confinement is incorporated via an effective potential that contains string-breaking effects. The charm or bottom and light quarks interact only via one-gluon exchange and, of course, are subject to the same confining potential-for a detailed review of the model we refer the reader to Refs. [30].
The quark-quark potential in the AL1 model contains a chromoelectric part made of a Coulomb-plus-linear interaction together with a chromomagnetic spin-spin term described by a regularized Breit-Fermi interaction with a smearing parameter that depends on the reduced mass of the interacting quarks. Further details of the AL1 model are given in Ref. [31].
The tetraquark wave function is taken to be a sum over all allowed channels with welldefined symmetry properties [37,38]: where χ csf κ are orthonormalized color-spin-flavor vectors and R κ (r 1 , r 2 , r 3 ) is the radial part of the wave function of the κ−th channel. In order to get the appropriate symmetry properties in configuration space, R κ (r 1 , r 2 , r 3 ) is expressed as the sum of four components, where w(k, r) = ±1. Finally, each R r κ (r 1 , r 2 , r 3 ) is expanded in terms of n generalized Gaussians where s 1 (r), · · · , s 3 (r) are equal to ±1 and a i κ , · · · , f i κ are variational parameters. The latter are determined by minimizing the intrinsic energy of the tetraquark -see Ref. [38] for further details about the wave function and the minimization procedure.
The tetraquark will be stable under the strong interaction if its total energy, E T QQ , lies below all allowed two-meson thresholds. Thus, one can define the difference between the mass of the tetraquark, E T QQ , and that of the lowest two-meson threshold, E(M 1 , M 2 ), namely: where E(M 1 , M 2 ) is the sum of the masses of the mesons M 1 and M 2 . When ∆E T QQ < 0, all fall-apart decays are forbidden and, therefore, a strong-interaction stable state is warranted.
When ∆E T QQ ≥ 0 one is simply dealing with a state in the continuum. Another quantity of interest is the root-mean-square (r.m.s.) radius of the tetraquark, RMS T QQ , given by [10]: The Gaussian nature of the radial functions R r κ (r 1 , r 2 , r 3 ) allows one to obtain an analytical expression for the overlap function F T QQ (T ). Substituting Eq. (8) into Eq. (4), the Wigner function ρ W int (r 1 , r 2 , r 3 ; k 1 , k 2 , k 3 ) can be written as where r ± i = r i ± r ′ i /2 and The overlap function F T QQ (T ) is obtained by performing the eight-dimensional integral over the variables r i , r ′ i and k i⊥ . The integrals can be done analytically, most easily using Cartesian coordinates, since all of them are of the form, with a real and b real or complex.
As we have discussed in the introduction, we will concentrate our discussion on the most promising exotic tetraquark candidate, the isoscalar doubly-bottom axial-vector tetraquark T bb , about whose existence there exists a broad theoretical agreement [4][5][6][7][8][9][10][11][12][13]. In the lowestlying tetraquark configuration all four-quarks are in a relative S wave. Thus, the tetraquark shows a separate dynamics for the compact heavy quark in a color antitriplet (see Fig.   8 of Ref. [4] and Table II of Ref. [37]), and therefore due to Fermi statistics spin 1, and the light antiquarks bound to a color triplet to obtain a total color singlet. To satisfy the Pauli principle, the flavor-antisymmetric light-antiquark pair must have spin 0 while the flavor-symmetric has spin 1. The one-gluon exchange is much more attractive for the good antidiquark, a color triplet with spin and isospin 0. Thus, the total spin and parity of the unavoidable T bb tetraquark are J P = 1 + and its isospin would I = 0 [6,7].

IV. RESULTS
We present results for N T QQ obtained with the χCQM and the AL1 models. We also investigate the effects of a finite-temperature partial chiral symmetry restoration on N T QQ using the χCQM.  Table I displays results for the masses, E T QQ , r.m.s. radii, RMS T QQ , and the binding energies, ∆E T QQ , of the (I)J P = (0)1 + T bb and T cc states in vacuum and for the selected temperatures. We note that the lowest two-meson threshold for T bb and T cc corresponds tō BB * and DD * in relative S−wave, respectively. It can be seen that the tetraquarks are

Vacuum
In-medium compact structures instead of molecular ones, with r.m.s. radii much smaller than 1 fm that remain almost constant with T . With respect to their binding energies, the temperature affects the stability of the tetraquarks making them more stable as T increases. This is mainly due to a larger threshold energy, as can be inferred from Fig. 2(a) of Ref. [39], where the temperature dependence of the meson masses has been evaluated. This latter feature is very important: even when the tetraquark masses are almost T −independent, they become more stable as T increases, indicating that chiral symmetry restoration has a larger impact on the masses of D andB mesons than on the tetraquarks, T QQ . Clearly, the improved stability of the tetraquarks at finite T is a welcome feature for their formation in the matter produced in a heavy-ion collision; because, as we discuss further ahead, once they have been formed, the probability to be destroyed by subsequent interactions with other hadrons (mainly pions) of the medium is diminished.
Next, we present results for the tetraquark yields, N T QQ , given by Eq. (5). For this purpose, it is necessary to specify the values of the hadronization temperature T , the volume V , and the quark numbers N i (i = q, b, c) (the latter being understood as being per unit of rapidity at midrapidity). We use the values given in Ref. [29] which are suitable for Pb+Pb collisions at the LHC energies of √ s N N = 2.76 TeV and 5.02 TeV: T = 156 MeV, V = 5380 fm 3 , N q = 700, and N b and N c are given in Table II. Note that the heavy quarks are produced by hard scatterings at the early stage of the collisions and as such are   Table II. It is important to emphasize that in a comparison of the production cross section of T bb states to that of doubly-bottom baryons, Ref. [23] finds that the latter is 2.4 the former. This result represents an excellent discovery potential of T bb tetraquarks in the near future in a dedicated search at the LCH at √ s N N = 13 TeV, and points to the necessity of obtaining predictions for the yields in that range of energies.  which is almost one order of magnitude smaller. As has been discussed above, Ref. [34] finds essentially the same numbers of Ref. [26]. In all cases, it is clear that the yield from the coalescence model for compact multiquark states is smaller than that for the usual quark configurations as a results of the suppression owing to the coalescence of additional quarks.
We have also obtained results with the AL1 quark model. That model predicts a compact T bb bound state with the same r.m.s. radius predicted by the χCQM, 0.22 fm, although with a binding energy 30% smaller. For the tetraquark yields, the model predicts: for √ s N N = 2.76 TeV, and N T bb = 3.3 × 10 −8 for √ s N N = 5.02 TeV. Although some very small differences can be observed between the yields predicted by the two realistic models, χCQM and AL1, the order of magnitude is the same, being smaller than the simplistic approximation of considering a single Gaussian for the tetraquark wave function. It is reassuring that the results are stable with respect to the difference in the binding energy When including effects due to partial chiral symmetry restoration, Table III reveals that the yields increase roughly by a factor two. Such a modest influence is a consequence of the small effect of partial chiral symmetry restoration on the r.m.s. radii. Although those effects also modify the masses m q entering the Boltzmann distribution, they essentially cancel out in the expression for N T QQ , see Eq. (1). However, as already mentioned, the important feature of the partial chiral symmetry restoration is the improved stability of the tetraquark in-medium, due to the larger threshold for two-meson decays.
The role played by hadronic effects, i.e., changes occurred in the production rate due to the interaction with other particles during the expansion of the medium, was discussed for the case of T cc states in Ref. [34]. The authors conclude that these hadronic effects are negligible for the case of compact states. We recall that the T bb wave function, Eq. (8) [38] for detailed discussions. Therefore, one can safely assume that the abundance of T bb calculated at the QGP phase will not change significantly during the expansion of the hadronic matter as a result of absorption by other hadrons in the medium. The signal for the formation of such states would be through the detection of their weak-decay products with several Cabibbo allowed two-and threebody decay channels [6,7]: , that offer enormous discovery potential as they do not contain identical quarks or antiquarks, which will induce a spin-statistic suppression. Recent flavor SU(3) relations based on a chromomagnetic model [42] confirm the adequacy of these channels to search for doubly heavy tetraquark states at the LHCb and Belle II experiments.

V. CONCLUSIONS AND PERSPECTIVES
The question on whether QQqq ≡ T QQ hadrons can be experimentally observed is of great contemporary interest. Observation of such states will be of help in our quest to understand the structure of the newly observed states with quark compositions beyond the traditional quark-antiquark and three-quark configurations. A major challenge in such a program is the lack of experimental information on the production of these exotic hadronic systems invacuum and in-medium. With this motivation, we have investigated the production of T bb and T cc tetraquarks in relativistic heavy-ion collisions in central Pb+Pb collisions at the LHC in the framework of the quark coalescence model. To our knowledge, this is the first study in which a four-quark wave function obtained by solving exactly the four-body problem using realistic constituent models was used to calculate the hadron Wigner function. In addition, the effects of a partial restoration of chiral symmetry on the coalescence probability have been investigated. We have found that the order of magnitude of the predictions for the tetraquark yields is not modified when using different realistic constituent models, either the χCQM or the AL1. However, the obtained production yields are typically one order of magnitude smaller than previous estimations based on simplified wave functions for the tetraquarks.
Our results indicate that the (I)(J) P = (0)1 + T bb tetraquark is a compact state. It becomes more stable in-medium, when effects of a partial restoration of chiral symmetry are taken into account, leading to an increase in the production yields by a factor roughly equal to two. The improved in-medium stability of the tetraquarks implies a smaller probability to be destroyed by subsequent interactions with other hadrons (mainly pions) of the medium. Therefore, the number of produced tetraquarks is given essentially by that calculated at the QGP phase, which essentially depends on the structure of the state. In short, our results also suggest that measuring the T bb tetraquark from heavy-ion collisions would inform us about the nature of its structure.