Doubly heavy tetraquark states in the constituent quark model using diffusion Monte Carlo method

We use the diffusion Monte Carlo method to calculate the doubly heavy tetraquark $T_{cc}$ system in two kinds of constituent quark models, the pure constituent quark model AL1/AP1 and the chiral constituent quark model. When the discrete configurations are complete and no spatial clustering is preseted, the AL1/AP1 model gives an energy of $T_{cc}$ close to the $DD^*$ threshold, and the chiral constituent quark model yields a deeply bound state. We further calculate all doubly heavy tetraquark systems with $J^P=0^+,1^+,2^+$, and provide the binding energies of systems with bound states. The $I(J^P)=0(0^+)$ $bc\bar{n}\bar{n}$, $0(1^+)$ $bb\bar{n}\bar{n}$, $0(1^+)$ $bc\bar{n}\bar{n}$, $\frac{1}{2}(1^+)$ $bb\bar{s}\bar{n}$ systems have bound states in all three models. Since the DMC method has almost no restriction on the spatial part, the resulting bound states have greater binding energies than those obtained in previous works.

A minimal constituent quark model used to describe multiquark states typically incorporates the one-gluon-exchange interaction and the confinement effect.Some quark models may also encompass additional interactions, such as one-boson exchange (OBE) interactions.They include pseudoscalar meson exchange interactions stemming from the spontaneous break-ing of chiral symmetry, vector meson exchange interactions, and scalar meson exchange interactions.Based on the inclusion or exclusion of the OBE interaction, constituent quark models can be broadly categorized into two types: the pure constituent quark model (PCQM) without OBE and the chiral constituent quark model (χCQM) with OBE.One of the objectives of this study is to investigate the difference between these two types of quark model interactions concerning their predictions of the doubly heavy tetraquark states.
With the quark interactions, one can solve the four-body Schrödinger equation based on the variational method via the basis expansions.Apparently, the precision of the solution depends on the choice of trial functions or basis functions.In the practical calculation, the diquark-antidiquark structure and the dimeson structure are two widely used structures to investigate tetraquark states.The diquark-antidiquark structure has a color wave function where two quarks (antiquarks) combine first and then form a color singlet together, namely [(qq) 3c (q q) 3c ] 1c and [(qq) 6c (q q)6 c ] 1c .In the spatial part, one can choose the corresponding Jacobi coordinates to depict the correlation as shown in Fig. 1, where the diquark and antidiquark form two clusters separately and then form the tetraquark state.The dimeson structure usually restricts the color wave function to be [(q q) 1c (q q) 1c ] 1c .It should be noticed that even if the [(q q) 8c (q q) 8c ] 1c is not included, the color bases with two dimeson structures by exchanging (anti)quarks is complete, though they are not orthogonal.In principle, one could choose different structures for wave functions of discrete quantum numbers and spatial wave functions.For example, the spatial wave functions are the dimeson structure, while allowing the color part to contain both [(q q) 1c (q q) 1c ] 1c and [(q q) 8c (q q) 8c ] 1c .If it is such a case, we will specify it explicitly.
In previous investigations of the ccnn system with (I, J) = (0, 1), it has been noted that the choice of different structural configurations yields disparate outcomes.For a PCQM, Yang et al. obtained a shallow bound state under the dimeson structure, while they found no bound state under the diquark-antidiquark structure [12].For the χCQM, they obtained a shallow bound state under the dimeson structure and a deeply bound state under the diquark-antidiquark structure.Similar discrepancy was also observed in other χCQM Two structures for the tetraquark system.(a) diquarkantidiquark structure.(b) dimeson structure.
works [16,47,50,52].In other words, to accurately determine the ground state energies of the multiquark states using a variational method-based approach, one must either conceive a good conjecture about the clustering behavior of the wave functions or employ very general trial functions.Alternatively, the diffusion Monte Carlo (DMC) method is a promising technique to precisely and efficiently ascertain ground state energies, without the need for a priori assignment of clustering behaviors.
DMC has been well used in the simulations of the molecular physics [70], solid physics [71], and nuclear physics [72].In hadronic physics, the DMC method has been used in quark models in several works [73][74][75][76][77][78][79].In the DMC formalism, the distribution of the so-called walkers is used to represent the spatial wave function, which gradually evolves to the ground state over time in principle.The DMC approach does not select a specific basis for the spatial wave function, resulting in a general wave function space.This approach could be more suitable for exploring the discrepancy caused by predetermined clustering structures.
This paper is arranged as follows.In Sec.II, the diffusion Monte Carlo method is introduced.In Sec.III, the Hamiltonians with different types of potential models are presented.In Sec.IV, the construction of the wave function for all doubly heavy tetraquark systems with J P = 0 + , 1 + , 2 + is presented.In Sec.V, we give the numerical results for the T cc system in different models, and the bound states of all other doubly heavy tetraquark systems with J P = 0 + , 1 + , 2 + .Finally, a brief discussion and summary are given in Sec.VI.

