1 Introduction

Until the current millennium, the spectrum of detected hadrons was limited to systems that fit simply into the patterns typical of constituent-quark models [1, 2], i.e. quark–antiquark \((q {\bar{q}})\) mesons and three-quark (qqq) baryons. Notwithstanding this, Refs.  [1, 2] also raised the possibility of complicated hadrons, e.g. \(qq {\bar{q}}{\bar{q}}\) and \(q{\bar{q}} qqq\). A few years later, Lichtenberg suggested that “exotic mesons might be realized by replacing the quark and antiquark of a normal meson with a diquark and antidiquark, respectively” [3, 4]. Investigations on the phenomenology of \(qq {\bar{q}} {\bar{q}}\) mesons were conducted in the 1970s, when Jaffe proposed a four quark interpretation for light scalar mesons [5, 6] and Iwasaki the possible existence of hidden-charm tetraquarks [7]. Later, a similar study related to multiquark systems indicated the possible existence of baryon–antibaryon mesons [8]. Potential models to study \(qq{\bar{q}} {\bar{q}}\) systems were presented in Refs. [9, 10], and fully-heavy tetraquarks were hypothesized in Ref. [11].

Today, a large amount of data, obtained at both \(e^+ e^-\) and hadron colliders, has provided evidence for the possible existence of such exotic hadrons. In particular, we need to mention a recent LHCb study [12], where the \(J/\psi \)-pair invariant mass spectrum was studied by using pp collision data [12]. A narrow structure, dubbed the X(6900), which matches the lineshape of a resonance, and a broad one, next to the di-\(J/\psi \) mass threshold, were found. The global significance of either the broad or the X(6900) structure was determined to be larger than five standard deviations. It was stated that “The X(6900) structure could originate from a hadron state consisting of four charm quarks, \(T_{cc {\bar{c}} {\bar{c}}}\), predicted in various tetraquark models.” [12].

Herein we focus on those exotic systems which may be considered tetraquarks, viz. systems with meson-like quantum numbers that can be built using two valence quarks and two valence antiquarks. Some tetraquark candidates cannot be described using typical constituent quark models [13,14,15,16,17,18,19,20,21,22] because they carry electric charge; hence, cannot simply be \(q {\bar{q}}\) systems. Consequently, they are good candidates for: hidden-charm/bottom tetraquarks [5, 9, 23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41],or molecular systems constituted from a pair of charm/bottom mesons [42,43,44,45,46,47]; or hadro-charmonia [48,49,50,51].

Owing to the large masses of the valence degrees of freedom, the possible existence of fully-heavy \(QQ {\bar{Q}} {\bar{Q}}\) states (\(Q=c,b\)) and similar, mixed systems (\(c{\bar{c}} b{\bar{b}}\), \({\bar{c}} {\bar{c}} b b \sim c c {\bar{b}} {\bar{b}}\)) can reasonably be explored using nonrelativistic tools for QCD phenomenology and theory. Here, in contrast to systems involving light-quarks, for which both light-meson and gluon exchange may play a role in tetraquark formation, binding in fully-heavy systems is very probably dominated by gluon-exchange forces because the typical gluon mass-scale (\(m_g \sim 0.5\,\)GeV [52]) is much lighter than that of any necessarily-heavy meson that could be exchanged between two subsystems within the tetraquark composite. It is thus natural to suppose that the favoured structural configuration for a fully-heavy tetraquark state is diquark\(+\)antidiquark. Fully-heavy tetraquark states have been studied by some theorists since the 1980s [10, 11, 53,54,55], even though those investigations did not receive much attention due to the lack of experimental data.

We compute the spectra of \(cc{\bar{c}} {\bar{c}}\), \(c{\bar{c}} b{\bar{b}}\), \({\bar{c}} {\bar{c}} b b \sim c c {\bar{b}} {\bar{b}}\), \(b b {\bar{b}} {\bar{b}}\) tetraquarks from the diquark\(+\)antidiquark perspective, using a potential model characterised by linear confinement and one-gluon exchange. Using the same approach and assuming isospin symmetry, we also calculate the ground-state masses of similarly viewed \(b q {\bar{b}} {\bar{q}}\), \(b b {\bar{q}} {\bar{q}}\) systems (\(q=u,s\)).

Finally, as a test of the diquark model approach, we compute the masses of fully-heavy baryons in the diquark model. Our results may be compared soon to the forthcoming experimental data for fully-heavy baryons.

2 Relativized diquark model

