A non-relativistic model for the $[cc][\bar{c}\bar{c}]$ tetraquark

We use a non-relativistic model to study the spectroscopy of a tetraquark composed of $[cc][\bar{c}\bar{c}]$ in a diquark-antidiquark configuration. By numerically solving the Schr\"{o}dinger equation with a Cornell-inspired potential, we separate the four-body problem into three two-body problems. Spin-dependent terms (spin-spin, spin-orbit and tensor) are used to describe the splitting structure of the $c\bar{c}$ spectrum and are also extended to the interaction between diquarks. Recent experimental data on charmonium states are used to fix the parameters of the model and a satisfactory description of the spectrum is obtained. We find that the spin-dependent interaction is sizable in the diquark-antidiquark system, despite the heavy diquark mass, and also that the diquark has a finite size if treated in the same way as the $c\bar{c}$ systems. We find that the lowest $S$-wave $T_{4c}$ tetraquarks might be below their thresholds of spontaneous dissociation into low-lying charmonium pairs, while orbital and radial excitations would be mostly above the corresponding charmonium pair thresholds. Finally, we repeat the calculations without the confining part of the potential and obtain bound diquarks and bound tetraquarks. This might be relevant to the study of exotic charmonium in the quark-gluon plasma. The $T_{4c}$ states could be investigated in the forthcoming experiments at the LHC and Belle II.


