Nonrelativistic Treatment of Fully-Heavy Tetraquarks as Diquark-Antidiquark States

The goal of the present work is to obtain a reliable estimate of the masses of the ground and radially excited states of fully-heavy tetraquark systems. In order to do this, we use a nonrelativistic model of tetraquarks which are assumed to be compact and consist of diquark-antidiquark pairs. This nonrelativistic model is composed of Hulthen potential, a linear confining potential and spin-spin interaction. We computed ground, first, and second radially excited $cc\bar{c}\bar{c}$ and $bb\bar{b}\bar{b}$ tetraquark masses. It was found that predicted masses of ground states of $cc\bar{c}\bar{c}$ and $bb\bar{b}\bar{b}$ tetraquarks are significantly higher than the thresholds of the fall-apart decays to the lowest allowed two-meson states. These states should be broad and are thus difficult to observe experimentally. First radially excited states are considerably lower than their corresponding (2S-2S) two-meson thresholds. We hope that our study may be helpful to the experimental search for ground and excited $cc\bar{c}\bar{c}$ and $bb\bar{b}\bar{b}$ tetraquark states.


I. INTRODUCTION
Quark model describes ordinary mesons as (qq) systems and baryons as (qqq) systems in terms of quarks q and antiquarksq. In addition to quark model, the existence of multiquark states such as tetraquarks (qqqq or qqqq), pentaquarks (qqqqq), and structures with more quarks was proposed decades ago [1][2][3][4]. These new structures present quantum numbers, masses, flavours, decay channels and widths in the experiments which cannot be fitted into conventional quark model. Now they are called exotic hadrons. Concerning tetraquarks, the first discovery of exotic states was made in 2003 by Belle Collaboration in the charmonium sector [5]. This state was named as X(3872) (now referred to name as χ c1 (3872)) with a mass of (3872.0 ± 0.6 MeV). Its extraordinary properties, such as isospin violation [6] and radiative decays [7,8] still make its structure unsolved. Many new exotic hadron candidates have been observed by decaying into final states of a charm and anticharm quarks. These states are called XY Z states. The X states are neutral and have positive parity quantum number and generally seen in J/ψ + pions decays. Typical example is the observation of χ c1 (3872). The Y states are neutral and have negative parity quantum number, and seen in e + e − annihilation with or without initial state radiation. Y (4260) is an example of this family [9]. The Z states are mostly charged or neutral, have typically positive parity and decay into J/ψ + π, h c (1P ) + π, χ c (1P ) + π. The charged Z c (3900) state is the famous example of Z family [10,11]. Observation of pentaquark states in 2019 made exotic hadrons more interesting [12].
The physics of exotic hadrons is a thorough piece of research which involves both short and long distance behaviors of QCD. At one side, increasing values of radius r are needed when considering the spectroscopy. This is the place where nonperturbative effects are at the stage. * hmutuk@omu.edu.tr On the other side hard processes, such as decays, occur at short distances, i.e., short values of radius r. This is the place where perturbative effects are at the stage. There are two different regimes which make theoretical predictions difficult. Many theoretical and phenomenological models are being studied to understand and interpret these exotic states such as lattice QCD, dynamically generated resonances, QCD sum rules, coupled channel effects and nonrelativistic effective field theories (see Ref. [13] and references therein).
Generally two pictures are taken into account in the approaches mentioned above: molecular and compact tetraquark pictures. Hadronic molecules are loosely bounded systems together by the exchange of pions and other light mesons. This scenario has received a lot of interest due to the masses of several XY Z hadrons are very close to the related meson-antimeson thresholds. In the case of χ c1 (3872) it has been suggested that if it has a binding energy of less than 200 keV with respect to the D 0D * 0 , according to the R = 1/ √ 2µE B , where µ is the reduced mass of the two-hadron system and E B is the binding energy, it would be at least as large as 10 fm [14]. From this perspective, hadronic molecules can be seen as extended objects. Compact tetraquarks are bound states of color nonsinglet diquark-antidiquarks, tightly bound by gluons, very much along the same lines as colored quark-antiquark pairs are bound into colorneutral mesons [15]. A diquark is a bound quark-quark (qq) pair, whereas an antidiquark is a bound antiquarkantiquark (qq) pair. These pairs are colored, i.e., have non-zero color charges and can have colorless combinations which turns out to be an ansatz of tetraquark paradigm. Tetraquark configurations are not ruled out by QCD. Indeed this context opens a new window of compact hadrons in QCD, which are even more numeorus than the conventional quark-antiquark mesons.
A very recent study of LHCb present a J/ψ-pair invariant mass spectrum by using pp collision data [16]. A narrow structure around 6.9 GeV/c 2 matching the lineshape of a resonance and a broad structure just above twice the J/ψ mass are observed. The deviation of the data from nonresonant J/ψ-pair production is above five standard deviations in the mass region between 6.2 and 7.4 GeV/c 2 , covering predicted masses of states composed of four charm quarks [16]. This energy range lie well above the experimentally known range for charmonium which is in the range of 3 -4.5 GeV and makes fullycharm tetraquark state cccc very interesting. The reason for this is that, the energy range of the XY Z states are in the same mass range of the conventional charmonium states and this can yield a confusion on these structures [17]. The observation of a possible cccc structure triggered many theoretical studies [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37].
In the side of bbbb structures, there is no observation up to now. An experimental study claimed the existence of a full-bottom tetraquark states bbbb, with a global significance of 3.6 σ and a mass around 18.4 GeV, almost 500 MeV below the threshold of ΥΥ [38]. However, LHCb Collaboration presented an intriguing analysis looking for the exotic bbbb tetraquark in the Υ(1S)µ + µ − final state and announced that no observation was made [39]. Prior to and after these experimental studies, the possible existence of bbbb state was investigated on the theoretical basis [40][41][42][43][44][45][46][47][48].
Motivated with this charm sector observation and a possible open window for fully-bottomed tetraquark states, in the present work we will use a nonrelativistic model to study tetraquarks as composed of diquarks and antidiquarks, which interact much like ordinary quarkonia. Owing to the masses of the valence degree of freedoms, fully-heavy four quark states can be investigated by nonrelativistic approach. We calculate mass spectra of fully-heavy tetraquark systems assuming that they are composed of doubly heavy quark (QQ) and antidiquark (QQ). Choosing this configuration is not just an assumption. Firstly, in light quark systems, the binding mechanism is maintained by the light-meson (such as pion, π) and gluon exchange. The binding mechanism in fully-heavy systems is probably dominated by the gluon-exchange forces since the typical gluon mass scale is m g ∼ 0.5 GeV. This mass value is much lighter than the possible force carriers of heavy-mesons that could be exchanged between the heavy diquark (QQ) and antidiquark (QQ) of tetraquark structure. Secondly, there is some evidence of diquark clustering in baryons. In 2017, the LHCb Collaboration reported the observation of a doubly charmed and doubly charged baryon, Ξ ++ cc a ccu state where the charm diquark may play a role in the structure, has lead to further attention on heavy-quark systems as the description of exotic hadrons [49]. This paper is organized as follows: In Section II, we describe nonrelativistic potential model of this work. In Section III, the masses of diquark/antidiquark and tetraquark systems are calculated. Detailed comparisons of diquark and tetraquark masses with previous studies within different approaches are given. Section IV is reserved for conclusion and summary of the obtained results.

