Tetraquark bound states in constituent quark models: benchmark test calculations

We investigate the tetraquark bound states that are manifestly exotic using three distinct few-body methods: Gaussian Expansion Method (GEM), Resonating Group Method (RGM), and Diffusion Monte Carlo (DMC). We refer to manifestly exotic states that do not involve a mixture with the conventional mesons through the creation and annihilation of $n\bar{n}$, where $n=u, d$. Our calculations are conducted with two types of quark models: the pure constituent quark model featuring one-gluon-exchange interactions and confinement interactions, and the chiral constituent quark model, supplemented by extra one-boson-exchange interactions. This study represents a comprehensive benchmark test of various few-body methods and quark models. Our findings reveal the superiority of GEM over RGM and DMC methods based on present implements for the tetraquark bound states. Additionally, we observe a tendency for the chiral quark model to overestimate the binding energies. We systematically explore the fully, triply, doubly, and singly heavy tetraquark states with $J^P=0^+,1^+,2^+$, encompassing over 150 states in total. We successfully identify several bound states, including $[cc\bar{n}\bar{n}]_{J^{P}=1^{+}}^{I=0}$, $[bb\bar{n}\bar{n}]_{J^{P}=1^{+}}^{I=0}$, $[bc\bar{n}\bar{n}]_{J^{P}=0^{+},1^{+},2^{+}}^{I=0}$, $[bs\bar{n}\bar{n}]_{J^{P}=0^{+},1^{+}}^{I=0}$, $[cs\bar{n}\bar{n}]_{J^{P}=0^{+}}^{I=0}$, and $[bb\bar{n}\bar{s}]_{J^{P}=1^{+}}$, all found to be bound states below the dimeson thresholds.


