The Doubly-heavy Tetraquarks ($qq'\bar{Q}\bar{Q'}$) in a Constituent Quark Model with a Complete Set of Harmonic Oscillator Bases

We have improved our previous variational method based constituent quark model by introducing a complete set of 3-dimensional harmonic oscillator bases as the spatial part of the total wave function. To assess the validity of our approach, we compared the binding energy, thus calculated with the exact value for the hydrogen model. After fitting to the masses of the ground state hadrons, we apply our new method to analyze the doubly-heavy tetraquark states $qq'\bar{Q}\bar{Q'}$ and compared the result for the binding energies with that from other works. We also calculated the ground state masses of $T_{sc} (ud\bar{s}\bar{c})$ and $T_{sb} (ud\bar{s}\bar{b})$ with $(I,S) = (0,1), (0,2)$. We found that $T_{bb} (ud\bar{b}\bar{b})$ and $us\bar{b}\bar{b}$, both with $(I,S) = (0,1)$, are stable against the two lowest threshold meson states with binding energies $-145$ MeV and $-42$ MeV, respectively. We further found that $T_{cb} (ud\bar{c}\bar{b})$ is near the lowest threshold. The spatial sizes for the tetraquarks are also discussed.


I. INTRODUCTION
Since the observation of X(3872) [1] and several exotic hadron candidates that followed, the structure of these particles and other potential flavor exotic configurations have become a central theme of study. Theoretical approaches on these topics range from direct lattice calculation [2], sum rule approach [3], effective models and constituent quark based models [4]. To better understand the data and to point to future searches, it becomes necessary to improve any simple minded model calculations to full details.
Quark model based on color and spin interaction was successful in describing the mass differences between hadrons within a multiplet [5]. A constituent quark model with a more realistic potential was proposed in Ref. [6] in the early 1980s, which gives a unified description of meson and baryon spectra. The potential in BCN model [6] was applied to investigate the baryon spectra in detail [7], and to various tetraquark states [8]. In Ref. [9], the authors introduced four different types of interquark potentials, which improves the simultaneous fit to both the meson and the baryon spectra.
In addition to building up the realistic potentials, there have been many works towards improving the accuracy of the spatial part of the total wave function. Among them, in Refs. [8,10], the harmonic oscillator bases were applied to construct the spatial part of the wave function for the baryon structure. In particular, Ref. [8] discussed the validity of the harmonic oscillator bases within the light baryons, composed of only u, d, and s quarks. They also extended their work to the tetraquarks in Refs. [8][9][10]. In Refs. [10,11], the authors built the spatial function using the hyper spherical coordinates for the tetraquark systems. In Ref. [12], using the spatial function that is made up of multiple Gaussians which allow for internal * sungsiknoh@yonsei.ac.kr † diracdelta@yonsei.ac.kr ‡ suhoung@yonsei.ac.kr angular momentum and satisfy the permutation symmetries restricted by the Pauli principle, the authors studied the properties of T cc and T bb , and confirmed the stability for the latter. Subsequently, in Ref. [13], the authors performed a similar calculations but with a more sophisticated bases of multiple Gaussians.
In this work, we have introduced a complete set of 3-dimensional harmonic oscillator bases with a rescaling factor that can be flexibly used for better convergence compared to the harmonic oscillator bases in the work [10], and applied them to construct the spatial part of the total wave function in a constituent quark model with hyperfine potential given in Eq. (1). In doing so, our spatial bases more efficiently describes the ground state wave functions not only for the meson and baryon structures, but also for the tetraquark structure.
In the following section, we first introduce the constituent quark model and compare the fittings in the present work to those in our previous work [14]. In Section III, each part of the total wave function is introduced, emphasizing the classificiation of the harmonic oscillator bases with the notion of quanta of the harmonic oscillator bases. In Section IV, the newly obtained numerical results are presented and compared with those in the previous work [14]. Furthermore, we show the relative distances between the quark pairs together with a pictorial description of the spatial sizes for the tetraquarks. We also discuss the tetraquarks as a baryon structure. Discussion and summary are given in Section V. In Appendix, to assess the validity of our approach, we compare the result obtained with our method for the hydrogen model with the exact solution. We also show the details in constructing the spatial wave function with the complete set of harmonic oscillator bases for the meson as well as the baryon structure. As a special case, in Appendix C, we present the method for constructing the bases for the proton.

