Exotic resonances of fully-heavy tetraquarks in a lattice-QCD insipired quark model

Fully-heavy tetraquark states, i.e. $cc\bar{c}\bar{c}$, $bb\bar{b}\bar{b}$, $bb\bar{c}\bar{c}$ ($cc\bar{b}\bar{b}$), $cb\bar{c}\bar{c}$, $cb\bar{b}\bar{b}$, and $cb\bar{c}\bar{b}$, are systematically investigated by means of a non-relativistic quark model based on lattice-QCD studies of the two-body $Q\bar{Q}$ interaction, which exhibits a spin-independent Cornell potential along with a spin-spin term. The four-body problem is solved using the Gaussian expansion method; additionally, the so-called complex scaling technique is employed so that bound, resonance, and scattering states can be treated on the same footing. Moreover, a complete set of four-body configurations, including meson-meson, diquark-antidiquark, and K-type configurations, as well as their couplings, are considered for spin-parity quantum numbers $J^{P(C)}=0^{+(+)}$, $1^{+(\pm)}$, and $2^{+(+)}$ in the $S$-wave channel. Several narrow resonances, with two-meson strong decay widths less than 30 MeV, are found in all of the tetraquark systems studied. Particularly, the fully-charm resonances recently reported by the LHCb Collaboration, at the energy range between 6.2 and 7.2 GeV in the di-$J/\psi$ invariant spectrum, can be well identified in our calculation. Focusing on the fully-bottom tetraquark spectrum, resonances with masses between 18.9 and 19.6 GeV are found. For the remaining charm-bottom cases, the masses are obtained within a energy region from 9.8 GeV to 16.4 GeV. All these predicted resonances can be further examined in future experiments.


I. INTRODUCTION
Many efforts have been made in the past twenty years in order to understand exotic tetra-, penta-and even hexa-quark systems, constituted mostly by heavy quarks. For example, in 2015, the hidden-charm pentaquark states P c (4380), P c (4312), P c (4440) and P c (4457) were discovered by the LHCb collaboration [1,2], adding them to a large number of unconventional heavy mesons called jointly XYZ states, e.g. the X(3872) discovered by the Belle Collaboration in 2003 [3], the Y (4260) discovered by the BaBar Collaboration in 2005 [4], Z c (3900) discovered by the BESIII Collaboration in 2013 [5], etc.
The fully-heavy tetraquark states QQQQ (Q = c, b) have recently attracted much attention. In 2017, the CMS Collaboration reported a benchmark measurement of Υ(1S)-pair production in pp collisions at √ s=8 TeV [6]. A preliminary analysis of the CMS data shows an excess at 18.4 GeV in the Υ(1S) + − decay channels [7][8][9]. This excess, if confirmed by future experiments, may indicate a fully-bottom tetraquark state (bbbb). Besides, a significant peak at 18.2 GeV was observed in Cu+Au collisions at RHIC [10] but the LHCb and CMS collaborations [11,12] were not able to confirm it from the Υ(1S)µ + µ − invariant mass spectrum. Recently, in the di-J/ψ invariant mass spectrum, a narrow peak at 6.9 GeV, a broad one between 6.2 and 6.8 GeV, and a hint for a possible structure around 7.2 GeV were reported by the LHCb collaboration, which could indicate the existence of fully-charm tetraquarks [13]. It is then expected that more investigations on the existence of fully-heavy tetraquark states will be performed in the future in, for instance, the LHCb experiment at CERN.
The search of fully-heavy tetraquark states is valuable from the experimental side but also from the theoretical view point [14][15][16]. The debates on the existence of QQQQ (Q = c, b) states were quite intense before the LHCb's findings in the di-J/ψ invariant mass spectrum. A bbbb bound state was supported by various models [17][18][19][20], QCD sum rules studies [21,22], and a naive diffusion Monte Carlo calculation [23]. The study of its decay properties concluded that a fully-bottom tetraquark state with J P = 2 + is possible [24]. As for the potential cccc state, it was suggested that its mass is located between 5 and 6 GeV through various phenomenological models [17,[25][26][27] and the Bethe-Salpeter equation approach [28]. Additionally, narrow bbbc and bcbc tetraquark states were predicted [29,30]. On the contrary, there were counter examples on the existence of QQQQ tetraquark states. For example, it was shown that bound states of cccc and bbbb systems are impossible through model investigations of Refs. [29][30][31][32][33][34] but also in preliminary lattice-QCD computations [35].
The interpretation of the newly reported structures in the di-J/ψ invariant mass spectrum have been carried out by various theoretical approaches. In a non-relativistic model with a compact diquark-antidiquark configuration, masses of the S-wave fully-charm tetraquark states are predicted to be between 5.96 and 6.32 GeV [36]. In a potential model including the linear-confinement and the one-gluon-exchange contributions, fully-charm tetraquark masses are predicted at 6.5 GeV and 6.9 GeV, which can be identified with the observed structures as S-and P -wave fully-charm tetraquarks [37]. This conclusion differs from a study with a dynamical diquark model, where the two resonances around 6.7 GeV and 6.9 GeV are identified with the 1P and 2S multiplets [38]. Further-arXiv:2104.08814v1 [hep-ph] 18 Apr 2021 more, the narrow and broad structures observed in the di-J/ψ invariant mass spectrum are explained as radial excitations of the fully-charm state within the QCD sum rules approach [39], the string-junction picture [40], and the extended relativized quark model [41]. On the other hand, based on the perturbative QCD approach [42], it was predicted that both the narrow and broad structures prefer molecular-like configurations. On the contrary, a holography inspired model suggests a compact picture for the structure around 6.9 GeV [43]. Concerning the quantum numbers of the di-J/ψ structures, the spin-parity of the X(6900) resonance is suggested to be J P = 2 + within the constituent quark model of Ref. [44] and the relativistic quark model of Ref. [45], while it is predicted to be 0 + or 1 + by an effective potential model [46] and 0 + by a diffusion Monte Carlo calculation [47].
In addition, a unitary coupled-channels approach seems to describe well di-J/ψ invariant mass spectrum finding hints of a near-threshold molecular state X(6200) with quantum numbers J P C = 0 ++ or 2 ++ [48]. Meanwhile, an ab initio perturbative QCD investigation suggested that there exists another state near the peak at 6.9 GeV and thus the nature of X(6900) can be uncovered by calculating the cross section of this undiscovered state [49]. On the other hand, strong decay properties of the fully-charm tetraquark states were investigated in Ref. [50], predicting that the broad structure is a S-wave state with either J P C = 0 ++ or 2 ++ and the narrow peak can be identified as a P -wave state with quantum numbers J P C = 0 −+ or 1 −+ . The production of the fullycharm tetraquark state at 6.9 GeV was also studied by a pp annihilation process [51], a dynamic simulation [52], and a model-independent fragmentation mechanism [53].
In this work, we explore the possibility of having bound, resonance and scattering states of fully-heavy quark systems, viz. QQQQ (Q = c, b), with spin-parity J P (C) = 0 +(+) , 1 +(±) and 2 +(+) in the S-wave channel. 1 We employ a non-relativistic quark model for the twobody interaction between heavy quarks according to the Lattice-QCD study of Ref. [55]. The four-body problem is solved by using the Gaussian expansion method [56], which has been demonstrated to be as accurate as a Faddeev calculation. The complete set of four-body configurations: meson-meson, diquark-antidiquark, and K-type arrangements, as well as their couplings, are considered in the calculation. Meanwhile, a coupled-channels calculation which treats the bound, resonance, and scattering states on the same footing is performed by employing the complex scaling method [57][58][59][60][61] according to the so-called ABC theorem [57,58]. This tool has been already used in previous studies of nuclear [62][63][64][65] and hadron [66][67][68][69] physics.
This manuscript is organized as follows. The theoretical framework is presented in Sec. II, including the lattice-QCD inspired potential model and the general structure of the four-body wave function. Section III is devoted to the analysis and discussion of the obtained results. And we summarize our theoretical work in Sec. IV.