I. INTRODUCTION
Since the discovery of X(3872) [1], many heavyquarkonium-like states have been observed in experiments.These states are challenging to be accommodated within the quark-antiquark meson spectrum predicted by quark models, as exemplified in [2,3].They are considered as candidates of the tetraquark states (for recent reviews, see [4][5][6][7][8][9][10]).However, apart from states with the exotic quantum numbers [11][12][13], most of the heavy-quarkonium-like states may be a mixture of the tetraquark states and quark-antiquark states, influenced by the unquenched dynamics such as the creation and annihilation of the light q q pairs (q = u, d, s) [14][15][16].
In the past three years, a series of exotic hadron states composed of at least four (anti)quarks have been observed.The LHCb collaboration first discovered the X(6900) with a quark composition of cccc [17].Subsequently, the CMS [18] and ATLAS [19] collaborations confirmed the existence of the X(6900) state and reported additional candidates for the fully charmed tetraquark states.In 2020, the BESIII collaboration reported the Z cs (3985) state with the minimal quark content ccsu in the recoil-mass spectra of K + in the process e + e − → K + (D − s D * 0 + D * − s D 0 ) [20].Later, the LHCb also reported the state Z cs (4000) with the same quark contents but with a slightly higher mass and larger width.The LHCb collaboration also reported spin-0 and spin-1 states in the invariant mass spectrum of D − K + in the decays B + → D + D − K + [21,22].These states were named as T cs0 (2900) 0 and T cs1 (2900) 0 according to the new naming convention [23].The two states are candidates of the csū d tetraquarks.In addition, two charmedstrange tetraquark states, T cs0 (2900) ++ and T cs0 (2900) 0 , with the minimal quark compositions of csu d and csūd, were observed by LHCb [24,25].The former is the first doubly charged tetraquark state observed in experiments, and the latter is likely its neutral partner.Furthermore, the doubly heavy tetraquark states, anticipated for about forty years [26][27][28][29], were also discovered.The T cc (3875) + state, composed of ccū d, was observed by the LHCb collaboration [30,31].All the aforementioned states consist of four (anti)quarks if we neglect the unquenched effect of the heavy quark-antiquark pairs.Undoubtedly, we are rapidly entering the era of the "genuine" multiquark states.These multiquark states have ignited heated theoretical discussions [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47].
Historically, various quark models have made various predictions regarding tetraquark states.In this study, we focus exclusively on the nonrelativistic quark potential models.We do not consider quark models that parameterize the matrix elements without considering the spatial wave function [48][49][50][51][52], or incorporate relativistic effects [53][54][55] in our analysis.In Fig. 1, we present the case of the ccnn state (n = u, d) with I(J P ) = 0(1 + ) as an example to illustrate the predictions of various quark models.It is evident that these results exhibit significant divergence.Some calculations anticipate deeply bound states located below the DD * threshold with binding energies up to 300 MeV, while others suggest loosely bound states.Additionally, some calculations place the state above the DD * threshold, rendering it unstable.These significant variations can be attributed to differences in the potential models and the methodologies to solve the few-body problems.
In this study, we conduct benchmark calculations aimed at assessing the performance of various few-body methods to analyze the tetraquark bound states.Similar benchmark tests have previously been carried out for the four-nucleon systems [67], in which seven different few-body methods were cross-referenced to compute the four-nucleon systems employing the same nucleon-nucleon (NN) interactions.For the tetraquark states, we will employ three distinct few-body methods, namely the Gaussian expansion method (GEM) [68], the resonating group method (RGM) [69], and the diffusion Monte Carlo method (DMC) [70][71][72][73].We check their con- FIG. 1. Theoretical predictions for the masses relative to the DD * threshold of the ccnn tetraquark state with I(J P ) = 0(1 + ) [28,54,[56][57][58][59][60][61][62][63][64][65][66], where theoretical uncertainties have been neglected.The red pentagram and blue cross correspond to the experimental result and theoretical predictions, respectively.
sistency by comparing their results for the tetraquark bound states using the same quark potential models.Another objective of this study is to explore the distinctions between different quark potential models for the tetraquark states.Unlike the nucleon systems where the NN scattering phase shifts constrain the interaction to several high-precision nucleon forces see [74][75][76], the interactions among quarks exhibit greater variability.Basically, there are two types of constituent quark models: those featuring one-gluon-exchange (OGE) interaction combined with the confinement interaction, and those encompassing both of these interactions alongside an additional one-boson-exchange (OBE) interaction.The latter category is referred to as the chiral constituent quark models (χCQM), owing to the inclusion of the pseudoscalar meson exchange interactions arising from the spontaneous breaking of chiral symmetry.In this work, we designate the former quark models without OBE as the pure constituent quark models (PCQM).The debate over which type of quark models is superior has persisted for many years, yet without definitive conclusions.For instance, it has been argued that the two types of quark models yield qualitatively consistent baryonbaryon scattering results [77].In this study, we select the AL1 and AP1 models [57,78] as the representative examples of PCQM, while we employ a quark model proposed by the group at Salamanca University [79,80] as an example of a chiral quark model.
After testing the reliability of the quark potential models and few-body methods, we proceed to predict the tetraquark states located below the strong decay threshold.We focus on the tetraquark states that have no admixture with the conventional mesons due to the creation and annihilation of the nn pairs, where n = u, d.It is important to note that we loosen the constraint and assume the unquenched effect of the ss pairs is suppressed.In our investigation, we examine systems of QQ Q Q, QQ Qq, QQq q, and Qq q q with J P = 0 + , 1 + , and 2 + , encompassing a total of over 150 systems.Here, Q = c, b and q = u, d, s.
The paper is organized as follows.In Sec.II, we present an introduction to three different quark models and three distinct few-body methods.In Sec.III, we provide a detailed ex-ploration of the possible bound solutions, including the fully, triply, doubly, and singly tetraquark states.In this section, we offer a comprehensive comparison of the results obtained from different models and different few-body methods, alongside a direct comparison with the lattice QCD results.In Sec.IV, we summarize our findings and assessment of various quark models and few-body methods.We also list the final candidates of the tetraquark bound states against strong decays.

A. Constituent quark models
In this work, we only focus on the nonrelativistic quark models with the Hamiltonian, where the kinetic energy of the center of mass T CM is subtracted.We only consider the pairwise interaction.One can find the effect of other kinds of interactions in Refs.[72,81].
A minimal quark model consists of the OGE interaction and confinement interaction.In this work, we choose the quark models proposed in Refs.[57,78], where the power coefficient is 1 for AL1 model and 2/3 for AP1 model.The λ c i is the Gell-Mann matrix for the SU(3) color symmetry, and σ i is the Pauli matrix in the spin space.m i is the quark mass.κ, κ ′ and λ are the coupling constants for the Coulomb interaction, Gaussian hyperfine interaction and confinement interaction, respectively.r 0 and Λ are the typical scale of the hyperfine interaction and overall shift parameter, respectively.All the above parameters were determined by the meson and baryon spectra.One can find their specific values in Ref. [57].For the minimal quark model, one can find different choices of parameter sets in Refs.[82] For the χCQM, we choose the quark model proposed by the group of University of Salamanca (SLM for short) as an example.The idea of this model can be traced back to the Ref. [69] to depict the NN system.The specific form was set up in Ref. [79] and the parameters were redetermined in Ref. [80] by fitting the meson spectrum.The interactions read, where the λ c i is the Gell-Mann matrix for the SU(3) color symmetry and σ i is the Pauli matrix in the spin space.m i is the quark mass.α s , a c are the coupling constants for the OGE interaction and confinement interaction, respectively.Here, the Yukawa-type hyperfine interaction is chosen with a typical scale r 0 .∆ is the overall shift parameter.For the confinement potential, the color screening effect is included, which becomes a linear interaction at the short distance and a constant at the long distance.In addition to the OGE and confinement interaction, the pseudoscalar meson exchange is included considering the spontaneous breaking of the chiral symmetry [83].Meanwhile, the meson-exchange interaction is extended to the scalar-meson-exchange interaction to mimic the two-pion-exchange interaction.All the above OBE interactions read, with where λ a is the Gell-Mann matrix in the SU(3) flavor symmetry.Y (x) = e −x /x is the Yukawa function.The coupling constant g ch is determined by the experimental N N π vertex [69].θ P is the mixing angle to introduce the physical η rather than the one in the octet of the SU(3) symmetry limit.m χ are the experimental masses for the π, K and η mesons.The m σ is determined via the PCAC relation The cutoff Λ χ and Λ σ are determined by fitting the meson spectrum.In this work, we use the parameter values in Ref. [80].It is worthwhile to mention that the vector-meson-exchange interactions are also incorporated in some chiral quark mdoels [84,85].
With the three different quark models, we present the theoretical ground meson spectra in Fig. I.One can see their numerical results agree with the experimental results well.For simplicity, we assume there is no mixing effect between η(nn) with I = 0 and η(ss), which are irrelevant to our tetraquark bound state predictions.

B. Gaussian expansion method
The first few-body method used to solve the tetraquark systems is the Gaussian expansion method [68].Namely, we ex-pand the spatial wave function of r using the following basis, where the r n is taken in geometric progression, r n = r 0 a n−1 .Y lm is the spherical harmonics representing the angular momentum.The basis functions are not orthogonal but could be approximately complete if a large range of r n were taken.It has been proved that the choice of the basis can embed both long-and short-range correlations simultaneously [68].
For a four-body system omitting the motion of the center of mass, there are three independent coordinates.As depicted in Fig. 2, various sets of Jacobi coordinates can be chosen.In principle, different choices of the set of Jacobi coordinate will give the same results if the basis functions are complete.One can choose either set of Jacobi coordinates and construct the basis with the total angular momentum J by combining the spatial angular momenta regarding three coordinates and the spin wave functions.To make the basis function complete, the orbital excited basis functions should be incorporated.However, it takes great pains to handle the angular momentum in GEM, although the strategy has been invented [68].Alternatively, we only use the l = 0 spatial wave functions in Eq. ( 6) but include different Jacobi coordinates to consider the different spatial correlations.In our calculation, we include three different sets of Jacobi coordinates as shown in Fig. 2. For each coordinate, we choose six basis functions with r n in geometric progression.In order to make the basis functions more efficient, we choose different r 0 and r max , r 0 = 0.1 fm, r max = 2 fm q − q or q − q r 0 = 0.1 fm, r max = 2 fm (qq) − (q q) r 0 = 0.1 fm, r max = 1 fm q − q r 0 = 0.1 fm, r max = 5 fm (q q) − (q q) , where we take a large r max for the spatial wave functions between two q q clustering to depict the possible molecular solutions.In general, there are four extra K-type coordinates as mentioned in Refs.[33,68].We have verified that the current selection without them already yields very precise results.

Diquark-antidiquark Di-meson Di-meson
Jacobi coordinates used in the GEM of this work.
In addition to the spatial wave functions, we also have different options to construct the discrete wave function.For the color wave functions, one could use either of the following color wave functions, color-I:  a For simplicity, we assume there is no mixing effects between η(nn) with I = 0 and η(ss), which are irrelevant to our tetraquark bound state predictions. color-II: color-III: The color-I is the diquark-antidiquark basis in the color space.
Color-II represents the dimeson basis, where two basis are not orthogonal but complete.The color-III is orthogonal and complete.For the spin-wave functions, one can choose one of the following basis functions, spin-I: S 12 = 0, 1; S 34 = 0, 1 spin-II: S 13 = 0, 1; S 24 = 0, 1 spin-III: where the S-wave orbital angular momentum is assumed.The spin-I is the diquark-antidiquark basis and spin-II and spin-III are two dimeson basis.Each of them is complete.In addition to above coupling modes in color and spin, one can also construct the discrete basis in a sequence of combining three (anti)quarks first and then with the forth (anti)quark.Similarly, one can find possible options of the complete flavor wave functions.
In our calculations, the wave functions of tetraquark states are expressed as the direct product of flavor wave function χ f , color-spin wave function ψ cs , and spatial wave function ψ: Here, A denotes the antisymmetrization operator, representing the exchange of identical quarks.For instance, in the case of bbnn states, the antisymmetrization operator is defined as , whereas for bcnn states, it becomes A 34 = (1 − P 34 ), where P ij permutes the i-th and j-th (anti)quarks.Our approach involves constructing basis wave functions with fixed quantum numbers, followed by the application of the antisymmetrization operator.It is important to note that antisymmetrization introduces additional constraints, potentially reducing the basis space.In other words, independent basis functions may become linearly dependent after antisymmetrization.To address this, we employ the algorithm outlined in Appendix A to automatically eliminate redundant bases in our calculations.
In our calculation, we test five different choices of wave functions numerically, with where χ f , χ s , χ c and ψ represent flavor, spin, color and spatial functions, respectively.Their superscript shows the specific choice of the wave functions in Eqs. ( 8

C. Resonating group method
In this work, we choose the formalism of RGM in momentum space [69].In the RGM formalism, the wave functions of the tetraquark states are formulated as, where ψ M 1 and ψ M 2 are spatial wave functions of two mesons, with the p 1 and p 2 the relative momentum of the quark and antiquark inside two mesons, respectively.Accordingly, the discrete wave functions χ M 1M 2 csf are also the dimeson type ones.In our calculation, we obtain the meson wave functions from the GEM.The ψ 12 is the relative wave functions of the two mesons with the corresponding relative momentum P .We could act ψ M 1 ψ M 2 χ M 1M 2 csf on the Schrodinger equation from the left, ĤΨ(P , p 1 , p 2 ) = EΨ(P , p 1 , p 2 ).We obtain an equation of ψ 12 , where V D and K Ex are the kernels stemming from the direct diagrams and exchange diagrams in Fig. 3. µ M 1M 2 is the reduced mass of the two clusters.In our calculation, the explicit V D and K Ex can be derived from the meson wave functions.We can solve Eq. ( 18) by performing the partial wave expansion and dicretizing the magnitude of the P and P ′ .In our calculation, we choose the coupled-channel formalism and sum over all the possible two ground meson states with the S-wave relative angular momentum in the wave function of Eq. ( 17).The Eq. ( 18) becomes the coupled-channel integral equations accordingly.One can find details about RGM in momentum space in Ref. [69].
It should be noticed that the diquark-antidiquark spatial wave functions are absent in the trail functions of RGM.If we use the notations of Sec.II B, the wave functions is similar to The trial functions are not as general as Ψ A used in GEM.Meanwhile, the meson wave functions are determined which correspond to the free mesons.The distortion effect of the meson wave functions within the tetraquark bound states is also neglected.Thus, from the variational principle, we expect the RGM will give a higher solutions than those from GEM.

Direct diagram Exchange diagram
FIG. 3. The direct and exchange diagrams in RGM.

D. Diffusion Monte Carlo method
Unlike the two previous methods based on the basis expansion (essentially the variational method), the DMC is a kind of projection Monte Carlo method.One can find the detailed formalism in Refs.[72,73].To make this paper self-contained, we introduce it briefly.The imaginary time Schördinger equation reads, with where E R is a shift parameter of the energy.Φ i are the eigenstates with the eigenvalue E i .One can see if we take the E R to approach the ground state energy, all Φ i except the Φ 0 will be suppressed exponentially after a long-time evolution.
The DMC is implemented by sampling the wave function with walkers.The distribution of walkers represents the wave function.The imaginary time Schrödinger equation is actually a diffusion equation with the source and sink.As shown in Fig 4, one can start from a Ψ int which is not orthogonal to the ground state wave function.For every small time step, the walkers will perform a random walk (diffusion process) and experience the death or birth (branch process).In the branch process, one walker could be replicated for several times or be deleted.After a long-time evolution, the distribution of the  (15).The energies are with respect to the lowest relevant thresholds (in units of MeV)."NB" represents that there is no bound solution." →Excited" labels the excited bound state solutions.

Systems
Thresh.walkers will approach the ground state wave function.For a practical calculation, the importance sampling is adopted where in addition to the diffusion and branch process, there is an extra drift process (see [72,73] for details).

AL1-GEM AP1-GEM SLM-GEM
FIG. 4. Illustration of the DMC method [87], where the wave function is sampled by the walkers.
As shown in Sec.II B, one needs to use either very general or very proper trial functions in the variational method-based approaches to get accurate solutions of a few-body problem.
In other words, one should assign a priori clustering behavior to solve the tetraquark states efficiently.However, DMC could get the ground state energy without the pre-assignment of the clustering behaviors.In principle, the wave function space allowed by the DMC calculation could be very general.The correct cluster behaviors could be obtained automatically after a long-time evolution.It has been proved that in molecular physics [88], solid physics [89], and nuclear physics [90], the DMC method is very efficient and precise.In hadronic physics, the DMC method has been used in quark models in several works [70][71][72][91][92][93][94].However, some advantages of the DMC have not been realized in the past calculations.For example, in Ref. [70], the authors performed the calculations for the fully tetraquark states via DMC and obtained many solutions above the corresponding the dimeson thresholds.In principle, one should get the dimeson thresholds if there are no bound state solutions.
Compared with the electron systems in the molecular physics and solid physics, and the nucleon systems in the nuclear physics, there are two distinct features in the multiquark systems.First, there are complicated color structures for the multiquark systems, which result in complicated discrete wave functions or more coupled channels.Meanwhile, the confinement effect makes the dimeson thresholds the only meaningful thresholds.The four quark threshold and the triquark-quark thresholds are meaningless.To suit the calculations of the multiquark systems, one has to adjust the implements of the DMC method.
In Ref. [72], we improved the implements of the DMC and finally got the dimeson thresholds for the fully tetraquark systems by taking more coupled channels into considerations.In Ref. [73], we adopted the same strategies to calculate the possible tetraquark bound states of the doubly heavy tetraquark states.In this work, we will further compare the results from the DMC with those from the variatonal methods, GEM and RGM.

III. NUMERICAL RESULTS
We adopt the GEM and RGM to calculate the fully, triply, doubly and singly heavy tetraquark states.We also compare the results from GEM and RGM with those from the DMC for the doubly heavy tetraquark states.In our tetraquark calculations, we focus on the difference between the tetraquark masses and the lowest dimeson thresholds in the same quark potential models rather than the absolute values.In Figs. 5, 6 and 7, we present the possible bound solutions, where the theoretical results are shifted to align the relevant dimeson thresholds to the experimental ones.

A. Fully heavy and triply heavy tetraquark
We calculate the J P = 0 + , 1 + , 2 + fully heavy tetraquark states with the following quark contents, where C = ± represents that both states with even and odd Cparities are investigated.The three quark models with different few-body methods yield consistent results.There do not exist the bound state solutions for these systems.These results agree with those qualitatively in Refs.[82,[95][96][97][98] where the AL1 and SLM were used, respectively.In lattice QCD simulations, it was shown that there are no bound states of J P C = 0 ++ , 1 +− and 2 ++ bb bb bound states below the noninteracting dimeson thresholds [99].In Ref. [100], the lattice QCD simulations disfavor the existence of the stable spin-0 bbcc state.These lattice QCD results are also consistent with our findings.
Our results indicate that there is no bound solution below the relevant dimeson thresholds, which is consistent with results in Ref. [101].The bbcq tetraquarks were also investigated in lattice QCD [100,102], where the existence of the bound solutions is inconclusive.It was shown that there is a spin-1 bbcs state below the threshold in Ref. [100], where the finite volume effect was not considered.To pin down its existence, more data is needed to handle the finite volume effect.
We first compare two methods based on variational methods, GEM and RGM.As we expected, the GEM gives the lower bound states than the RGM.For example, for the J P = 1 + [ccnn] I=0 state (the candidate of the experiment T cc (3875) + state), the GEM yields a bound state solution in the AL1 model while the RGM calculation indicates no bound sate.For the same state in the SLM model, GEM gives a very deep bound state with the binding energy about 200 MeV while the binding energy in RGM is about 4 MeV.These disparities arise from the fact that the trial functions or basis functions of GEM are more general than those of RGM where only the dimeson-type wave functions are used.In order to testify the above statement, we also employ the GEM with only the dimeson-type functions, see (19), for the J P = 1 + states.We compare the results with those from the RGM and GEM with general wave functions in Fig. 6.One can see the GEM results with only the dimeson-type trial functions are very similar to those from the RGM.For example, the bound solution of [ccnn] I=0 in AL1 model disappears and the result in SLM model becomes a loosely bound state.We can find the GEM with only the dimeson-type wave functions still yields slightly lower solutions than the RGM.This is because in the implement of the RGM, the meson wave functions is constrained to be the same as those of the free mesons.However, in the GEM-dimeson scheme, the meson wave functions are also determined by the variational parameters, where the possible distortion effect of the meson wave functions is included.
It should be noticed that the constrained trial functions in RGM could change the results qualitatively.One will miss the bound solution of J P = 1 + [ccnn] I=0 state in AL1 model.Meanwhile, one could identify the J P = 1 + [bbnn] I=0 state as the molecular states rather than the compact tetraquark quark states via RGM.Thus, the GEM with the general basis functions is superior than the RGM.
We can also compare the bound state solutions from DMC in Fig 5 and Table III with those from GEM and RGM.It is evident that the majority of DMC results exhibit lower energies compared to the RGM outcomes.Qualitatively, the DMC results are consistent with those obtained through GEM, as exemplified by the presence of deep bound states for the J P = 1 + [ccnn] I=0 configuration in the SLM model.In essence, the DMC method, devoid of any prior constraints on the clustering behavior of the wave functions, naturally provides the roughly proper clustering behavior in the present implementation.A more comprehensive discussion is available in Ref. [73].When compared to GEM results, the masses obtained through DMC generally are higher.For instance, in the case of the J P = 1 + [ccnn] I=0 configuration in the AL1 model, DMC only yields the dimeson thresholds, whereas GEM predicts the existence of a bound state with a binding energy of approximately 14 MeV.It should be noted that the statistical uncertainties in the current DMC calculations are estimated to be at the order of 1 MeV [73].Consequently, it is reasonable to conclude that the observed differences stem from systemic uncertainties within the DMC method that have yet to be fully comprehended.One plausible explanation could be that the importance functions employed in Ref. [73] may not be optimized for addressing the coupled-channel complexities inherent to the tetraquark systems.In principle, the antisym-  15), GEM with only the dimeson-type wave function in Eq. ( 19) and RGM.The relevant theoretical dimeson thresholds are aligned to the physical ones.Based the above comparisons, it is evident that the GEM is superior than the RGM and DMC method.With this in mind, we will utilize the results obtained through the GEM to conduct a comprehensive comparison across the three distinct quark potential models.
The experimental T cc (3875) + state is the candidate of J P = 1 + [ccnn] I=0 tetraquark state.Notably, both the PCQM and χCQM yield the bound state solutions, which can be viewed as the predictions made by these quark models prior to the experimental observation.In the original work of AL1 and AP1 models [57], the existence of the bound states was not conclusively established due to the limitations in computational resources at the time.Meanwhile, the results from AL1 and AP1 are more consistent with the experimental result of T cc (3875) + , a very loosely bound state.Conversely, the SLM model suggests the presence of a compact tetraquark state well below the DD * threshold, with a substantial binding energy of approximately 200 MeV.In Ref. [103], employing the same SLM, a loosely bound state was achieved through the RGM.The result stems from the constrained basis wave function in RGM.Before the experimental observations, investigations of the J P = 1 + [ccnn] I=0 tetraquark state were undertaken in lattice QCD simulations in Refs.[100,104,105], but the existence of this state remained inconclusive.Subsequent to the experimental observations, lattice QCD simulations based on Lüscher's method [106] and the potential method (HAL QCD method) [107] reported virtual states in this channel.
In all three quark models, we consistently obtain very deep bound state of J P = 1 + [bbnn] I=0 , which is the possible heavy quark flavor symmetry partner of T cc state.The three quark models indicate that both the ground states and the first excited states of the system are the bound states.It is worth noticing that the binding energy obtained in SLM is still much larger than those from AL1 and AP1 models.This state has also been extensively investigated in lattice QCD simulations [100,[108][109][110], which consistently establish the existence of the bound states with the binding energies about 100-200 MeV.Additionally, the static potentials from lattice QCD [111][112][113][114] also indicate the existence of the deeply bound T bb states.Furthermore, these three distinct models also predict the existence of J P = 1 + [bcnn] I=0 bound states.It is worthwhile to notice that the result from SLM is significantly deeper than those from the AP1 and AL1 models.For this state, the lattice QCD simulations have not provided a consistent conclusion [102,105,115,116].
In our calculations, there is no isovector bound state for the QQnn systems.However, in Ref. [103], a J P = 1 + [bbnn] I=1 bound states were obtained in SLM within the RGM framework.Apparently, this result conflicts with our calculation.
The absence of an isovector bound state in I=1 system is widely acknowledged and is a consensus reflected in the majority of pertinent publications, such as Ref. [61,117,118].
Another state worth attention is the J = 1 + bbns state.We obtain the bound state below the Bs B * threshold with the binding energy 30-70 MeV.In lattice QCD simulations, the existence of this bound state was also implied [100,108,115,119].
In addition to the above bound states existing consistently in three different quark models, the SLM model also predicts a bcns bound state, which is absent in AL1 and AP1 models.
2. J P = 0 + and J P = 2 + For the doubly heavy tetraquark states with J P = 0 + and 2 + , the PCQM and χCQM both predict the [bcnn] I=0 bound states.The main difference is the J P = 0 + state from SLM is much lower than those from AL1 and AP1.Additionally, the SLM model predicts the extra bcns states for J P = 0 + and 2 + .The AP1 model predicts the extra isovector tensor bbnn bound state.

C. Singly heavy tetraquark states
For the singly heavy quark states, we investigate the J P = 0 + , 1 + and 2 + system with the following quark contents, [bnsn] I=1 , [bsnn] I=0,1 , bssn , bnss , where the [Qnsn] I=0 and [Qnnn] I=1/2 are not considered because we only focus on the manifestly exotic states.We use the GEM and RGM to solve the spectra in three models.The results are presented in Table IV and Fig. 7.
Comparing the results from GEM and RGM, it is apparent that the results from RGM could be biased due to the constrained basis functions.We still employ the GEM to compare the results from different quark potential models.
One can see that three quark models all predict the J P = 0 + [csnn] I=0 , J P = 0 + , 1 + [bsnn] I=0 bound states.For all the above results, the predictions from SLM tend to be much deeper.It should be noticed that the J P = 0 + [csnn] I=0 is irrelevant to the experimental T cs0 (2900) state which is a resonance state close to the D * K * threshold.Additionally, the SLM model predicts several extra bound states, J P = 1 + [csnn] I=0 state, J P = 2 + cssn state, J P = 2 + [bsnn] I=0 state, and J P = 2 + bssn state.
In Ref. [120], the authors investigated the spin-0 and spin-1 csnn and cnsn states in SLM and obtained no bound solution, which conflicts with our J P = 0 + , 1 + [csnn] I=0 bound states.In their calculations, several virtual states are found.7. Bound states of the Qq q q systems from GEM and RGM in three quark models.The relevant theoretical dimeson thresholds are aligned to the physical ones.TABLE IV.Bound states of the Qq q q systems.The notations are the same as those in Table III

IV. DISCUSSIONS AND SUMMARY
In this study, we conduct benchmark test calculations to investigate the tetraquark bound states that are manifestly exotic across three distinct quark models (AL1, AP1, and SLM) and employ three different few-body methods (GEM, RGM, and DMC).In the GEM calculations, we find the results is insensitive to the choice of the discrete wave functions once they are complete.However, when it comes to the spatial wave functions, the inclusion of both diquark-antidiquark and dimeson types is imperative for obtaining precise solutions.We use the GEM with the general wave functions as the standard to compare with the results from the RGM and DMC.
Our results show that the GEM is superior than RGM and DMC by obtaining the more lower ground state energies and detecting extra bound solutions.Because GEM is a method based on the variational principle, the lower ground state solutions mean more precise results.While RGM can identify the molecular-type states, it tends to produce biased results for the deeply bound states due to its constraints on the basis functions, particularly the dimeson-type wave functions.DMC can identify both the compact tetraquark states and loosely bound molecular states without a priori assumption about the clustering behaviors of the wave functions.However, when compared to GEM, the current implementation of DMC falls slightly short in terms of precision.This could be attributed to the possibility that the importance functions used in Ref. [73] may not have been optimized to suit the multiquark systems.Up to now, the advantages of the DMC method, as demonstrated in atomic and molecular physics as well as nuclear physics, has not been fully exploited in the multiquark systems.DMC remains a promising approach distinct from the variational method.The DMC method circumvents the challenges associated with the exponentially growing basis set with number of particles and the intricate integrals tied to few-body forces in the variational approach.Its advantages are likely to become apparent in the case of multiquark states with a large number of quarks [92] and in systems featuring flux-tube few-body potentials [72].It holds potential for further development and refinement.
In Table III, we also consider the potential excited states [bbnn] I=0 J P =1 + .The basis expansion method can be extended to encompass these excited states, although the Ritz theorem in quantum mechanics textbooks typically emphasizes the variational method's role in establishing an upper limit for ground states.The expectation value of the Hamiltonian remains stationary in the vicinity of its discrete eigenvalues.Consequently, with an increase in the dimensions of the basis space, the eigenvalues within this space converge towards the exact solutions of the Hamiltonian, encompassing both ground and excited states.Further details can be found in Ref. [121].Alternatively, the Diffusion Monte Carlo (DMC) method can be extended to identify excited states as well (see [122]) We use the results from GEM to compare the results from three quark models belonging to two types, PCQM (AL1 and AP1) with the minimal OGE interaction and confinement interaction, and χCQM (SLM) with both two interactions and extra OBE interaction.Our analysis reveals that the SLM models tend to yield more deeper bound states or provide extra bound state solutions than the AL1 and AP1 models.
We perform the calculations for over 150 tetraquark states for the fully, triply, doubly and singly heavy tetraquark systems with J P = 0 + , 1 + and 2 + .We summarize the bound states existing in all three quark potential models in Table V.We find there are no fully heavy and triply heavy tetraquark bound states.For the doubly heavy tetraquark systems, we find the [ccnn] I=0 J P =1 + (candidate state of the experimental T cc (3875) + state), and its heavy quark flavor symmetry partners [bbnn] I=0 J P =1 + and [bcnn] I=0 J P =1 + are all bound states.In addition, the [bcnn] I=0 J P =0 + ,2 + and [bbns] J P =1 + are also bound states.It is worthwhile to mention that the [bbns] J P =1 + bound state was also supported by the lattice QCD simulations.For the singly heavy systems, we find the [bsnn] I=0 J P =0 + ,1 + and [csnn] I=0 J P =0 + bound states.We hope these stable states against the strong decays may be searched for in the experiments.
We find the SLM quark potential models tend to give more and deeper bound states and the RGM tends to underestimate the bind energy.However, it is still irrational to discard them.On the one hand, we still have the room to refine the parameters in SLM to fit the experiment results.Meanwhile, the SLM was originally proposed to depict the NN scattering phase [69] via RGM with presuming di-baryon clustering behaviors.If one adopts a general trial wave functions of the six quarks, one perhaps obtains quite different solutions of SLM (e.g.deep bound states) instead of the NN scattering states or deuteron.In other words, the SLM model and RGM were used in combination at the birth of this quark potential model.Perhaps, they should still be used in combination.The combination of the SLM and RGM provides a loosely bound solutions for [ccnn] I=0 J=1 , which is consistent with the experimental T cc (3875) + state.If that is the case, it is still reasonable to use the SLM model when it is assumed in advance that the target state is a molecular state.
In fact, we do have some hints that the mixing effect between the molecular configurations and the diquarkantidiquark configurations is suppressed.For example, in the flux-tube models, one can model the complicated dynamics of the sea quarks and gluons as the flux-tubes.There could be the dimeson-type flux-tube and the diquark-antidiquarktype (butterfly-type) flux-tube.When considering the possible mixture of the dimeson constructions and diquark-antidiquark constructions, in addition to the valence quark wave functions, one has to consider the flux-tube wave functions, which represent the dynamics of the sea quarks and gluons.The small overlap of the flux-tube wave functions for the different configurations could suppress their mixing effect.In the weak mixing limit, we could get the molecular states without the effect from the diquark-antidiquark configurations.For these states, it is reasonable to use the dimeson basis functions to expand the wave functions like RGM.One can find similar discussion in Ref. [73,81,123].One can also assume the SLM is a quark potential model only working for the molecular configurations.In this way, the combination of the RGM and SLM becomes reasonable.
In our benchmark test calculations, we unveil discrepancies among the various quark potential models and few-body methods on the market.Furthermore, when we compare our results with those in employing the same quark potential models, we continue to encounter inconsistencies.While some of these inconsistencies may be attributed to limitations in computational precision in earlier years, others remain unexplained.As we are rapidly entering the era of the "genuine" multiquark states, it becomes increasingly vital to conduct additional benchmark tests of quark model calculations, particularly when involving different research groups.In the variational method based on basis expansion, the final step involves solving the generalized eigenvalue problem: Here, H and N represent the Hamiltonian matrix and the overlap matrix, respectively.λ and v correspond to the eigenvalue and eigenvector, and |i⟩ and |j⟩ denote the basis states, which may not be orthogonal.In a more general context, basis states may exhibit dependencies, either in a rigorous sense, such as complete basis states becoming linearly dependent after antisymmetrization, or in a less strict sense, where almost "parallel" basis states are considered linearly dependent, taking into account the machine precision truncation error.The presence of (nearly) dependent bases can lead to ill-conditioned matrices, breaking down algorithms designed to solve the general eigenvalue problem.Therefore, for robust results, it is crucial to eliminate (nearly) redundant bases.
To this end, various methods can be employed, including the Gram-Schmidt process.In our calculations, we opt for the diagonalization of the overlap matrix N, In practical applications, one can establish a tolerance for the eigenvalues of N and eliminate basis states associated with very small eigenvalues.This approach effectively removes (nearly) redundant bases, contributing to the robustness of the algorithm.

80 FIG. 5 . 8 FIG. 6 .
FIG.5.Bound states of the QQq q systems from GEM, RGM and DMC in three quark models.The relevant theoretical dimeson thresholds are aligned to the physical ones.

T
TNT † = D, where the transformation matrix are decomposed to two blocks.The total number of bases is denoted as n = a + b, with only a of them being linearly independent.The matrix D is diagonal, and it is apparent that N is semi-positive-definite, ensuring that the diagonal elements of D are all positive.The matrix N becomes positive-definite if and only if the basis states |i⟩ are linearly independent.To eliminate the (nearly) redundant basis vectors, we introduce a set of new bases |β⟩ with a number of a, αi ⟨i|j⟩T † jβ = d α δ αβ , (A4) where d α is the elements of matrix D. One can further make the set of bases normalized, | β⟩ = |β⟩ d β .One can get the matrix elements of the Hamiltonian under the orthogonal and normalized bases,

TABLE I .
[86] spectra of the ground state mesons from three different quark models (in units of MeV).The "Exp."represents the experimental results[86]as a comparison.
A is the most general basis.Ψ B and Ψ C are basis wave functions with general spatial wave function but dimeson and diquark-antidiquark discrete wave functions, respectively.The discrete wave functions of Ψ D and Ψ E are general but with dimeson and diquarkantidiquark spatial wave functions, respectively.We use [ccnn] I=0 J P =1 + , [bbnn] I=0 J P =1 + and [bcnn] I=0 J P =2 + as three examples to test different choices of the basis wave functions.The results are presented in II.It should be noticed that the basis wave functions which yield a lower ground state are more precise according to the variational principle.Our numerical results show that using different discrete basis functions makes little difference once they are complete.Meanwhile, including both the dimeson and diquark-antidiquark spatial functions is very important.Otherwise, one can obtain biased results.For example, neglecting the dimeson spatial wave functions makes the bounded [bcnn] I=0 J P =2 + states unbound in the result of Ψ E .The absence of the diquarkantidiquark spatial wave functions in Ψ D makes the deeply bound [ccnn] I=0 J P =1 + state in SLM model a loosely bound state.In our final results, we choose the most general Ψ A as our basis functions of GEM.

TABLE II .
Comparisons of the GEM results in different basis functions in Eqs.

TABLE III .
Bound states of the QQq q systems.The energies are with respect to the lowest relevant thresholds (in units of MeV)."NB" represents that there is no bound solution."..." labels systems that are not investigated in the literature." →Excited" labels the excited bound state solutions.of the identical Fermions and the associated sign problem should be addressed through proper choices of the importance functions, thereby opening up avenues for future improvements. metrization .

TABLE V .
Final results of the bound states that exist in all three quark potential models.