II. FORMALISM
We use a nonrelativistic Hamiltonian for the constituent quarks, which is the same as in our previous work [14].
where m i is the quark mass and λ c i /2 is the SU(3) color operator for the i-th quark. The internal quark potentials V C ij and V CS ij are taken to be the same forms as in our previous work [14]: Here where r ij = |r i − r j | is the relative distance between the i and j quarks, and σ i is the spin operator. The parameters appearing in Eqs.
(2)-(5) are determined by fitting them to the experimental ground state masses of the hadrons listed in Tables I, II. The fitting has been done mostly within 0.7% error as can be seen in Tables I, II. In particular, we focused on trying to minimize the error on the doubly-charmed baryon, Ξ cc , and on the baryons Λ, p and also on the mesons D, D * , B, B * , B s which comprise the lowest thresholds for the doubly-heavy tetraquarks of interest in this work.
In comparison to our previous work [14], where only hadrons with at least one c or b were fit, here we fit all hadrons including those with only light quarks. The number of fitting hadrons is increased by almost a factor of 2 compared to that in our previous work [14]. The better fit is a consequence of using the complete set of harmonic oscillator bases approach. Using the harmonic oscillator bases, the new fitting parameters are as follows.

III. WAVE FUNCTION
Here, we will study the tetraquark states with the total orbital angular momentum equal to zero in the constituent quark model. However, it should be noted that the p-wave or even higher orbital states between the quarks can play a crucial role in lowering the total hadron energy and contribute to the total wave function. This can be successfully done by introducing the complete set of 3-dimensional harmonic oscillator bases. The total wave function of the Hamiltonian consists of the spatial, color, spin, and the flavor parts of bases. We adopt the harmonic oscillator bases as the spatial part of the wave function. The other parts of the wave function are the same as in the work [14]. We will thus discuss the spatial part in detail, while the other parts will be mentioned briefly.

A. Jacobi Coordinates Sets
To set up the spatial function, we first set the Jacobi coordinates, representing the relative positions of all the quarks in the tetraquark configuration. The Jacobi coordinates sets in each configuration can be written as follows.
• Coordinates Set 1 For symmetry reason, we take the coordinates set 1 as our reference, and use the transformations between the above sets of Jacobi coordinates to calculate relevant matrix elements involving two quarks.

B. Color, Spin Bases and Flavor
The most stable state for the doubly-heavy tetraquarks can be found in the spin 1 channel [14]. The colosr-spin (CS) space for the spin 1 tetraquark system is spanned by the six CS bases due to the fact that the tetraquark configuration in the total spin 1 state can be described by two color bases, and three spin bases. In the configuration of the Jacobi coordinates set 1, the six CS bases are as follows [14].
where the superscript indicates the color state, and the subscript indicates the spin state for the subparticle systems in the tetraquark structure.
For the falvor part, we will consider the isospin 0 for T cc (udcc), T bb (udbb), T cb (udcb), T sc (udsc), T sb (udsb), and the isospin 1/2 for usbb. The basis of the Hamiltonian is determined to satisfy the symmetry constraint due to the Pauli principle, and constructed by combining the CS basis with the spatial part. The permutation symmetries for the CS bases and for the flavor part are summarized in Table III.