II. POTENTIAL MODEL
The fundamental assumption in the quark potential model is that, if one integrate out the gluon fields in the QCD action, he/she can hopefully obtain an effective Hamiltonian with suitable potentials which describes the physics of hadrons fairly good. The quark masses in the resulting effective Hamiltonian might not be the same quark masses in the QCD Lagrangian, that is why they are named as constituent quark masses.
The use of quark potential models to describe the energy spectra of mesonic and baryonic systems gave reliable results. Phenomenological models with a simple relativistic kinetic energy term plus a scalar potential term which incorporates the so called linear confinement plus a term related to short distance which incorporates color-Coulomb interaction stemmed from QCD, give good results and descriptions of the observed spectra of both heavy and light quark mesons and baryons.
In the limit b → 0, the Hulthen potential approaches from above the Coulomb-like potential In other words, the Hulthen potential behaves like the Coulomb potential as r → 0 but decreases exponentially in the asymptotic region when r is sufficiently large. V C is the confining part of the potential with where c is a constant. V 0 is also a constant which would act as a zero-point energy. The parameters for this potential are given in the  To see whether or not V (r) (Eq. 1) coincides with the properties of Cornell potential exhibit, we can plot them. At first step, it can be shown that for large r, the Hulthen potential approaches to zero faster than the Coulomb-like potential, V (r) = − 4 3 αs r with α s = 0.5202 [17]. This can be seen in Figure 1. The choice of Hulthen potential could be explained as follows. As can be seen from Figure 1, the contribution of large r value of Hulthen potential is less than the Coulomb potential as r increases. This is due to the nature of short range potentials. Furthermore, the general theory of scattering is not immediately applicable to the Coulomb potential case because it decreases too slowly as the distance increases, which can be seen in Figure 1. The Coulomb potential falls off so slowly that it continues to influence the particles even as they move apart [66]. In order to apply the general theory of scattering to the Coulomb potential case, some modifications of the Coulomb potential have been considered. For example replacing Coulomb potential by the Yukawa potential and subsequently making the exponential factor approach unity is one example of this modification. Hulthen potential is a good approximation of Yukawa potential, V (r) = −V 0 e −αr r . Furthermore, Hulthen potential combined with a linear term, ∼ r n , can be a theoretical playground for the different forms of linear part of full potentials.
In addition to above discussion, the spectrum generated by the Hulthen potential could be investigated. Since this is done in literature excessively, it will be good to mention the well-known operator inequality. It was shown in [67] that for h b ≤ κ. This operator inequality between the Coulomb and Hulthen potentials is shared by the corresponding Hamiltonians and thus carries over to the entailed (ground state) energy eigenvalues.
On the other side, the plots of Cornell potential V (r) = − 4 3 αs r + br with α s = 0.5202 and b = 0.1463 GeV 2 [17], and potential model of this work are given in Figure 2. It can be seen from Figure 2 that in the range of 0 < r < 1 fm, both potentials have similar behaviour. In literature, all the phenomenological potentials have almost similar behaviour in the range of confinement, r < 1 fm. This is the characteristic interval of heavy quarkonium systems, charmonium cc and bottomonium bb. At its present status, it is not possible to obtain an exact interquark potential in the whole range of distances from the first principle of QCD.
A nonrelativistic approach with static potentials can be reasonable under the condition that the kinetic energy is much less than the rest masses of the constituents, which is usually the case when considering heavy quark bound states. Another advantage of assuming tetraquarks are composed from the heavy diquark (QQ) and heavy antidiquark (QQ) is that it significantly simplifies calculations, since the four-body problem of a four quark structure is very complicated. In the context of diquark-antidiquark configuration, the problem turns out to be a two two-body problem. In a two-body problem which involves a central potential, it is legitimate to work in the center-of-mass frame (CM), where one can use spherical coordinates to separate the radial and angular parts of the wave function. In terms of the reduced mass, let µ ≡ m 1 m 2 /(m 1 + m 2 ) where m 1 and m 2 are the constituent masses of quark 1 and quark 2, respectively. Since we are dealing equal masses of quarks, m ≡ m 1 = m 2 yields reduced mass as µ = m/2. Hence, the time-independent radial Schrödinger equation can be written as with the orbital quantum number L and the energy eigenvalue E. In order to account full spectra and splittings in energy between states with different quantum numbers, one can include spin-dependent terms to the potential. Due to the two nature of energies (high distances and short distances), it is not clear how the spin dependent forces change with distance. Assuming that spin dependent forces stem from short distance potential, observed hadron masses can be explained within well accuracy. For more about the nature of these spin dependent forces see Refs. [68,69]. The spin-dependent parts of the potential are sensitive to the Lorentz structure of the interquark potential. Hence, one has to split the total potential as the sum of Lorentz vector V v (r) and Lorentz scalar V s (r). Based on the Breit-Fermi interaction for one-gluon-exchange [70][71][72], the following spin dependent terms can be added to the potential. For equal masses m = m 1 = m 2 , with m being the constituent mass of the two-body problem (charm quark, bottom quark, or diquark) spin-spin interaction can be defined as This interaction is related to the S−wave splittings. This term has no effect for = 0 states. The expectation value of the operator for the spin-spin interaction can be calculated in terms of the spin quantum numbers by using where S, S 1 , and S 2 is the total spin, the spin of quark 1, and the spin of quark 2, respectively. The contribution of spin-dependent terms can be calculated by writing total potential Eq. (1) as sum of Lorentz vector structure and Lorentz scalar structure: The Lorentz structures are well known, for example, of Cornell potential but it is not clear in the Eq. (1). To account these Lorentz structures, it was assumed that the terms in Eq. (1) have partly vector and partly scalar [73]. In our formalism, as an alternative we use a different manner. The spin-spin correction to the nonrelativistic potential can be obtained as [74]: If the spin-spin interaction was treated as a first-order perturbation without the Gaussian smearing, as it is the case in this work, it would be proportional to modulus of the wavefunction at the origin, |Ψ(0)| 2 . Since spin-spin interaction only occurs in S−wave, only S−wave states (i.e., orbital angular momentum = 0) have non-zero value of the wavefunction at the origin. Therefore for S−wave state we have [70]: |R n, (0)| 2 can be obtained directly from the numerical calculations and is related to the radial potential as where µ is the reduced mass, µ = m 1 m 2 /(m 1 + m 2 ). Then it is possible to replace |Ψ(0)| 2 in Eq. (10) with |R n, (0)| 2 .
In QCD, the coupling constant α s is not actually a "constant". It changes according to the energy scale of each bound state. Therefore it changes with respect to energy scale and that is why it is called "running" coupling constant. We have used the nonrelativistic model of Ref. [65] where α s is constant and has different values for fitting quarkonium spectra. Choosing α s constant is a common approach in many of the nonrelativistic quark potential models.
By solving Schrödinger equation, one can obtain the energy eigenvalue E. The result will depend on the number of nodes of the wave function n (or principal quantum number N = n + 1) and the orbital angular momentum number . The mass calculated in this case is called spin-averaged mass. If one include spin-dependent terms in the potential, the result will also depend on the total spin S of the constituents spins S 1 and S 2 . The Schrödinger equation has no analytical solution for the potential given in this work. So we solve it numerically inspired from Ref. [75].
After having solved Schrödinger equation, the mass M of a particular state can be written as where m q is the mass of corresponding quark/antiquark and m d is the mass of corresponding diquark/antidiquark. We have considered the diquarks and antidiquarks as constituents of the tetraquarks to predict the tetraquark masses. Therefore, we first calculate the diquark masses and then tetraquark masses.

