A Non-Relativistic Model of Tetraquarks and Predictions for Their Masses from Fits to Charmed and Bottom Meson Data

We investigate a non-relativistic model of tetraquarks, which are assumed to be compact and to consist of diquark-antidiquark pairs. We fit, for the first time, basically all currently known values for the measured masses of 45 mesons, including both charmed and bottom mesons, to the model and predict masses of tetraquarks as well as diquarks. In particular, we find masses of four axial-vector diquarks, i.e. $qc$, $cc$, $qb$, and $bb$, where $q = u,d$, and 24 ground-state tetraquarks, including both heavy-light tetraquarks ($qc\overline{qc}$ and $qb\overline{qb}$) and heavy tetraquarks ($cc\overline{cc}$ and $bb\overline{bb}$). In general, our results for the masses of $qb\overline{qb}$, $cc\overline{cc}$, and $bb\overline{bb}$ are largely comparable with other reported results, whereas our results for the masses of $qc\overline{qc}$ are slightly larger than what has been found earlier. Finally, we identify some of the obtained predictions for masses of tetraquarks with masses of experimental tetraquark candidates, and especially, we find that $\psi(4660)$, $Z_b(10610)$, and $Z_b(10650)$ could be described by the model.


I. INTRODUCTION
The concept of hadrons was introduced in 1962 by Okun [1] and developed into the quark model in 1964 independently by Gell-Mann [2] and Zweig [3,4], describing ordinary mesons (qq) and baryons (qqq) in terms of quarks q and antiquarksq. In addition to the quark model, the possible existence of exotic hadrons, such as tetraquarks (qqqq or qqqq) and pentaquarks (qqqqq), consisting of four or more quarks was proposed in Gell-Mann's seminal work [2], but it was not until the beginning of the 21st century when the first claimed observations of exotic hadrons were made [5]. Today, a large amount of data, obtained at both electron-positron and hadron colliders, has provided evidence for the possible existence of such exotic hadrons. Concerning tetraquarks, the first discovery was made in 2003 by the Belle collaboration that observed a resonance peak at (3872.0 ± 0.6) MeV [6], which was named X(3872) (and now sometimes referred to as χ c1 (3872) [7]), and then confirmed by several other experiments (see e.g. the reviews [8,9] and references therein). Many proposed exotic hadrons only appear in one decay mode, although X(3872) can be observed in several other decay modes as was discovered by the BaBar [10], CDF [11], and D∅ [12] collaborations. Later, the ATLAS, CMS, and LHCb collaborations were able to contribute with a massive amount of data on the electrically-charge neutral X(3872) and its current mass is determined to be (3871.69 ± 0.17) MeV [7,13]. It is the most studied exotic hadron, but its nature is still fairly unknown. It has similar properties to the charmonium state cc and * lundhamm@kth.se † tohlsson@kth.se was first believed to be an undiscovered excited state of cc, but a closer investigation of the decay modes X(3872) → J/ψ π + π − and X(3872) → J/ψ ω shows violation of isospin [14,15], which is unusual for cc. If X(3872) is an exotic hadron, then a common description is that it contains two quarks and two antiquarks forming ucuc. However, it is still an open problem how it is bound together. Since the discovery of X(3872), many new exotic hadron candidates have been claimed to be observed with final states of a pair of heavy quarks and a pair of light antiquarks, which are labeled as X, Y , and Z states by experimental collaborations and collectively referred to as XY Z states [16]. Examples of candidates for XY Z states are Z c (3900) [17,18], Z c (4025) [19,20] (now known as X(4020) ± [7]), Z b (10610) [21], and Z b (10650) [21]. The dynamics of the XY Z systems involves both short and long distance behaviors of QCD, which make theoretical predictions difficult. Hence, many competing phenomenological models currently exist for such states, including lattice QCD, compact tetraquark states, molecular states, QCD sum rules, coupled-channel effects, dynamically generated resonances, and non-relativistic effective field theories (see Ref. [22] and references therein). Many models view the exact nature of the inner structure of tetraquarks to be compact and to consist of socalled diquark-antidiquark pairs [23]. A diquark is a bound quark-quark pair, whereas an antidiquark is a bound antiquark-antiquark pair. These pairs are not by themselves colorless, but are proposed in the context of tetraquarks to form colorless combinations. Exotic hadrons of such pairs are thus not ruled out by QCD, but cannot be accommodated within the naive quark model. Modeling of tetraquarks containing only heavy quarks is therefore of special interest and easier to study theoretically, since several assumptions can be justified.
Recently, the LHCb collaboration reported the observation of a doubly-charmed and doubly-charged baryon Ξ ++ cc [24], which has lead to further attention on heavyquark systems as the description of exotic hadrons. Many tetraquark candidates are not possible to describe within quark models, since they have electric charge, and therefore cannot be charmonium or bottomonium, but are potential candidates for hidden-charm or hidden-bottom tetraquarks, molecular systems of charmed or bottom mesons or hadroquarkonia (see Ref. [25] and references therein).
In this work, we will study a non-relativistic model describing tetraquarks as composed of diquarks and antidiquarks, which interact much like ordinary quarkonia. Performing numerical fits to masses of mesons, masses of tetraquarks and the underlying diquarks will be predicted. For the first time, we will use both charmed and bottom mesons in the same fit and data of in total 45 charmed and bottom mesons (e.g. charmonium, bottomonium, D mesons, and B mesons) will be considered. We will predict masses of 24 tetraquark states, which is more than what has previously been performed in the literature.
This work is organized as follows. In Sec. II, we present the non-relativistic diquark-antidiquark model describing tetraquarks that can predict their masses and describe the numerical fitting procedure of meson data. Then, in Sec. III, we perform numerical fits of this model and state the results of the fits, including predicted masses of diquarks and tetraquarks. We will also present a thorough discussion on the results obtained and comparisons to other works, both theoretical and experimental. Finally, in Sec. IV, we summarize our main results and state our conclusions.