C. Harmonic Oscillator Bases as The Spatial Function
Generalizing the method used in Appendix A, and B for the meson and baryon structures, we construct the complete set of harmonic oscillator bases for the tetraquarks. Since there are three internal orbital angular momenta l 1 , l 2 , and l 3 , there are three ways to combine them depending on the order of addition. Choosing to combine l 1 and l 2 first, the spatial function can be constructed as follows.
where the coefficient C(l 1 , m 1 , l 2 , m 2 ; L 1,2 , m 1,2 ) is the Clebsch-Gordan(CG) coefficient for the decomposition of the subtotal angular state |L 1,2 , m 1,2 in terms of |l 1 , m 1 |l 2 , m 2 , and the subtotal angular state of the subparticle corresponding to the quark pair (1, 2) is restricted to |L 1,2 = l 3 , m 1,2 = −m 3 to satisfy the total l = 0 state. R ni,li (x i ) has the same form as in Eq. (A.2), and Y mi li (θ i , φ i ) is the spherical harmonic function for the angular part of the i-th Jacobi coordinate x i . For the other two types of combinations, Eq. (11) is modified by adding l 1 , l 3 or l 2 , l 3 first.
Since there are a large number of harmonic oscillator bases in the tetraquark system, there has to be a criterion when adding the harmonic oscillator bases as a spatial function into the calculations. The magnitude of the eigenvalues of the Hamiltonian is mostly affected by the diagonal elements rather than the off-diagonal ones. Especially, except for the rest masses, the largest part in the diagonal element is the kinetic energy. In addition, non-zero kinetic energy contribution is mostly coming from the diagonal elements. Therefore, in classifying the harmonic oscillator bases, we define the quanta of the harmonic oscillator bases by the expectation value of the kientic energy of the diagonal component for the extreme case where all the constituent quark masses are identical. Introducing the center of mass frame, the kinetic energy denoted by T c becomes as follows. where for udcb, for udsc, for udsb, For the extreme case where all the masses of the constituent quarks are identical (m 1 = m 2 = m 3 = m 4 ≡ m), and thus also the variational parameters (a 1 = a 2 = a 3 ≡ a), the diagonal component of the kinetic energy reduces to Tc =h 2 c 2 a m 2n1 + l1 + 3 2 Therefore, one notes for this special cases that the kinetic energy is the same for all the possible combinations of the quantum numbers (n 1 , n 2 , n 3 , l 1 , l 2 , l 3 ) if the sum (Q = 2n 1 + 2n 2 + 2n 3 + l 1 + l 2 + l 3 ) is unchanged. One can now organize the harmonic oscillator bases according to the quanta Q appearing in Eq. (14). Then, for mesons, each of the spatial functions is categorized in the different quanta. For the tetraquarks, the spatial functions included in our calculations are classified as follows.
In terms of ψ spatial [n1,n2,n3,l1,l2,l3] without specifying the argu- As can be seen in the above classification, we have considered bases with internal angular momenta up to the l i = 3 states because the contribution from the higher angular momentum basis is not so significant. There are other types of harmonic oscillator bases constituting the complete set, such as ψ Spatial [0,0,0,1,1,1] . However, there is no transitional matrix element of the Hamiltonian between this type of bases and the other bases presented above. Thus, they do not contribute to the ground state masses of the tetraquarks. In doing so, we have observed the convergence behavior of the ground state mass, and have included the harmonic oscillator bases up to the 5th quanta in the calculations.
On the other hand, in order to satisfy the Pauli principle, it is necessary to consider the permutation symmetries of the spatial functions as well as the other parts of the wave function. In the Jacobi coordinates set 1 of our reference, it is obvious that T QQ configuration has symmetries under the permutations (12) and (34). T QQ ′ has symmetry under the permutation (12), while usbb has symmetry under the permutation (34) as summarized in Table III. The symmetry property of the spatial functions can be summarized for each tetraquark configuration in Table IV.
The masses and binding energies BT of the tetraquark states obtained with the fitting parameters in Eq. (6). The binding energy BT is defined by the difference between the tetraquark mass and the sum of the masses of the lowest threshold mesons, The values in the parentheses are the results in our previous work [14].
Note that there is antisymmetric part in some of the spatial bases due to the symmetry property of the harmonic oscillator bases as seen in Table. IV. Thus, contrary to the previous work [14] where only the fully symmetric spatial basis ψ Spatial [0,0,0,0,0,0] was used in the calculations, all the CS bases are used in the present work.

IV. NUMERICAL ANALYSIS
We use the total wave function discussed in the previous section, and perform a variational method to determine the ground state masses for the tetraquarks. The results are shown in Table V. We found that udbb and usbb are stable against the lowest strong decay threshold. In particular, one can see that the effect from the harmonic oscillator bases appears to lower the binding energies. As a result, T cb (udcb) is found to be below the lowest threshold. One of the major contributions to this is coming from the excited orbital states.
In this section, we first discuss a similar tendency of the ground state wave functions between the tetraquark and the meson structures. Then we will discuss the effects of the excited orbital states. Also, we investigate the spatial distibutions and the sizes for the tetraquarks. Finally, we compare our model with other works.