II. DIFFUSION MONTE CARLO METHOD
We adhere to the formalism presented in Ref. [78].The DMC algorithm can be derived from the imaginary-time Schrödinger equation (in natural units ℏ = c = 1) and its subsequent solution.
where R ≡ (r 1 , r 2 , ..., r m ) represents the positions of m particles.E R is the shift of energy, and Φ i (R) is the eigenfunctions.When E R is taken close to the ground state energy E 0 , the wave function Ψ(R, t) will approach Φ 0 (R) after a sufficiently long time evolution (as long as c 0 is not too small), and other components will be suppressed exponentially [80].The Eq. ( 1) has a solution in the form of path integral, where Theoretically, this form of the solution can be readily implemented using an algorithm.In this algorithm, a substantial number of walkers sample the wave function Ψ(R, t) and the distribution of walkers represents the amplitude of Ψ(R, t).However, this algorithm is usually unstable due to the drastic statistic fluctuation.We need to further use the importance sampling technique to reduce the fluctuation [81].In this approach, rather than directly sampling the wave function Ψ(R, t), we sample a newly defined function, denoted as where ψ T (R) is a time-independent importance function, which should be chosen as close to Φ 0 (R) as possible.In this work, we use the form following Ref.[74], Here a ij are adjustable constants and their values are set to minimize the fluctuation.The Schrödinger equation with importance sampling reads where is the drift force acting on particle i.The path integral solution of Eq. ( 8) is with The DMC algorithm employed to obtain the solution described above, along with the formalism utilized to address the coupled-channel, has been comprehensively elucidated in Ref. [78].

III. HAMILTONIAN
The nonrelativistic Hamiltonian of a four-quark system reads where m i and p i are the mass and momentum of quark i. T CM is the center-of-mass kinematic energy, which automatically vanishes in the evolution, because the system will tend to the lowest energy state.
In this study, we will employ two types of non-relativistic quark models: pure constituent quark model (PCQM) and chiral constituent quark model (χCQM).The primary distinction between these two types is that the χCQM incorporates the OBE between quarks, a feature not included in the PCQM.
For the PCQM, we adopt AL1 and AP1 quark models [8,82], which contain the one-gluon-exchange (OGE) interaction and the confinement effect.The potential reads where The power parameter p of the confinements are 1 and 2/3 for AL1 and AP1, respectively.λ c represents the SU(3) color Gell-Mann matrix.s i is the spin operator of quark i.The parameters of the two models are taken from Ref. [82].There are some similar models, such as Barnes, Godefrey and Swanson model [83].However, the parameters for the light quarks were not determined.
For the χCQM, we refer to the model introduced in Ref. [84,85] as the Salamanca model (SLM), which encompasses the OBE interaction, OGE interaction, and a confinement potential with a screening effect, In this work, only the central term is considered.
The OGE potential includes a Coulomb term and a spindependent color-magnetic term, where r 0 (µ ij ) = r0 µij , µ ij is the reduced mass of quark i and j. λ c represents the SU(3) color Gell-Mann matrix.The effective scale-dependent strong coupling constant α ij s is given by The confinement part is This screened form is introduced to simulate the string breaking effect [84].It exhibits a linear behavior at short distances and becomes constant at large distances.The meson exchange interactions only occur between the light quarks q = u, d, s.The OBE interactions read where x. r ij is the relative distance between quark i and j. σ i is the spin Pauli matrix for quark i. λ a represents the SU(3) flavor Gell-Mann matrix.
For interaction between q i qj , λ a i • λ a j should be replaced by λ a i • λ a j T .m π , m K , and m η are the pseudoscalar meson masses, and as suggested by Ref. [86].The angle θ P is introduced to obtain the physical η meson instead of the flavor octet state η 8 .The chiral coupling constant g ch can be determined from the πN N coupling constant through The parameters of SLM listed in Table I are taken from Ref. [85], which are fitted over all meson spectra.The calculated meson masses in AL1, AP1 and SLM models are shown in Table II.It can be seen that these potential models can describe the experimental meson spectra well.