II. MODEL AND FITTING PROCEDURE
In this section, the model of tetraquarks viewed as diquark-antidiquark systems is presented and the method used to prescribe some tetraquark states quantitative masses is derived. This is preformed by firstly considering a quark-antiquark system and describing the Hamiltonian of that system with an unperturbed one-gluon exchange (OGE) potential and a perturbation term taking the spin of the system into account. This gives rise to a model with four free parameters, which are then fitted to meson data. Secondly, the model is expanded to incorporate composite quark-quark systems, which are called diquarks (the antiquark-antiquark systems are called antidiquarks). Thirdly, with the masses of the diquarks determined, the initial stage of the model describing quark-antiquark systems is then used to describe the diquark-antidiquark systems, which are interpreted as bound states of tetraquarks. First, considering a model of the quark-antiquark system q1q2.
Second, extrapolating the model to also describe the quarkquark (or diquark) system q1q2. Third, modeling a tetraquark q1q2q1q2 in the same way as the quark-antiquark system, but with diquarks and antidiquarks as constituents.

A. Model Procedure
The modeling procedure can be outlined and summarized as follows: 1. Fitting a quark-antiquark model to meson data to obtain the parameters of the effective potential.
2. Using that set of parameters to determine the diquark and antidiquark masses by changing the color constant and the string tension of the potential.
3. Considering the diquarks and antidiquarks as constituents of the tetraquarks to predict the tetraquark masses, see Fig. 1 for a schematic overview of the modeling procedure.
We begin by considering the interaction between a quark and an antiquark. In quark bound state spectroscopy, a commonly used potential describing the unperturbed contribution is the so-called Cornell potential [26] V (r) = κα S r + br , where κ is a color factor and associated with the color structure of the system, α S the fine-structure constant of QCD, and b the string tension. The first term in Eq. (1), i.e. V V (r) ≡ κα S /r, is the Coulomb term and associated with the Lorentz vector structure. It arises from the OGE between the quarks. The second term in Eq. (1) is associated with the confinement of the system. A non-relativistic approach is legitimate under the condition that the kinetic energy is much less than the rest masses of the constituents, which is usually the case considering heavy-quark bound states. We formulate the Schrödinger equation in the center-of-mass frame. Using spherical coordinates, one can factorize the angular and radial parts of this Schrödinger equation. Now, 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. In the case that m ≡ m 1 = m 2 , it holds that µ = m/2. Thus, the time-independent radial Schrödinger equation can be written as with the orbital quantum number L and the energy eigenvalue E. Substituting ψ(r) ≡ r −1 ϕ(r), Eq. (2) transforms into Based on the Breit-Fermi Hamiltonian for OGE, one can include a spin-spin interaction on the form [27][28][29][30][31] (4) In this model, we incorporate the spin-spin interaction V S (r) in the unperturbed potential V (r) by replacing the Dirac delta function with a smeared Gaussian function, depending on the parameter σ, in the following way as performed in Ref. [32]. Now, Eq. (3) takes the simple form where the effective potential V eff (r) is given by taking into account the spin-spin interaction. Equation (6) can be solved numerically for the energy eigenvalue E and the reduced wavefunction ϕ(r). The mass M of the bound quark-antiquark system can then be expressed as Note that this model consists of five unknown free parameters, namely the masses m 1 and m 2 of the two constituents, the fine-structure constant α S of QCD, the string tension b, and the parameter σ of the spin-spin interaction.