A. Ground State of Tetraquarks
We first analyze the ground state of the tetraquarks by expanding the wave function in terms of the complete set of harmonic oscillator bases. Using the variational method, we find that the expansion coefficients of the bases rapidly converge to zero.
As discussed in Section III C, each of the harmonic oscillator bases for mesons is classified by different quanta. Then, as can be seen in Figure 1 for the D meson, the expansion coefficient monotonically decreases as desired when the quanta of the bases increases. This tendency can be seen also for the tetraquark states in Figures 2-5, respectively. In addition to the convergence of the expansion coefficients, one can observe the convergence of the ground state masses as the number of the harmonic oscillator bases and their quanta increase as seen in Figure 10 for the D meson. We then determined the ground state masses for the mesons, baryons, and the tetraquarks when the convergent values change by just few MeV for  the entire bases contained in the last quanta, which turns out to be the 5th quanta. On the other hand, the convergence behavior of the expansion coefficients for the tetraquarks is not monotonic. This is so because while the nonzero orbital bases are ordered later, their contributions to the tetraquark configurations are important due to the attraction coming from the dipole and quadrupole moments. Still the average value of the coefficients in each quanta monotonically decreases. Discussing the tetraquark state in some more detail, one finds that from the second quanta on, many of the largest coefficients correspond to the harmonic oscillator bases with l 1 = l 2 = 0, which implies that the contributions from this type of bases are important to obtain the exact ground state masses for the tetraquarks. This can be quantitatively seen by the mass changes when this type of bases are additively included in the calculations. For simplicity, considering only the dominant CS bases, we first evaluate the mass with only the spatial basis ψ Spatial  Table VI. The contribution from the basis ψ Spatial [0,0,0,1,1,0] appears to lower the tetraquark masses by 11 MeV in T bb and 34 MeV in T cc , respectively. The effect is larger in T cc because the magnitudes of the corresponding expansion coefficients in T cc are obviously larger than the others in the second quanta, while it is not so in T bb .
Also, considering only the dominant CS bases, the coefficient in the first quanta corresponding to the spatial basis ψ Spatial [0,0,0,0,0,0] is larger in T bb than that in T cc . The first coefficient in T bb is 0.94, which implies that there are less contributions from the other bases in the tetraquarks of heavier antiquarks than in the tetraquarks of lighter antiquarks. However, compared to the total changes in masses from the value only with the first basis in Table VI to the exact ground state masses in Table V, the contribution from the basis ψ Spatial [0,0,0,1,1,0] is still important in T bb as well.
Let us now discuss the coefficients corresponding to the rest types of excited orbital bases. For T cc in Figure 2 and T cb in Figure 4, the coefficients in T cb are not so small relative to those in T cc . Obviously in the second quanta, it is comparable to the others in T cb , while they are not so in T cc . As a result, as can be seen in Table VI, the mass is lowered by 1 MeV in T cc but it is 4 MeV in T cb when adding the bases ψ Spatial [0,0,0,1,0,1] and ψ Spatial [0,0,0,0,1,1] to the calculations. This is due to the symmetry breaking in the flavor part of antiquarks, relative to T cc structure.
A similar behavior appears also in the comparison between T bb and usbb depicted in Figures 3,5. In this case, the flavor symmetry is bronken in the quark part. Comparing the changes in masses between T bb and usbb, the contributions from the bases ψ Spatial are little larger in usbb. Therefore, all the types of the excited orbital states are necessary to obtain the exact ground state masses for the tetraquark states.

B. Spatial Size of Tetraquarks
It is also useful to investigate the spatial size of the tetraquarks, and the spatial distribution of the constituent quarks in the tetraquark structure. The relative distances between the quarks are listed in Table VII. The relative distance between the heavier quarks are in general closer than that of the lighter quarks [14]. This tendency is maintained in each tetraquark state as seen in Table VII. Looking into the relative distances, it is found that the relative distances except for (1, 2) and (3,4) pairs are all the same in T QQ structure. For T QQ ′ structure, the relative distances for the pairs (1, 3) and (2, 3) are the same, and also for the pairs (1, 4) and (2,4). Likewise, the relative distances for the pairs (1, 3) and (1,4) are the same as are those for the pairs (2, 3) Figure 2 but for T cb . The number of Hamiltonian bases for T cb is 501, but the spatial bases are the same with those for Tcc. and (2,4) in usbb. This is due to the flavor symmetry in each tetraquark structure, and can be simply evaluated through the permutation symmetry for the ground state wave function in each tetraquark state. Since the total wave function satisfies the Pauli principle in each tetraquark, if we denote the ground state wave function by |Ψ T etraquark G ≡ |ψ Spatial × |ψ CS , the permutation symmetries for each ground state are as follows.

