Mass spectra and wave functions of $T_{QQ\bar Q\bar Q}$ tetraquarks

The compact tetraquark states with fully heavy quark contents $QQ\bar Q\bar Q$ are studied as the bound states of the diquark-antidiquark within the Bethe-Salpeter framework. The (anti)diquark masses and form factors used are the same as we calculated the doubly heavy baryons in a previous work. Under the instantaneous approximation, the three-dimensional (Bethe-)Salpeter equation of the tetraquarks is derived and solved numerically to obtain the corresponding mass spectra and wave functions of the tetraquarks with $J^{PC}=0^{++}$, $1^{+-}$, and $2^{++}$. Our results show that the three ground states of $cc\bar c\bar c$ locate in the mass range of $6.4\sim6.5$GeV, and the $bb\bar b\bar b$ states in mass range of $19.2\sim19.3$GeV. The obtained relativistic wave functions naturally include the mixing effects from the possible $D$\,(or $G$) partial waves, and then can be further used to do precise calculations of the tetraquark decays. Based on the obtained results, the LHCb's observation $X(6900)$ is less likely to be the ground states of compact $cc\bar c\bar c$ tetraquarks but might be the first or second excited states. In addition, a widely used propagator-like form factor is also investigated and discussed.


Introduction
The quantum chromodynamics and the quark model allows not only the well-known traditional hadrons, such as the qq-type mesons and qqq-type baryons, but also the exotic tetraquark states and the pentaquark baryons [1,2]. About 50 years after the predictions of these exotic states, the LHCb Collaboration first detected the pentaquark baryons [3,4], and then in 2020 reported a new narrow structure labeled as X(6900) which covers the predicted masses of states composed of four charm quarks [5].
Although not the first hint of the tetraquark states, X(6900) causes great attention in hadron physics for its fully heavy quark contents. Inspired by this observation, the fourcharm states around 6.9 GeV have been investigated in several models or approaches [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20], such as, solving the two-or many-body time-independent Schrödinger equation [6][7][8][9][10], the chromomagnetic interaction models [11,12], the QCD sum rules combined with the Regge trajectories [13], the extended Godfrey and Isgur (GI) quark model [14], the relativistic quark model based on the quasipotential approach or effective Hamiltonian [15,16], and also the Bethe-Salpeter framework with different interaction kernels and approximation methods [17,18]. Also notice by using the heavy diquark limit M i → ∞, Ref. [18] is in fact dealing with a Schrödinger equation combined with the Regge trajectories to obtain the mass spectra. Notice most of these previous studies are based on the nonrelativistic Schrödinger equation or effective Hamiltonian methods while the relativistic effects and the possible S-D or D-G mixing effects are not included properly. In these previous studies, the ground states of the cccc structures are predicted to be 6.1 ∼ 6.5 GeV, and the X(6900) are usually tentatively identified as the radial or orbital excited states of tetraquark cccc in these recent studies, or interpreted as the coupled-channel or rescattering effects of two (or more) charmonia [21][22][23]. Other interpretations, such as the gluonic tetracharm [24] or light Higgs-like boson [25] are also proposed. So far, there is still no strong evidences whether these observed exotic hadrons are the genuine multiquark bound states or just the loosely bound molecules of the traditional mesons and baryons.
The fully heavy tetraquark states have several advantages in both experimental and theoretical researches. On the one hand, the mass of fully heavy tetraquarks locates far away from those of the traditional mesons and doubly heavy tetraquark QQqq and they can be clearly identified from the the known hadron spectra. On the other hand, the molecular states of two charmonia can not be bounded by the light boson exchange which makes the molecule configuration much more difficult to produce such states. Hence, the exotic states consisting of four heavy quarks are more likely to be the real compact tetraquarks.
In this work, we will try to deal with the mass spectra and wave functions of the fully heavy tetraquark states T QQQQ , namely, cccc or bbbb, within the framework of the Bethe-Salpeter equation. The fully heavy tetraquarks are assumed to be formed by the diquark QQ and the antidiquarkQQ. The diquark QQ is further assumed to be in the color3 configuration in order to produce the attractive force, and then similarly, the antidiquark QQ is in the color 3 state. This compact color-3 diquark picture has already been used in a previous work to study the doubly heavy baryons [26], where the obtained mass of Ξ ++ cc is just 20 MeV lower than the experimental measurements, and other predictions are also consistent with the recent theoretical researches especially the Lattice QCD results [27]. We will use the previously calculated mass spectra and form factors of the J P = 1 + cc and bb diquarks in this work. In addition, a propagator-like form factor is investigated and the corresponding cutoff dependence is also studied. The color3 QQ diquark and the color 3 QQ antidiquark can finally form a compact tetraquark in color singlet by the one-gluonexchange interaction. Based on the above analysis, the four-body tetraquark problem can be first reduced into two two-body bound problems of fermion (antifermon) system, which can be solved by calculating the original Bethe-Salpeter equation (BSE) [28,29]. Then we need to deal with a two-body problem of boson system which will be the main focus of this work.
The Bethe-Salpeter framework has great advantages in dealing with the two-body bound states for the relativistic interaction kernel and corresponding wave functions. The constructed relativistic wave functions are based on the good quantum number J P (C) rather than the nonrelativistic characteristics spin and orbital angular momentum 2S+1 L J . The BSE framework has been successfully used in the mass spectra of mesons [30][31][32], traditional baryons and pentaquarks [26,33], hadronic transitions, electro-weak decays, and etc [34][35][36][37][38][39][40]. In this work we would try to push the BSE framework further to study the fully heavy QQQQ system and develop a precise and systematic approach to describe the compact tetraquark states.
This manuscript is organized as: in section 2 the (Bethe-)Salpeter equation is derived in the perspective of the tetraquarks taken as the axial-vector diquark-antidiquark bound states; in section 3 the Salpeter wave functions of the tetraquarks with J P C = 0 ++ , 1 +− , and 2 ++ are constructed; then in section 4 the obtained mass spectra and wave functions are given and discussed; finally a brief summary is presented.