IV. WAVE FUNCTION CONSTRUCTION
The experimental T + cc (3875) [42,43] is a candidate of the ccnn (n=u and d) tetraquark state with quantum number (I)J P = (0)1 + .We sequentially label the (anti-)quarks of ccnn as 1, 2, 3, 4. The two heavy quarks, designated as 1 and 2, are indistinguishable, as are the two light quarks, labeled as 3 and 4. The construction of the wave function needs to fulfill the Pauli principle.For I = 0, the flavor wave function is symmetric for cc, while antisymmetric for nn.The remaining color-spin-spatial wave functions allowed by the Pauli principle are where S represents exchange symmetric, and A represents exchange antisymmetric.For the system with isospin I = 1, the flavor wave function for nn is symmetric.The remaining color-spin-spatial wave functions need to be changed accordingly |T In addition to the T cc states, we will calculate the ground states of all doubly heavy tetraquark systems with J P = 0 + , 1 + , 2 + .The two heavy quarks could be bb, cc and bc.The two light antiquarks could be nn, ss and sn.In Table III, we list them specifically and their corresponding dimeson thresholds.We use the same methods to construct the wave functions satisfying the Pauli principle.The spin-color configu-

System
I Thresholds ration channels included in our calculation are listed in Table IV.For each channel, the exchange symmetry of its spatial part should be determined according to the quark composition (and isospin) of the system.In contrast to other frameworks, such as those discussed in Refs.[88,89], the calculations do not presume the clustering behavior of quarks.

V. NUMERICAL RESULTS
In our simulation, we employ 1 × 10 4 walkers to implement the DMC algorithm.The ensemble undergoes 2 × 10 4 steps, each with a time increment of ∆t = 0.01GeV −1 , to ensure stability.The resulting energy is averaged over the last 15000 steps to mitigate fluctuations.To estimate statistical uncertainty, we consider correlations among adjacent steps.To achieve this, we divide the steps into blocks and calculate block averages.We observe that for the four-quark systems, the blocks become uncorrelated when the block size reaches 1000 steps.Consequently, we partition the last 15000 steps into 15 blocks, each consisting of 1000 steps.We then employ the Jackknife resampling method [90] to calculate the uncertainty, which is found to be less than 1 MeV.For more details about the statistical uncertainty analysis, refer to Ref. [78].