I. INTRODUCTION
The existence of multiquark states with four or more quarks has been proposed decades ago [1,2]. The early papers on tetraquark configurations were based on the MIT bag model with light quarks only. Later on, the tetraquark picture was extended to heavy quarks [3,4]. The interest in this subject was renewed in the past decade due to the experimental observation of states which are not combinations of three quarks (qqq) or of quark and antiquark (qq). These new states present quantum numbers, masses, decay channels and widths that cannot be explained with the conventional meson or baryon models (therefore they are called Exotics) [5][6][7][8][9][10][11]. Some of them were even found to be charged, what establishes unambiguously their exotic nature [12,13].
In the present work we focus on tetraquarks composed by a single flavor, charm quarks only, using a diquarkantidiquark picture [cc][cc], which we will call T 4c or "the all-charm tetraquark".
The first work on the all-charm tetraquark was published in 1975 by Iwasaki [14]. In a subsequent paper Chao studied the T 4c in the diquark-antidiquark picture with orbital excitations and its production in e + e − annihilation [15] including an interesting analysis of the possible decay channels. Later, in the eighties and nineties, several works with different approaches addressed the question of the existence of this cccc state [16][17][18][19][20]. In more recent years, after the discovery of the X(3872), a new series of theoretical works on the subject appeared [21][22][23][24][25][26][27][28][29][30][31].
On the experimental side, recent measurements of J/ψ pair production are very promising and might be the ideal starting point to search for the all-charm tetraquark. They have been studied in the LHC, by the LHCb [32,33], CMS [34] and ATLAS [35] collaborations. Double cc production has also been observed by the Belle collaboration [36]. Specially in Refs. [32,33] one can see that there is an enhancement in the differential production cross-section for J/ψ pairs between 6 and 8 GeV. Further investigating the invariant mass distribution in this energy range with high statistics would bring very useful information about the possible existence of the T 4c .
Most of the predictions for the T 4c mass lead to values around 6 GeV, therefore lie well above the charmonium experimentally known range (which is concentrated within 3 -4.5 GeV). This energy gap makes the all-charm tetraquark a special object in the sector of exotic multiquarks. The most discussed tetraquark candidates (the X, Y , Z states) are in the same mass range of conventional charmonium states and this can lead to confusion.
The absence of light quarks in the T 4c makes it unlikely to be a meson-meson molecule, since it is not easy to describe this binding in terms of pion exchange or light vector meson exchange. If it exists, the T 4c is bound by QCD forces and studying its spectrum will lead to a more complete understanding of QCD interactions. If it does not exist we have to understand why.
We will describe the T 4c as a two-body non-relativistic system, made of a cc diquark and acc antidiquark, which interact through a Cornell-like potential. We choose the diquark and antidiquark to be in the color antitriplet and triplet representations, respectively. The energy states are then obtained by solving the appropriate Schrödinger equation.
Why do we choose this model? We choose to work with the Cornell model because it is able to capture the essential aspects of the heavy quark-antiquark interactions. It has almost never been too wrong and when it was, there was something really new happening. Moreover, the quark-antiquark potential can be continuously improved [37] and its parameters can be adjusted so as to incorporate the most recent experimental information on the charmonium spectrum. Finally, we will study systems with angular momentum and all kinds of spin interactions. With more constituents, we may form systems with higher spin and total angular momentum. With the Cornell model (unlike in some other approaches) we can identify the individual contribution of each one of these interactions.
We choose to work with diquarks, not only because they simplify the calculations, but also because there is some evidence of diquark clustering in baryons. In the case of heavy diquarks the interaction has a stronger short distance component, in which the perturbative one gluon exchange may be attractive. In particular, the cc diquark became more interesting after the prediction of the T cc [38] and even more after the very recent discovery of the baryon Ξ ++ cc [39], a ccu state where the charm diquark may play a role.
In the literature we find some calculations which are very simple and strongly based on the existing empirical information, as in Ref. [30] and very sophisticated as the lattice calculations of Ref. [25] or the QCD sum rules calculations of [28]. Our model is at an intermediate level, being more precise than the estimates made in [30] and more transparent than the results found in [25,28], where it is very difficult to access the role of spin interactions. Ideally, all these approaches should converge and the origin of the remaining discrepancies should be well understood.
In the end of this work we will present a comparison with the results obtained in other approaches.
A pictorial representation of the all-charm tetraquark in the diquark-antiquark scheme of our model can be seen in Figure 1. One of the most common functional form of the zeroth-order potential, V (0) (r), employed in heavy quarkonium spectroscopy is the Coulomb plus Linear, where the Coulomb term arises from the one gluon exchange (OGE) associated with a Lorentz vector structure and the linear part responsible for confinement, which is usually associated with a Lorentz scalar structure. It is given by: where κ s , sometimes called "color factor", is related to the color configuration of the system (it can be negative or positive), α s is the QCD fine structure constant and b, sometimes called "string tension", is related to the strength of the confinement. One could also add a constant term, which would act as a zero-point energy. Usually, in heavy quark bound states the kinetic energy of the constituents is small if compared to their rest energy, hence a non-relativistic approach with static potentials can be a reasonable approximation. In two-body problems involving a central potential, it is convenient to work in the center-of-mass frame (CM), where one can use spherical coordinates to separate the radial and angular parts of the wavefunction, and the kinetic energy is written in terms of the reduced mass µ = (m 1 m 2 )/(m 1 + m 2 ). We start with the time independent Schrödinger equation: We first solve this radial equation to obtain the energy eigenstate and the wavefunction of each particular state. Next, spin-dependent terms are included as perturbative corrections. They account for the splitting between states with different quantum numbers. Based on the Breit-Fermi Hamiltonian for one gluon exchange [40][41][42][43], we introduce three spin-dependent terms: V SS (Spin-Spin), V LS (Spin-Orbit) and V T (Tensor). For equal masses m 1 = m 2 = m, they are given by: where the radial-dependent coefficients come from the vector V V and scalar V S parts of the potential in Eq. (1), where m is the constituent mass of the two-body problem (charm quark, or diquark). The second term in the spinorbit correction (proportional to the scalar contribution) is a Thomas precession, which follows from the assumption that the confining interaction comes from a Lorentz scalar structure. Notice that if we introduce a constant term V 0 in the potential it won't affect these radial coefficients, since only derivatives appear in them. In fact, adding a constant term only shifts the whole spectrum, forcing a change in the parameters such to reproduce the charmonium spectrum, without actual improvement in the quality of the fit. These spin-dependent terms are proportional to 1/m 2 , what justifies their treatment as first-order perturbation corrections in heavy quark bound states. The expectation value of their radial-dependent coefficients can be calculated using the wavefunction obtained with the solution of the Schrödinger equation.
This framework appears frequently in quarkonium spectroscopy, but a better agreement between predicted states and experimental data of cc mesons can be obtained by including the spin-spin interaction in the zeroth-order potential used in the Schrödinger equation (as done in Refs. [44][45][46][47]), with the artifact of replacing the Dirac delta by a gaussian function which introduces a new parameter σ. Then the spin-spin term becomes When the term S 1 · S 2 acts on the wavefunction it will generate a constant factor, so we still have a potential as function only of the r coordinate. The expectation value of the operator of the spin-spin interaction can be calculated in terms of the spin quantum numbers using the following relation, S 1 · S 2 = 1 2 S 2 − S 1 2 − S 2 2 , where S 1 and S 2 are the spins of particles 1 and 2 respectively, and S is the total spin in consideration.
The expectation value of the operator of the spin-orbit interaction can be calculated in terms of the quantum numbers of total angular momentum J (defined by the vector sum: J = L + S), total spin S and orbital angular momentum ℓ using the following relation, L · S = 1 2 J 2 − L 2 − S 2 . For S-wave states (ℓ = 0), the spin-orbit term L · S is always zero.
The tensor interaction demands a bit of algebra. For convenience, we redefine the tensor operator with an extra factor 12, which we remove from its radial coefficient in Eq. (8) The results for the diagonal matrix elements of the tensor operator between two spin 1/2 particles, like in the cc mesons, can be found in Refs. [41,48] and also (with more details) in Ref. [49]. The expectation value of the tensor is non-zero only for 1) ℓ = 0 and S = 1 (triplet), 2) J = ℓ, or J = ℓ − 1, or J = ℓ + 1.
After some manipulations of the spin operators, with the aid of some relations of spherical harmonics and the Pauli matrices with respective eigenvalues, we can obtain the following general result, which satisfies the above conditions (it always vanishes if ℓ = 0 or S = 0): for any of the allowed values of J and ℓ. For instance, for ℓ = 1 we have S 12 = − 2 5 , +2, −4, for J = 2, 1, 0, respectively. These results are valid for diagonal matrix elements. The tensor actually has non-vanishing non-diagonal matrix elements, but as a first-order perturbation correction they can be neglected. They would be important if the tensor operator were to be used as part of the potential, what would cause the mixing of the wavefunction itself, as in the deuteron [50,51].
Notice that in order to obtain these three general cases of non-vanishing diagonal matrix elements of the tensor operator for two spin 1/2 particles in Eq. (12), it is necessary to make use of a few relations that are valid only for Pauli matrices [48,49], like its eigenvalues and the anticommutation relation. Therefore we cannot use this result in the diquark-antidiquark tensor interaction (if we wish to treat it as a two-body problem), since the diquarks can have spin 0 or 1. This issue will be discussed later when we address the tetraquark interaction.
Regarding the wavefunction, we will consider only pure states where ℓ (orbital), S (total spin), and J (total angular momentum) are good quantum numbers. Then the wavefunction will be composed of a radial part, and an angular part which comes from the coupling of spherical harmonics and spin functions in a specific value of J.
Solving the eigenvalue equation (2), one can obtain the interaction energy E and the wavefunction y(r) of the two-body system under consideration, where both depend on the number of nodes of the wavefunction n (or principal quantum number N = n + 1), on the orbital angular momentum number ℓ, and in the case of the spin-spin correction included in V (0) , they will also depend on the total spin S and on the constituents spins S 1 and S 2 . Since the Schrödinger equation has no analytical solution for the potentials that are relevant here, we solve it numerically, using an improved version of the code published in Ref. [52].
An interesting quantity that can be used to check the validity of the non-relativistic approximation is the velocity of the constituents in each of the systems in consideration: the quark velocity inside the meson or the diquark velocity inside the tetraquark. As discussed in Ref. [53], the mean square velocity can be obtained from the kinetic energy, which can be calculated directly from the Hamiltonian, or using the Viral Theorem: where V (0) (r) is the effective zeroth-order potential placed in the Schrödinger equation. Both methods yield approximately the same result within the numerical precision employed. One interesting aspect of the non-relativistic approach is that, even though the charmonium system is not completely non-relativistic, a surprisingly good reproduction of its mass spectrum can be obtained. As discussed in Ref. [54], where a charmonium model is developed with completely relativistic energy and also with non-relativistic kinetic energy, a good agreement with the experimental data can be obtained with both methods, just by using a different set of parameters in the effective potential employed.
The value of the square modulus of the wavefunction at the origin, |Ψ(0)| 2 , is an important quantity. If the spinspin interaction was treated as a first-order perturbation without the gaussian smearing, it would be proportional to |Ψ(0)| 2 because of the Dirac delta. Decay widths can also be calculated using the wavefunction or its derivative at the origin. Only S-wave states (ℓ = 0) have non-zero value of the wavefunction at the origin. For states with orbital angular momentum (ℓ = 0), the centrifugal term in the Schrödinger equation creates a "centrifugal barrier", which makes the wavefunction at the origin vanish. Thus, for ℓ = 0 we will assume |Ψ(0)| 2 = 0 and for S-wave we have In fact, the important quantity is the square modulus of the radial wavefunction at the origin |R n,ℓ (0)| 2 , which can be obtained directly from the numerical calculations. In the literature of quarkonium models, we find the following formula (see Ref. [41] for a deduction) that relates the wavefunction at the origin |Ψ(0)| 2 to the radial potential V (0) (r): We have checked that the result obtained directly from the numerical method is compatible with the one using the formula above.
In more sophisticated quarkonium models, such as the relativized potential model of Ref. [44], the coupling constant α s is considered as a "running" parameter, that changes according to the energy scale of each bound state. However, we have chosen to adopt the non-relativistic model of Ref. [47], where α s is as a constant in the potential, which is also a common approach in many charmonium models.
The value of α s , the charm quark mass m c , the string tension b, and the gaussian parameter σ, will be obtained from a fit of the charmonium experimental data, and once the best set is found, they are kept fixed to generate the whole mass spectrum.