Tetraquarks as the bound states of diquark and antidiquark
Considering the exclusive principle, the diquark (antidiquark), consisting of two charm quarks (antiquarks) in the color3 (3) configuration, could only be in the J P = 1 + spin configuration in orbital ground states, since the flavor wave function is naturally symmetric and spatial wave function is also symmetric in ground states. In current work, we do not consider the excitation of the diquarks or antidiquarks. Now we try to deal with the BSE of the tetraquark states consisting of a 1 + diquark and a 1 + antidiquark, which could form three ground states with spin-parity configuration J P C = 0 ++ , 1 +− , and 2 ++ respectively. The Greeks (red) are used for the Lorentz indices. P, p 1 (k 1 ), p 2 (k 2 ) denote the momenta of the tetraquark state, constituent diquark, and the constituent antidiquark respectively; K αα ;β β denotes the effective interaction kernel between the diquark and the antidiquark based on the one-gluon exchange.

Bethe-Salpeter equation of two vector bosons
The Bethe-Salpeter equations of the bound states consisting of two vector (or axialvactor) constituents are schematically depicted in Fig. 1. The corresponding BSE are expressed as the four-dimensional integral of the inner relative momentum k, where Γ αβ (P, q, ξ) denotes the vertex of the two axialvector constituents; symbol P is used to denote the momentum of the tetraquark state, and we have P 2 = M 2 with M representing the tetraquark mass; symbol ξ indicates the polarization state, with ξ = ±2, 0 for the J P = 2 + tetraquark states, ξ = ±1, 0 for the 1 + ones, and ξ = 0 for the scalar ones; K αα ;β β (P, k, q) is the interaction kernel of the diquark pair based on the onegluon exchange. The constituent mass of the diquark (antidiquark) is represented by M 1 (2) . The constituent masses of the J P = 1 + cc and bb diquarks are obtained by solving the corresponding BSE [26]. The inner relative momenta q and k are defined as q = α 2 p 1 −α 1 p 2 and k = α 2 k 1 − α 1 k 2 respectively, with α i ≡ M i M 1 +M 2 . The effective propagator of the 1 + diquark reads, 1 +i is the usual scalar propagator. As usual, the Bethe-Salpeter wave function describing the tetraquark states can be defined as T αβ (P, q, ξ) = D αµ (p 2 )Γ µν (P, q, ξ)D βν (p 1 ). (2. 3) The symbols P and ξ in the BS wave function T αβ (P, q, ξ) and vertex Γ µν (P, q, ξ) will be omitted unless it is necessary to specify them. By using the definition of BS wave function, the BSE in Eq. (2.1) can also be rewritten as the integral of T αβ , where the inverse of the vector propagator reads, It can be easily checked that ϑ αβ (p i ) fulfills the condition 2 denotes the kinetic energy of the ith constituent.