FIG. 4. Same as
Then the relative distances in each tetraquark structure can be obtained as follows. For For usbb, Therefore, the relative distances reflect the symmetry property.
On the other hand, without loss of generality, one can place one of the quarks at the origin of the Cartesian system, and also the second quark can be located on one of the axes. Then the degrees of freedom for the tetraquark system is in general 7. By the symmetry property, the degrees of freedom reduces to 3 for T QQ or 4 for T QQ ′ and usbb. Now we can prove for T QQ structure that the three independent vectors, R (1,2) ≡ (r 1 −r 2 ), R (3,4) ≡ (r 3 −r 4 ), and R ′ ≡ 1 2 (r 1 + r 2 − r 3 − r 4 ), which correspond to the Jacobi coordinates x 1 , x 2 , and x 3 are orthogonal to each other.
These relations show that we can describe the positions of the four quarks in T cc or T bb with the three independent vectors as shown in Figures 6, 7. For T cb , R (1,2) and R (3,4) can be taken as two orthogonal vectors, but for the third vector However, R ′ = (r 1 + r 2 ) /2 − (m 3 r 3 + m 4 r 4 ) / (m 3 + m 4 ) and linearly independent from R (3,4) and thus span the linearly independent third direction as shown in Figure 8.
One can also show that the non-vanishing component along the R (3,4) direction for T cb and R (1,2) direction for usbb in the wave functions are described through the excited orbital states that produce the asymmetry in the configuration of the Jacobi coordinates set 1 as can be seen in Figures 8 and 9, respectively. Specifically, since the spatial bases with (l 1 , l 2 , l 3 ) = (0, 0, 0) are of even powers with respect to the Jacobi coordinates x 1 , x 2 , and x 3 , one can find that the spatial part integration of becomes zero if considering only the spatial bases with (l 1 , l 2 , l 3 ) = (0, 0, 0). This implies that the calculations only with the spatial basis ψ Spatial [0,0,0,0,0,0] lead to the T QQ like structure even for both T cb and usbb, which is far away from the real structures in nature. In addition, using the relations of the relative distances in Table VII, it is possible to specify the spatial positions of the quarks in the tetraquark structure as the points on the surface of a sphere. The results are depicted in Figures 6-9 with the center specified as R c in the same scales.