A. Charmonium
In order to make good estimates for diquarks and tetraquarks, we first study the spectrum of charmonium. In this case, the color factor κ s in Eq. (1) should be the one of a color singlet state, since for cc mesons we have |qq : 3 ⊗ 3 = 1 ⊕ 8 [55,56]. The result for the color singlet is κ s = −4/3 [49,55].
After having solved the Schrödinger equation, the mass of a particular state will be given by: The parity and charge conjugation quantum numbers of qq states are given by [55] P = (−1) ℓ+1 and C = (−1) ℓ+S repectively. Using the equation above we calculate the masses M calc i of the i states with well defined P and C, then we fit the experimentally measured masses M exp i and determine the parameters minimizing the χ 2 , defined as: Following Refs. [47,54] we choose w i = 1, which is equivalent to give the same statistical weight to all the states used as input. This way we ensure the resulting set of parameters will handle simultaneously the spin-spin splitting in S-wave, the spin-orbit and tensor splitting, which are specially important in P -wave, and the radial excitations as well.

B. Diquarks
In the study of tetraquarks, we shall treat the full four-body problem as three two-body problems. Repeating the steps described in the previous subsection, we first compute the mass spectrum of the diquark, then we do the same for the antidiquark and finally we solve once again the Schrödinger equation for a two body system composed by the diquark and antidiquark. The inspiration of this factorization is the color structure behind it.
A diquark is a cluster of two quarks which can form a bound state. This binding is caused by one gluon exchange between the quarks. In this interaction the factor κ s can be negative, then the short distance part (∝ 1/r) of the potential will be attractive. The SU (3) color symmetry of QCD implies that, when we combine two quarks in the fundamental (3) representation, we obtain: |qq : 3 ⊗ 3 = 3 ⊕ 6. Similarly, when we combine two antiquarks in the 3 representation, they can form an antidiquark in the 3 representation. Then the diquark and antidiquark can be combined according to |[qq] − [qq] : 3 ⊗ 3 = 1 ⊕ 8 and form a color singlet, for which the one gluon exchange potential is also attractive (see Refs. [56][57][58] ). The antitriplet state is attractive and yields a color factor κ s = −2/3, while the sextet is repulsive and yields a color factor κ s = +1/3 [49,55]. Therefore we will consider only diquarks in the antitriplet color state. We will use a diquark [cc] in the ground state, with no orbital nor radial excitations, such that we have the most compact diquark. We choose the attractive antitriplet color state, which is antisymmetric in the color wavefunction. Then, in order to respect the Pauli principle (the two quarks of same flavor are identical fermions), the diquark total spin S must be 1. In this way the total wavefunction of the diquark will be antisymmetric.
Notice that going from the color factor −4/3 (for quark-antiquark in the singlet color state) to the color factor −2/3 (for quark-quark in the antitriplet color state) is equivalent to introducing a factor 1/2, which one would expect to be a global factor since it comes from the color structure of the wavefunction. Because of that, it is very common to extend this factor 1/2 to the whole potential describing the quark-quark interaction. This rule is motivated by the interactions inside baryons, where two quarks can also be considered to form a color-antitriplet diquark, which then can interact with the third quark and form a color-singlet baryon. Since this seems to give satisfactory results in baryon spectroscopy, it has also been extended to diquarks inside tetraquarks. The general rule would be simply V qq = V qq /2. Many authors with different tetraquark models, for instance Refs. [59,60], also divide the confining part of the potential by 2 in order to adapt it to the diquark case. In our model, besides the change in the color factor, the string tension b, obtained from the fit of cc spectra, will be also divided by 2.
The calculation of the total mass of the diquark is completely analogous to the cc mesons, as in Eq. (16). The spin-dependent corrections are also analogous since we are still dealing with a two-body system composed of two spin 1/2 particles.