B. Color Structure
Hadrons are only stable when the colors of their constituent quarks sum up to zero, and thus, every naturally occurring hadron is a color singlet under the group symmetry SU (3). This means that a hadron only occurs if the product color state of the constituent quarks decomposes to an irreducible representation with dimension equal to one. Mesons consist of quarks in the color triplet state 3 and antiquarks in the color antitriplet state3, yielding product color states, which can be decomposed to the following irreducible representations: including a color singlet 1, and thus describing a naturally occurring hadron. In our modeling procedure, we consider the system of a diquark consisting of two quarks and an antidiquark consisting of two antiquarks in the triplet state, yielding a decomposition into a color singlet.
The difference in color structure between the quarkantiquark and quark-quark systems allows us to extend the model of the quark-antiquark system to also be valid considering a quark-quark system by only changing the color factor κ and the string tension b. The SU(3) color symmetry of QCD implies that the combination of a quark and an antiquark in the fundamental color representation can be reduced to |qq : 3 ⊗3 = 1 ⊕ 8, which gives the resulting color factor for the color singlet as κ = −4/3 for the quark-antiquark system. When combining two quarks in the fundamental color representation, it reduces to |qq : 3 ⊗ 3 =3 ⊕ 6, i.e. a color an-titriplet3 and a color sextet 6. Similarly, when combining two antiquarks, it reduces to a triplet 3 and an anti-sextet6. Furthermore, combining an antitriplet diquark and a triplet antidiquark yields |[qq]−[qq] : 3⊗3 = 1⊕8, thus forming a color singlet for which the Coulomb part of the potential is attractive. The antitriplet state is attractive and has a corresponding color factor of κ = −2/3, while the sextet state is repulsive and a color factor of κ = +1/3. Therefore, we only consider diquarks in the antitriplet state. Thus, the effect of changing from a quark-antiquark system with color factor κ = −4/3 to a diquark system with color factor κ = −2/3 is equivalent of introducing a factor of 1/2 in the Coulomb part of the potential for the quark-antiquark system. It is common to view this factor of 1/2 as a global factor, since it comes from the color structure of the wavefunction, thus also dividing the string tension b by a factor of 2. For further details, see Ref. [31]. We apply this change of the color factor when considering diquarks. Given the parameters of the potential, we obtain the mass of the corresponding diquark in a similar manner as when considering the quark-antiquark system, only changing the string tension b → b/2 and κ → κ/2, due to the change in the color structure of the system, and thus finding the energy eigenvalues of the diquark systems.

C. Fitting Procedure
The fitting procedure of the model is described as follows: A fit of the four parameters of the model to experimental data is performed by finding the parameters v ≡ (m, α S , b, σ), where m ≡ m 1 = m 2 , that minimizes the function where N data is the number of experimental data and M exp,i is the experimental mass of the corresponding mass M model,i (v), which is given in the model as a function of v. Each term in Eq. (9) is then weighted with w i for each mass. Following Ref. [33], we will only consider w i = 1, giving the same statistical significance to all states used as input. It should be noted that choosing w i = 1, the χ 2 function in Eq. (9) will be dimensionful. However, we will choose to present the values of this χ 2 function (as well as individual pulls) without units.