Diquark form factors and tetraquark interaction kernel
For the tetraquark states consisting of the 1 + diquark and antidiquark constituents, each vector constituent has internal structure and usually can not be regarded as the pointlike particle. The diquark (antidiquark) is usually described by the corresponding form factors. Generally speaking, the form factors have great effects on the energy splittings of the tetraquark bound states. Then the potential V (s) between the two constituents will be smeared by the two form factors. With these two form factors, the interaction kernels of tetraquarks can be expressed as, where f 1V and f 2V denote the vector form factor of the diquark and antidiquark, and in the fully heavy cccc or bbbb system, we have f 1V = f 2V = f V ; V (s) denotes the onegluon exchange potential with s ≡ (k − q) being the momentum of the exchanged gluon. Compared with the case in baryon problem [26], there are two form factors to describe the nonpointlike structures, which make it much more complicated to deal with the tetraquark system.
Under the instantaneous approximation, the one-gluon-exchange potential is assumed to be static, namely, V (s) ∼ V (s ⊥ ) with s ⊥ = s − P ·s M , which reads in the Coulomb gauge as [41][42][43][44] where 4 3 is the color factor in the color singlet; a 1(2) is introduced to avoid the divergence in small momentum transfer zone; the potential V Conf ( s ) describing the confinement effects is introduced phenomenologically, which is characterized by the the string constant λ and the factor a 2 . The potential used here is based on the famous Cornell potential [45,46], which behaves as the one-gluon exchange Coulomb-type potential at short distance and a linear growth confinement one at long distance, and then modified as the aforementioned one to incorporate the color screening effects [47,48] in the linear confinement potential. V 0 is a constant fixed by fitting to the meson data. The strong coupling constant α s has the following form, where Λ QCD is the scale of the strong interaction; N f is the active flavor number which is 3 for the cc interaction while 4 for the bb interaction; a = e is a regulator constant. For later convenience, we further split V ( s ) into two parts as Also since the spatial parts are suppressed by a factor of ( v c ) 2 , we only consider the dominate time component (ν = 0) in the one-gluon-exchange Lorentz structure, namely, where we used the abbreviations i ≡ 2(α i M + q P ) with definition q P ≡ P ·q M . Using the notations and approximations introduced above, the tetraquark interaction kernel behaves as where κ ≡ f 2 V V ( s ); the factor 1 2 has been stripped off for later convenience. Notice the kernel κ(s ⊥ ) here has no dependence on the time component of momentum transfer s.
The diquark form factors f V (s 2 ) can be calculated by using the BS wave functions of the diquark. The attractive diquark is in the color antitriplet with the corresponding color factor (− 2 3 ), which makes the bound interaction be half of the corresponding meson. Namely, in the rainbow-ladder truncation the effective interaction for the diquarks is reduced by a factor of 2 compared to that in the meson channel. Then after a charge conjugate transformation, a diquark fulfills the same Bethe-Salpeter equation with a meson system with only the interaction kernel halved and the parity flipped [49][50][51]. Then by solving the BSE of the J P = 1 − QQ meson system, we can obtain the corresponding mass spectra and wave functions for the doubly heavy diquarks cc and bb (see Ref. [26] for detailed calculations). The corresponding form factors f V describing the interaction between the diquark and a gluon have also been obtained in the previous work [26], which are showed graphically in Fig. 2(a) labeled as method I. On the other hand, a colored diquark is never observed in nature alone and the confinement in diquark is an open problem in quark model. Therefore, we also calculated the corresponding form factors without considering the confinement item V Conf (s) in the Cornell potential of the diquark interaction kernel, and the obtained results are showed in Fig. 2(b) labeled as method II. Both form factors will be used in the calculation of the tetraquarks. The form factors fall faster with s 2 in method II than in method I. Also notice the form factors used here are calculated from the relevant BS wave functions combined with the Mandelstam formalism, which makes this work a self-contained framework and is different from the parameterized treatment in Ref. [17]. Also notice the form factors in method I and II are calculated based on the on-shell diquark bound state, and then are simply generalized to the off-shell diquark propagator to describe the structure's non-pointlike effects. In order to investigate the influences of different form factors and make comparisons, we will also use another phenomenological form factors in this work, namely, where Λ c is a introduced regulator parameter. We notice that this propagator-like form factor is widely used in literature. Generally speaking, the parameter Λ c should be determined by fitting to the data. However, we found that this propagator-like function can not describe the obtained form factors well in a large range of the momentum transfer. Hence, in this work, we only fit in the small zone of s 2 , namely, about 1(3) ∼ 0 GeV 2 for cc (bb) diquark, and the determined Λ c is 1.7 and 2.8 GeV for the cc and bb diquarks respectively. Also we will let Λ c change from 0.1 GeV to 20 GeV to investigate the corresponding dependence.