C. Tetraquarks
The all-charm tetraquark will be treated as a two-body (cc -cc) system with m cc = mcc. The color factor should correspond to the color singlet, therefore we will use κ s = −4/3 and also the same parameters α s , b and σ obtained from the fit of the cc spectrum. The calculation of its total mass will also be analogous to the charmonium: In order to properly calculate the spin-dependent corrections we need to remember that the diquarks have spin 1. Then, for the coupling of spin 1 diquark and spin 1 antidiquark, we will have the total tetraquark spin S T = 0, 1, 2. Besides that, we will also allow radial and/or orbital excitations in the diquark-antidiquark system. In our non-relativistic approach, we use ordinary quantum mechanics to couple the total spin S T to the orbital angular momentum L T into the total angular momentum J T .
For the spin-spin and spin-orbit corrections, we can obtain the angular factors from the spin, orbital and total angular momentum quantum numbers. However, for the tensor correction we only have a general result (in terms of eigenvalues) for the interaction between two spin 1/2 particles, shown in Eq. (12). Then, for a proper treatment of the tensor interaction in the diquark-antidiquark system we will explicitly apply the tensor operator on the angular part of the tetraquark wavefunction, as we will describe below.
Let us focus on the spatial and spin components of the wavefunction. We factorize the radial wavefunction from the angular one that combines orbital angular momentum and spin, which are coupled using Clebsh-Gordan coefficients. We will use the indices 1 and 2 for the two quarks inside the diquark, 3 and 4 for the two antiquarks inside the antidiquark (see Fig. 2). To illustrate our procedure to treat tensor interactions, we present one specific example with total spin S T = 2, L T = 1 and J T = 2. S d and Sd will denote the total spin of the diquark and antidiquark, respectively. We write the possible couplings in a generic form |S, M S , where S is a total spin and M S is its zcomponent. The arrows denote the spins of each constituent, in the order 1, 2 for the diquark and 3, 4 for the antidiquark. As usual the arrow up denotes spin up: | 1 2 , 1 2 and the arrow down denotes spin down: | 1 2 , − 1 2 . We show it in terms of diquark and antidiquark spin basis, and also in terms of the two quarks and two antiquarks spin basis (each group of four arrows is always in the order "1234"). These wavefunctions were inspired by the ones presented in Refs. [45,58,61,62], and we generalized them to include orbital angular momentum between diquark and antidiquark. For the choices mentioned above the wavefunction reads: For L T = 1 we have seven combinations to get J T if we are considering spin 1 diquark and antiquark: one for S T = 0 (J T = 1), three for S T = 1 (J T = 0, 1, 2) and three for S T = 2 (J T = 1, 2, 3). We now explicitly apply the tensor operator on the above angular wavefunction and we note that within our approximations, it is equivalent to apply this operator directly on the pair diquark-antidiquark (in spin 1 basis) or consider a sum of four tensor interactions between each pair of quark-antiquark (spin 1/2 basis) as illustrated in Fig.  2, as it would be expected from the angular momentum algebra 1 . We have: Since the tetraquark is treated as a two-body system, the expectation value of the radial wavefunction between every qq pair is the same and can be factorized. In a four-body problem (using Jacobi coordinates for example) where all the four constituents are allowed to move and interact at the same time with each other, this would not be true. This type of approach can be found in other models of tetraquarks, for instance in Ref. [61,62]. Usually in this kind of approach only the ground state is considered, with no orbital excitations, and hence only the spin-spin interaction is relevant since the spin-orbit and tensor vanish for ℓ = 0. Besides, in order to tackle the four-body problem one needs to resort to a variational approximation with gaussian trial wavefunctions or similar methods, therefore there will always be a compromise between the precision of the numerical solution and the reliability of the assumptions. In order to deal with the generalization of the tensor interaction to the tetraquark case, we will rewrite the tensor in a form that allows us to recover the same results that we already know for the particular case of two spin 1/2 particles and that can also be used as a generalization to other cases, such as the interaction between two spin 1 diquarks. The operator S 12 in Eq. (10) is a "rank-2" tensor which can be written in terms of spin operators and spherical harmonics, as shown in the textbook [63]. An extensive discussion on this approach can be found in Ref. [49].
The following functional form does not use any particular relation or eigenvalues for spin 1/2 particles, only general properties of angular momentum elementary theory. One can write the unity vectorr in spherical coordinates and the spin operators in cartesian components. Then they can be rearranged into raising, lowering and z-component spin operators and spherical harmonics of ℓ = 2, and we can write: where With the expression above we can take the expectation value of the tensor operator in the angular wavefunctions, as in Eq. (19), and use the selection rules of the spherical harmonics to find the non-vanishing terms.
To close this subsection, we discuss the tetraquark quantum numbers, as in Refs. [64,65]. We can use the diquarkantidiquark basis to label the possible quantum numbers J P C of the tetraquark. We shall use the following notation: where S d is the total spin of the diquark, Sd is the total spin of the antidiquark, S T is the total spin of the tetraquark, assumed to come from the coupling S d ⊗ Sd, L T is the orbital angular momentum relative to the diquark-antidiquark system (in the two-body approximation), and J T is the total angular momentum of the tetraquark, assumed to come from the coupling S T ⊗ L T . The general formula for charge-conjugation and parity of the tetraquark is: Since we are interested in the T 4c tetraquark, where the diquarks are composed by two charm quarks with spin 1 in the antitriplet color configuration, for the S-wave states we have the following possibilities: Note that all the S-wave tetraquarks described above have positive parity. The introduction of the first orbital excitation will bring a factor (−1) both in parity and in charge conjugation. Then all the P -wave states (with L T = 1) will have odd parity and the opposite charge conjugation in comparison with the S-wave states. In the Table  I we list the J P C quantum numbers of the 10 possibilities which we consider for the S-wave and P -wave all-charm tetraquarks built with spin 1 diquarks (also in accordance with Refs. [66] and [15]).