A. Data Sets
The model will be numerically fit to five different data sets. First, a data set consisting entirely of charmonium mesons (in total 15 mesons). Second, a data set consisting entirely of bottomonium mesons (in total 15 mesons). Third, a data set consisting of D mesons (in total 8 mesons). Fourth, a data set of B mesons (in total 7 mesons). Fifth, a fit to both the charmonium and bottomonium meson data will be made (in total 30 mesons). A meson consisting of two charm quarks is a good candidate to fit to the model, since it has a relatively large constituent mass compared to light quarks, and therefore, a non-relativistic approach can be justified. Both charmonium and bottomonium are heavy mesons and well suited to the restrictions of this model. For reference, the data set of charmonium mesons is called I, the data set of bottomonium mesons II, the data set consisting of only D mesons III, the data set consisting of only B mesons IV, and finally, the data set containing both charmonium and bottomonium mesons V. In Tab. I, the data used are presented.

B. Numerical Fits and Results
In this subsection, the results of the fitted data sets, and subsequently, the resulting masses of different diquarks and tetraquarks are presented. The procedure can be divided into three main parts. First, fitting the model to each data set I-V to obtain five sets of parameter values for the free parameters of the model. Next, using the sets of parameter values obtained by fitting data sets I-IV to calculate the masses of different diquarks. In detail, the sets of parameter values obtained by fitting data sets I, II, III, and IV are used to calculate the masses of the cc, bb, qc, and qb diquarks, respectively, with q being either an up quark (u) or a down quark (d). Finally, using the calculated diquark masses to calculate the masses of different tetraquarks. The set of parameter values used for this computation is the one obtained by fitting data set V to the model.
The number of free parameters when fitting the model to data sets I and II is four, since the masses of the constituent quarks for those data sets are equal, i.e. m = m 1 = m 2 . When fitting the model to data sets III-V, where q = u, d). Note that we use a spectroscopic notation, where N denotes the principal quantum number, S the total spin quantum number, L the orbital quantum number, and J the total angular momentum quantum number.
we use the values for the constituent masses of the charm and bottom quarks obtained in the fits to data sets I and II, which means that the number of free parameters is three. Also, when considering data sets III and IV, we use the value 0.323 GeV as the constituent mass of an up quark or a down quark, which is taken from Ref. [34] (see also p. 1 in the review "59. Quark Masses" [7]).
In practice, we are solving Eq. (6) in the eigenbasis of the spin operators S, S 1 , and S 2 , thus effectively replacing the product S 1 · S 2 by where S, S 1 , and S 2 is the total spin, the spin of quark 1, and the spin of quark 2, respectively. However, note that this modeling procedure is able to split the masses of states with equal principal (N ), orbital (L), and spin (S) quantum numbers, but not different total angular momentum (J) quantum numbers, i.e. the model is independent of J. Solving the Schrödinger equation numerically is performed by assuming Dirichlet boundary conditions at r = 0 and r = r 0 when using Eq. (6). The value of the parameter r 0 is chosen so that the energy eigenvalue E is independent of r 0 up to five significant digits. This approach was inspired by the method described in Ref. [35]. Next, the minimization of Eq. (9) is initially carried out by performing a random search with n = 100 000 points in the parameter space spanned by the parameters v = (m, α S , b, σ). The conditions on the parameters are chosen to be 0.05 ≤ α S ≤ 0.70 , 0.01 GeV 2 ≤ b ≤ 0.40 GeV 2 , 0.05 GeV ≤ σ ≤ 1.50 GeV as well as 1.00 GeV ≤ m ≤ 2.00 GeV for data set I and 4.00 GeV ≤ m ≤ 5.00 GeV for data set II. Furthermore, for data sets III-V, the values of m obtained for data sets I and II are used as input values. After the initial random search, an iterative adaptive method, using the same technique but with narrower conditions on the parameters and significantly smaller number of points, is performed to optimize the coarse point found during the initial random search in order to obtain the (local) bestfit point that minimizes the χ 2 function in Eq. (9). In Tab. II, the resulting values of the χ 2 function for the five data sets I-V are presented as well as each pull of each meson is listed, and in Tab. III, the resulting parameter values for m, α S , b, and σ when fitting the model to the respective data sets are given.