BSE under the instantaneous approximation
Generally speaking, solving the four dimensional BSE 2.1 is not an easy computational problem. However, under the instantaneous approximation introduced above, following Salpeter's procedures to cope with the fermion-antifermion system [29], we can reduce the four-dimensional BSE of the tetraquarks into the three-dimensional integral equation.
Inserting the instantaneous kernel 2.10 into the BSE 2.1, the tetraquark vertex can now be further expressed as Γ αβ (q) = 1 2 Θ αβ (q ⊥ ), where the three-dimensional vertex Θ αβ is defined as where the Salpeter wave function is defined by absorbing the integration over the time component as usual, The Salpeter wave function ϕ αβ (k ⊥ ) is only explicitly dependent on the three-dimensional momentum k ⊥ . Following the standard procedure, performing the contour integral over q P on both sides of Eq. (2.3), we obtain the three-dimensional Salpeter equation (SE), (2.14) Also we can define the positive and negative energy wave functions as and we have ϕ αβ = ϕ + αβ + ϕ − αβ . In the weak binding condition M ∼ (w 1 + w 2 ), ϕ + αβ ϕ − αβ , and the positive energy wave function ϕ + αβ (q ⊥ ) dominates. The SE can be further rewritten as the following simple Shrödinger type (2.16) The obtained three-dimensional BSE, namely, Eq. (2.16), indicates that the mass of the tetraquark state consists of two parts, the kinetic energy and the potential energy. Also we notice Eq. (2.16) is in fact the integral equation of the Salpeter wave function ϕ αβ (q ⊥ ), and M 2 behaves as the eigenvalue of the corresponding Salpeter wave function. By solving this eigenvalue equation, we can obtain the mass spectra and wave functions of the corresponding tetraquarks.
The normalization condition of the Bethe-Salpeter wave function for the tetraquark state is generally expressed as, where the conjugate wave function is defined asT αβ ≡ γ 0 T † αβ γ 0 ; the integral kernel reads, Both the inverses of the propagators and the interaction kernel K are dependent on P 0 and q P . Inserting the inverses of the propagators, the integration involving the propagators behaves as, (2.17) While the partial differential of the kernel gives and then we obtain normalization related to the interaction kernel as, where we have used the Salpter equation (2.14); w q ≡ α 2 w 1 + α 1 w 2 . On the other hand, the vertex Θ αβ can also be expressed by the wave function as, Putting the two part together and inserting the equation above, we obtain the normalization of the Salpeter wave function, where the abbreviation N 0 ≡ 2w 1 w 2 (w 1 +w 2 ) + wqM 2 (w 1 +w 2 ) 2 − w q is used. By the dimensional analysis, we can conclude that the dimension of the Salpeter wave function ϕ αβ is (−2) in units of mass.
It should be noted that, the obtained three-dimensional (Bethe-)Salpeter equation 2.16 and the corresponding normalization condition 2.21 are universal to any fully heavy tetraquark states consisting of two (axial)vector constituents, and do not depend on the specific properties of the total angular momentum J or parity. These obtained equations are applicable to the Salpeter wave functions with J P C = 0 ++ , 1 +− , 2 ++ or any other possible spin-parity. This is different from the approaches adopted in Ref. [17], where the Salpeter equations are coupled with three certain wave functions.The obtained equations and relevant results can also be further extended to deal with the molecular states consisting of two vector mesons.