A. Pure Constituent Quark Model
The calculation is performed using the AL1 and AP1 potentials.We first give the results of the (I)J P = (0)1 + system, which is the candidate of the experimental T cc state.A complete set of six spin-color-spatial configurations |T 0 1 ⟩, |T 0 2 ⟩, ..., |T 0 6 ⟩ is included.The results are shown in Table V.The binding energy is zero and -1 MeV for the AL1 and AP1 potentials respectively.With the current level of uncertainty (less than 1 MeV), it is still inconclusive whether they are below threshold or not.But the results indicate that if they exist, the binding energies are likely to be very small.
We further calculate the ground states of all doubly heavy tetraquark systems with J P = 0 + , 1 + , 2 + .The systems with bound state solutions in AL1 and AP1 models are shown in Table V.One can identify the unbound system by comparing the Tables V and III.There could be resonance solutions for these systems, which will not be investigated in this work.
Comparing our results with those obtained in the same potential models using the variational approach in Ref. [8], we can find that most of our masses are lower than theirs, and more bound states are obtained.In Ref. [18], Meng et al. adopted a tuned AP1 model, using Gaussian expansion method (GEM) [91] and obtained a binding energy of −23 MeV.Due to the differences in parameters, only a qualitative comparison can be made.Our results are generally consistent with theirs.In Ref. [12], Yang et al. adopted a quark model which is very similar to the AL1 and AP1 models.
Their results showed that the (I, J) = (0, 1) ccnn system is unbound under the diquark-antidiquark structure and bound with a binding energy of −1.8 MeV under the dimeson structure.If the mixing of these two structures is taken into account, a deeper bound state of −23.7 MeV is obtained.For the (I, J) = (0, 1) bbnn system, the diquark-antidiquark and dimeson structures give the binding energies of −120.9MeV and −11.5 MeV respectively, and their mixing gives a deeper −160.1 MeV binding energy.Our results are basically consistent with their mixing picture with a loosely bound ccnn state and a deeply bound bbnn state.This indicates that DMC is effective even without the priori assignment of the clustering behaviors.

B. Chiral Constituent Quark Model
As before, we first give the result of the (I)J P = (0)1 + T cc system.The calculation is performed using the SLM potential.A complete set of six spin-color-spatial configurations |T 0 1 ⟩, |T 0 2 ⟩, ..., |T 0 6 ⟩ is included.The result is shown in Table V.It can be seen that the SLM yields a deeply bound T cc with a binding energy of −156 MeV.This indicates that once the spin-color configurations are complete and the spatial wave function has no presumed clustering, the χCQM will give exactly a deeply bound state instead of a shallow one near the threshold.
Using the SLM, we further calculate the ground states of all doubly heavy tetraquark systems with J P = 0 + , 1 + , 2 + .Their corresponding dimeson thresholds are listed in Table III.The bound states solutions are shown in Table V.By comparing Tables III and V, one can identify the unbound systems.
In Ref. [92], Ortega et al. used the same SLM model and got a very shallow bound 0(1 + ) ccnn state as listed in Table V.This is perhaps due to the Resonating Group Method (RGM) they used.In this method, the meson wave functions are calculated in advance and input inside directly, where only the S-wave dimeson structure is considered.The higher partial wave component inside the q q pair and the distortion effect of the meson inside the tetraquark systems are neglected.This significantly restricts the wave function space and shifts a deeply bound state to the vicinity of the threshold.In the DMC calculations in this work, the spatial wave function is unclustered with no assumption of either diquark-antidiquark or dimeson structure, and all walkers diffuse freely in space.Therefore we can naturally obtain a deeply bound state.Similarly, their T bb bound state with quantum number I(J P ) = 0(1 + ) is also shallower than ours.It should be noted that the corresponding threshold of the 0(0 + ) bbnn system actually cannot be B B, which is forbidden because of the Bose symmetry. 1 In Ref. [52], Deng et al. used a similar χCQM.In the dimeson structure with only the 1 c ×1 c color configuration, a bound state with a binding energy of only −0.34 MeV was found.Even after adding the 8 c × 8 c color configuration, it remains a shallow one.For comparison, we also perform a calculation using the same model and obtain a deep binding energy of −110 MeV.The difference arises from the dimeson S-wave constraint in their variational approach.So essentially, it is for the same reason as described earlier about RGM.The same situation happens in Ref. [12].Two shallow bound states with binding energies of −0.6 and −0.2 MeV are obtained using two χCQMs under the dimeson structure.
In Ref. [52], the results of the diquark-antidiquark structure are also provided.They included two spin-color-space configurations |T 0 1 ⟩ and |T 0 2 ⟩, and the resulting binding energy is −60 MeV, shallower than our −110 MeV.Comparing with the configurations |T 0 1 ⟩, |T 0 2 ⟩, ..., |T 0 6 ⟩ included in our calculation, they only considered the two channels corresponding to the spatial diquark-antidiquark S-wave wavefunction.However, the ground state does have higher partial wave components as explained earlier, and the remaining four channels with higher partial wave components should be included.In Ref. [12], Yang et al. also gives the results under the diquark-antidiquark structure and under the mixing of the dimeson and diquark-antidiquark structures.Deeply bound states of −142.4MeV and −202.7 MeV are obtained respectively.The mixing one is deeper because of the larger basis function space. 1 The B meson is a boson with I(J P ) = 1 2 (0 − ).The total wave function of two identical B mesons has to be exchange symmetric.When constructing a I(J P ) = 0(0 + ) B B total wave function, the spin part 0s ⊗ 0s → 0s is symmetric.The isospin part 1 2 I ⊗ 1 2 I → 0 I is antisymmetric.To satisfy the exchange symmetry of two boson B B, the remaining spatial wave function has to be exchange antisymmetric, i.e., the orbital angular momentum quantum number l is odd.However, the parity P becomes (−1) l = −1.Therefore the 0(0 + ) B B system is forbidden.
There are some other works using the GEM under the similar χCQM.Their results are qualitatively consistent with ours.In Ref. [10], Vijande et al. used the same form of the potential as SLM, but with different parameters.They considered the spatial diquark-antidiquark structure and obtained a −129 MeV ccnn bound state and a −341 MeV bbnn bound state, both with (S, I) = (1, 0).In Ref. [17], Tan et al. included both the dimeson and diquark-antidiquark structures and obtained deeply bound ccnn and bbnn states.