C. Tetraquarks in a Baryon Structure
The structure (QQqq) is expected to be similar to that of an antibaryonQ 1qq with the heavy antiquarkQ 1 replaced by the diquark (QQ) [10]. However, in qqQ 1 structure, Q 1 is in fact composed of two heavy antiquarks so that it can be of either color [3] or [6]. For an antiquark pair of color [6], it cannot be regarded as a point particle because, at short distance, the Coulomb potential gives strong repulsion due to the fact that the color matrix element λ c i λ c j is 4/3 for the color [6] while it is −8/3 for the color [3]. Thus, in terms of color structure, this means that the (q 1 q 2 )3 ⊗(q 3q4 ) 3 channel, which is the only color configuration in qqQ 1 , dominates the (q 1 q 2 ) 6 ⊗ (q 3q4 )6 channel in T QQ ′ tetraquark structure. This assumption can be tested by computing the antidiquark (QQ ′ ) mass in the color [3] state, and then putting it into the evaluation of the baryon-like mass (q-q ′ -Q 1 ) where Q 1 replaces the antidiquark (QQ ′ ). An isolated diquark is in principle an ill defined concept as its nonzero color makes it a gauge dependent quantity so that one can add any other gauge dependent gluon field configuration to change its mass. What that means is that inside a tetraquark or baryon we are free to include any fractional amount of the interaction between the diquark and other color source to define the diquark mass. For example, the division of the interaction terms between light quark and heavy quark pairs in a tetraquark given in Table XII is in principle arbitrary. Furthermore, the different spin structure between a heavy quark and diquark, which can be of either spin 1 or 0, could induce different spin interaction for the two different configurations depending on the quantum numbers. Thus, we define the diquark mass in a constituent quark model to be the sum of the masses, the interactions between the two antiquarks, and their relative kinetic term all in an isolated color [3] and spin 0 or 1 configuration. We can then calculate the hypothetical baryon-like mass assuming a (q-q ′ -Q 1 ) structure and call it the tetraquark in a baryon structure, where the subscript 1 denotes the spin of the antidiquark. As shown in Table IX, the tetraquark mass is almost reproduced in the case of T bb and the difference becomes larger in the other tetraquark states. The dif- ference is related to how much the color (q 1 q 2 )3 ⊗ (q 3q4 ) 3 channel dominates over the (q 1 q 2 ) 6 ⊗ (q 3q4 )6 channel. As shown in Figures 2-5, the contribution of the dominant CS basis in the first quanta is much larger than the others in each tetraquark state, tendency of which is more apparent in T bb than the other tetraquarks. Moreover, looking at the sizes of the antidiquarks in Table X, the corresponding relative distances in the tetraquark structure in Table VII are larger. However, it becomes larger only by 0.006 fm in T bb while it is 0.149 fm in T cc . This implies that treating the antiquark pair such as (bb) in T bb as an isolated diquark seems a better approximation when the mass becomes heavier. Now we look at the size of (ud) pair in various configurations in Table XI. The size of (ud) pair in the diquark configuration becomes smaller in the baryons Λ c and Λ b due to the interaction with the heavy quark. Comparing the size of (ud) pair in Λ b and Λ c , one finds that the size of (ud) pair becomes smaller when the heavy quark is closer to the (ud) pair. This tendency can be seen also in the baron-like (q-q ′ -Q 1 ) structure. However, in the tetraquark configuration as can be seen in Table VII, the size of (ud) pair becomes smaller in T bb than that in T cc even though the relative distance (1,2)-(3,4) is smaller in T cc , which shows an opposite tendency to that in the baryon structure. TABLE XI. The relative distances between the quarks in various configuration. The size of an isolated ud diquark in the color [3] and spin 0 state is given in Column 2. The distances in the baryon-like (q-q ′ -Q1) structure are given in Columns 3-5. The distances in the baryons Λc and Λ b are given in Columns 6-7. Considering the tetraquarks T cc and T bb as a baryonlike (q-q-Q 1 ) structure, the size of (ud) pair in the baryonlike T cc (T bb ) is close to that in the baryon Λ c (Λ b ). However, as can be seen in Table VII, the size of (ud) pair in T cc is in fact 0.830 fm, which is not so close to that in the baryon-like T cc , and is even larger than the size of (ud) pair in the diquark configuration, while the size of (ud) pair in T bb is close to that in the baryon-like T bb . The analysis so far points out that T bb indeed has a similar structure as a baryon, while T cc does not.
As discussed in the work [14], the total mass obtained with the Hamiltonian in Eq. (1) can be divided into 'light quark', 'heavy quark', 'CS' parts as given in Table XII. The constant −D term in the Hamiltonian is divided into each quark by multiplying a factor of 1/2. The relative kinetic energy involving p 3 , which corresponds to the relative coordinate x 3 connecting the quark and the antiquark pairs can be divided according to the relative contributions depending on the mass of either the quark pair or the antiquark pair.
Considering only the spatial basis ψ Spatial [0,0,0,0,0,0] , the relative distance between the light quarks in T cc is the same with that in T bb [14]. Comparing the results between the full calculation and the simplified calculation denoted by 1-Basis in Table XII, the hyperfine interaction V CS (i, j) in T cc becomes much stronger than that in T bb . This can be understood from the change in V CS (1, 2), which implies that the attraction coming from V CS (1, 2) partially spreads into the attraction coming from V CS (i, j) in T cc . As a result, for T bb as shown in Table XII, the attraction from the hyperfine part betweenb and u (or d) quarks is just -5.7 MeV, while that in T cc is -69.4 MeV. This is one of the major reasons that the distance between the light quark pair and the heavy antiquark pair is closer in T cc than that in T bb as seen in Table VII. This closer distance in T cc causes the slightly larger size of (ud) pair in T cc than in T bb .  (udbb) and Tcc(udcc) masses from the present work. (i, j) denotes the i and j quarks, where i, j = 1, 2 label the light quarks, and 3, 4 are for the heavy antiquarks in each configuration.
V C (i, j) and V CS (i, j) covers pairs (i, j) except for (1,2) and (3,4) pairs. D is separately added and not included in V C (i, j). mQ is the heavy quark mass, and m ′ i are defined in Eq. (13) for each configuration. pi is the relative momentum corresponding to the i-th Jacobi coordinate xi. 1-Basis is the result with only one spatial basis ψ Spatial [0,0,0,0,0,0] and the corresponding dominant CS basis.