BS wave functions of the tetraquark states
A 1 + diquark and a 1 + antidiquark can form a boson with J P = 0 + , 1 + , or 2 + in the ground state, namely, 1 × 1 = 0 + 1 + 2. For the tetraquark state, the spatial parity is expressed as P = (−1) l , with l stands for the quantum number of the angular momentum. A tetraquark state containing two c quarks and twoc quarks also occupies the definite Cparity. The charge conjugate parity (C-parity) is expressed as C = (−1) l+s with s stands for the spin of the particle. In order to simplify the expressions, from now on, we will always use the abbreviation x α ≡ q ⊥α | q | . Since throughout this manuscript we work in the momentum space, this abbreviation is supposed not to cause confusion.
According to Lorentz condition, spin-parity and also considering the relativistic covariance, the Salpeter wave function of tetraquark states with J P = 0 + can be generally expressed as, whereP α = Pα M ; the radial wave function g i (| q |) (i = 1, 2) just depends on | q | explicitly. It is clear to see that g 1 corresponds to the S-wave component, and g 2 contributes to both the S and D partial waves (see Ref. [26] for a detailed expression in terms of the spherical harmonics Y m l ). By inserting the two wave functions into Eq. (2.21), the normalization condition of the 0 ++ Salpeter wave function is obtained, Similarly, the J P C = 1 +− Salpeter wave function can be constructed by using the the antisymmetric Levi-Civita tensor αβµν as where αβeP = αβµν e µP ν ; e α (ξ) is the polarization vector with ξ (ξ = 0, ± 1) denoting the possible polarization states, and fulfills the following conditions It is clear to see that h 1 and h 2 parts in Eq. (3.3) represent the S and D-wave components, respectively. . The J P C = 1 +− tetraquark wave function is antisymmetric under the interchange of the two free Lorentz index, which is different from the 0 ++ case. The 1 +− normalization is finally expressed as, Notice in a similar work [17] only part of the S-wave components in the wave functions of J P C = 0 ++ and 1 +− are included, while our results show that the D-wave components, namely, the g 2 and h 2 items, and the the possible S-D mixing effects also play important roles especially in the excited states (see the obtained wave functions in Fig. 5).
Notice in the construction of the Salpeter wave functions, we work based on the good quantum number J P C while not the nonrelativistic characteristics, such as spin S or orbital angular momentum L. Then the contributions from the D or G-wave components are naturally included in the 0 ++ , 1 +− , and 2 ++ Salpeter wave functions, and are determined by the dynamics of the Bethe-Salpeter equation while not the man-made mixing effects. This is one of the advantages of the relativistic Bethe-Salpeter methods. The relativistic effects are naturally included in both the adoption of the BSE and the construction of the wave functions though the instantaneous approximation partly destroyed the covariance. Inserting these constructed Salpeter wave functions into the three-dimensional Salpeter equation 2.16 and solving this eigenvalue problem numerically, we can obtain the corresponding mass spectra and wave functions, which are presented in following section.