III. NUMERICAL ANALYSIS AND DISCUSSION
To study tetraquark systems in terms of diquark states, we will elaborate four-particle system as two two-body systems, i.e., diquark-antidiquark system. We assumed that the interaction between the diquarks and antidiquarks is to be effectively the same for ordinary quarkonia. In this picture, a diquark and antidiquark interact as a whole which means interactions between quarks from a diquark system (QQ) with antiquarks from an antidiquark system (QQ) are not considered. Since quark and antiquark masses are same one does not have to calculate antidiquark masses separately. We first calculate diquark masses (antidiquark masses at the same time) and use these diquark (antidiquark) masses to calculate tetraquark masses in diquark-antidiquark picture. The motivation for this factorization is the color structure of diquarks.
Hadrons can exist only when their total color charge of constituent quarks are zero. Technically, this means that every naturally occuring hadron is a color singlet under the group symmetry SU(3). A diquark is composed of two quarks (qq) whereas conventional quark-antiquark (also called quarkonium) is composed of (qq). The difference between quark-quark systems and quark-antiquark systems is mainly due to the color structure.
Regarding the Cornell potential, V (r) = κ αs r + br, where κ is the color factor, α s is the QCD fine structure constant and b is string tension and related to the strength of the confinement, the color factor change in the diquark configuration. The SU(3) color symmetry of QCD implies that, when we combine a quark and an antiquark in the fundamental color representation, we obtain |qq : 3 3 = 1 8. This representation gives the color factor for the color singlet as κ = −4/3 of the quark-antiquark system. When we combine two quarks in the fundamental color representation, it reduces to |qq : 3 3 =3 6, a color antitriplet3 and a color sextet 6. In a similar way, combining two antiquarks reduces to |qq :3 3 = 3 6 , a triplet 3 and antisextet6. Accordingly, combining an antitriplet diquark and a triplet antidiquark reduce to | [qq] − [qq] : 3 3 = 1 8, and form a color singlet for which the one-gluon exchange potential is attractive. The antitriplet state has a color factor κ = −2/3 which is attractive whereas the sextet state has a color factor κ = +1/3 which is repulsive. Therefore we will only consider diquarks in the antitriplet color state. This conclusion was achieved, for example by Ref. [76], in which it is shown that for single-flavor tetraquarks, only the antitriplet diquarks can build pure states.
This difference in color structure of the quarkantiquark and quark-quark systems make possible to extend the quark-antiquark model to be valid in the model of quark-quark system, just by changing the color factor κ and the string tension b. Changing color factor κ = −4/3 (for quark-antiquark system in color singlet state) to κ = −2/3 (quark-quark system in the antitriplet color state) is equivalent of introducing a factor of 1/2 in the Coulomb part of the Cornell potential for the conventional quark-antiquark system. This factor should be taken as a global factor since it comes from the color structure of the wave function. Therefore the string tension should also be divided by a factor of 2. The general rule for diquark potential from quark-antiquark potential is making V qq = V qq /2. This conclusion was done in different tetraquark models [17,[77][78][79]. So we also divide our potential (Eq. 1) by a factor of 2 for obtaining diquark spectra.

A. Diquarks
We now present our numerical calculations for masses of diquarks composed of charm (c) and bottom (b) quarks, which are also equivalent for antidiquark cases in this work. We assume that information in the spindependent interaction described in Eq. (7) is inherited when changing from quark-antiquark system to quarkquark system. In other words, spin-spin interaction is encoded in the diquarks. Similar assumptions were made in [17,79,80] and gave reliable results. In the present work, we choose the attractive color antitriplet state which is antisymmetric in the color wavefunction. In order to maintain Pauli exclusion principle, the diquark total spin must be 1, S = 1. So the total wavefunction of the diquark will be antisymmetric.
Based on the previous arguments, the results for the ground and first radially excited diquark masses are presented in Table II.   TABLE II. Results for the four diquark masses for cc and bb. Corresponding diquark masses are same as cc and bb, respectively. All results are in MeV unit.

Diquark
N 2S+1 LJ Mass cc 1 3 S1 3114 2 3 S1 3443 bb 1 3 S1 9792 2 3 S1 10011 Comparison of ground and first radially excited states for the cc diquark masses of earlier works are shown in Table III. As can be seen from Table III, masses for ground and first radially excited cc diquark states are in good agreement. We note that Ref. [17] used a Cornell potential with nonrelativistic framework and Ref. [81] used a nonrelativistic QCD motivated potential. We also compare ground states masses of cc diquark with the available theoretical studies in Table IV. IV. Comparison of cc diquark masses of this work and results from other works. All diquarks are considered to be in the ground state N 2S+1 LJ = 1 3 S1. All results are in MeV.

Reference
Mass This work 3114 [17] 3133 [26] 3226 [41] 3130 [43] 3204 [79] 3128 [81] 3130 [82] 3329 [83] 3144 [84] 3510 As can be seen from Table IV, cc diquark masses lie in the range of 3.1 − 3.5 GeV. Differences are due to the models of the references. In the present work, diquark masses depend on the parameters of the potential model, which were given in Section II. Compared with the values of the reference works, the results deviate with at most 400 MeV. Our result is in good agreement with the results in Refs. [17,41,79,81] which used nonrelativistic quark model. Ref. [26] used diquark-antidiquark picture in the framework of the relativistic quark model based on the quasipotential approach. In their model, quark-quark and diquark-antidiquark interactions in the potential are constructed similar to the mesons and baryons. Their cc diquark mass is approximately 100 MeV heavier than our result. In Ref. [43], diquark masses are calculated by a framework based on existing empirical information about mesons and baryons. The predicted value is approximately 100 MeV higher than our result. In Ref. [82], diquark masses are calculated by taking into account spinspin, spin-orbit, and tensor interactions in quark model. The predicted value of Ref. [82] is higher by about 200 MeV compared to our result. In Ref. [83], Schwinger-Dyson and Bethe-Salpeter equations are solved in order to obtain diquark masses by vector × vector interaction model (denoted as CI in reference paper) in which the obtained mass agree well compared to our result. A QCD Sum Rule study was done in Ref. [84] yielding a mass value which is consistently bigger by about 400 MeV than our result. Ref. [85] made an exploratory study with doubly heavy baryon which is composed of a heavy diquark and a light quark and calculated diquark masses in the Bethe-Salpeter (BS) formalism. They used two different parameters sets with confining parameter κ and coupling strength α s . They also took care of heavy quark limit. Mod. Ia refers m Q → ∞, i.e., taking the heavy quark limit while Mod. Ib refers m Q → finite, i.e., without heavy quark limit with appropriate values of κ and α s . The heavy quark limits are the same but parameters are different for Mod. IIa and Mod. IIb. Our results are at most 300 MeV lower than the results of Ref. [85].
In Table V, we compare bb diquark masses of ground and first radially excited states of different works. It can be seen from Table V that ground and first excited bb diquark masses are in good agreement. We note that Ref. [86] used the QCD potential of Buchmüller and Tye under nonrelativistic quark model framework.
As happened before in the cc diquark case, we compare ground state masses of bb diquark with the available theoretical studies in Table VI. The bb diquark masses lie in the range of 8.6 − 10.0 GeV. This gap is bigger with respect to cc diquark mass which was around 400 MeV. The reason for this could be the constituent mass of b quark, which is roughly three or four times c quark mass in constituent quark models. Nonrelativistic quark model results agree well with our result [17,41,79,81,86]. Ref. [87] calculated bbbb tetraquark mass in a nonrelativistic effective field theory and in a relativized diquark model characterized by one-gluon-exchange (OGE) plus a confining potential. In order to do this, they estimated bb diquark mass by binding the OGE plus a confining potential. Their result agree well with our result.
The differences of our results in cc and bb diquark masses apart from nonrelativistic formalisms could be as a result of the relativistic effects, spin-orbit and tensor interactions, model parameters or method uncertainties (there exist a typical 10% uncertainty in the results of QCD Sum Rules).

B. Tetraquarks
In tetraquark systems, we treated fully-heavy tetraquark as a two body (QQ −QQ) system with equal masses, m QQ = mQQ. We also treated them as axialvector states composed of diquarks and antidiquarks and we assumed that the interaction between the diquarks and antidiquarks mimickes the ordinary quarkonia interaction. It should be also noted that, tetraquarks made up of quarks of the same flavor can only be made up of axial-vector diquark-antidiquark pairs because of Pauli principle.
One of the important advantages of quark model is that, it allows to study systems with various combinations of principal quantum number, orbital angular momentum, spin and total angular momentum. Coupling of spin 1 diquark and spin 1 antidiquark produces tetraquark of total spin 0, 1 and 2 which we denote as S T = 0, 1, 2. We transferred quantum mechanic couplings into this nonrelativistic approach where total spin S T and orbital angular momentum L T couple into the total angular momentum J T . Obtaining charge and parity quantum numbers of tetraquark states is quite different than the conventional method. We obtain charge and parity quantum numbers of tetraquarks as explained in [17,80,88]. We will use the following notation where S d is the total spin of the diquark, Sd is the total spin of the antidiquark, S T is the total spin of the tetraquark, L T is the orbital angular momentum relative to the diquark-antidiquark system (in the twobody scheme), J T is the total angular momentum of the tetraquark comes from the coupling of S T and L T . Then, charge and parity quantum numbers can be written as Since we use diquarks/antidiquarks with spin 1 in the antitriplet color configuration, the resulting tetraquark states in S−wave can have the following possibilities: Considering diquarks and antidiquarks as consituents of tetraquarks, masses of the tetraquarks can be calculated using the diquark masses with N 2S+1 L J = 1 3 S in Table II for cc and bb diquarks. We also show the values of ∆ = M tetra − M threshold , where M tetra is the tetraquark mass and M threshold is the mass of its lowest meson-meson threshold. A negative ∆ means that the tetraquark lies below the threshold of the fall-apart decay into two mesons and thus should be a narrow state. Besides that a state with small positive ∆ value could also be observed as a resonance since its decay rate will be suppressed by the phase space. All other states with large positive ∆ values are expected to be broad and difficult to observe in the experiments. The results are presented in Table VII.
We see from this table that the predicted masses of ground and first radially excited cccc and bbbb tetraquarks lie significantly higher than the corresponding lowest (1S-1S) meson-meson thresholds. However masses of first radially excited cccc and bbbb tetraquarks are below than their corresponding (2S-2S) meson-meson thresholds. 1S cccc tetraquarks are located at ∼ 6.3 GeV and the mass splitting is about 30 MeV. J P C = 0 ++ 1S cccc state is about ∼ 300 MeV above than η c η c mass threshold while ∼ 100 MeV above than J/ψJ/ψ mass threshold. This suggests that J P C = 0 ++ cccc state might easily decay into η c η c and J/ψJ/ψ final states through quark rearrangements, respectively. J P C = 1 +− 1S cccc state lies about 270 MeV above than the η c J/ψ mass threshold. It might easily decay into η c J/ψ final state through quark rearrangements. J P C = 2 ++ 1S cccc tetraquark is about 200 MeV above than the mass threshold of J/ψJ/ψ, which suggests that it might easily decay into J/ψJ/ψ final state through quark rearrangements. The decays of 1S cccc tetraquarks are not suppressed dynamically or kinematically. Therefore these states should be broad and hard to be observed in experimental studies.
The previous discussion can be easily adapted to 1S bbbb tetraquarks. The J P C = 0 ++ , 1 +− , and J P C = 2 ++ 1S bbbb tetraquarks are located at ∼ 19.6 GeV and the mass splitting is about 10 MeV. This result is expected since spin-spin interaction term depend on mass as ∼ 1/m 2 . As the mass increases the contribution of spin-spin interaction term becomes small. J P C = 0 ++ 1S bbbb tetraquark is about ∼ 800 MeV above than η b η b mass threshold and ∼ 700 MeV above than ΥΥ mass threshold, respectively. This significant amount suggest that it might easily decay into η b η b and ΥΥ final states through quark rearrangements, respectively. J P C = 1 +− 1S bbbb state lies about 800 MeV above than the η b Υ mass threshold. This significant amount suggests that it might easily decay into η b Υ final state through quark rearrangements. J P C = 2 ++ 1S bbbb tetraquark is about 700 MeV above than the mass threshold of ΥΥ, and might easily decay into ΥΥ final state through quark rearrangements. The decays of 1S bbbb tetraquarks are not suppressed dynamically or kinematically. Therefore these states should be broad and difficult to be observed in experimental studies. This conclusion about 1S bbbb tetraquarks agree with current experimental data of LHCb [39] and CMS [89] Collaborations. No narrow beautiful tetraquarks in the Υ(1S)-pair production were observed in these experiments. Indeed, a lattice QCD study was performed to investigate the existence of a bbbb tetraquark bound state TABLE VII. Masses of ground, first, and second radially excited fully-heavy cccc and bbbb tetraquark states. E th is the threshold mass of two heavy QQ mesons. All results are given in MeV.
The obtained mass values of 2S cccc and bbbb tetraquarks are significantly above than the fall-apart decays to the lowest allowed (1S-1S) two-mesons. These states should be broad and difficult to be observed experimentally.
On the other hand, the predicted masses of radially excited 2S cccc and bbbb tetraquarks are lower than their corresponding thresholds. The J P C = 0 ++ , 1 +− , and J P C = 2 ++ 2S cccc tetraquarks are located at ∼ 6.6 GeV and mass splitting is about 30 MeV. The J P C = 0 ++ , 1 +− , and J P C = 2 ++ 2S bbbb tetraquarks are located at ∼ 19.8 GeV and mass splitting is about 10 MeV. Very recently, the LCHb Collaboration observed a broad structure at around 6.2-6.8 GeV which is very close to twice the J/ψ mass threshold, a narrow structure around 6.9 GeV, and a vague structure around 7.2 GeV [16]. We have several candidates of excited 2S cccc tetraquarks in the 6.2-6.8 mass range for the broad structure observed by LHCb Collaboration [16]. Narrow structure around 6.9 GeV might be orbital excitations of 2S or 3S cccc tetraquarks. A recent study about LHCb of J/ψ pairs was performed by using coupled channel methods in two channels J/ψJ/ψ and ψ(2S)J/ψ and in three channels J/ψJ/ψ, ψ(2S)J/ψ, and ψ(3770)J/ψ [93]. They found that in both channels, there is a hint of an existence of a near-threshold in the J/ψJ/ψ system with J P C = 0 ++ or J P C = 2 ++ quantum numbers having a mass around 6.2 GeV which they refer this state as X(6200). Our results of J P C = 0 ++ and J P C = 2 ++ cccc tetraquarks are ∼ 100 MeV above than their prediction. In Ref. [94], it was suggested that the structure around 7.2 GeV which was named as X(7200) can be related to χ c1 (3872) (formerly known as X(3872)) via heavy anti-quark di-quark symmetry. Orbital angular momentum addition to our calculation for 3S cccc tetraquarks could yield a candidate for X(7200).
The obtained mass values of 2S cccc and bbbb tetraquarks are significantly above than the thresholds of decays to the lowest (1S-1S) two-mesons. These states should be broad and are difficult to observe experimentally.
We compare our predictions for the masses of fullyheavy tetraquarks with the results of previous calculations. At first step we compare results of some works which considered excited fully-heavy tetraquarks. The results can be seen in Table VIII.
In Table VIII, Ref. [95] calculated the masses and wave functions of the exotic hadron states cccc and bbbb solving the four-body Schrödinger equation by including Cornell potential and spin-spin interaction in vacuum and strongly interacting matter in quark model. Ref. [96] used a relativized diquark model Hamiltonian built with relativistic kinetic energies, a one gluon exchange potential and linear confinement, in order to obtain masses of cccc and bbbb, and some other systems. Our results are compatible compared to reference studies. The reason of deviances in the results should be the different models used in the studies. At the second step, we compare our results for the ground states of fully-heavy tetraquarks of available studies. The results can be seen in Table IX. Ref. [97] used a nonrelativistic quark model within a potential model by including the linear confining potential, Coulomb potential, and spin-spin interactions in diquark-antidiquark picture. Ref. [98] studied the mass spectra of the S−wave fully-heavy tetraquarks in two nonrelativistic quark models. Color-magnetic interaction model (CMIM), a traditional constituent quark model (CQM) and a multiquark color flux-tube model (MCFTM) were used to investigate the fullyheavy tetraquarks in diquark-antidiquark scenario [99]. Ref. [100] used nonrelativistic chiral quark model using the Gaussian expansion method for ground state of bbbb tetraquark. Refs. [46,101] made QCD Sum Rule studies on fully-heavy tetraquarks. In another work, only fullycharmed tetraquarks were considered using a parameterized Hamiltonian with a sufficiently large but finite oscillator basis [102]. Ref. [103] calculated tetraquark masses in a constituent quark model, where the Cornell-like potential and one-gluon exchange spin-spin coupling were taken into account.
It can be seen from Tables VIII and IX that there are significant disagreements between the masses of fully-heavy cccc and bbbb tetraquarks. Especially from Table VIII, in the cccc sector, masses of J P C = 0 ++ , 1 +− , and J P C = 2 ++ lie in the range of 5.8 − 7.0 GeV, 6.0 − 6.8 GeV and 6.1 − 6.9 GeV, respectively. In the bbbb sector, masses of J P C = 0 ++ , 1 +− , and J P C = 2 ++ lie in the range of 18.7 − 20.0 GeV, 18.3 − 20.2 GeV and 18.3 − 20.2 GeV, respectively. Refs. [17,29,43,46,79,96,101] predict ground state fully-heavy tetraquark masses below or slightly above the corresponding two meson threshold masses which means these states could be observed in the experiments. On the other side, present model of this work and other approaches predict ground state fully-heavy tetraquark masses significantly above the corresponding thresholds, and thus it is difficult to observe these states in the experiments. 2S fully-heavy tetraquark masses are significantly above then the corresponding lowest (1S-1S) two-meson threshold masses. Concerning excited 2S fully-heavy tetraquark masses, our results are significantly below than the corresponding two-meson thresholds and there can be an open window for observing these structures in the experiments.
It would be good to analyze the response of the obtained results to minor variations of the numerical input of the potential model. The potential model has eight parameters. Therefore we divided into two group. In the first group we have changed mass m and α s values and fixed the rest of parameters. In the second group we have changed coupling strength h, range b, c and V 0 values and fixed the rest of parameters. The parameters of group 1 is shown in Table X.   Comparison of calculated fully-heavy tetraquark cccc and cccc bbbb masses for ground (1S) and radial excited states (2S, 3S) with quantum numbers J P C = 0 ++ , 1 +− and 2 ++ . All results are in MeV. State  1S  2S  3S  1S  2S  3S  1S  2S  3S  cccc  This work As can be seen from Tables XI and XII, the responses in the numerical input changes of group 1 did not significantly effect the results. The parameters of group 2 are given in Table XIII.   Table XIV for cccc tetraquark and in Table XV for bbbb tetraquark. As can be seen from Tables XIV and XV, the responses in the numerical input changes of group 2 did not significantly effect the results.
We believe that, the results of this work should serve as a benchmark for other studies, which are necessary to identify and develop theoretical models. This effort is important since predictions with respect to different models or approaches should be in well agreement and depend mildly to the corresponding model parameters.

IV. SUMMARY AND CONCLUDING REMARKS
In this work, adapting a perspective in which tetraquarks are assumed to be compact and consist of axial-vector diquarks and antidiquarks, masses of fullyheavy tetraquark states were calculated. We used a nonrelativistic potential quark model with the linear confining potential, Hulthen potential, and spin-spin interactions. We predicted ground, first and second radially excited mass spectra of fully-heavy tetraquarks, cccc and bbbb. In order to do this, we adopted a model which was formulated for quarkonia systems and made use of this model for diquark and tetraquark systems.
We have factorized four-body system into two twobody systems: quarks form diquark and antiquarks form antidiquarks, then diquark and antidiquark form tetraquark. By doing this, the complex structure of fourbody problem became more accessible. This flowchart allows us to study both ground and excited tetraquark states. Furthermore, various spin configurations were conducted easily in this manner.
At first step, we obtained diquark/antidiquark masses. Since diquark masses equal antidiquark masses, in the second step we computed tetraquark masses. In comparison with other models, the obtained cc and bb diquark masses are in good agreement with earlier results. This situation can be seen in Tables IV and VI. Indeed, radial excited states of cc and bb diquark masses agree well with available results in the literature which can be seen in Tables III and V. After obtaining diquark masses, we computed ground, first, and second radially excited cccc and bbbb tetraquark masses. The results are presented in Table VII. It was found that both ground states of cccc and bbbb tetraquarks masses are considerably higher than their corresponding lowest (1S-1S) twomeson decay threshold masses. Therefore they might easily decay into the lowest two-meson allowed states and these decays proceed through quark rearrangements. Such decays are not suppressed dynamically or kinematically. These states should be broad and are difficult to observe experimentally. We observed that first excited states of fully-charmed and fully-bottomed tetraquarks are lower than their corresponding (2S-2S) two-meson thresholds. These tetraquarks might be observed by obtaining enough center of mass energy and luminosity in the experiments. We also observed that the masses of 2S cccc and bbbb tetraquarks are significantly higher than the corresponding lowest (1S-1S) two-meson threshold masses. First and second radially excited masses of fully-heavy tetraquarks are compatible with the reference studies given in Table VIII. The comparison of ground state masses of cccc and bbbb tetraquarks are presented in Table IX. As argued before, the results deviate up to 1.0 GeV for cccc tetraquark and 2.0 GeV for bbbb tetraquark which are considered as high differences. Significant deviations between diquark and tetraquark masses may present a deeper perspective. In quark model, quarks are point-like objects. The dynamics, interactions and bound states for two and three quarks systems are described well. The diquarkantidiquark models require additional assumptions. In this perspective most diquark-antidiquark models assume diquarks and antidiquarks are point-like objects so that the interactions are effectively the same as for ordinary quarkonia systems. This was argued, for example in Ref. [26], in which diquarks have a structure (size) and this effect the results. In addition to this, chromomagnetic interaction inside the diquarks or between diquarks are hypothesized [80,104].
To sum up, we believe that our predictions for fully-charmed and fully-bottomed ground and excited states may help to experimental studies for a possible tetraquark signal in corresponding channels.