D. Comparison with Other Models
We now compare our present results with other works in Ref. [15] and Ref. [8]. In a simplified constituent quark model [15], the masses for T cc and T bb are given as where m b c,b,q are the constituent quark masses of c, b, and light quark q inside a baryon, and B(cc) (B(bb)) is the binding between the charm(bottom) quarks, which can be understood as coming from the extra attraction between two c(b) quarks due to that shorter interquark distance as compared to two light quarks: this attraction can be estimated by studying the quark attractions inside Λ c and Ξ cc [14]. a's are multiplicative constants for the CS interaction. Here, a QQ /(m b Q ) 2 is the CS interaction between the two heavy quarks Q corresponding to V CS (3,4) in our model, while −3a/(m b q ) 2 is that between the light quarks q denoted by V CS (1, 2) in our model. Then, treating B(cc) (B(bb)) as part of the two charm(bottom) quark system, the energy in the simplified model in Eq. (21) can be divided into the charm(bottom) quark, light quark, and CS interaction parts. The subtotal values in Table XII can be regarded as the constituent quark masses in the simplified quark model in Eq. (21). The additional attraction for the heavy quark pairs with respect to the light quark pair in the simple model [15] can be seen from V C (3, 4) being much more attractive than V C (1, 2). While the importance of the large additional attraction for heavier quarks denoted by B(QQ) remains a valid point, detailed mass and size differences change other parameters such as the effective quark masses in Eq. (21) so that a simple parametrization formula given in Eq. (21) for the tetraquark masses might become problematic. Furthermore, while Eq. (21) does not allow for the color spin interaction between light and heavy quarks, our full calculations show the presence of such terms as given in Table XII, which becomes −5.7 and −69.4 MeV for T bb and T cc , respectively. Also, the discrepancy becomes larger for T bb suggesting that nonlinear quark mass dependence becomes also important. Therefore, care should be taken when simple parametrizations that work for normal hadrons are generalized to more complicated configurations.
There is another work [8] using the complete set of harmonic oscillator bases, which is similar to ours but different in part as we introduced the rescaled form of the harmonic oscillator bases. The hamiltonian in their model is as follows. with Comparing with our results in Table XIII, the two independent models give almost the same binding energies for each tetraquark state. However, there is a large difference in the mass of T cc , which is due to the larger obtained masses for the D and D * mesons in their model, which were found to be 1891 MeV and 2021 MeV, respectively. Compared to the total mass of the D and D * mesons in our model, the difference is 52 MeV, which approximately accounts for the difference in the mass of T cc between the two models. Also, in their model, the lowest threshold for usbb is B B * s , while it is B * B s in our model as in the experiment. In the experiment, the total mass of the B and B * s mesons is slightly larger than that of the B * and B s mesons by approximately 2.7 MeV, while the corresponding difference in our model calculation is 4.7 MeV.

V. SUMMARY AND DISCUSSION
We have improved our constituent quark model by introducing the complete set of 3-dimensional harmonic oscillator bases. We also assessed the validity by comparing the ground state wave function for the meson structure with the exact solution of the hydrogen model. The effect turns out to lower the binding energies for the tetraquark systems. In particular, the harmonic oscillator bases for the excited orbital states play a crucial role in obtaining the exact ground state for the tetraquark systems. We have also successfully fitted the parameters in the Hamiltonian in Eq. (1) to most of the observed mesons and baryons allowed in our model. The results are summarized in Table V. Also, from the relative distances between the quarks given in Table VII, we have described the spatial distributions of the quarks in the tetraquark structure in Figures 6-9. By comparing with our earlier work shown in Table V, and the discussions in Sections IV B, IV C, we have found that a simple Gaussian spatial function fails to provide precise information on the stability and the structure for the tetraquarks so that the detailed treatment as presented in this work should be performed. Also, by comparing with other work, fitting procedure is important to evaluate the exact values for the masses. One can conclude that while simplified constituent quark model based on universal constants that does not depend on specific configuration and/or simple models based on universal diquarks are intuitively important to identify possible attractive configurations, a detailed full constituent quark model calculation is needed to assess the stability and existence of a compact exotic multiquark configuration.

ACKNOWLEDGMENTS
This work was supported by Samsung Science and Technology Foundation under Project Number SSTF-BA1901-04, and by the Korea National Research Foundation under the grant number 2017R1D1A1B03028419(NRF).