III. RESULTS
In this section we present the results of the calculations with the formalism outlined in the previous section. We use the following notation in our tables: the principal quantum number is N (N = 1 for the ground state, N = 2 for the first radial excitation and so on), ℓ is the orbital angular momentum, S is the total spin and J the total angular momentum. In spectroscopy notation the states are usually labeled by N 2S+1 ℓ J , with ℓ = 0, 1, 2, 3, . . . → S, P, D, F, . . . , for example 1 3 S 1 for J/ψ.

A. Charmonium
In order to make good estimates of diquark and tetraquark properties, we first study the spectrum of the conventional charmonium states to observe how well we can fit the experimental data. In our model we considered the zeroth-order potential of the form Coulomb plus linear plus smeared spin-spin interactions. We separate spin triplet (S=1) and spin singlet (S=0) before solving the Schrödinger equation. Using κ s = −4/3, S 1 = S 2 = 1/2 and S = 0 or S = 1, we replace the operator S 1 · S 2 by the constant [S(S + 1) − S 1 (S 1 + 1) − S 2 (S 2 + 1)]/2 and we find the wavefunction y(r) and the eigenvalue E. In Fig. 3 we show the zeroth-order potential for total spin 0 or 1. Later the spin-orbit and tensor corrections are included, splitting orbitally-excited states. We performed a fit with experimental values from the PDG [67]. The four parameters were allowed to vary in the following range: 1.1 < m c < 1.9 GeV, 0.1 < α s < 0.7, 0.050 < b < 0.450 ; GeV 2 , 0.7 < σ < 1.3 GeV. The results are also very similar to the ones from Refs. [46,47], which were obtained with the fit of 11 cc states with equal statistical weight. We have included two more, h c (1P ) and χ c2 (2P ), in a total of 13 states as input, obtaining the following values: m c = 1.4622 GeV, α s = 0.5202, b = 0.1463 GeV 2 , σ = 1.0831 GeV (26) Several fits with different number of input states and alternative models were tested in Ref. [49]. There is one particular alternative case worth mentioning. In this case, we considered the spin-spin interaction as first-order perturbation, proportional to the wavefunction at the origin, with the radial coefficient given by Eq. (6) (without the gaussian smearing), and also removed the Thomas precession term from the spin-orbit interaction, which is proportional to the string tension b on Eq. 7. In this way the spin-dependent corrections come exclusively from the Breit-Fermi Hamiltonian describing one gluon exchange, as in Ref. [43]. In this scheme it was possible to fit very accurately the 6 ground states 1S and 1P : , with the parameter set m c = 1.2819 GeV, α s = 0.3289 and b = 0.2150 GeV 2 . This set is appealing since the mass of the charm quark is exactly the PDG value [67] obtained in the M S scheme: 1.28 ± 0.03 GeV, and the coupling constant α s is also smaller, favoring the assumption of the perturbative regime of QCD. However, for radial excitations, specially above the DD threshold, this scheme does not work very well and hence we restrict ourselves to present the results obtained with the model that gives the best agreement with the whole experimental data set, since we believe this might yield better predictions to higher new charmonium states and also for the diquark and tetraquark. In Table II we present the wavefunction properties. Notice that the inclusion of the spin-spin interaction in the zeroth-order potential creates a small difference between the wavefunction of spin singlet and spin triplet. The spin 0 states receive a negative contribution from this interaction in the potential, what causes the short-distance region of the potential (small r coordinate) to be "more negative" generating states with smaller root mean square radius, higher value of the wavefunction at the origin and higher quark velocity.
The spin-dependent interactions are very important in charmonium spectroscopy because they can test the QCD dynamics in the heavy quark context, lying between the perturbative and the non-perturbative regime. Particularly interesting is the role of the spin-spin interaction in orbitally-excited states. It is convenient to define the spin-average mass of a multiplet (spin here means J), also know as "center-of-weight" or "center-of-gravity" (c.o.g.): For P -wave ground state, for example, we have:  Interestingly, in the spin average mass the spin-orbit and tensor corrections cancel each other and hence if the spin-spin correction is zero in the orbitally-excited singlet state (1 1 P 1 for instance), its mass should be equal to this spin average. However, the spin-spin correction is zero for orbitally-excited states only if it is treated as a first-order perturbation proportional to the wavefunction at the origin. In our model, where we include the gaussian term non-perturbatively, there will be a small difference.
In Table III we present the results for the masses including the spin interactions. Note that the contribution of the spin-spin interaction to orbitally-excited states is not zero, specially in P -wave, even tough the wave function at the origin is still compatible with zero. Because of that the spin singlet in orbitally-excited states is slightly different from the spin-average (c.o.g.). The experimental measurements of 1P states suggest that they should be very close (see Table IV for experimental values). As pointed in Ref. [68], a precise measurement of the difference between the c.o.g. of the 1 3 P J states and the singlet 1 1 P 1 can provide useful information about the spin-dependent interactions in heavy quarks. Actually, the prediction for h c (1 1 P 1 ) is already close to the experimental value and even more so if one considers its mass as the spin-average of the 3 P J states (as done in Ref. [47] for the calculations where its mass was required). Also, the inclusion of the recently measured χ c2 (2P ) [69,70] did not affect much the resulting set, even though the prediction for its mass is a little higher than the experimental value.  [47] the hc(1P ) is taken as the spin-average (c.o.g.) of the P -wave states, which is in better agreement with experimental data.
In the 2014 edition of the PDG [71] the X(3915) was assigned to the 2 3 P 0 cc state, the χ c0 (2P ), but due to many reasons [72] it has been removed from this position. The X(3915) has still the status of exotic resonance. A discussion about its nature (and also about the χ c2 (2P ) state) can be found in Ref. [73]. A recent example of the X(3915) interpreted as a diquark-antidiquark tetraquark [cs][cs] can be found in Ref. [74]. In Ref. [75] an analysis of Belle [69] and BaBar [70] data showed some evidence of the "real" χ c0 (2P ) indicating that its mass could be around 3837.6±11.5 MeV, which is in better agreement with quarkonium models. Recently, the Belle Collaboration found a candidate for the χ c0 (2P ) in the data of e + e − → J/ψDD [76] with a mass of 3862 +26+40 −32−13 MeV and a width of 201 +154+88 −67−82 MeV. Finally, in Table IV we compare the results of the model with the experimental data, which are illustrated in the mass spectrum presented in Fig. 4. We can see that the agreement with the experimental data is satisfactory.