II. THEORETICAL FRAMEWORK
The kinetic motion for heavy quarks (charm and bottom) can be treated in a non-relativistic way, and their interactions can be described via a potential model. Here, we employ a potential model inspired by the Lattice-QCD investigation of Ref. [55], viz. the interaction between a heavy quark and a heavy antiquark can be well approximated by the spin-independent Cornell potential along with a spin-spin term. For the four-body system QQQQ (Q = c, b), the Hamiltonian results to be where M Q is the mass of the heavy quarks and antiquarks. The two-body interaction potential can be written as which includes the coulomb, linear-confining and spinspin interactions. The color-dependence of the interaction is encoded in the SU(3) Gell-Mann matrices, λ a i (a = 1, 2, ..., 8). Since the two-body interaction depends only on the distance between the two quarks/antiquarks labeled by i and j, |r ij | = |r i − r j |, the four-body problem can be factorized into a center-of-mass motion and relative motions. In the following, we subtract the centerof-mass motion and focus on the relative motions of the four-body system. The model parameters α, β, γ, and σ can be determined via a calculation of the mass spectrum of the S-wave QQ mesons. Table I lists the values of the model parameters, which are the ones collected in Ref. [70] except the spin-spin coefficient of the charmbottom case, β cb , which is determined herein. The calculated masses of the S-wave QQ mesons along with their experimental values are listed in Table II. The theoretical results are consistent with the experimental data, providing a solid ground to study possible bound and resonance states in the four-body system QQQQ. Figure 1 shows a complete set of configurations for the fully-heavy tetraquark system QQQQ (Q = c, b). Figs. 1(a) and 1 (b) correspond to the mesonmeson structure and diquark-antidiquark one, respectively. Figs. 1(c)-(f) are the so-called K-type configurations. In general, all these configurations are coupled to each other; on the other hand, it was shown [71,72] that if all possible excited states are included, it is enough to consider only the color-singlet channel. Since it is difficult to include all excited states and we are interested in the low-lying ones, it is more convenient to consider all possible configurations and the couplings among them. The price to pay following this procedure is that the overlap matrix becomes singular durinf the process of diagonalization because the exchange of identical particles. This can be fixed by diagonalizing the overlap matrix and drop the eigenvectors which correspond to vanishing eigenvalues. Then, a new Hamiltonian matrix can be generated using the rest of eigenvectors and the energy levels of the four-body system can be obtained from the reconstructed Hamiltonian matrix.
A. Color, flavor, and spin structure There are four fundamental degrees of freedom for heavy quarks, i.e. color, flavor, spin, and spatial degrees of freedom.
The color-space of a tetraquark system should be spanned by 12 basis functions. In our description, they correspond to the complete configurations shown in Fig. 1 and each one gives two basis wave functions or color channels. That is to say, the meson-meson configuration, All possible configurations in the fully-heavy tetraquarks. Panel (a) is the meson-meson structure, panel (b) is the diquark-antidiquark one and the other panels are of K-type. (Q = c, b) Fig. 1(a), gives two color channels, ϕ 1 and ϕ 2 , corresponding to the product of two color-singlets (1 c ⊗ 1 c ) and the product of two color-octets (8 c ⊗ 8 c ). They are explicitly given by (rr +ḡg +bb) × (rr +ḡg +bb), (3brrb + 3ḡrrg + 3bgḡb + 3ḡbbg + 3rgḡr The diquark-antidiquark configuration, Fig. 1 (b), gives again two color basis, ϕ 3 and ϕ 4 , corresponding to the product of a color triplet and a color antitriplet (3 c ⊗3 c ) and the product of a color sextet and a color antisextet (6 c ⊗6 c ). They are explicitly given by (rrḡg −ḡrrg +ḡgrr −rgḡr +rrbb −brrb +bbrr −rbbr +ḡgbb −bgḡb +bbḡg −ḡbbg), (2rrrr + 2ḡgḡg + 2bbbb +rrḡg +ḡrrg +ḡgrr +rgḡr +rrbb +brrb +bbrr +rbbr +ḡgbb +bgḡb +bbḡg +ḡbbg). Each K-type configuration, Figs. 1(c)-(f), corresponds with two color channels. For convenience, we denote the K-type configurations from (c) to (f) in Fig. 1 as K 1 , K 2 , K 3 , and K 4 , respectively. The corresponding color wave functions have been constructed in Ref. [67]. The color wave functions ϕ 5 and ϕ 6 , corresponding to the K 1 configuration [ Fig. 1(c)], are given by The color wave functions ϕ 7 and ϕ 8 from the K 2 configuration [ Fig. 1(d)] are similar. We have The color wave functions ϕ 9 and ϕ 10 from the K 3 configuration [ Fig. 1(e)] can be shown to be ϕ 9 = 1 2 √ 6 (rbbr +rrbb +ḡbbg +ḡgbb +rgḡr +rrḡg +bbḡg +bgḡb +ḡgrr +ḡrrg +bbrr +brrb) + 1 √ 6 (rrrr +ḡgḡg +bbbb), And, finally, the color wave functions ϕ 11 and ϕ 12 from the K 4 configuration [ Fig. 1 The flavor structure of the QQQQ (Q = c, b) tetraquark system is trivial and will be neglected. The total spin quantum number can be S = 0, 1, and 2. Since there is no spin-orbital interaction in the Hamiltonian, the third component of the total spin (M S ) can be set to be equal to S without loss of generality. The four-body spin wave functions, ξ α S , for the configurations shown in Fig. 1 can be listed as The indices u 1 , ..., u 4 , v 1 , ..., v 6 , and w 1 number the spin channels corresponding to the configurations shown in  Table III. These spin wave functions are obtained by considering the cou-pling of sub-clusters with the SU(2) algebra. The necessary bases are given by

B. Computational method
We solve the four-body problem by means of an exact and high-efficiency numerical approach, the Gaussian expansion method [56]. Three relative orbital motions of the four-body system are all expanded by using the Gaussian bases, where the widths ν n are taken as the sizes of the geometric progression 2 and the normalization factor N nl is not shown explicitly. The general spatial wave function of the four-body system can be formally expressed as where ρ, λ, and R are the internal Jacobi coordinates for the four-body configurations in Fig. 1. For the mesonmeson configuration, Fig. 1(a), they are given by For the diquark-antdiquark configuration, Fig. 1 (b), they read For the K-type configurations in Figs. 1(c)-(f), we have where the indices i, j, k, l should be properly assigned according to different K-type configurations. Furthermore, with these coordinates the center-of-mass motion can be completely eliminated. 2 The details of the Gaussian parameters can be found in Refs. [73].
The complete basis set of the four-body wave functions that fulfill the Pauli principle can be expressed as where J is the quantum number of the total angular momentum. Here the antisymmetrizing operator A is added considering the exchange of identical particles. Particularly, for the fully-charm, -bottom andcbcb tetraquark states, it is explicitly given by Meanwhile, for thecccb andbbbc tetraquark states, it reads as and A=1 for theccbb case. This antisymmetrizing procedure is necessary since the basis set is constructed from sub-clusters, i.e. meson-meson, diquark-antidiquark, and K-type structures. The Hamiltonian matrix can be constructed using this complete basis set of the four-body wave functions.
We employ the complex scaling method (CSM) [57][58][59][60][61][62][63] according to the ABC theorem [57,58] in order to search for possible resonant states of the fully-heavy tetraquark system. Using the CSM, the resonance energy (its position and width) is obtained as a stable eigenvalue of the complex scaled Schrödinger euqation: where H(θ) is obtained from the original Hamiltonian by transforming the four-body Jacobi coordinates with respect to a common complex scaling angle θ: According to the ABC theorem [57,58], the eigenvalues of Eq. (19) can be separated into three groups: (i) A bound state corresponds to a pole on the real axis whose energy is smaller than the meson-meson threshold for a specific J P (C) state. It remains unchanged under the complex scaling transformation and such energy is identified as the mass of the four-body bound state.
(ii) The discretized continuum states, associated with the cuts, are rotated downward by an angle of 2θ with respect the real axis.
(iii) A resonant state corresponds to a complex pole, E res , that is independent of the angle θ. Moreover, it is isolated from the discretized continuum spectrum, lying along the 2θ-rotated line when the relation tan 2θ > −Im(E res )/Re(E res ) is satisfied. The mass M and width Γ of the four-body resonance are given by Predicted resonances of the fully-charm tetraquark system from the complete coupled-channel calculation with CSM. Their masses and widths are summarized in the third and fourth column, respectively, in units of MeV.

III. RESULTS AND DISCUSSION
We present herein our results for the low-lying S-wave states (L = 0 and hence J = S) of the fully-heavy tetraquark systems QQQQ (Q = c, b). They are organized as cccc, bbbb, cbcc, cbbb, bbcc and cbcb tetraquark states. Since the wave function of the tetraquark system is constructed at the fundamental quark level, some forbidden di-meson channels at hadron level, such as the di-J/ψ structure J/ψJ/ψ in the 1 +(−) state, are included; however, they will exhibit self-consistently scattering nature in our numerical calculations. These forbidden dimeson channels, interpreted as thresholds, are actually useful for us to identify possible tetraquark states.
A. Fully-charm system cccc We find no bound state for all S-wave channels J P (C) = 0 +(+) , 1 +(−) , and 2 +(+) of fully-charm system cccc. However, using the CSM, we identify several narrow resonances. Table IV summarizes our findings, the results are discussed in the following.
The J P (C) = 0 +(+) channel: To seek for possible bound states of the fully-charm tetraquark system, we first consider the real case (θ = 0 • ), shown in Table V. Therein, the lowest energies of all channels from a fully-uncoupled calculation (without considering any channel-channel coupling), the lowest energies of all structures from a partially-coupled calculation (considering only the couplings among the channels belonging to each structure), and the lowest energy from a fullycoupled calculation (considering all channel-channel couplings), are listed. The energy (mass) spectrum from the fully-coupled calculation is also shown in Fig. 2 (green circles on the real axis).
The meson-meson configuration includes two chan- V. Calculation of the J P (C) = 0 +(+) state of the fully-charm system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first two meson-meson structures give the theoretical results of the noninteracting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV. Then in the fully-uncoupled calculation, bound state is unavailable in these channels. Particularly, the masses of the hidden-color channels η c η c and J/ψJ/ψ are both around 6.4 GeV, which is also very close to those of the diquark-antidiquark channels. Meanwhile, the masses of the K-type channels, which are all above the di-meson thresholds, generally range from 6.2 to 6.4 GeV.
Complex energy spectrum of the fully-charm tetraquark system with J P (C) = 0 +(+) from the complete coupled-channel calculation with CSM. The parameter θ varies from 0 • to 6 • .
In the partially-coupled calculations, which are performed for each kind of configuration, the calculated masses are still all above the lowest threshold value, 5936 MeV for η c η c . From Table V we find that, the masses of the excited states for the diquark-antidiquark and Ktype configurations range from 6.15 to 6.35 GeV. Since we are treating a tetraquark system in which all quarks and antiquarks have the same flavor, no interaction in the Hamiltonian (1) can distinguish between the configurations K 1 and K 2 , as well as between K 3 and K 4 ; Thus, the results are degenerate and this is also true for other J P (C) states.
Finally, in the fully-coupled calculation, the lowest mass of the spectrum is 5936 MeV, which is still slightly above the η c η c threshold. Therefore, we can conclude that no bound state exists with J P (C) = 0 +(+) . However, as it will be shown next, several resonances are found in the complex-range and complete coupled-channel calculation. Figure 2 presents the complex energy spectrum with the angle θ ranging from 0 • to 6 • . In the mass range of 5.9 − 7.4 GeV, there are five meson-meson thresh- The calculated complex eigenenergies are generally aligned along the threshold lines, which are rotated downward making an angle of 2θ with the real axis, indicating the nature of the scattering states. Being cautious due to the calculation noise in the highly excited region, we can identify five fixed poles that are independent of θ, as circled in Fig. 2. Apparently, they can be identified as the resonance states of the fully-charm system cccc.
The first two poles, with masses and widths (6640 MeV, 0.66 MeV) and (6762 MeV, 1.14 MeV), are just above the η c (1S)η c (2S) threshold lines. Therefore, they TABLE VI. Calculation of the J P (C) = 1 +(−) state of the fully-charm system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first two meson-meson structures give the theoretical results of the noninteracting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.

Structure
Channel Complex energy spectrum of the fully-charm tetraquark system with J P (C) = 1 +(−) from the complete coupled-channel calculation with CSM. The parameter θ varies from 0 • to 6 • . and 6987 MeV are both quite close to the narrow structure around 6.9 GeV claimed by the LHCb collaboration [13] and hence they could be candidates for this. Meanwhile, the remaining η c (1S)η c (2S) and J/ψ(1S)ψ(2S) resonances are consistent with the broad structure between 6.2 and 6.8 GeV, and the possible peak shown around 7.2 GeV [13]. Their predicted small widths indicate that these resonances are stable against two-meson strong decays, indicating that they could be confirmed in future experiments.
The J P (C) = 1 +(−) channel: There are 23 channels in this case as listed in Table VI. We have four meson-meson channels, including the color-singlet and hidden-color states of η c J/ψ and J/ψJ/ψ, one diquarkantidiquark channel [cc] 6c [cc]6 c , and 18 K-type ones. A bound state is still unavailable in each single-channel calculation. In particular, the theoretical masses for the two lowest di-meson states, η c J/ψ and J/ψJ/ψ, in the color-singlet channel, are 6070 MeV and 6204 MeV, respectively. Moreover, the coupled-channel effect is too weak to acquire binding energy within these two channels. In the partially-coupled calculation, the masses of K 1 and K 2 configurations are remarkably low due to the strong coupled-channel effect. However, these masses (6271 MeV) are still above the theoretical thresholds of η c J/ψ and J/ψJ/ψ. In the fully-coupled calculation, the lowest mass of the spectrum is 6070 MeV, which is still slightly above the η c J/ψ threshold. Therefore, we conclude that no bound state exists with J P (C) = 1 +(−) .
Within the CSM, the complete coupled-channel calculation for J P (C) = 1 +(−) cccc tetraquark system has been performed and the complex energy spectrum is shown in Fig. 3. While most of the complex eigenenergies decline with increasing angle θ indicating that they are scatter-TABLE VII. Calculation of the J P (C) = 2 +(+) state of the fully-charm system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The value in the parenthesis after the first meson-meson structure gives the theoretical result of the noninteracting meson-meson threshold. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.  The J P (C) = 2 +(+) channel: We have two dimeson channels, J/ψJ/ψ and J/ψ 8 J/ψ 8 , one diquarkantidiquark channel, [cc] 6c [cc]6 c , and six K-type ones, as listed in Table VII. Again, no bound state is found in this highest total-spin case by all kinds of computations, including the partially-coupled and fully-coupled ones. The lowest mass of the spectrum from the fully-coupled calculation is 6204 MeV, which is slightly above the theoretical threshold of di-J/ψ.
Using the CSM, we obtain the complex energy spectrum from a complete coupled-channel calculation in Complex energy spectrum of the fully-charm tetraquark system with J P (C) = 2 +(+) from the complete coupled-channel calculation with CSM. The parameter θ varies from 0 • to 6 • .  Since the mass of this resonance is close to 6.9 GeV, the first radial excitation state of di-J/ψ, i.e. J/ψ(1S)ψ(2S), in J P (C) = 2 +(+) state is also a possible candidate for the X(6900).

B. Fully-bottom system bbbb
We do not find bound states for all S-wave cases: J P (C) = 0 +(+) , 1 +(−) , and 2 +(+) , of the fully-bottom system bbbb. However, performing a complete coupledchannel calculation within the CSM formalism, several narrow resonances are found and they are summarized in Table VIII. The results are discussed in the following.
The J P (C) = 0 +(+) channel: Table IX lists the four di-meson channels, including the color-singlet and hidden-color configurations of η b η b and ΥΥ, the two diquark-antidiquark channels, [bb] 3c [bb]3 c and [bb] 6c [bb]6 c , and twelve K-type allowed structures. For the colorsinglet channels of di-η b and di-Υ configurations, the lowest masses are just their corresponding theoretical thresholds, 18802 MeV and 18926 MeV, respectively. For other exotic structures, most of their masses are above 19.1 GeV, while two channels in the K 1 and K 2 configurations possess masses around 19 GeV, close to the threshold of the ΥΥ channel. When a partially-coupled calculation is performed for each structure, masses are ranging from 18.8 to 19.1 GeV. The lowest mass of the spectrum from the fully-coupled calculation is 18802 MeV, which is slightly above the threshold of di-η b . Therefore, no bound state exists for J P (C) = 0 +(+) .
We now look for possible resonances by performing a complete coupled-channel calculation within CSM. The complex energy spectrum is presented in Fig. 5, with energy ranging from 18.8 to 19.5 GeV. There are four scattering states, and Υ(1S)Υ(2S); with varying angle θ, the complex eigenenergies are generally aligned along the cut lines. In addition to the scattering states, four resonance poles are found whose structure and mass can be iden- Although the fully-bottom tetraquark state has not yet been reported experimentally, our theoretical findings could be valuable for future experimental studies. In fact, the Υ(1S)Υ(1S)(19019) resonance is also supported TABLE IX. Calculation of the J P (C) = 0 +(+) state of the fully-bottom system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first two meson-meson structures give the theoretical results of the noninteracting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.

Structure
Channel in other theoretical works as, for instance, Ref. [52]. The J P (C) = 1 +(−) channel: there are 23 channels listed in Table X, including four di-meson channels with η b Υ and ΥΥ di-meson configurations in the singletand hidden-color states, one diquark-antidiquark channel, [bb] 6c [bb]6 c , and 18 K-type ones. First of all, the lowest channel η b Υ is unbound, with a mass of 18864 MeV. The other color-singlet di-meson channel, ΥΥ, is also unbound with the lowest mass 18926 MeV. The masses of the remaining configurations populate densely a small energy range, 19.1 − 19.2 GeV; except for two channels from K 3 and K 4 structures whose masses are located at ∼19.6 GeV. The lowest mass of the spectrum from the TABLE X. Calculation of the J P (C) = 1 +(−) state of the fully-bottom system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first two meson-meson structures give the theoretical results of the noninteracting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupledchannel calculation. All results are in units of MeV.

Structure
Channel fully-coupled calculation is 18864 MeV, which is slightly above the threshold of η b Υ. Therefore, no bound state exists for J P (C) = 1 +(−) . Again, we look for possible resonances by performing a complete coupled-channel calculation with CSM. The complex energy spectrum is shown in Fig. 6. Apart from four continuum states, η b (1S)Υ(1S), Υ(1S)Υ(1S), The J P (C) = 2 +(+) channel: we have 9 channels in this case, including two di-meson channels, ΥΥ configurations in the singlet-and hidden-color states, one diquark-antidiquark channel [bb] 6c [bb]6 c , and six K-type ones. The calculated lowest mass for ΥΥ is 18926 MeV, which is equal to its theoretical threshold. The masses of the other channels are all above 19.1 GeV, and there is an approximate degeneracy among the channels [bb] 6c [bb]6 c , K 3 , and K 4 . The lowest mass of the spectrum from the fully-coupled calculation is 18926 MeV, which is slightly above the threshold of ΥΥ. Therefore, no bound state exists for J P (C) = 2 +(+) . Fig. 7 shows the complex energy spectrum for J P (C) = 2 +(+) , obtained by using CSM. We find two resonance poles above the Υ(1S)Υ(2S) threshold. They can be identified as Υ(1S)Υ(2S)(19473) and Υ(1S)Υ(2S)(19633), with their widths 0.42 MeV and 0.40 MeV, respectively.

C. Charm-bottom systemcccb
For the charm-bottom systemcccb, narrow resonances are found in J P = 0 + , 1 + , and 2 + states, respectively. Table XII summarizes the results and the details are discussed in the following.
The J P = 0 + channel: Table XIII lists the lowest masses of each four configurations ofcccb tetraquark TABLE XI. Calculation of the J P (C) = 2 +(+) state of the fully-bottom system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The value in the parenthesis after the first meson-meson structure gives the theoretical result of the noninteracting meson-meson threshold. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV. in 0 +(+) states. In the real-range calculation, θ = 0 • , the two dimeson color-singlet channels η c B c and J/ψB * c are just scattering states, and their hidden-color cases present masses about 9.6 GeV. The two diquarkantidiquark structures are almost degenerate at 9.67 GeV. As for the other K-type configurations, which are all excited channels, the calculated single-channel masses are generally in a mass region of 9.45 − 9.69 GeV. Then, in the coupled-channel investigations on each configuration, bound states are still not found, the lowest mass is just the theoretical threshold value of η c B c at 9243 MeV, and the other excited configurations are generally at 9.5 GeV. This fact remains in a complete coupled-channels calculation, whose lowest mass at 9243 MeV is shown in the last row of Table XIII. In a further step, when a complex-scaling study is considered in a fully-coupled calculation, one narrow resonance state is found. Fig. 8 presents the distributions of the complex energies ofcccb tetraquarks with rotated angle θ varied from 0 • to 6 • . In a mass region from 9.2 GeV to 10.2 GeV, the scattering nature of η c B c and J/ψB * c both in the ground and first radial excitation states are well established. One can also see in Fig. 8 a stable resonance pole, whose mass and width are 9813 MeV and 4.02 MeV, respectively. This resonance is below the threshold lines of η c (1S)B c (2S), and the dominant channel, therefore, can be identified as the J/ψ(1S)B * c (1S) state. The J P = 1 + channel: There are in all 30 channels listed in Table XIV to be considered for thecccb tetraquark in 1 + state. In particular, six channels devote to the meson-meson structures, three channels belong to the diquark-antidiquark ones, and 21 channels correspond to the K-type configurations. When they are considered as isolated channels or when a coupled-channels calculation in each configurations is performed, no bound state is found. The lowest mass, 9317 MeV, is the theoretical threshold value of η c B * c ; the other two dimeson channels, J/ψB c and J/ψB * c , are at 9377 MeV and 9451 MeV, respectively; the exotic configurations, i.e. hiddencolor, diquark-antidiquark and K-types structures, are generally located in a energy range from 9.62 GeV to TABLE XIII. Calculation of the J P = 0 + state of thecccb system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first two meson-meson structures give the theoretical results of the noninteracting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.

Structure
Channel Lowest mass of the fully-coupled result: 9243 9.68 GeV; and the lowest mass in each of the conficuration's coupled-channels calculation is around 9.60 GeV.
However, several resonances are obtained when the complex scaling method is employed in a complete coupled-channels calculation. Figure 9 presents our results; in particular, the top panel shows the distributions of η c B * c , J/ψB c and J/ψB * c continuum states, their ground and first radial excitation states basically align along the corresponding threshold lines. Nevertheless, two stable resonance poles are circled therein and they can be identified as J/ψ(1S)B * c (1S)(9823) and The J P = 2 + channel: We have ten channels listed in Table XV. There is only one channel for the dimeson structure, J/ψB * c , another one for the diquarkantidiquark case, [cb] 3c [cc]3 c , and the rest belong to the K-type configurations. We obtain a lowest mass of 9451 MeV, which is just the theoretical threshold value of J/ψB * c . Generally, the unconventional structures have masses which lie within the range 9.65−9.71 GeV. Moreover, no bound state is found when coupled-channel effects are taken into account in each tetraquark configuration separately. Figure 10 shows the distribution of the complex energies obtained for thecccb tetraquark system. The nature as scattering states of the J/ψB * c system in ground and first radial excitation are well shown. However, from 9.45 GeV to 10.65 GeV energy region, a narrow resonance pole is obtained at 10.56 GeV, and its width is 6.61 MeV. Due to the location in the complex plane, we can identify it as a ψ(2S)B * c (1S)(10563) resonance. We find some resonance structures with quantum numbers J P = 1 + and 2 + states for the charm-bottom system bbbc. Table XVI summarizes our findings, the details of the calculation are discussed in the following.
The J P = 0 + channel: There are 20 channels listed in Table XVII. In particular, there are two meson-meson channels, η bBc and ΥB * c , two diquark-antidiaquark ones, T h e f i r s t r a d i a l e x c i t a t i o n s t a t e s . [cb] 6c [bb]6 c and [cb] 3c [bb]3 c , and 14 channels of K-type configurations. Firstly, single-channel calculations deliver a lowest mass of 15676 MeV for the η bBc , confirming its nature as a scattering state. The color-singlet channel of ΥB * c configuration is also a scattering state with a mass of 15812 MeV. The remaining excited channels are located within 15.92−16.48 GeV. Although the coupledchannels effect helps in pushing down the lowest masses of each configuration, they are still above the threshold values of η bBc and ΥB * c . Moreover, when a complete coupled-bases calculation is performed, the coupling is still too weak to have a bound state.
A complex range analysis on the fully coupled-channels calculation is done and the results are presented in Fig. 11. One can notice that, within a mass region of 15.  XV. Calculation of the J P = 2 + state of thecccb system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The value in the parenthesis after the first meson-meson structure gives the theoretical result of the noninteracting meson-meson threshold. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV. The J P = 1 + channel: There are 30 channels when considering thebbbc tetraquark system in J P = 1 + . Table XVIII shows the masses of the color-singlet channels of dimeson structures, η bB * c (15750 MeV), ΥB c (15738 MeV) and ΥB * c (15812 MeV); these results indicate the nature of scattering states. As for the other excited channels, their calculated masses generally locate in the range between 15.95 GeV and 16.44 GeV. In the coupledchannels studies of each configuration, the lowest mass of color singlet di-meson structures remains to be 15738 IG. 10. Complex energy spectrum of thecccb system with J P = 2 + from the complete coupled-channel calculation with CSM. The parameter θ varies from 0 • to 6 • .
IG. 11. Complex energy spectrum of thebbbc system with J P = 0 + from the complete coupled-channel calculation with CSM. The parameter θ varies from 0 • to 6 • .
MeV, and the other exotic cases are about 16.0 GeV. Meanwhile, last row of Table XVIII indicate that the lowest channel, ΥB c , is still unbound even in a complete coupled-channels calculation.
In order to find possible resonance states, as those obtained in the previous tetraquark systems, a CSM is employed in the fully coupled-channels computation, and the results are shown in Fig. 12. Within a energy gap of 15.7−16.5 GeV, the three ground states, η bB * c , ΥB c and ΥB * c , along with their first radial excitations, are well shown in the top panel. Most of the complex-energy dots are well align along the corresponding cut lines, reflecting the nature of scattering states. However, one can TABLE XVII. Calculation of the J P = 0 + state of thebbbc system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first two meson-meson structures give the theoretical results of the non-interacting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.

Structure
Channel notice that there is a stable resonance pole circled by an orange line, the calculated mass and width are 16031 MeV and 24.61 MeV, respectively. It can be identified as a Υ(1S)B * c (1S) resonance. In the bottom panel, we further enlarge the energy region between 16.3 GeV and 16.5 GeV. One can see that the five radial excited states are clearly shown as scattering states. However, one narrow pole, which is independent of the rotated angle θ, is obtained. Thus, it can be identified as a Υ(2S)B c (1S) resonance state with a mass and width 16303 MeV and 1.80 MeV, respectively. The J P = 2 + channel: Twelve channels are under our investigation in the highest spin state ofbbbc tetraquark system, and the results are listed in Table XIX. First of all, in each single channels calculations, the lowest mass is 15812 MeV for the color-singlet channel of ΥB * c state, and this is also the theoretical threshold value. The other excited channels masses are all above 16.0 GeV. The fact of unbound state ofbbbc tetraquark system holds for different kinds of coupled-channels cal-

E. Charm-bottom systemcbcb
For the charm-bottom systemcbcb, some resonance structures are found with quantum numbers J P (C) = 0 +(+) , 1 +(−) and 2 +(+) . Table XX summarizes the results, and the details are discussed in the following. XIX. Calculation of the J P = 2 + state of thebbbc system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The value in the parenthesis after the first meson-meson structure gives the theoretical result of the noninteracting meson-meson threshold. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.
The J P (C) = 0 +(+) channel: There are 18 channels, as listed in Table XXI IG. 13. Complex energy spectrum of thebbbc system with J P = 2 + from the complete coupled-channel calculation with CSM. The parameter θ varies from 0 • to 6 • . diquark-antidiquark and K-types configurations, the obtained single-channels masses are generally located in the energy range 12.71−12.89 GeV. When the coupledchannels calculations are considered, no bound state is found. Figure 14 shows our results when a complete coupledchannels calculation using the CSM formalism is performed.
In an energy region of 12.  IG. 14. Complex energy spectrum of thecbcb system with J P (C) = 0 +(+) from the complete coupled-channel calculation with CSM. The parameter θ varies from 0 • to 6 • .
The J P (C) = 1 +(−) channel: we have 23 channels as listed in Table XXII. They involve two di-meson con-TABLE XXI. Calculation of the J P (C) = 0 +(+) state of thē cbcb system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first two meson-meson structures give the theoretical results of the noninteracting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.

Structure
Channel figurations, B c B * c and B * c B * c in the singlet-as well as in the hidden-color wave functions, one diquark-antidiquark channel [bb] 3c [cc]3 c , and 18 K-type structures. Firstly, the meson-meson configurations in color-singlet channels are unbound and their calculated masses are 12624 MeV and 12698 MeV, respectively. The two hidden-color channels are characterize by a mass of ∼12.9 GeV, this value is also compatible for the diquark-antidiquark case. The K-type tetraquarks are at ∼12.87 GeV. Additionally, we do not find any bound state in coupled-channels calculations, which imply couplings in each certain tetraquarj configuration and a complete one.
By employing the CSM in a complete coupled-channel  XXII. Calculation of the J P (C) = 1 +(−) state of thecbcb system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first two meson-meson structures give the theoretical results of the noninteracting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.

Structure
Channel IG. 16. Complex energy spectrum of thecbcb system with J P (C) = 2 +(+) from the complete coupled-channel calculation with CSM. The parameter θ varies from 0 • to 6 • .  ( Fig. 1(e)) and K 4 -type ( Fig. 1(f)) configurations. Hence, only K 1 -and K 3 -types are employed in the following Tables without losing any generality. The J P (C) = 0 +(+) channel: There are 16 channels listed in Table XXV. In particular, the η c η b and J/ψΥ di-meson channels have masses at 12369 MeV and 12565 MeV, and so they are both scattering states. Their corresponding hidden-color channels are almost degenerate with masses close to 12.85 GeV. The diquark-antidiquark and K-type configurations are all located in the energy range 12.77−12.89 GeV, except the K 1 channel with a mass at 12.60 GeV. No bound state is available in any of the three possible calculations: single-channel, configuration-type channels calculation and full-channel computation.
Let us now focus on Fig. 17, which presents the distribution of complex eigenenergies of J P (C) = 0 +(+)c cbb tetraquarks in a complete coupled-channels calculation by CSM, in a mass region 12.3−13.5 GeV. The scattering TABLE XXV. Calculation of the J P (C) = 0 +(+) state of theccbb system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first two meson-meson structures give the theoretical results of the noninteracting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.  The J P (C) = 1 +(±) channel: Table XXVI shows the 22 channels available in this case. The three mesonmeson structures of η c Υ, J/ψη b and J/ψΥ in colorsinglet channels are all unbound, and their calculated masses are 12431 MeV, 12503 MeV and 12565 MeV, respectively. Their hidden-color channels are around 12.86 GeV. As for the other excited channels in diquarkantidiquark and K-type structures, they lie in the energy range 12.75−12.90 GeV. The coupled-masses of each cer- Additionally, an energy zoom shown in the bottom panel of Fig. 18 gives a more clear picture on the scattering nature for J/ψ(1S)η b (2S), J/ψ(1S)Υ(2S), η c (2S)Υ(1S), ψ(2S)η b (1S) and ψ(2S)Υ(1S) and the fact that no resonance state is found.
The J P (C) = 2 +(+) channel: 8 channels characterize this case and they are listed in Table XXVII: one meson-meson structure, J/ψΥ, two diquark-antidiquark structures, [cb] 3c [cb]3 c and [cb] 6c [cb]6 c , and four K-type configurations. The lowest computed mass, 12565 MeV, is at the theoretical threshold of J/ψΥ. The other excited channels locate within 12.77−12.93 GeV. Meanwhile, the coupled-channels effect is still too weak to form a bound state.

IV. SUMMARY
The fully-heavy tetraquark system QQQQ (Q = c, b) in its S-wave state, viz. the cccc, bbbb, cbcc, cbbb, bbcc and TABLE XXVI. Calculation of the J P (C) = 1 +(±) state of theccbb system in the real case (θ = 0 • ). The first and second columns show all possible structures and channels. They are obtained by the necessary bases in spin and color degrees of freedom summarized in the third column. The values in the parentheses after the first three meson-meson structures give the theoretical results of the noninteracting meson-meson thresholds. The lowest energies of all channels without considering any channel-channel coupling are listed in the fourth column. The last column shows the lowest energies in a calculation taking into account the couplings among the channels belonging to each structure (configuration). The last row gives the lowest energy from the complete coupled-channel calculation. All results are in units of MeV.

Structure
Channel cbcb tetraquark lowest-lying states with quantum numbers J P (C) = 0 +(+) , 1 +(±) , and 2 +(+) , have been systemically investigated by using a potential model which mimics the lattice-QCD findings about the interaction between a heavy quark and antiquark pair. That is to say, the quark model involves the so-called Cornell potential supplemented by a spin-spin dependent interaction. Furthermore, the four-body problem is solved by a Gaussian expansion method where each relative motion is expanded by a Gaussian basis whose sizes are in geometric progression. The complex scaling method is used to identify bound, resonance and scattering states in the complex energy plane. And all four-body configurations, including meson-meson, diquark-antidiquark, and K-type structures, as well as the couplings among themselves, are comprehensively considered.
Several narrow resonances have been determined in any of the all-heavy tetraquark systems investigated. Tables IV, VIII, XII, XVI, XX and XXIV summarize our findings. Amongst all the results presented, the following are of particular interest: IG. 19. Complex energy spectrum of theccbb system with J P (C) = 2 +(+) from the complete coupled-channel calculation with CSM. The parameter θ varies from 0 • to 6 • .
prove the phenomenological potential model applied to multiquark systems.