A. Harmonic Oscillator Bases in Mesons
To construct the spatial function for the meson structure, we solve the Schrödinger equation for the 3dimensional symmetric harmonic oscillator. In the spher-ical coordinates system, the wave function can be separated into the radial part and the angular part such as ψ(r, θ, φ) = R(r)Y m l (θ, φ). The solution of the angular part is known to be the spherical harmonics. For the radial part of the equation, the solution is obtained in terms of the associated Laguerre polynomials. The orthonormalized radial part wave function is known as follows.
R n,l (r) = 2 Γ(n + 1) where L l+ 1 2 n (r 2 ) is the associated Laguerre polynomials. For our purpose of introducing the harmonic oscillator bases to our model, we have modifed Eq. (A.1) by rescaling the radial distance r to √ 2a x, where x is the magnitude of the Jacobi coordinate x, connecting the quark and the antiquark in the meson structure, and a is the variational parameter corresponding to the coordinate x. Then the spatial part of the total wave function is constructed by combining the spherical harmonics as follows.
where R n,l (x) is the rescaled radial part wave function.
Here, the quantum numbers n, l indicate the principal quantum number, and the orbital angular momentum, respectively. For mesons in the l = 0 state, each of the harmonic oscillator bases should be of l = 0, and they compose the spatial part of the wave function in the meson structure.
On the other hand, the permutation symmetry of the spatial bases depends on the power of the Jacobi coordinate |x|, which is contained only in the radial part of the spatial function. From Eq. (A.1), it is recognized that the permutation symmetry is determined by the angular momentum quantum number l. In our case, for all the bases in the calculations, the angular momentum is l = 0. Therefore, for the mesons in the ground state, all the spatial bases are symmetric under the permutation (12).
To assess the validity of the harmonic oscillator bases approach, we compared the meson structure to the hydrogen atom in hadron picture. To do this, we considered only the kinetic energy and the Coulomb potential in the Hamiltonian of Eq. (1). Then the Hamiltonian reduces to the following form with the kinetic energy in the center of mass frame. where m ′ = (2m 1 m 2 )/(m 1 + m 2 ), and p x is the relative momentum between the quark and the antiquark. We take m 1 = m c , m 2 = m q , and the values in Eq. (6) for the parameters m c , m q , and κ. For our purpose of the comparison to our model, we also modify the exact solution of the hydrogen atom as that of the relative coordinate x. Then, in hadron picture, the radial part of the ground state wave function for the hydrogen, R Hydrogen where µ is the reduced mass of the system in MeV unit. In comparing, we choose the D meson as a target, then the reduced mass µ in Eq. (A.4) becomes mumc mu+mc . The results are shown in Figure 10. It is obvious from Figure 10 that the more harmonic oscillator bases are included, the more precisely it describes the actual ground state of the meson structure. We have included the har- monic oscillator bases up to n = 4, which is enough to obtain the convergent values for the ground state masses of the tetraquarks of interest in this work.

B. Harmonic Oscillator Bases in Baryons
In baryons, there arise additional degrees of freedom due to the second Jacobi coordinate x 2 . Furthermore, there are contributions from combinations of nonzero internal relative orbital angular momenta, satisfying zero total orbital angular momentum of the ground state baryon structure. Using a similar method as in the meson structure, we construct the complete set of orthonormalized harmonic oscillator bases.
The following are the Jacobi coordinates for the baryon structure: we choose the coordinates set 1 as our reference.
• Coordinates Set 1 As in the mesons, the permutation symmetry of the spatial bases is determined by the angular momentum quantum numbers l 1 and l 2 . As can be seen in the coordinates set 1, x 1 is antisymmetric under the permutation (12), while x 2 is symmetric. Therefore, the bases with even number of l 1 is symmetric, while the bases with odd number of l 1 is antisymmetric under the permutation (12).
In the case of the proton where all the constituent quarks are identical, the total wave function should be fully antisymmetric. Thus, we need to specify the symmetries for the permutations (13) and (23) as well. However, it is not clear in the coordinates set 1. Therefore, it is instructive to investigate the method of constructing the spatial bases for the proton. In such a way above, it is possible to relate the harmonic oscillator bases to each of the Young tableaux for the permutation group S 3 . To fully construct the total wave function for the proton, it is necessary to construct the remaining parts of the wave function. For the color basis, it is always fully antisymmetric for the baryons. Thus, we focus on the spin and the isospin parts, and ψ Spatial ×ψ Isospin ×ψ Spin part should be fully symmetric. Since both the spin and the isospin bases are constructed by SU (2), we obtain the isospin-spin basis from the inner product of the flavor SU(2) and the spin SU (2). For the proton,