VI. DISCUSSION AND SUMMARY
We perform a diffusion Monte Carlo (DMC) calculation of the I(J P ) = 0(1 + ) T cc system and other doubly heavy tetraquark systems with J P = 0 + , 1 + , 2 + .We use two kinds of potential models, the PCQM (specifically AL1 and AP1) that contain the OGE and confinement interactions, as well as the χCQM (specifically SLM) containing the OGE, confinement, and OBE potentials.
We include a complete set of six spin-color configurations to perform the DMC evolution with no spatial clustering assumed.For the PCQM, we are not certain if there are bound T cc state due to the limitation of uncertainty.But if it exists, the binding energy will be small.We further calculate the ground states of all doubly heavy tetraquark systems with J P = 0 + , 1 + , 2 + and give the binding energies of the systems with bound state solutions.Among them, the I(J P ) = 0(0 + ) bcnn, 0(1 + ) bbnn, 0(1 + ) bcnn, 1  2 (1 + ) bbsn, and 0(2 + ) bcnn systems are found to have bound states.Our results are basically consistent with the mixing picture ones using the GEM [12], indicating that the DMC method can effectively avoid the difference in results caused by clustering in GEM.
For the χCQM case, the results show that when the colorspin-isospin configurations are complete and no spatial clustering is assumed, the SLM yields a deeply bound state with a binding energy of −156 MeV.This result is qualitatively consistent with the works using diquark-antidiquark structure or the mixing of dimeson and diquark-antidiquark structure in GEM.The appearance of the shallow bound states in some previous works are due to the restricted basis space.While in the DMC calculations in this work, the spatial wave function is unconstrained with no assumption of diquark-antidiquark or dimeson structure.Therefore we can obtain a deeply bound state that faithfully presents the solution of the χCQM.
The inability to directly obtain the experimentally observed shallow bound state T + cc (3875) under the χCQM may be attributed to several reasons.The first possible one is that the parameters of the models are obtained by fitting the hadron spectrum and the hadron-hadron scattering data.So extending them directly to the tetraquark systems may no longer be applicable.The original motivation of the SLM is aimed to depict the NN scattering phase [94] via RGM with presuming di-baryon clustering behaviors.If one incorporates the more general correlation of the six quarks, one perhaps obtains quite different solutions (e.g.deep bound states) instead of the NN scattering states or deuteron.The debate over whether the PCQM or χCQM is more suitable for depicting NN scattering phase shifts has persisted for many years, as exemplified in [95], but without arriving at a definitive conclusion.The preceding analysis, grounded on a single experimental result involving T cc , does not imply the superiority of the PCQM over the χCQM.We still have substantial room for refining the parameters of a chiral quark model that can effectively describe the tetraquark states.
Secondly, the calculations at hadron level or under the dimeson structure usually result in shallow bound states, which indicates that there may be other mechanisms that favor the dimeson structure rather than all structures being equally likely to participate in mixing.One can understand it through the Born-Oppenheimer approximation.The potential in constituent quark model should be regarded as the Born-Oppenheimer (BO) energy in a scenario with valence quarks as slow degrees of freedom, and sea quarks and gluons as the fast degrees of freedom.Under the BO approximation, one can first fix the positions of the valence quarks and calculate the energy of the systems, so-called BO energy, which should depend on the position of valence quarks.In the second step, the BO energy is treated as the potential of valence quarks, which is the interaction appearing in the constituent quark models.For valence quarks with different color configurations, the complicated dynamics of sea quarks and gluons are quite different, which are simplified as the flux-tubes with different topological structures.For example the di-meson configuration in Fig. 1 should be |1⟩⟩ = |[(q q) 1c (q q) 1c ] 1c ⟩|FT 1−1 ⟩, where |[(q q) 1c (q q) 1c ] 1c ⟩ and |FT 1−1 ⟩ represent the wave functions of valence quarks (slow degrees of freedom) and flux-tube (fast degrees of freedom) respectively.The notation |...⟩⟩ is introduced to discriminate with the naive valence quark states.In this picture, the matrix elements of |1⟩⟩, |1 ′ ⟩⟩, |3⟩⟩ and |6⟩⟩ will depend on both the valence quark wave function and flux-tube wave functions, where |1 ′ ⟩⟩ represent another dimeson state by exchanging quarks.It is possible that the |1⟩⟩ and |1 ′ ⟩⟩ are favored in the tetraquark states considering the dynamics of the fast degrees of freedom [96].In this way, it is reasonable to incorporate the dimeson configurations in tetraquark states exclusively.One can find a similar picture in Ref. [56].
The recent experimental observation of tetraquark states provides us with a valuable opportunity to critically examine the popular quark models on the market and further advance them, where the diffusion Monte Carlo method without the preassignment clustering shall play a pivotal role.