We assume that the tetraquark states are colour-antitriplet (\({\bar{3}}_c\)) diquark + colour-triplet (\(3_c\)) antidiquark (\(D {\bar{D}}\)) states. Furthermore, the constituent D, \({\bar{D}}\) are each treated as being inert against internal spatial excitations [56,57,58,59]. This should be a fair approximation for fully-heavy systems owing to the suppression of quark exchange between the diquark subclusters in this case [60]. Consequently, dynamics within the \(D {\bar{D}}\) system can be described by a single relative coordinate \(\mathbf{r}_\mathrm{rel}\), with conjugate momentum \(\mathbf{q}_\mathrm{rel}\).

To describe the internal dynamics of a \(D_a{\bar{D}}_b\) system, we choose the Hamiltonian constrained elsewhere [33, 35]:

$$\begin{aligned} {\mathcal {H}}^\mathrm{REL}&= T + V(\mathbf{r}_\mathrm{rel})\,, \end{aligned}$$
(1a)
$$\begin{aligned} T&= \sqrt{{{\mathbf {q}}}_\mathrm{rel}^2+m_{D_a}^2} + \sqrt{{{\mathbf {q}}}_\mathrm{rel}^2+m_{{\bar{D}}_b}^2}, \end{aligned}$$
(1b)

with the interaction being the sum of a linear-confinement term and a one-gluon exchange (OGE) potential: [13, 35, 61, 62],

$$\begin{aligned} \begin{array}{rcl} V(r_\mathrm{rel}) &{} = &{} - \frac{3}{4} \Bigg [\beta r_\mathrm{rel} + G(r_\mathrm{rel}) + \frac{2 \mathbf{S}_{D_a} \cdot \mathbf{S}_{{\bar{D}}_b}}{3 m_{D_a} m_{{{\bar{D}}}_b}} \; \nabla ^2 G(r_\mathrm{rel}) \\ &{} &{}- \frac{1}{3 m_{D_a} m_{{\bar{D}}_b}} \left( 3 \mathbf{S}_{D_a} \cdot {{\hat{r}}}_\mathrm{rel} \; \mathbf{S}_{{\bar{D}}_b} \cdot {{\hat{r}}}_\mathrm{rel} - \mathbf{S}_{D_a} \cdot \mathbf{S}_{{\bar{D}}_b}\right) \\ &{} &{} \times \left( \frac{\partial ^2}{\partial r_\mathrm{rel}^2} - \frac{1}{r_\mathrm{rel}} \frac{\partial }{\partial r_\mathrm{rel}}\right) G(r_\mathrm{rel}) + \Delta E \Bigg ] \mathbf{F}_{D_a} \cdot \mathbf{F}_{{\bar{D}}_b}; \end{array}\nonumber \\ \end{aligned}$$
(2)

it is worth to note that this is the same as [35, Eq. (7a)]; here, we have only included the product of color matrices, \(\mathbf{F}_{D_a}\) and \(\mathbf{F}_{{\bar{D}}_b}\), explicitly.Footnote 1 The Coulomb-like piece is [13, 62]

$$\begin{aligned} G(r_\mathrm{rel}) = - \frac{4 \alpha _\mathrm{s}(r_\mathrm{rel})}{3 r_\mathrm{rel}} = - \sum _k \frac{4 \alpha _k}{3 r_\mathrm{rel}} \text{ Erf }(\tau _{ D_a {\bar{D}}_b \, k} r_\mathrm{rel}). \end{aligned}$$
(3)

Here, Erf is the error function and [13, 62]:

$$\begin{aligned} \tau _{ D_a {\bar{D}}_b\, k}&= \frac{\gamma _k \sigma _{D_a {\bar{D}}_b}}{\sqrt{\sigma _{D_a {\bar{D}}_b}^2+\gamma _k^2}}, \end{aligned}$$
(4a)
$$\begin{aligned} \sigma _{D_a {\bar{D}}_b}^2&= \frac{1}{2} \sigma _0^2 \left[ 1 + \left( \frac{4 m_{D_a} m_{{\bar{D}}_b}}{(m_{D_a} + m_{{\bar{D}}_b})^2}\right) ^4\right] \end{aligned}$$
(4b)
$$\begin{aligned}&\quad + s^2 \left( \frac{2 m_{D_a} m_{{\bar{D}}_b}}{m_{D_a} + m_{{\bar{D}}_b}}\right) ^2. \end{aligned}$$
(4c)
Table 1 Parameters specifying the Hamiltonian in Eq. (1). Here: \(n = u\) or d; and the superscripts “sc” and “av” indicate scalar and axial-vector diquarks, respectively. \(M_{n,s,c,b}\) are the valence quark masses