B. Diquarks
We present now our calculations for heavy diquarks composed by two charm quarks cc (which are equivalent for antidiquarkscc in our framework). We use the model for charmonium discussed in the previous subsection, except that due to the different color structure, the color factor now is κ s = −2/3, which corresponds to the attractive antitriplet color state, and the string tension b will be half of the one obtained for the cc charmonium mesons. We will adopt the parameter set obtained by fitting this model to 13 cc states.
In Tables V and VI we present the results for the diquark wavefunctions and masses, respectively. For completeness we also show diquarks in 1P , 2S and 2P states. Because of the restrictions due to the Pauli exclusion principle the possibilities are much less numerous. Also, since the P -wave introduces a −1 factor in the parity, the antisymmetric restriction in the wavefunction implies that their total spin S should be 0 if they are in the antitriplet color state.   In Table VII we show a few results of other works about cc diquarks. Due to differences in the models and presentation in each reference, we show only the information that can be compared to our results. In particular, we select only the ones that correspond to the (attractive) antitriplet-color configuration. As it can be seen, the 1S diquark is very similar in all the models, with a mass around 3.1 GeV.

C. Tetraquarks
As discussed above, the diquark-antidiquark tetraquark is treated as a two-body system. The diquarks masses were presented in the previous subsection and the parameter set was obtained from the fit of the charmonium data. The tetraquark spectrum is calculated replacing the charm quark mass by the diquark mass m cc .
We now present the spectrum of the all-charm tetraquark considering the ground states 1S and the first orbital excitations 1P (relative to the diquark-antidiquark system), including all the possible combinations of total spin and total angular momentum. We also include the radial excitations 2S and 2P , in a total of 20 T 4c states built with two cc diquarks, each one of them being in an antitriplet color state and spin 1 (1 3 S 1 ). These 20 states were built considering the coupling of the total spin of the tetraquark S T (composed by the coupling of the total spins of the diquark S d and antidiquark Sd) with the relative orbital angular momentum L T between diquark and antidiquark, resulting in a total angular momentum J T of the tetraquark, in analogy to the cc charmonium spectrum. The corresponding parity and charge-conjugation quantum numbers of each combination are compiled in Table I. In our model the spin-spin interaction is treated non-perturbatively. In mesons and diquarks we had only two possibilities of total spin when combining two spin 1/2 particles S = 0, 1. Now, since we consider spin 1 diquark and antidiquark, we have three possibilities of total spin S T = 0, 1, 2, and therefore three different zeroth-order potentials, consequently three wavefunctions for each N L T state, as presented in Table VIII. The splitting structure from the perturbative corrections (spin-orbit and tensor) also has more possibilities, as presented in Table IX with the masses of the 20 T 4c states. In Figure 5 we show the mass spectrum.
From Table VIII we can observe that the tetraquark is very compact, in fact, its r 2 1/2 is even smaller than the ground state diquark. This should be taken as a warning against the point-like approximation, very common in the literature involving diquarks. However, since the color structure is assumed to be the leading mechanism responsible for the binding, due to the one gluon exchange between diquark and antidiquark, the two-body factorization should still be quite reliable. As suggested in Ref. [15], the two-body approximation is better for orbitally-excited states, such as the P -wave considered here, since the centrifugal barrier would suppress overlap at the origin. As we can see in Table VIII, the compactness of the T 4c is also reflected in the value of the wavefunction at the origin for the 1S states, which is very large.    Table IX we see that the compactness of 1S states is mainly caused by the Coulomb interaction. This suggests that the one gluon exchange is indeed the dominant mechanism responsible for the very strong binding between diquark and antidiquark, which causes the energy eigenvalue E to be negative. This also implies that the spin-spin interaction is strong. In this case we must have in mind that the factors coming from S 1 · S 2 are larger for the coupling of two spin 1 than for two spin 1/2 particles. It is interesting to see that even though the spin-dependent terms are now suppressed by a factor 1/m 2 cc and one would naturally expect them to be smaller when compared to the corresponding terms in cc mesons, the color interaction brings diquark and antidiquark so close that the suppression due to this factor is overwhelmed by the huge superposition at the origin of the system. The confinement term on the other hand, increases its contribution as radial or orbital excitations are included, as in cc mesons.
In Fig. 5 we can see that the masses of the 20 states are concentrated in the range between 6 and 7 GeV. Among the 1S states, the lowest one, with J P C = 0 ++ lies very close to the η c pair threshold. Within our uncertainties (both from the choice of parameters as well as the assumption of the diquark-antidiquark structure), we cannot say that this state is below or above such threshold. If it is above, it could be seen as a narrow state in the η c η c invariant mass. If it is below, then it would be stable against the rearrangement in cc pairs and other mechanisms would be necessary. Several possibilities are discussed in Ref. [30], with special attention to the 0 ++ lowest state, such as T 4c → DD through cc → g → qq . On the other hand, in Ref. [15] several decay possibilities of the orbitally-excited states are also discussed and branching fractions are estimated. It is interesting to see that even our estimates for the excited states with N L T = 2S, 2P are below the threshold of decay into doubly-charmed baryon pair due to light quark pair creation, cccc → (ccq) + (ccq), which is above 7 MeV.
The second lowest state, with quantum numbers J P C = 1 +− could rearrange itself into η c J/ψ. However, this state seems to be more than 50 MeV below this two-meson threshold, therefore it should be stable. The highest 1S state, with quantum numbers J P C = 2 ++ is also more than 50 MeV below the corresponding J/ψ pair threshold. It could still decay into η c pair in D-wave, but this mechanism should be suppressed.
In order to be consistent with our cc results, in Fig. 5 the two-meson thresholds are shown using the values of charmonium masses obtained with our model, which were compared to the experimental values in Table IV. In Table X we compare all the 1S and 1P T 4c states with the corresponding lowest S-wave two-meson thresholds. We see that while the 1S states lie close or below their thresholds, the orbitally-excited ones are close or above the corresponding thresholds. Therefore it would be interesting to search for these states in the two-meson invariant mass distributions, since some of them could show up as narrow peaks just around the threshold, like the 1 −− state (from 1 5 P 1 configuration) in the η c (1S) h c (1P ) invariant mass at 6.50 GeV, the 2 −− (from 1 5 P 2 configuration) in the J/ψ(1S) χ c1 (1P ) invariant mass at 6.60 GeV, and the 3 −− state (from 1 5 P 3 configuration) in the J/ψ(1S) χ c2 (1P ) invariant mass at 6.65 GeV.
One of these orbitally-excited state is of particular interest since it presents exotic quantum numbers that cannot be obtained as a simple cc system: the 1 −+ (from the 1 3 P 1 configuration). This state could be searched for in the η c (1S) χ c1 (1P ) invariant mass. However, it might be quite broad since our predictions show that it is about 80 MeV above its two-meson threshold.  Next, we comment the results of other works which also investigate the existence and properties of this state composed of four charm quarks. Some of them also consider the sextet structure of the diquarks (which can also lead to a color singlet tetraquark). In the following tables we present a compilation of the main results.  First, we show the results of Ref. [2] in Tables XI and XII. In this work a variational method with Gaussian trial wavefunctions was employed to study all-heavy tetraquarks, using a four-body coordinate system. The interactions were described with a potential due to the exchange of color octets in two-body forces. Two potentials were used: model I is a Cornell-type (Coulomb plus linear) and the model II is of the form A + Br β . Also a version of the MIT Bag model was used with the Born-Oppenheimer approximation. Both color structures were considered,3 − 3 and 6 −6. S-wave and P -wave were considered with both potentials, and spin shifts were calculated with the Cornell-like potential.
In Table XIII we compile the results of Refs. [66,78], where the T 4c production was studied. The estimates for the T 4c are very similar to those presented in this work. The authors used the diquark results of Ref. [77], where the cc diquark was calculated as a baryon constituent (we also compared these diquark results with ours). The same strategy of dividing the problem in two-body problems was used, but only S-wave states were calculated, and the spin-spin splitting was considered between each spin 1/2 constituent pair, using the wavefunction at the origin of the diquark or of the charmonium, depending on the interacting pair. It is interesting to see that the 0 ++ state is very close to our result, and the 1 +− is also below the η c J/ψ threshold. However, the 2 ++ is about 20 MeV above the J/ψ J/ψ threshold, indicating that this state could be seen in the J/ψ J/ψ invariant mass.  [66,78]. In Table XIV we compare our results for the S-wave T 4c with those of the recent diquark-antidiquark works: the one with antitriplet diquarks [66,78], and those with the color-magnetic model [27] and with QCD Sum Rules [28].  [66,78] Ref. [27] Ref. [ In Table XV we compare our results with the contribution of each term used to calculate the 0 ++ T 4c in Ref. [30], which was based in meson and baryon masses. The constant V 0 is obtained as twice the constant term S obtained from the fit of baryon and meson masses (which is added only into baryon masses, related to the QCD string junction, as discussed in that reference). Remember that in our model the spin-spin interaction is contained in the energy eigenvalue. At last, in Table XVI we compare our results for the P -wave T 4c with the old diquark-antidiquark predictions of Chao [15], the recent diquark-antidiquark predictions of QCD Sum Rules [28] and with lattice results [24].
Ref. [15] Ref. [28] Ref. [24] 1 −− 1 1 P 1 6.5771 6.55 − 6.82 6.83 − 6.84 6.420 1 −− 1 5 P 1 6.4954 6.39 The use of the Cornell potential allows us to study the charmonium spectrum without the confining interactions, which can be easily "switched off" by choosing the string tension to be zero. We can thus repeat all our calculations and check whether we find bound diquark states and also a bound T 4c . We have done these calculations and we find both diquark and tetraquark bound states. The obtained diquark and T 4c ground states have masses equal to m cc = 2881.4 MeV and T 4c = 5.3 − 5.4 GeV (for the lowest 1S states), respectively, as shown in the Appendix. These results can have applications in the context of relativistic heavy ion collisions, where a deconfined medium is formed (the quark gluon plasma, QGP). Our results suggest that the T 4c can be formed and perhaps survive in the QGP phase.

IV. CONCLUSION
In this work we have first updated the Cornell model, (a very well known and accepted model for the charmonium) obtaining a satisfactory reproduction of the charmonium spectrum, including the most recently measured states. We have then extended this model to study the all-charm tetraquark (cccc).
We explored a diquark-antidiquark configuration, including P -wave tetraquarks and we extended the spin-dependent interactions between diquarks, including a consistent strategy to deal with the tensor interaction between two objects of spin 1. The fact that our model is relatively simple compared to the four-body approach and to the relativistic models allows us to study many of the T 4c properties with clarity, specially the role of the spin interactions.
We were able to study the behavior of the all-charm tetraquark when radial and orbital excitations are included, investigating the contribution of the one gluon exchange, the confinement and the spin-dependent interactions, providing for the first time detailed results which elucidate the dynamics of the diquark-antidiquark structure.
The inclusion of one orbital excitation in the tetraquark also leads to a significant increase in the possibilities of quantum numbers and to the prediction of the exotic state with J P C = 1 −+ . The orbitally-excited cc states, having specific masses and quantum numbers, have different decay channels, which may be investigated experimentally.
Our model is simple and instructive, specially in what concerns spin interactions. For the lowest T 4c states our predictions are compatible with those made with other approaches. For the higher states, in particular those with orbital excitations, we make novel predictions which can be tested. In this region our predictions are more reliable, since the diquark-antidiquark spatial separation is bigger.
Our results and the others found in the recent literature on the T 4c tetraquark, taken together, should encourage a careful experimental search of these states at LHCb and Belle II. XVII: Results for lowest T4c states (1S) using ground state (1 3 S1) diquarks. Both diquark and tetraquark are calculated without the linear confinement term. Parameters are mcc = 2881.4 MeV, αs = 0.5202, b = 0, σ = 1.0831 GeV.