a
D1D)P /(D * 0 D * )P a DD * (D1D)P /(D * 0 D * )P There is no relevant S-wave thresholds of two ground state mesons in this quantum numbers considering the exchange symmetry of two identical bosons.The subscript represents the relative angular momentum is P-wave.

TABLE II .
[87]masses (in MeV) of mesons in AL1, AP1 and SLM models, compared with the experimental results taken from Ref.[87].The experimental results are averaged over the isospin multiples.

TABLE III .
The doubly heavy tetraquark states explored in this study and their associated dimeson thresholds.

TABLE V .
Doubly heavy tetraquark systems with bound state solutions in AL1, AP1, and SLM models.I is the isospin.Eth , E and ∆E are the corresponding lowest dimeson threshold, system energy, and binding energy (in MeV) respectively in our calculations."NB" represents no bound solution.We also present the ∆E in literature using the same interactions."..." represents it was not calculated in the corresponding work.

TABLE VI .
The proportion of configurations of bound states in AL1 model.I is the isospin.E and ∆E are the system energy and binding energy (in MeV).

TABLE VII .
The proportion of configurations of bound states in SLM model.I is the isospin.E and ∆E are the system energy and binding energy (in MeV).

TABLE VIII .
The rms radii for bound states in AL1 model.I is the isospin.E and ∆E are the system energy and binding energy (in MeV).» ⟨r 2 ij ⟩ is the rms radius between i, j (in fm).C and M represent the compact tetraquark and molecular type respectively.

TABLE IX .
The rms radii for bound states in SLM model.I is the isospin.E and ∆E are the system energy and binding energy (in MeV).» ⟨r 2 ij ⟩ is the rms radius between i, j (in fm).C and M represent the compact tetraquark and molecular type respectively.