The parameters defining our Hamiltonian are listed in Table 1. The strength of the linear confining interaction, \(\beta \), and the value of the constant, \(\Delta E\), in Eq. (2) are taken from [35, Table I]; and in Eqs. (3), (4), the values of the parameters \(\alpha _k\) and \(\gamma _k\) (\(k = 1,2,3\)), \(\sigma _0\) and s are drawn from Refs.  [13, 62]. The valence quark masses are extracted from [13, Table II]. This leaves the diquark masses; they are determined by “binding” a \((q_1 q_2)^\mathrm{sc,av}\) system, \(\{q_1, q_2 = n, s, c, b\}\), \(n=u=d\), by means of a OGE potential; see [13] and [63, Table 1]. Hence, the results we subsequently report are parameter-free predictions.

If we hypothesize the existence of a dynamical baryon–meson supersymmetry in hadrons [64, 65], the Hamiltonian of Eq. (1) can also be used to calculate the spectrum of baryons in the quark–diquark approximation. In the diquark model, baryons and tetraquarks are described as quark–diquark and diquark–antidiquark configurations, respectively, whose internal dynamics is described in terms of a single relative coordinate; see e.g. Refs. [56,57,58]. Because of these premises: (I) the solution of either the diquark–antidiquark or quark–diquark eigenvalue problems is substantially the same; (II) not only the Hamiltonian of Eq. (1) is unchanged, but also the values of the diquark model parameters are expected to remain the same both in the baryon and tetraquark sectors, provided that one substitutes the valence antidiquark, \({\bar{D}}\), with a quark, q, with \(q = n, s, c\) or b.

3 Fully-heavy baryon and tetraquark states

3.1 \(bb{\bar{q}} {\bar{q}}\) and \(b q{\bar{b}} {\bar{q}}\) ground-state masses

As an exploratory exercise, we first compute the masses of \(J=0^{++}\) heavy-light tetraquarks – \(b q {\bar{b}} {\bar{q}}\), \(b b {\bar{q}} {\bar{q}}\) systems (\(q=n,s\)).

Using the Hamiltonian specified by Eq. (1) and the parameters in Table 1, we obtain the following ground-state masses:

$$\begin{aligned} M_{bn{\bar{b}} {\bar{n}}}^\mathrm{gs} = \left\{ \begin{array}{rl} 10.29\,(08) \text{ GeV } &{} (\text{ sc--sc } \text{ configuration}) \\ 10.12\,(11) \text{ GeV } &{} (\text{ av--av } \text{ configuration}) \end{array} \right. \end{aligned}$$
(5)

and

$$\begin{aligned} M_{bs{\bar{b}} {\bar{s}}}^\mathrm{gs} = \left\{ \begin{array}{rl} 10.52\,(08) \text{ GeV } &{} (\text{ sc--sc } \text{ configuration}) \\ 10.35\,(10) \text{ GeV } &{} (\text{ av--av } \text{ configuration}) \end{array} \right. , \end{aligned}$$
(6)

where the energies of the two possible \(D {\bar{D}}\) configurations are both shown, viz. scalar-scalar and axial-vector–axial-vector. Evidently, when combining \({\bar{3}}_c\) and \(3_c\) constituents, the OGE colour-hyperfine interaction favours a lighter av–av combination. This is because the spin-spin interaction in Eq. (2) is attractive. The nature of our uncertainty estimate is discussed in Sect. 3.3 (first procedure).

Turning now to \(bb{\bar{q}} {\bar{q}}\) systems, we obtain

$$\begin{aligned} M_{bb{\bar{n}} {\bar{n}}}^\mathrm{gs}&= 10.31\,(17)\, \text{ GeV }, \end{aligned}$$
(7a)
$$\begin{aligned} M_{bb{\bar{s}} {\bar{s}}}^\mathrm{gs}&= 10.53\,(16) \,\text{ GeV }. \end{aligned}$$
(7b)

(Owing to Pauli statistics, only av–av configurations are allowed in these cases.)

In Sect. 3.3 (second procedure), we suggest an alternative prescription to estimate the theoretical uncertainties on our results. If we make use of it, we get for the ground-state (av–av configuration): \(M_{bn{\bar{b}} {\bar{n}}}^\mathrm{gs} = 10.12\,(19)\) GeV, \(M_{bs{\bar{b}} {\bar{s}}}^\mathrm{gs} = 10.35\,(19)\) GeV, \(M_{bb{\bar{n}} {\bar{n}}}^\mathrm{gs} = 10.31\,(22)\) GeV, and \(M_{bb{\bar{s}} {\bar{s}}}^\mathrm{gs} = 10.53\,(21)\) GeV. These errors are slightly larger than those in Eqs. (57), even though the order of magnitude remains the same.

Table 2 Row 1: experimental meson-meson thresholds. They are \(2 \eta _c\), \(2 \eta _b\), \(\eta _b \eta _c\) and \(B_c {\bar{B}}_c\) in the \(cc{\bar{c}}{\bar{c}}\), \(bb{\bar{b}}{\bar{b}}\), \(bc {\bar{b}} {\bar{c}}\) and \(bb {\bar{c}} {\bar{c}}\) cases, respectively. Row 2: our computed results for the masses of \(J=0^{++}\) ground-state tetraquark systems. For comparison, the other rows list values obtained elsewhere. Masses are in MeV. All entries describe av–av configurations, except those marked by an asterisk, which are sc–sc. The entry highlighted by \(^\dagger \) was obtained previously [33]

3.2 \(cc{\bar{c}} {\bar{c}}\), \(bb{\bar{b}} {\bar{b}}\), \(b c{\bar{b}} {\bar{c}}\), \(b b {\bar{c}} {\bar{c}}\) ground-states

In \(QQ{\bar{Q}} {\bar{Q}}\) systems treated as colour triplet-antitriplet pair configurations, fermion statistics also precludes a role for scalar diquarks. Consequently, the ground-state \(cc{\bar{c}} {\bar{c}}\) is an av–av combination; and using Eq. (1) we find

$$\begin{aligned} M_{cc{\bar{c}} {\bar{c}}}^\mathrm{gs} = 5.88\,(17)\,\mathrm{GeV}. \end{aligned}$$
(8)

Using the same framework, the calculated mass of the analogous \(b b{\bar{b}} {\bar{b}}\) system is

$$\begin{aligned} M_{bb{\bar{b}} {\bar{b}}}^\mathrm{gs} = 18.75(07)\,\mathrm{GeV}. \end{aligned}$$
(9)

Table 2 lists our prediction for the mass of these systems alongside a sample of values obtained elsewhere [33, 66,67,68,69,70,71,72,73,74]. We also need to mention: Ref. [53], where the authors stated that for any confining potential there is no state below the threshold; Ref. [75], where the authors studied tetraquark states in a potential model calculation. They used a kind of individual-particle wave function, with a removal of the center-of-mass energy, \(H_\mathrm{cm} - \frac{3}{2} \hbar \Omega \); see [75, Eq. (1)]. \(\Omega \) was determined from \(c {\bar{c}}\) states and, as the four-body system is more diluted than a quark–antiquark one, there is a risk of over-correction; Ref. [76], where the authors investigated 4-quark states in a potential model with a chromo-electric interaction.

Considering the \(b c{\bar{b}} {\bar{c}}\) case, both sc–sc and av–av configurations can exist; and we find

$$\begin{aligned} M_{bc{\bar{b}} {\bar{c}}}^\mathrm{gs} = \left\{ \begin{array}{rl} 12.52\,(08) \text{ GeV } &{} (\text{ sc--sc } \text{ configuration}) \\ 12.37\,(09) \text{ GeV } &{} (\text{ av--av } \text{ configuration}) \end{array} \right. . \end{aligned}$$
(10)

One can also imagine \(J=0^{++}\) \(bb {\bar{c}} {\bar{c}}\) (\({\bar{b}} {\bar{b}} c c\)) configurations. In this case, only the av–av \((bb)_{{\bar{3}}_c}({\bar{c}}{\bar{c}})_{3_c}\) configuration is possible and its ground-state mass is

$$\begin{aligned} M_{bb{\bar{c}} {\bar{c}}}^\mathrm{gs} = 12.45\,(11)\, \text{ GeV } . \end{aligned}$$
(11)

Analogously to what is done at the end of Sect. 3.1, we make use of the alternative prescription from Sect. 3.3 (second procedure) to give a second estimate of the theoretical uncertainties. We get: \(M_{cc{\bar{c}} {\bar{c}}}^\mathrm{gs} = 5.88\,(17)\) GeV, \(M_{bb{\bar{b}} {\bar{b}}}^\mathrm{gs} = 18.75 (23)\) GeV, \(M_{bc{\bar{b}} {\bar{c}}}^\mathrm{gs} = 12.37\,(20)\) GeV, \(M_{bb{\bar{c}} {\bar{c}}}^\mathrm{gs} = 12.45\,(21)\) GeV. Again, even though the theoretical uncertainties resulting from the application of this second prescription are larger than those from Eqs. (811), they are still of the same order of magnitude.

3.3 Estimate of model uncertainty

In the following, we discuss two prescriptions to estimate the theoretical errors on our predictions of Eqs. (511).

First procedure In order to provide an estimate of our theoretical uncertainty, we also compute tetraquark masses using the Relativized Quark Model (RQM) Hamiltonian introduced in Ref.  [13]. Only a few obvious changes are necessary because this Hamiltonian was also constructed to bind a colour triplet-antitriplet pair into a colour-singlet system. The theoretical error on our results is then given by the difference between the predictions for the tetraquark masses obtained by using: I) the tetraquark model Hamiltonian of Eq. (1), with the values of the model parameters reported in Table 1; II) the RQM Hamiltonian [13], where the masses of the quarks are substituted with those of the diquarks, extracted from Table 1 herein.

When forming a S-wave system from two axial-vector constituents, the only contribution from spin-dependent interactions in the RQM is that produced by the contact term, \(V_\mathrm{cont}\):

$$\begin{aligned}&\left\langle S_1' \; S_2' \; S' \right| V_\mathrm{cont}(\mathbf{r}) \left| S_1 \; S_2 \; S \right\rangle = \left\langle 1 \; 1 \; 0 \right| V_\mathrm{cont}(\mathbf{r}) \left| 1 \; 1 \; 0 \right\rangle \nonumber \\&\quad \propto \frac{1}{2} \left( \mathbf{S}^2 - \mathbf{S}_1^2 - \mathbf{S}_2^2 \right) = -2. \end{aligned}$$
(12)

Contrarily, in the case of tensor, \(V_\mathrm{tens}\), and spin-orbit, \(V_\mathrm{so}\), interactions, one obtains the matrix elements

$$\begin{aligned}&\left\langle S' \; L' \; J' \right| V_\mathrm{tens}(\mathbf{r}) \left| S \; L \; J \right\rangle = \left\langle 0 \; 0 \; 0 \right| V_\mathrm{tens}(\mathbf{r}) \left| 0 \; 0 \; 0 \right\rangle \nonumber \\&\quad \propto \left\langle L' \right| Y^{(2)} \left| L \right\rangle \propto \left( \begin{array}{ccc} 0 &{}\quad 2 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \end{array}\right) = 0, \end{aligned}$$
(13)

where \(Y^{(2)}\) is a \(L=2\) spherical harmonic [77], and

$$\begin{aligned}&\left\langle S' \; L' \; J' \right| V_\mathrm{so}(\mathbf{r}) \left| S \; L \; J \right\rangle = \left\langle 0 0 \; 0 \right| V_\mathrm{so}(\mathbf{r}) \left| 0 \; 0 \; 0 \right\rangle \nonumber \\&\quad \propto \sqrt{L (L+1) (2L+1)} = 0. \end{aligned}$$
(14)

The smearing function coefficient employed in Ref.  [13], \(\sigma _{C_1 C_2}\), with \(C_{1,2}\) denoting the constituents, is the same as that we use, given by Eq. (4).

As an illustrative example, consider the fully-b \(J=0\) tetraquark. Our prediction for the ground-state mass is reported in Eq. (9). Using the RQM Hamiltonian and a computed value of \(\sigma _{(bb)({\bar{b}} {\bar{b}})} = 77.9\) fm\(^{-1}\), one finds

$$\begin{aligned} E^\mathrm{gs, RQM}_{bb{\bar{b}} {\bar{b}}} = 18822\,\mathrm{MeV}. \end{aligned}$$
(15)

Our mass prediction cannot be judged more accurate than the difference between this result and that in Eq. (9), viz. 74 MeV. We therefore list this value as the uncertainty in Eq. (9).

Second procedure Below, we discuss an alternative prescription to estimate the errors on our diquark model results. In this second scheme, the errors are extracted by first providing estimates of the uncertainties on the model parameters in Table 1 and then by propagating those theoretical errors on the tetraquark model results. We start by evaluating the theoretical errors on the diquark masses.

The diquark masses are computed by binding a quark–quark pair via the Relativized QM Hamiltonian [13, 63]. The errors on both the quarkonium and diquark masses can be estimated by comparing the RQM predictions of Refs. [13, 78,79,80,81,82] to the existing experimental data for heavy and light quarkonia [83] and extracting the average mass deviations between them. In the case of heavy and heavy-light mesons, we estimate the deviation to be of the order of \(0.5\%\). We thus get: \(M_{cc}^\mathrm{av} = 3329 \pm 17\) MeV, \(M_{bn}^\mathrm{sc} = 5451 \pm 27\) MeV, \(M_{bn}^\mathrm{av} = 5465 \pm 27\) MeV, \(M_{bs}^\mathrm{sc} = 5572 \pm 28\) MeV, \(M_{bs}^\mathrm{av} = 5585 \pm 28\) MeV, \(M_{bc}^\mathrm{sc} = 6599 \pm 33\) MeV, \(M_{bc}^\mathrm{av} = 6611 \pm 33\) MeV and \(M_{bb}^\mathrm{av} = 9845 \pm 49\) MeV. In the case of light diquarks, \(M_{nn}^\mathrm{av}\) and \(M_{ss}^\mathrm{av}\), the previous procedure may underestimate the theoretical uncertainties. We thus assign the \(M_{nn}^\mathrm{av}\) and \(M_{ss}^\mathrm{av}\) masses an error of 17 MeV, which is the same as that on \(M_{cc}^\mathrm{av}\). We have: \(M_{nn}^\mathrm{av} = 840 \pm 17\) MeV and \(M_{ss}^\mathrm{av} = 1136 \pm 17\) MeV.

The error on the linear confining potential parameter, \(\beta \), is roughly estimated as the difference between the value used herein and in Refs. [33, 35] and that from Ref. [13]. We get: \(\beta = 3.90 \pm 0.72\) fm\(^{-2}\). We do the same for the constant \(\Delta E\): \(\Delta E = - 370 \pm 117\) MeV. The other model parameters, \(\alpha _i\) and \(\gamma _i\) (with \(i = 1, 2, 3\)), \(\sigma _0\) and s, assume the same values both in the RQM [13] and the diquark model of Refs. [33, 35]. Therefore, in first approximation, we assume that their values have no theoretical error. By doing this, we somehow underestimate the theoretical error in this second procedure.

Table 3 Spectra obtained by solving the eigenvalue problem defined by Eq. (1). Following Sect. 3.3 (first procedure), the model uncertainty in each result is \(\lesssim 2\,\)%. The states are labelled thus: N is the radial quantum number (\(N=1\) is the ground state); \(S_\mathrm{D}\), \(S_{\bar{\mathrm{D}}}\) are the spin of the diquark and antidiquark, respectively, coupled to the total spin of the meson, S; the latter is coupled to the orbital angular momentum, L, to get the total angular momentum of the tetraquark, J. Degenerate states are orthogonal combinations of diquark+antidiquark spin vectors [35]

Finally, we can use the model parameters of Table 1 with their theoretical errors, see above, to calculate the tetraquark masses with their theoretical uncertainties. For example, in the fully-b tetraquark case of Eq. (9) we have: \(M_{bb{\bar{b}} {\bar{b}}}^\mathrm{gs} = 18.75\,(23)\) GeV.

Table 4 As Table 3, but for \(b c {\bar{b}} {\bar{c}}\) states

3.4 Complete tetraquark spectra

The Hamiltonian in Eq. (1) predicts a rich spectrum; and in Tables 3, 4 we report the lightest states in the spectra of \(cc{\bar{c}} {\bar{c}}\), \(b b {\bar{b}} {\bar{b}}\), \(b c{\bar{b}} {\bar{c}}\) and \(b b {\bar{c}} {\bar{c}}\) systems. These results should serve as useful benchmarks for other analyses, which are necessary in order to identify model-dependent artefacts and develop a perspective on those predictions which might only be weakly sensitive to model details. Moreover, given that the decay modes of fully-b \(J=0^{++}\) tetraquarks may be difficult to access experimentally [84,85,86], our predictions for orbitally-excited and \(J\ne 0\) tetraquarks may serve useful in guiding new experimental searches for fully-heavy four-quark states. Our predictions for the masses of the fully-charm states can be compared to the recent LHCb experimental results [12]. In particular, it can be noticed that we have several candidates in the 6.2–7.4 GeV energy range with \(0^{++}\), \(1^{++}\) and \(2^{++}\) quantum numbers.

3.5 \(\eta _c\), \(\eta _b\) and \(B_c\) ground-state masses

The Relativized QM [13] and Relativized Diquark Model [33, 35] Hamiltonians can also be used to calculate the masses of fully-heavy mesons, including the \(\eta _c\), the \(\eta _b\) and the \(B_c\) states. These estimates can provide further informations regarding the theoretical uncertainties in our approach.

Starting from the \(\eta _c\), we have

$$\begin{aligned} M_{\eta _c} = 2.91\,(06)\,\mathrm{GeV}. \end{aligned}$$
(16)

The Hamiltonian that we use here is that of Eq. (1) where, according to hadron supersymmetry, we can substitute the valence diquark/antidiquarks with antiquarks/quarks; the model parameters, including the valence quark masses, are reported in Table 1. The theoretical error, 0.06 GeV, is extracted according to the procedure discussed in Sect. 3.3 (first procedure).

Analogously, for the \(\eta _b\) and \(B_c\) we have

$$\begin{aligned} M_{\eta _b} = 9.33\,(07)\,\mathrm{GeV} \end{aligned}$$
(17)

and

$$\begin{aligned} M_{B_c} = 6.18\,(08)\,\mathrm{GeV}, \end{aligned}$$
(18)

respectively.

It is worth noting that the theoretical errors on the meson ground-state masses of Eqs. (1618), \({{\mathcal {O}}}\) (100 MeV), are of the same order of magnitude as those on the tetraquark energies of Eqs. (511). Therefore, we can state that the theoretical error on our results for both tetraquarks (Sects. 3.13.4) and ordinary mesons (Sect. 3.5) is \({{\mathcal {O}}}\) (100 MeV).

Table 5 \(\Omega _{ccc}\) and \(\Omega _{bbb}\) baryon spectra (up to the 2nd radial excitation), obtained by solving the eigenvalue problems of the Relativized QM (RQM) [13, Eqs. (2)] and relativized Diquark Model (RDM) [Eq. (1) herein] Hamiltonians. The values of the model parameters are extracted from [13, Table II] and Table 1, respectively. In the diquark model, the quark–diquark states are labelled as \(N(S_q,S_D)S, L; J^P\). Here, N is the radial quantum number; \(S_q\), \(S_D\) are the spin of the quark and diquark, respectively, coupled to the total spin of the baryon, S; the latter is coupled to the orbital angular momentum, L, to get the total angular momentum of the baryon, \(J^P\), with parity P. For \(J^P = \frac{1}{2}^+\) states, we have \(N(\frac{1}{2},1)\frac{3}{2}, 2; \frac{1}{2}^+\); for \(J^P = \frac{3}{2}^+\) states, we have \(N(\frac{1}{2},1)\frac{3}{2}, 0; \frac{3}{2}^+\), where \(N = 1, 2, 3\). Our results are compared with those of lattice QCD calculations, Ref. [87] for the \(\Omega _{ccc}\) sector and Ref. [88] for the \(\Omega _{bbb}\) sector

3.6 Fully-heavy baryon masses

Here, we extend the diquark model results of the previous sections and Refs. [33, 35] to the study of fully-heavy baryons.

We calculate the ground-state masses of the \(\Omega _{ccb}\), \(\Omega _{bbc}\), \(\Omega _{ccc}\) and \(\Omega _{bbb}\) fully-heavy baryons and also those of their radial excitations in the quark–diquark picture. We do that by making use of both the Hamiltonian of Eq. (1), with the model parameters reported in Table 1, and that of the relativized QM of Ref. [13] (where the diquark mass values are extracted from Table 1 herein); analogously to the \(Q {\bar{Q}}\) meson case of Sect. 3.5, we also need to substitute a valence antidiquark with a quark. \(\Omega _{QQQ}\) baryons are made up of a heavy quark, \(Q = c\) or b, and an axial-vector diquark, \(\{Q,Q\}\). Scalar-diquark configurations, [QQ], are only permitted in the \(\Omega _{ccb}\) and \(\Omega _{bbc}\) cases.

Starting from the \(\Omega _{ccc}\) and \(\Omega _{bbb}\) ground-states, characterized by \(J^P = \frac{3}{2}^+\) quantum numbers, we have

$$\begin{aligned} M_{ccc}^\mathrm{gs} = 4.66\,(08)\,\mathrm{GeV} \end{aligned}$$
(19)

and

$$\begin{aligned} M_{bbb}^\mathrm{gs} = 14.16\,(09)\,\mathrm{GeV}; \end{aligned}$$
(20)

these states are described as S-wave quark-axial-vector diquark configurations. The errors on these results are estimated by using the procedure outlined in Sect. 3.3 (first procedure). As an example, we discuss the calculation of the \(\Omega _{ccc}\) mass with its theoretical uncertainty. If we make use of the Relativized Diquark Model (RDM) [33, 35], we get \(M_{ccc}^\mathrm{gs} = 4655\) MeV; on the contrary, if we calculate the mass of the \(\Omega _{ccc}\) ground-state in the Relativized QM (RQM) [13], we obtain \(M_{ccc}^\mathrm{gs} = 4732\) MeV. Our result of Eq. (19) is thus given by the RDM prediction, 4.66 GeV, with the theoretical error being the difference between the previous RDM and the RQM results. Our predictions for the \(\frac{1}{2}^+\) and \(\frac{3}{2}^+\) \(\Omega _{ccc}\) and \(\Omega _{bbb}\) and their radial excitations are reported in Table 5. It is worth to note that our RQM and RDM results differ by \({\mathcal {O}}\) (100 MeV); however, the splittings among radial excitations in the two approaches are very similar. For example, those between the first and second \(\frac{3}{2}^+\) radial excitations are 304 MeV (RDM) and 334 MeV (RQM).

In the case of \(\Omega _{ccb}\) and \(\Omega _{bbc}\) ground-states, we get

$$\begin{aligned} M_{ccb}^\mathrm{gs} = 7.72\,(07)\,\mathrm{GeV} \end{aligned}$$
(21)

and

$$\begin{aligned} M_{bbc}^\mathrm{gs} = 11.01\,(12)\,\mathrm{GeV}; \end{aligned}$$
(22)

these states are described as S-wave quark-axial-vector diquark configurations, \(\{c,c\}b\) and \(\{b,b\}c\), with \(J^P = \frac{1}{2}^+\) quantum numbers. For their radial excitations, we get

$$\begin{aligned} M_{ccb}^\mathrm{re} = 8.26\,(13)\,\mathrm{GeV} \end{aligned}$$
(23)

and

$$\begin{aligned} M_{bbc}^\mathrm{re} = 11.53\,(16)\,\mathrm{GeV}. \end{aligned}$$
(24)

Our results can be compared to those of previous studies [87,88,89,90,91,92,93,94]. The recent experimental observation of doubly-charmed baryons [95, 96] will pave the way, hopefully within a few years, for that of fully-heavy three-quark systems. Once experimental results for the fully-heavy baryons are released, the quality of the match between our predictions and the data will give indications on the validity of the diquark model approach and the possible existence of a baryon–meson dynamical supersymmetry.

4 Summary and perspective

Adopting a perspective in which tetraquarks are viewed as states made up of elementary colour-antitriplet diquarks and colour-triplet antidiquarks and baryons as two-body quark–diquark configurations and using a well-constrained model Hamiltonian, built with relativistic kinetic energies, a one-gluon exchange potential and linear confinement (Sect. 2), we computed the masses of \(b {\bar{b}} q {\bar{q}}\), \(bb {\bar{q}} {\bar{q}}\), \(cc{\bar{c}} {\bar{c}}\), \(b b {\bar{b}} {\bar{b}}\), \(b c{\bar{b}} {\bar{c}}\), \(b b {\bar{c}} {\bar{c}}\) states and also fully-heavy QQQ baryons (\(Q = c\) or b) (Sect. 3). The eigenvalue problems were solved using a numerical variational procedure in concert with harmonic-oscillator trial wave functions. In each channel, we compared our prediction for the mass of the \(J^{PC}=0^{++}\) tetraquark ground-state (including its theoretical error) with those of previous studies. Our results for the masses of \(1^{++}\), \(2^{++}\), \(0^{++}\) four-quark excited states, and so on, may be of help to the experimentalists in their search for a tetraquark signal in those channels which do not involve the decay of the \(0^{++}\) ground-state tetraquark. Our predictions for QQQ states will be of use to the experimentalists as the recent experimental observation of doubly-charmed baryons [95, 96] will pave the way, hopefully within a few years, for that of fully-heavy three-quark systems. Once experimental results for the fully-heavy baryons are released, the quality of the match between our predictions and the data will give indications on the validity of the diquark model approach and the possible emergence of a baryon–meson dynamical supersymmetry.

While there is a large number of studies on the spectroscopy of tetraquark states, there are few investigations regarding their production [97,98,99,100,101,102,103] and principal decay modes [33, 69, 85, 86, 104,105,106,107,108,109]. Some progresses have been made. In this regard, one can notice that there is a big difference between resonances that can spontaneously dissociate into two mesons, and bound states that decay weakly (e.g., \(bb {\bar{c}} {\bar{c}}\)) [107] or by internal annihilation (e.g., \(bb {\bar{b}} {\bar{b}}\)) [86]. We also observe that while the \(\Upsilon \) is so narrow, \({\mathcal {O}}\) (50 keV), because of the Zweig rule, the fully-b ground-state tetraquark case is different; see Ref. [86]. Indeed, in this second instance one can make use of a Fierz transformation to re-write the ground-state tetraquark wave function as a linear combination of \((b {\bar{b}})_1 (b {\bar{b}})_1\) and \((b {\bar{b}})_8 (b {\bar{b}})_8\) components [86, Eq. (7)], where the subscripts 1 and 8 indicate color-singlet and octet components, respectively. Then, a color-octet quark–antiquark component, \((b {\bar{b}})_8\), can couple to a gluon (or, in other words, annihilate into one), giving rise to the possibility of hadronic decays characterized by a relatively large, \({\mathcal {O}}\) (10 MeV), decay width [86].