Diquarks
Given the best-fit values for the free parameters of the model, we find the diquark masses by calculating the energy eigenvalues, changing κ → κ/2 and b → b/2 in order to compensate for the change in color structure of the quark-quark system (compared to the color structure of the quark-antiquark system, see the discussion in Subsec. II B). The sets of parameter values obtained when fitting the model to data sets I, II, III, and IV are used to calculate the masses of the cc, bb, qc, and qb diquarks, respectively. We consider only diquarks in the ground state N 2S+1 L J = 1 3 S 1 , which are known as axial-vector diquarks and named good diquarks by Jaffe [5]. In Tab. IV, the results for the four diquark masses are presented.

Tetraquarks
For tetraquarks, we consider them to be composites of (axial-vector) diquarks and antidiquarks and the interaction between the diquarks and the antidiquarks is assumed to be effectively the same as for ordinary quarkonia. Thus, the parameter set obtained when fitting data set V to the model is used in the effective potential for all tetraquarks. However, when considering cccc and bbbb tetraquarks, we also compute the tetraquark masses with the parameter sets found by fitting the model to data sets I and II, respectively. Since the diquarks and antidiquark are in the antitriplet and triplet color states, respectively, the color structure of tetraquarks has identical color structure as the mesons, and subsequently, the same color factor κ = −4/3. In addition, the same string tension b is used for tetraquarks as for the mesons. Thus, considering diquarks and antidiquarks as consituents of tetraquarks, the tetraquark masses can be calculated using the diquark masses in Tab. IV. In Tab. V, the results for the masses of 24 tetraquark states are presented.