Numerical results and discussions
Before giving the numerical results of the mass spectra and wave functions, we specify the numerical values of the parameters first. The model parameters used in this work are kept the same with that we applied in previous meson and baryon calculations [26,31,[36][37][38][39][40] The free parameter V 0 plays a role in shifting the mass spectra and is determined by the spin-weighted average methods. By using the Clebsch-Gordan coefficients, the J P = 0 + tetraquark formed by the two 1 + diquark and antidiquarks can be decomposed as, where label 1, 2 denote the quarks, and 3, 4 denote the antiquarks; |(12) 1 (34) 1 means quark-1 and quark-2 are in the spin-1 state, while the antiquark-3 and antiquark-4 are also in the spin-1 state; then other notations are also implied. Then the parameter V 0 for 0 + cccc will be expressed as V 0(cccc) = 1 4 V 0(J/ψ) + 3 4 V 0(ηc) . By similar analysis, the J P = 1 + tetraquark is decomposed as, and then the corresponding parameter is determined as V 0(cccc) = 1 2 V 0(J/ψ) + V 0(ηc) . The V 0 for 2 + cccc is totally decided by the V 0(J/ψ) since any two quarks (anitquarks) inside are in the spin-1 state. Finally, the obtained parameters V 0 for 0 + , 1 + , and 2 + (cccc) are −0.304, −0.276, and −0.221 GeV respectively; for bbbb are −0.207, −0.192, and −0.162 GeV, respectively.
Mass (GeV) 6  The obtained mass spectra of the tetraquarks cccc and bbbb are showed in Fig. 3, where method I denotes the results with considering the confinement item V Conf in calculating the diquark masses and the corresponding form factors, while method II not. The results labeled as method III represents the ones where the diquark form factors are assumed to be propagator-like type in Eq. (2.11) with Λ c = 1.7 GeV for cc diquark and 2.8 GeV for bb diquark respectively. In order to see the dependence on the regulator Λ c in method III, we show the variation of the mass spectra along with Λ c in Fig. 4 by changing the values of Λ c from 0.1 GeV to 20 GeV. With the increases of Λ c , f V becomes more and more flat, which means the diquark is more and more similar to a pointlike particle. Also it should be pointed out that the regulator parameter Λ c can not be taken too large which may cause the instability of the instantaneous BSE in high momentum zone. The plots reveal that the mass spectra in ground states are sensitive to the diquark form factors. The obtained mass spectra show that the ground cccc tetraquarks locate in the range 6.4 ∼ 6.5 GeV when considering V Conf in diquarks. The mass splittings for the 1 +− and 2 ++ to 0 ++ states are about 35 MeV and 100 MeV respectively. When the confinement item is not included in the diquarks, the corresponding masses then locate about 40 MeV lower. The first excited states are always about 400 MeV higher than their corresponding ground states. A comparison of   Tab. I  are 6.36, 6.43, and 6.48 GeV for the 0 ++ , 1 +− , and 2 ++ cccc in ground states respectively; and the corresponding standard derivations are 0.13, 0.10 and 0.09 GeV respectively. All the results listed in Tab. I are above the threshold of the lowest quarkonium pair η c η c . Hence these three ground states are expected to be broad, since all of them can decay to a pair of quarkonia η c η c or J/ψJ/ψ through the (anti)quark rearrangements. These kind of decays are favored both dynamically and kinematically. From Tab. I we can conclude that the obtained masses of the ground cccc states are usually much lower (about 400 ∼ 500 MeV) than the X(6900) observed by the LHCb collaboration. The observed X(6900) is less likely to be the ground state of the compact tetraquark cccc states, but might be the first or second radial excited states. However, more detailed information is needed to investigate the inner structure of X(6900). We also notice that the masses of the cccc in ground states is near to the X(6900) with Λ c ∼ 1 GeV in method III.
The obtained bbbb masses in ground states are in the range 19.2 ∼ 19.3 GeV, which are higher than the ΥΥ threshold 1.892 GeV and η b η b threshold 1.880 GeV but lower than the χ b0 χ b0 threshold [55]. Also we notice that M bb in method II is about 80 MeV lower than that in method I, while the obtained mass spectra of bbbb are about 20 MeV lower in ground states. The comparison of our predictions with other researches is collected in Tab. II, from which we can see that the theoretical masses of bbbb in ground states locate in a large range of 18.75 ∼ 19.35 GeV in researches. The mean values of the ground bbbb tetraquarks in Tab. II are about 19.14, 19.20, and 19.22 GeV for the 0 ++ , 1 +− , and 2 ++ , respectively, where the corresponding standard derivations are 0.20, 0.18 and 0.16 GeV respectively. Except the 3 results predicted in Ref. [16] and the 0 ++ state in Ref. [12], all other results listed in Tab. II are higher than the ΥΥ threshold. Therefore, the three states can decay to a pair of quarkonia η b η b or ΥΥ through the quark rearrangements, and hence are expected to be broad.
Tab. II: Comparison of the bbbb masses in ground states with J P C = 0 ++ , 1 +− , and 2 ++ in units of GeV.  5 shows the obtained Salpeter wave functions for the first four (cccc) states obtained by the method I. Notice that the D partial wave gives important contribution in the third and the forth states, and the possible D or G-wave mixing are included naturally. The obtained wave functions show the rich information of the inner structure, and can be further used to do precise calculations on the decays, magnetic moments, or other properties of the tetraquarks.

Summary
In this work, we study the compact tetraquark states with fully heavy quark contents QQQQ based on the Bethe-Salpeter equation. The compact tetraquark is taken as the bound state of the diquark-antidiquark where the (anti)diquark form factors are calculated with and without considering the confinement item respectively. In addition, a propagator-like form factor is also used to calculate the tetrequark states and the corresponding parameter dependence is also investigated. Under the instantaneous approximation, the three-dimensional (Bethe-)Salpeter equation of the tetraquarks are derived, and the Salpeter wave functions with J P C = 0 ++ , 1 +− , and 2 ++ are then constructed and solved numerically to obtain the corresponding mass spectra of the tetraquarks. The compact tetraquarks as the bound states of the diquarks and antidiquarks follows a natural step of our previous work to deal with doubly heavy baryons. Therefore some results and approaches, such as the relevant interaction kernel, diquark constituent masses and form factors, the reduction of BSE, etc., can be directly adopted or extended to deal with compact tetraquarks, which allows us to take a systematic, relativistic, and unified framework to treat the heavy hadron systems, including the mesons, diquarks, baryons, and tetraquarks. Also notice all the parameters used to deal with the tetraquarks have already cccc tetraquarks but might be the first or second radially excited states. The obtained relativistic wave functions are based on the good quantum numbers and can naturally include the mixing effects from the possible D(or G)-wave components.