C. Motivation of Parameters and Comparison with Other Works
Similar models to the model presented in this work have been proposed in Refs. [31,36]. In Ref. [31], the authors were using the same model as in this work, but also taking into account perturbation in the spin-orbit and tensor interactions, although considering only fullycharmed diquarks and tetraquarks (i.e. cc and cccc), and thus, their results can be compared to the results for the set of parameters fitted to data set I. The models are identical for those states, where the perturbation energy is zero. In Ref. [36], the same authors considered X(3872) (also known as χ c1 (3872) [7]) under the hypothesis that its constituents are consisting of a diquark qc and an antidiquark qc and their parameter values could therefore be compared to the parameter values found by fitting data set III. They fit the model in order to investigate if Z c (4430) could be an exited state of X(3872). Furthermore, we will compare the masses of diquarks and tetraquarks calculated in this model with those presented in Refs. [25,31,[37][38][39][40][41][42][43][44][45][46][47][48][49][50]. In Tabs. VI, VII, and VIII, we 1.887 · 10 −2 5.485 · 10 −3 2.182 · 10 −2 1.696 · 10 −4 1.694 · 10 −1 display the different comparisons.

D. Discussion on Results
A thorough discussion on the results obtained in this work is in order. In Tab. II, the pulls and the values of the χ 2 function from the five fits to data sets I-V are presented. Comparing the values of the χ 2 function among the fits, we observe that the value for data set V is the largest with an order of magnitude of 10 −1 , the values for data sets I and III with orders of magnitude of 10 −2 , the value for data set II with an order of magnitude 10 −3 , and finally, the value for data set IV is the smallest with an order of magnitude of 10 −4 . The discrepancy in the values of the χ 2 function between data sets IV and V could be a consequence of the much larger variation of the masses in data set V compared to the variation of the masses in data set IV. Also, one could expect that, when fitting the model to data set V, the value of the χ 2 function would be of the same order of magnitude as the ones when fitting the model to data sets I and II, since data set V consists of quarkonia, which are well suited for this model. Nevertheless, comparing the pull values obtained for data set V, we note that almost all charmonium mesons yield positive pull values and almost all The constituent diquarks are assumed to be in the ground state N 2S+1 LJ = 1 3 S1, which are to be found in Tab. IV. All tetraquarks are considered in the states N 2S+1 L = 1 1 S, N 2S+1 L = 1 3 S, and N 2S+1 L = 1 5 S. The label "Data sets n1 + n2" indicates that the input values for the diquark masses are adopted from data set n1 and the input values for the other parameters are taken from data set n2.
bottomonium mesons yield negative pull values, implying a skewed adjustment of the model to this data set. The smallest absolute value of the pulls from this data set is 3.144 · 10 −3 for χ b0 (1P ), whereas the largest absolute value of the pulls is 1.598 · 10 −1 for χ c0 (1P ). In general, the deviation in pull values is difficult to explain. It could originate from the fitting procedure being not suitable to assign the same parameter values for both charmonium and bottomonium mesons or simply by the inclusion of more data points contributing to the total value of the χ 2 function. Overall, data set IV fits the model the best and data set V the worst. In Tab. IV, the predicted values for the masses of the diquarks are given, and in Tab. VII, a comparison with other works is presented. In this work, the masses of the diquarks are dependent on the parameters of the effective potential obtained from fitting data sets I-IV to the model. Note that the parameters α S and σ obtained in fitting data set II-V sometimes assume the upperend values of the intervals in which they are allowed to vary, which could imply the existence of more suitable parameter sets for these data sets if the intervals constraining the parameters would be enlarged. However, compared with the values for the diquark masses of the different works presented in Tab. VII, they deviate with at most about 250 MeV and are generally in excellent agreement with the results in Refs. [31,44], which values are also predicted in the framework of nonrelativistic quark models. In Ref. [37], the diquark masses are studied by means of the so-called Schwinger-Dyson and Bethe-Salpeter equations, which take into account the kinetic energy as well as splittings in the spin-spin, spin-orbit, and tensor interactions. The predicted values for the diquark masses in this work are consistently smaller by about 100 MeV compared to the values in Ref. [37]. Relativistic models, such as the ones presented in Refs. [41-43, 45, 51], and models based on QCD sum rules, such as the one in Ref. [38], all predict larger diquark masses. The differences could be a consequence of the introduction of more and updated data in this work or relativistic effects may play a significant role, since such are not taken into account in this work. In Tab. V, the resulting mass spectrum for the ground states of the tetraquarks are presented. An overall feature is that lighter tetraquarks have a larger spread in the energy eigenvalues than heavier ones, giving a larger relative difference in the masses among the states for lighter tetraquarks compared with heavier ones. In Tab. VIII, our predicted values for tetraquark masses and the corresponding ones from other works are shown. Regarding qcqc tetraquarks, the results obtained in Refs. [40,41,48,51] all predict smaller masses for all states. In our model, the masses of qcqc tetraquarks are sensitive to the parameters used in the effective potential, which means that a possible explanation for this deviation could be the skew fit of the model to data set V. Also, a non-relativistic framework may not be suitable when considering heavylight tetraquark systems, since relativistic effects play a significant role in such systems. In general, the predicted masses for cccc tetraquarks are in good agreement with the values in Refs. [25,31,39,44,46,47,50]. However, the 1 1 S state differs by about 1 GeV in comparison to Ref. [49]. Furthermore, the predicted masses for qbqb tetraquarks are in excellent agreement with the values in Refs. [25,41,48], and the relative deviation among the masses of different tetraquark states is overall small for this type of tetraquarks. Concerning bbbb tetraquarks, the predicted masses are in very good agreement with Refs. [25,43,44], but consistently smaller by about 0.5 GeV-1.0 GeV compared to the values in Refs. [39,49,50]. For the heavy tetraquarks (i.e. cccc and bbbb), the predicted masses obtained in Refs. [39,49,50] are consistently and significantly larger than those obtained in this work. In Ref. [39], the color-magnetic interaction is adopted to calculate the masses, and in Refs. [49,50], a similar model to the one used in this work is considered, but the variational principle is applied when solving the Schrödinger equation. This difference in the modeling approach could be the reason for the differences in the results.

E. Comparison with Experimental Results
Considering experimental results, there are about ten candidates of tetraquarks that are listed in the particle listings of the Particle Data Group [7]. These experimental tetraquark candidates are χ c1 (3872) [ which are both potential qbqb tetraquarks. There exists a classification of tetraquark states based on "good" diquarks (N 2S+1 L J = 1 3 S 1 ) such that [25,51,52] Diquark M [MeV] Ref. [37] Ref. [31] Ref. [41] Ref. [42] Refs. [43,51] Ref. [44] Ref. [45] Ref. [38] qc Ref. [40]  investigate. Furthermore, if one considers spin-1 diquarks and antidiquarks (i.e. axial-vector diquarks), then one has only three possibilities for the total spin S = 0, 1, 2 of the tetraquarks, i.e. three wavefunctions for each state N 2S+1 L: N 1 S, N 3 S, N 5 S, N 1 P , N 3 P , N 5 P , N 1 D, N 3 D, and N 5 D, which are nine possibilities [31]. Note that the total angular momentum quantum number J is dropped from the spectroscopic states of the tetraquarks, since our model is independent of J. Therefore, we should compute the following eight interesting states 1 1 S, 1 3 S, 1 5 S, 1 1 P , 1 5 P , 1 1 D, 1 3 D, and 1 5 D (i.e. not 1 3 P included), which are all ground states (N = 1), and compare the experimental values for the masses of tetraquarks with our theoretically predicted values of the masses using the allowed ground states for each tetraquark candidate. In comparing our theoretical predications in Tab. V with the experimental values for the tetraquark masses (cf. Ref. [7]), we find the following agreement within 100 MeV for qcqc tetraquarks Thus, it seems that the most likely tetraquark candidate to be described with our model is ψ(4660) as either a 1 1 P state of mass 4582 MeV or a 1 5 P state of mass 4591 MeV (see Tab. V). Furthermore, the qcqc tetraquark candidates Z c (3900), X(3915), and ψ(4360) as well as the qbqb tetraquark candidates Z b (10610) and Z b (10650) could be described by our model. Unfortunately, the most studied tetraquark candidate χ c1 (3872) cannot be accommodated with any state in our model due to the lowest state having the mass value 4076 MeV (see Tabs. V and VIII). Finally, the cccc and bbbb tetraquarks are interesting objects to study in the sector of exotic hadrons. The results obtained in this work suggest that the mass of the fullycharmed tetraquark could be about 5960 MeV or above in its ground state, whereas the mass of the fully-bottom tetraquark could be as large as 18720 MeV (see Tabs. V and VIII).

IV. SUMMARY AND CONCLUSIONS
We have investigated a model of tetraquarks, assumed to be compact and to consist of diquark-antidiquark pairs, in a non-relativistic framework and predicted mass spectra for the qcqc, cccc, qbqb, and bbbb tetraquarks. Considering tetraquarks as bound states of axial-vector diquarks and antidiquarks, a simple model originally formulated for quarkonia has been adopted and used to calculate and predict the masses of different tetraquark states. For the first time, a total number of 45 mesons, including both charm and bottom quarks, and the most recent corresponding data on the masses of these mesons [7] have been used to fit the free parameters of the model. Particularly, we have found predictions for four axial-vector diquark masses, and subsequently, a total number of 24 tetraquark masses, which are all presented in Tabs. IV and V. In comparison with other nonrelativistic models, our results for the cccc, qbqb, and bbbb tetraquarks are shown to be in excellent agreement with earlier results presented in the literature. However, considering qcqc tetraquarks, our results deviate slightly from earlier results and the predicted masses of these tetraquarks are consistently larger than the ones found in the literature. For the masses of heavy-light tetraquark states, i.e. qcqc and qbqb, we have been able to identify some of these states with experimentally proposed tetraquark candidates. One such identification includes the ψ(4660) tetraquark candidate, which can be proposed to be the qcqc tetraquark in either the state 1 1 P or the state 1 5 P . For qbqb tetraquarks, the tetraquark candidates Z b (10610) and Z b (10650) could both be identified with the state 1 3 S. Concerning the heavy tetraquark states, i.e. cccc and bbbb, the model predicts the mass of the fully-charmed tetraquark to be 5960 MeV and the mass of the fully-bottom tetraquark to be 18720 MeV, both values correspond to values obtained for their respective ground states. Finally, our model could also be used to predict masses for other potential tetraquark states for which no experimental data exist today.