All-charm tetraquark in front form dynamics

We study all-charm tetraquarks in the front form of Hamiltonian dynamics using the many-body basis function approach known as basis light-front quantization. The model Hamiltonian contains transverse and longitudinal confining potentials and a one-gluon-exchange effective potential. We calculate masses of two-charm-two-anticharm states focusing on the lowest state. We also calculate two-quark and four-quark estimates of meson-meson breakup threshold. The results suggest that the lowest two-charm-two-anticharm state is not a tightly bound tetraquark. We discuss implications of the cluster decomposition principle for theories formulated on the light front and present our treatment of identical particles together with color-singlet restrictions on the space of quantum states.


I. INTRODUCTION
Even though four-quark states, called tetraquarks, have been studied for a long time (see Refs. [1][2][3][4] for some early studies), the stability of tetraquarks is still debated. One of the basic questions is whether there exist four-quark states whose masses are smaller than the sum of masses of two mesons, into which the tetraquark could potentially decay through rearrangement of quarks. Because ab initio calculations in QCD are challenging, researchers make use of various strategies and approaches to estimate the masses of tetraquarks and their results are often in conflict with each other . The goal of our paper is to initiate studies of tetraquarks within the framework of Front Form of Hamiltonian dynamics [41] and Basis Light-Front Quantization (BLFQ) [42], an approach whose ultimate goal is to achieve ab initio calculations in QCD. Therefore, our study is focused on the development of the approach as much as on providing a preliminary answer to the main tetraquark problemwhether or not four heavy quarks can form a bound state.
We choose to study heavy quarks (charm quarks) because for heavy quarks one expects that the proper, QCD-based, theoretical description can be simplified. Asymptotic freedom, which is believed to be relevant for heavy quarks, allows for perturbative expansion of the QCD Hamiltonian and produces some confidence that the simple Hamiltonian with confining and one-gluon-exchange potentials that we use shares important features with the full QCD Hamiltonian. Due to asymptotic freedom and quark masses much larger than the strong interactions scale Λ QCD , charm quarks are expected to be relatively slow in comparison with the speed of light, hence, additional pairs of heavy charm quarks cannot be easily produced and should not contribute significantly to the tetraquark dynamics. Tetraquarks of any kind are an interesting topic of study because they are exotic, i.e., they are neither mesons nor baryons, therefore, they provide opportunities to test and extend our understanding of hadron physics beyond the boundary of fairly well-established meson and baryon physics.
Finally, studies of all-heavy tetraquarks received recently additional motivation in the form of first experimental identification of all-charm tetraquark resonance X(6900) [43]. The discovery of a doubly charm tetraquark is also worth noting [44].
BLFQ has already been used with success to study various mesons and baryons [45][46][47][48][49][50][51][52][53][54] as well as in QED, see for example Ref. [55]. However, most of those studies involve only one Fock sector, with recently appearing extensions [56]. Questions like "how does confinement work?" cannot be fully answered by studying quark-antiquark or three-quark systems alone, even if one uses phenomenologically successful confining potentials. If one is to believe that gluon strings are formed in a Hamiltonian approach to QCD (as seems to be the case for Lattice QCD), then one is necessarily forced to explicitly include many-gluon sectors in addition to the leading "valence" Fock sector. Furthermore, breaking of those strings requires Fock sectors with additional quark-antiquark pairs. The strength of BLFQ stems from the fact that, in principle, it can handle many Fock sectors, each of which can contain many particles, in a straightforward manner.
The QCD Fock space is rich in structure and even with the help of supercomputers the calculations are challenging, because the dimensionality of required spaces of states grows quickly with the addition of new Fock sectors. The QQQQ sector is one of the natural next targets after the QQ and QQQ sectors.
Another important challenge resides in how to renormalize divergent interactions of QCD.
The eventual success of the approach will probably require an adoption of effective interactions calculated from QCD using, for example, Renormalization Group Procedure for Effective Particles (RGPEP) [57]. The Hamiltonian of bare, pointlike quarks and gluons leads to the problem of overlapping divergences [58]. RGPEP by defining effective, finite-size particles can tame singular interactions and reduce the number of Fock sectors necessary to obtain satisfactory results. Effective Hamiltonians computed using the closely related Similarity Renormalization Group [59] (see also Ref. [60]) have been successfully used in combination with many-body methods in ab initio calculations in nuclear physics, see for example [61][62][63]. However, a relativistic quantum field theory such as QCD is much more complicated than the non-relativistic nuclear many-body problem of interacting nucleons.
Since we choose to deal with only charm quarks and antiquarks we take into account antisymmetrization of identical particles. This is also the first system treated within BLFQ where the question about color-dependence of the confining potential needs to be addressed because there are two color-singlet combinations in the QQQQ sector whereas both the QQ and the QQQ sectors admit only one color singlet each. We adopt the commonly-used assumption that the confining potential depends on color in exactly the same manner as one-gluon-exchange interactions depend on color. We also add a color-independent term in the longitudinal direction. Without this added term we find some spurious, unphysical solutions with negative mass squared.
In Sec. II we present our model many-body Hamiltonian and derive Schrödinger-like equations for three cases -describing one meson, a tetraquark, and two mesons. The twomeson system allows us to discuss the cluster decomposition principle on the Light Front.
Section III is devoted to a description of the main elements of the computational framework of BLFQ. Our results for masses in the three mentioned cases and a discussion about whether all-charm tetraquarks are stable against dissociation are given in Sec. IV. Section V concludes the paper. Color factors between color-singlet states are given in the Appendix, where we describe the procedure that takes into account Pauli exclusion principle and allows us to work with color singlets only.

A. Front Form of Hamiltonian dynamics
Before we introduce our model Hamiltonian we mention a few aspects of the framework we use that are important in the context of our long-term goal of ab initio calculations in QCD. The Front Form of Hamiltonian dynamics [41] has two important advantages with respect to other Hamiltonian approaches. One of them is the fact that particles cannot be created from free vacuum in a way that they can be created, for example, in the Instant Form of Hamiltonian dynamics. Front-Form theories conserve total longitudinal momentum of particles taking part in the interaction, where the longitudinal momentum of a particle is defined as p + = p 0 + p 3 . In Hamiltonian approaches particles are on mass shell, hence, p 0 ≥ |p 3 | and p + cannot be negative. At the same time vacuum should have p + = 0, therefore, all particles created from vacuum should have exactly p + = 0. Since for massive particles p + → 0 means energy diverging to infinity one should regularize the theory and remove p + = 0 states, called zero modes. However, it is also known that one cannot simply discard those states and zero modes have to be taken into account in some way. Even though it is an open question about exactly what way zero modes need to be included we still gained something: the difference between the free vacuum and the interacting vacuum can only be contained in the singular point p + = 0. Therefore, to a large extent one can separate zeromodes from p + > 0 region, where most of the usual dynamics happen [similar in form to the Schrödinger equation or quark model Hamiltonians, see Eqs. (23), (28) and (36)]. This is in contradistinction from the Instant Form, in which particles of arbitrary momenta can be created from the free vacuum making the interacting vacuum a complicated state upon which one-, two-and many-particle states are to be built.
Another advantage of the Front Form of Hamiltonian dynamics is the fact that one can freely boost particles and wave functions can be decomposed into products of total and relative motion factors. This is very fortunate because one can use exactly the same wave functions that describe the internal structure of a hadron regardless of how fast the hadron is moving in the laboratory frame. Hence, the Front Form is uniquely suited to describe highenergy processes and offers practical advantages for building a Poincaré-covariant quantum theory in a Hamiltonian approach.
In the Front Form of Hamiltonian dynamics the Hamiltonian is P − = P 0 − P 3 . The momentum operators are P + = P 0 +P 3 , which is the longitudinal momentum, and transverse momenta P 1 and P 2 . We denote two-dimensional transverse vectors with a bold font, e.g., P = (P 1 , P 2 ). The evolution of quantum states is given by the analog of the Schrödinger equation, which in the stationary version is P − |Ψ = E |Ψ , where E is the eigenvalue of operator P − . One can also study the closely related eigenvalue equation, where the eigenvalue M 2 is the invariant mass squared of the eigenstate |Ψ . The eigenvalue M 2 depends only on relative motion of constituents and not their absolute motion. Since we work with P µ P µ instead of P − it is convenient for us to call H = P µ P µ the Hamiltonian.
It is sometimes referred to as "light cone Hamiltonian" [64]. Therefore, In the Front Form of Hamiltonian dynamics operators P + and P are kinematic while P − is dynamic. In other words, P − contains interactions, while P + and P are the same regardless of what interactions are present in the theory.

B. Hamiltonian
The model Hamiltonian that we use to study four-quark systems is where H kinetic , H transverse , H longitudinal and H OGE stand for kinetic term, transverse confining potential term, longitudinal confining potential term and one-gluon-exchange (OGE) term, respectively. The kinetic energy Hamiltonian is where P − 0 stands for the noninteracting, kinetic part of P − . The momentum operators are where p − 1 = (m 2 + p 2 1 )/p + 1 , where m is the quark mass and b 1 and d 1 are annihilation operators of quark and antiquark with label 1, respectively. Moreover, where c 1 and σ 1 are the color and the light-front helicity of particle 1, respectively. The where U conf,⊥ and U conf,z are the interaction kernels that depend on momenta and helicities of particles 1, 2, 1 , and 2 . The momentum conservation Dirac delta is The color dependence is encoded in BD OGE and BD CI , where t a ij stands for χ † c i T a χ c j , where T a = 1 2 λ a , with λ a a Gell-Mann matrix (a = 1, 2, . . . , 8) and χ c = [δ c,1 , δ c,2 , δ c,3 ] T is a three dimensional vector while c = 1, 2, 3 is the color quantum number. In other words, t a ij is half of the matrix element of matrix λ a in the c i th row and c j th column. The color dependence of BD OGE is the same as the color dependence of the one gluon exchange, hence, the subscript "OGE." On the other hand, BD CI is diagonal in color, color independent, hence, the subscript "CI." Both BD OGE and BD CI have three terms each that describe pair-wise interactions in quark-quark, quark-antiquark, and antiquark-antiquark pairs. The factor 1/2 that multiplies quark-quark as well as antiquark-antiquark terms is present because the two quarks, or the two antiquarks, that interact are indistinguishable.
Finally, a is a constant between 0 and 1, and C F = (N 2 c − 1)/(2N c ) = 4/3 is the value of quadratic Casimir operator in a fundamental representation of SU (N c ), N c = 3. We choose a = 0.85, therefore, in our Hamiltonian 85% of longitudinal confining strength in a meson comes from the OGE-like term and 15% comes from the color-independent term. See below for more detailed discussion.
The kernels are where κ is the interaction strength parameter, x 12 = p + 1 /(p + 1 + p + 2 ) is the longitudinal momentum fraction of particle 1 with respect to 2 and x 21 = 1 − x 12 is the longitudinal momentum fraction of particle 2 with respect to particle 1. Relative transverse momentum is k 12 = x 21 p 1 − x 12 p 2 . Moreover, and Objects with subscript 1 2 are defined in the same way as objects with subscript 12, except that 1 is replaced with 1 and 2 is replaced with 2 .
The confining potential is determined by the anti-de Sitter (AdS)/QCD holography [65] and its transverse part reproduces the AdS/QCD harmonic oscillator in the QQ sector. In appropriate momentum variables [66], in the QQ sector, the longitudinal and transverse terms complement each other and form a three-dimensional, rotationally invariant harmonic oscillator, see Eq. (26). The potentials in the QQ sector are naturally extended to other sectors through Eqs. (9) and (10), which act in all sectors. The extension, however, is not unique. For example, the factor p + 1 + p + 2 could be replaced with the total P + . Moreover, BD OGE and BD CI evaluate to the same expression between states in the QQ sector up to a factor of C F . Their combination, as in Eq. (10), gives the result that is independent of a in the QQ sector. Our choice of the confining potential was obtained after a study of several variants and searching for acceptable spectral behavior of the solutions.
We found that removing color independent part or replacing p + 1 +p + 2 with P + leads to the appearance of unphysical solutions with negative mass squared. While in general tachyonlike states can be a sign of unstable equilibrium in a linear approximation of a field theory, see Ref. [67], our approach is nonperturbative and we are dealing with model Hamiltonians. Therefore, we regard the candidate model Hamiltonians with such tachyonic solutions as unphysical. The properties of those states are very far from properties expected of bound tetraquark states. For example, the dominant components of wave functions of those nonphysical states reveal very fast motion of quarks with respect to each other making them more like highly excited, high momentum scale states than like states characterized by low relative momenta appropriate to our model. a = 0.85 is the largest value of a that guarantees no negative M 2 states appear up to K = 50 for N max = 6 (see Sec. III). It is worth noting that two-and more-gluon exchange potentials are in general mixtures of OGE-like and color independent parts. Hence, our confining potential appears reasonable, apart from the fact that our CI potential confines at large distances. However, the states which should be affected the most by this confinement are the excited states while we focus mainly on the ground state.

The Hamiltonian term of the one gluon exchange interaction is
The kernel of the OGE term is, where D is the energy denominator, with µ being a fictitious gluon mass. We use the same spinors as those in Ref. [55] and u 1 γ µ u 1 ū 2 γ µ u 2 can be found in Table I therein. The fictitious gluon mass µ is introduced to regulate the Coulomb singularity: if we take p 1 = p 1 and which is in the denominator of Eq. (20), becomes zero, unless µ = 0. This singularity is integrable if momenta are continuous, however, in BLFQ we discretize longitudinal momenta and the singularity has to be somehow regulated. Even though diagonal matrix elements of the discretized version of H OGE diverge as µ → 0 the eigenvalues and eigenvectors approach a finite limit.
The Hamiltonian of Eq.
(3) provides a unified description of QQ and QQQQ systems. In fact, one could apply this Hamiltonian in sectors with arbitrary number of heavy quarks and antiquarks. We use it in three separate calculations for three purposes. In all three cases we restrict the space of states to color singlets, which can be achieved since H conserves color.
Details are provided in the Appendix. Firstly, we solve the QQ eigenvalue problem and, by fitting the numerical spectrum to the experimental spectrum of charmonium, we fix the free parameters of the Hamiltonian: quark mass m, confining potential strength parameter κ and OGE coupling constant g. Secondly, we solve the QQQQ eigenvalue problem to find the four-quark ground state mass. Thirdly, we solve the QQQQ eigenvalue problem with some interactions turned off. The interactions that are kept allow one quark to form a meson with one antiquark and the other quark to form a meson with the other antiquark. There is no interaction between the two mesons and we restrict the space of states to the states in which both mesons are color singlets separately. This way we can numerically estimate the two-meson threshold, which can be different than the sum of masses of two mesons obtained in the QQ calculation due to finite basis. Below we briefly present the three cases.

C. Eigenvalue equation for mesons
The Hamiltonian can have many eigenvectors of various forms. States that describe a single meson with fixed momenta P + M and P M are of the form, The "front-form energy" of the meson is where M is the mass of the meson.
One can simplify the form of this equation considerably by changing variables from k 12 and x 12 to q 12 and q z 12 (collectively denoted q 12 ) introduced in Eqs. (16) and (17). Moreover, we assume that the meson is a color singlet state. Therefore, We get, This equation looks much like nonrelativistic Schrödinger equation in momentum space. The Laplacian acting on the wave function is equivalent to a rotationally symmetric harmonic oscillator potential and the OGE potential is written in a generic form. It is worth noting that the same confining potential can be derived using RGPEP with a gluon mass ansatz [68,69].
Our OGE potential is different from the Coulomb plus Breit-Fermi of Ref. [68] and is taken instead from Ref. [55]. The choice was dictated by the availability of software implementation of the latter potential. Similarly, instead of the longitudinal potential given by Eq. (15) we could have chosen a kernel that would give us ∂ x x(1 − x)∂ x potential of Ref. [45,46]. In the limit of relative momenta vanishing with respect to quark masses the two potentials become equal, hence, both should be suitable for phenomenology. It is sufficient for our purposes to select one longitudinal confining potential and one OGE potential and work with them.

D. Eigenvalue equation for tetraquarks
Tetraquark states have a form very similar to meson states, This state has fixed momenta P + T and P T , while The eigenvalue equation where x i = p + i /P + T and 4 i=1 x i = 1, while k i is a transverse momentum of particle i in a rest frame of the bound state where 4 i=1 k i = 0. The harmonic oscillatorŨ ij for quarkantiquark interaction is given in Eq. (24) with 1, 2, 1 , and 2 replaced by i, j, i , and j , respectively. For quark-quark and antiquark-antiquark it is, W OGE (i, j; i , j ) ψ is different depending on i and j. For the quark-quark interaction, i.e., i = 1 and j = 2, For the antiquark-antiquark interaction, i = 3 and j = 4, For quark-antiquark interactions, i = 1 or 2 and j = 3 or 4, where ψ i j = ψ T (1 23 4), ψ T (1 234 ), ψ T (12 3 4), and ψ T (12 34 ) for ij = 13, 14, 23, and 24, respectively. The interaction kernels are antisymmetrized as a result of having identical

E. Eigenvalue equation for two mesons
To describe two separate mesons, A and B, we choose, Meson A has momentum components P + A and P A , while meson B has momentum components P + B and P B . By placing the two mesons far enough from each other we can make the total interaction between them to be arbitrarily small. We simulate this situation by turning off all interactions except those between particles 1 and 3, which form meson A, and between particle 2 and 4, which form meson B. Moreover, two identical quarks contained in two separated mesons are in practice distinguishable. Therefore, in this section we treat all particles as distinguishable. Since there are no interactions between the two mesons, we expect that in the general eigenvalue equation, the eigenvalue M 2 can be written as the invariant mass of two mesons with mass M A and M B , The relative transverse momentum between mesons is where Note, that the relative transverse kinetic energy between the mesons in the eigenvalue, canceled with the transverse kinetic energy between mesons in H kinetic , which is fixed by the choice of state |ψ AB . We separate Eq. (36) into two, (24). Using the same kind of substitution as in Sec. II C, we get two eigenvalue equations, We could formally restore the decomposition principle by replacing p + 1 + p + 2 in Eqs. (9) and (10) with P + , but this would lead to the appearance of spurious states as described in Sec. II B. We prioritize the acceptable spectrum over exact conservation of the decomposition principle, since the former is more important in practice, while the latter can be approximately restored. Since charm quarks are heavy, the two-meson system and tetraquark can be considered as nonrelativistic. Therefore, x A ≈ x B ≈ 1/2 and one can partly restore the cluster decomposition principle by rescaling κ in the QQQQ sector. In other words, in the QQQQ sector we use κ T = 2 1/4 κ instead of κ. This guarantees that As opposed to the confining potential, the OGE potential fully obeys the cluster decomposition principle. This is to be expected because it can be derived from QCD in perturbation theory. In fact, all two-body potentials in QCD have the same generic form of Eq. (19) (apart from the color factors, which may differ). It is important that U OGE (1, 2; 1 , 2 ) depends only on relative momenta between particles 1 and 2 and between 1 and 2 , and not on, e.g., momentum fractions x 1 or x 2 , which depend on the total P + of the system. Therefore, Eq. (19) illustrates the general form of two-body operators that admit the cluster decomposition principle in the sense described here. A more general treatment of relativistic theories obeying cluster separability can be found in Ref. [70].
In the eigenvalue equations the cluster decomposition principle manifests itself by the presence of 1/x A and 1/x B factors in Eq. (36) and the 1/(x i + x j ) factor in Eq. (28) that multiply interaction terms that depend only on relative momenta within the interacting pair with no trace of the total P + . Note that momentum conservation in i j (p + i + p + j )δ ij.i j fixes p i + p j and one is left with an integral over relative momenta x i j and k i j . This leads to discretization of the longitudinal momenta. We apply antiperiodic boundary condition for the quark field, which means that quark longitudinal momenta can only take values, where k is called the longitudinal quantum number and it is a positive half-integer. In sectors with many particles the total longitudinal momentum is by definition P + = 2π L K, where K = i k i is the sum of longitudinal quantum numbers of all particles. In Sec. II P + denoted the momentum operator, from now on P + means the eigenvalue of the operator P + and we keep it fixed (we use only eigenstates of the operator P + with eigenvalue P + ). For a given particle i, the longitudinal momentum fraction x i is, The longitudinal continuum limit is L, K → ∞ while keeping P + fixed. None of the quantities that we calculate depend on the exact values of P + and L due to Front-Form boost invariance.
For transverse momenta we introduce the harmonic oscillator basis [42]. We define new creation and annihilation operators, Note that the operators B i and D i depend on discrete quantum numbers n i , m i , k i , σ i and c i , while plane-wave operators b i and d i depend on continuum transverse momentum p i = √ x i q, discretized longitudinal momentum p + i = 2πk i /L (or equivalently on k i ), and on spin and color σ i and c i . Operators B i and D i are normalized to unity, that is, The basis wave functions are, where L |m| n are the associated Laguerre polynomials, q = (q 1 ) 2 + (q 2 ) 2 , ϕ = arg q, and b is a selectable positive parameter of dimension of P. The principal quantum number n is a non-negative integer, while m can be an arbitrary integer. The choice of the harmonicoscillator wave functions is compatible with our choice of the transverse confining potential and is important for the factorization of the center-of-mass motion, which we describe in detail later in this section.
In practice one has to truncate the many-particle basis in the transverse direction by limiting the allowed radial numbers n i and angular numbers m i by a cutoff in the number of oscillator quanta in each basis state, Removing this truncation is equivalent with taking the limit N max → ∞. In addition, we require our multi-particle basis state to have total angular momentum projection, where σ i = ± 1 2 is the fermion light-front helicity. Throughout this article we limit our attention to M J = 0 states. 1 It is also worth mentioning that the truncation of the basis breaks the cluster decomposition principle. For example, if we consider our two-meson example from Sec. II E and if the quantum numbers of particles forming meson A already almost saturate Eq. (49), then the particles of meson B will be restriced to a much smaller space of states than the particles of meson A. The opposite situation is also possible and included in the truncated basis. Therefore, one meson can influence the other through the truncation, even if there are no interactions between them. Moreover, each of the mesons in the QQQQ sector is subject to a different truncation than the one meson in the QQ sector and meson masses in the QQ and in the QQQQ sectors can differ slightly, but the difference should vanish as the basis size is increased.
It is straightforward to rewrite the Hamiltonian presented in Sec. II using new operators B and D. One has to additionally discretize the longitudinal momenta by following a simple prescription, 4πδ(p is only a matter of computing matrix elements of H and diagonalizing the obtained matrix to obtain eigenstates of H and their masses. Computation of matrix elements between states containing two quarks and two antiquarks is not much more complicated than the analogous computation between states containing only one quark and one antiquark because no particle nor any pair of particles is distinguished. One, obviously, has to calculate terms for all six pairs of particles instead of just one, and interactions between identical particles must be property antisymmetrized. Using a basis in relative momenta of Jacobi type, for example, would require us to use different formulas for different pairs of interacting particles. It should be evident that the addition of more particles in our calculation (including gluons) would be straightforward. Admittedly, this comes at a cost of larger matrices (effectively one more particle per Fock sector compared to Jacobi coordinates), but the larger matrices are also more sparse which aids applications on modern computers, while the simplicity makes the software development more reliable. Probably the most important complication is introduced by restricting our space of states to only color-singlet states. This important, but rather technical topic is described in more detail in the Appendix. Similar basis spaces restricted to include only color singlets have been implemented for a BLFQ treatment of glueballs with Fock spaces having up to six gluons [71].
Since BLFQ implements states using single-particle transverse motion instead of relative where P is the transverse momentum operator and R is the transverse center-of-mass position operator, Eigenvalues of H CM are n · 2b 2 λ CM , where n is a non-negative integer. n = 0 corresponds to the ground state of center-of-mass motion and n ≥ 1 correspond to excited states of centerof-mass motion. In a typical scenario among eigenstates of H with the lowest eigenvalues there will be states with the same relative motion but different center-of-mass motion. To keep only the eigenstates with the ground-state center-of-mass motion we diagonalize H + H CM instead of H. Since H and H CM commute they have the same eigenvectors, while the eigenvalues of the sum will be the sum of the eigenvalues of H and H CM . Therefore, states with excited-state center-of-mass motion will be shifted up by a multiple of 2b 2 λ CM .
Choosing λ CM sufficiently large and positive, all states with excited center-of-mass motion will have eigenvalues larger than the eigenvalues of the limited number of states that we obtain numerically. We use λ CM = 50 in our calculations.

IV. ANALYSIS OF BINDING ENERGY
To address the question whether there exist cccc states that cannot break up into two charmonia we need to know the mass of the lowest tetraquark state and the value of the two-charmonium threshold taking into account the implications of the truncated basis space for the subsystems. We obtain estimates of both by numerically diagonalizing truncated matrices of our model Hamiltonian obtained using BLFQ. We therefore solve three problems which correspond to three eigenvalue equations presented in Secs. II C, II D and II E. The discrete spectra of truncated Hamiltonians should look more and more like the spectrum of the untruncated, infinite Hamiltonian as N max → ∞ and K → ∞.
The Hamiltonian matrix in the sector with one meson is used to fix free parameters of the model, m, κ and α = g 2 /(4π). The gluon mass µ = 10 MeV and the basis parameter b in the meson calculation is fixed to be equal to κ. We fit the lowest eight states in the spectrum of charmonium. The root mean square difference between fitted and experimental masses is 31 MeV. The parameters are given in Table I, while Table II lists   Tetraquark masses are calculated for N max = 6, 8, 10, 12 and for K = 6, 10, 14, 18.
We calculated three sets of tetraquark masses: masses with all interactions turned on, cf.
Eq. (28) Moreover, as already mentioned in Sec. II F, instead of κ we use κ T = 2 1/4 κ for the confining strength parameter.
One of the sources of systematic errors of the framework we adopt originates from the fact that a pair of particles in cccc system has a minimal nonzero kinetic energy with respect to the other two particles. The minimal kinetic energy should approach zero from above as N max approaches infinity, but may be of importance for finite N max . This artifact of a finite basis is called "kinetic energy penalty" in Ref. [6]. Here we estimate it in the following way, where N 2 and K 2 are fixed by conditions N 1 + N 2 = N max and K 1 + K 2 = K. M 2 free cccc and M 2 free cc are tetraquark and meson ground-state masses squared, respectively, computed with all interactions turned off. We use masses squared instead of masses because they, and not the masses, are the eigenvalues of our Hamiltonian. Moreover, the two two-quark masses, M 2 free cc (N 1 , K 1 ) and M 2 free cc (N 2 , K 2 ), need to be combined according to Eq. (35) to get the invariant mass of the full state. To get the minimal invariant mass we put k AB = 0 and minimize over all possible values of N 1 , N 2 and K 1 , K 2 into which N max and K can be partitioned. Hence, ∆M 2 , being the difference between the actual tetraquark mass and the minimal possible mass of two separate two-quark subsystems, is a measure of minimal k AB between the two subsystems. Table III lists the values we obtain. Note, that ∆M 2 does not depend on K for the choice of Ks we made. It should stay the same for all K = 2 (mod 4).
We correct the actual eigenvalues of the truncated Hamiltonians by subtracting the kinetic energy penalty,  We introduce three estimates of the threshold with which we compare our numerical tetraquark masses. One estimate uses the same idea behind the second term in Eq. (53) but with full meson masses that include interactions, where the minimum, just like in Eq. (53), is taken over all possible values of N 1 and K 1 .
Threshold T 1 gives an unexpectedly poor estimate. It is substantially smaller than twice our fitted numerical mass of η c . The reason seems to be overestimation of the OGE potential for small values of K because the minima of T 1 tend to be reached at the minimal K 1 = 1, while turning off OGE potential makes the minima to appear for K 1 = K 2 = K/2. In fact, one naively expects the minimum in the definition of T 1 to be reached for K 1 = K/2, because it implies x A = x B = 1/2, which means zero relative longitudinal momentum between the two mesons (as long as they have equal masses). Moreover, the actual M 2 full cc turns out to be negative for some K = 1 cases, which is unacceptable. Therefore, we define another estimate of the threshold, for which both N max and K are equally partitioned among N 1 , N 2 and K 1 , K 2 , i.e., This estimate seems to be more reasonable and it is in rough agreement with a third estimate of the threshold provided below.
By turning off interactions between particles that do not belong to the same meson we can compute numerically the invariant mass of two mesons occupying almost the same finite basis -we have to make identical particles distinguishable because otherwise one would not be able to consistently turn off, for example, an interaction between 1 and 4 and at the same time keep interaction between 1 and 3 turned on. Therefore, we define, where M 2 two-meson is the ground state mass in the aforementioned calculation of two-meson system in a tetraquark calculation. Results for threshold estimates and tetraquark masses are summarized in Table IV and plotted in Fig. 1. to go down in the limit N max → ∞ (provided we do not refit our meson masses), but we expect that the shift should be much smaller than the shift due to K → ∞ extrapolation.
Our gluon mass introduces additional shift upwards on the order of the value of µ, i.e., 10 MeV. All tetraquark masses lie substantialy above all threshold estimates, including the extrapolations. These results indicate that the lowest cccc eigenstate of our model Hamiltonian is not bound with respect to breakup into two separated mesons. It could be a resonant state. However, such a conclusion would require additional confirmation in the form of decay analysis.

V. CONCLUSION
We have done the first, to our knowledge, study of all-heavy tetraquark states using a Hamiltonian in the Front Form of dynamics, where all quarks are treated individually, color degrees of freedom are unconstrained (apart from the restriction to global color singlets) and antisymmetrizations due to identical particles are taken into account.
We note, however, that our confining potential breaks the cluster decomposition principle, but the breaking should be rather small for a nonrelativistic system like all-charm tetraquark. Attempts to restore it exactly lead to unphysical states with negative mass squared. Therefore, our confining potential should be regarded as an approximate effective potential with a limited range of applicability.
Even without the negative M 2 problem, confining long-range forces lead to problematic long-range van der Waals forces [72,73]. Such long-range forces are unlikely to be present in QCD. A more probable picture would involve effective, massive gluons to be the source of confining forces. They may or may not form strings, but in any case a force mediated by such gluons would be short-ranged.
All our estimates for the cccc ground state mass turn out to be substantially higher than the estimates we made for the lowest threshold for breakup into two cc mesons. Therefore, in our model, the ground state tetraquark is unstable against dissociation into two charmonia. It is also worth noting that the results for tetraquark masses seem to be much more reliable than the threshold estimates that we obtain, as can be seen in Fig. 2  Ref. [42] for detailed examples of numbers of color singlets in sectors with more particles.
One has to, however, invest extra effort in the evaluation of the matrix elements of the Hamiltonian.
The first step, which needs to be done in any case, is to define a space of states with arbitrary color that takes into account that some particles are identical. One can use states Since quarks are fermions this can be done using a relation of strict order. We say that 1 > 2 if and only if k 1 > k 2 , or k 1 = k 2 and n 1 > n 2 , or k 1 = k 2 and n 1 = n 2 and m 1 > m 2 , or k 1 = k 2 and n 1 = n 2 and m 1 = m 2 and σ 1 > σ 2 , or k 1 = k 2 and n 1 = n 2 and m 1 = m 2 and σ 1 = σ 2 and c 1 > c 2 . We define our basis to contain only such states |1234 for which 1 > 2 and 3 > 4. Taking another such state |1 2 3 4 with 1 > 2 and 3 > 4 we have In the second step we need to find color-singlet states, which are defined as the kernel of the quadratic Casimir operator C 2 = 8 a=1T aT a , wherê where 12 is the sum over all quantum numbers of particles 1 and 2. We omit the part for gluons, because we do not have gluons in our model. The color operatorsT a do not change any of the momentum and spin quantum numbers, hence, C 2 is diagonal in momentum and spin. Therefore, we can separately diagonalize C 2 in subspaces of fixed momentum and spin quantum numbers. Note that our relation i > j compares colors c i and c j in the very end, only if all other quantum numbers turned out to be the same.
There are four different kinds of subspaces. To classify them it is convenient to introduce another two relations. We say that i ≈ j if all quantum numbers of i and j except color are the same (colors can be arbitrary). We say that i j if i > j and not i ≈ j. In other words, either k i > k j or k i = k j and n i > n j or k i = k j and n i = n j and m i > m j or k i = k j and n i = n j and m i = m j and σ i > σ j . Hence, is just like > except it does not take into account color. We can now easily classify the four cases of color spaces.
where kets on the right hand sides are denoted by colors |c 1 c 2 c 3 c 4 and we omit momentum and spin quantum numbers, which are the same for each ket. Instead of 1, 2, 3, colors are called r, g, b, respectively, for quarks andr,ḡ,b, respectively, for antiquarks. For completeness, r < g < b andr <ḡ <b. Note that |1234, S is symmetric for either 1 ↔ 2 or 3 ↔ 4, while |1234, A 1 is antisymmetric for either 1 ↔ 2 or 3 ↔ 4. The third and last step is to calculate the common factors in Hamiltonian matrix elements between states with different color-singlet configurations that arise due to color and antisymmetrization. We summarize the results. In general we need to evaluate matrix elements of the following types of operators, where V 5 ;5 and V 5 ,6 ;5,6 depend on all quantum numbers except color, while C qq 5 ,6 ;5,6 , Cqq 5 ,6 ;5,6 and C qq 5 ,6 ;5,6 depend only on color. There are two types of interactions: color-independent ones, for which, C qq 5 ,6 ;5,6 = Cqq 5 ,6 ;5,6 = C qq 5 ,6 ;5,6 = δ c 5 ,c 5 δ c 6 ,c 6 , (A. 13) and OGE-like interactions, for which, Note, that C qq 5 ,6 ;5,6 = Cqq 5 ,6 ;5,6 . We label color singlets with capital letters I and J that can equal to S, A 1 , A 2 , A 3 or A 4 . The general matrix elements are, where e i stands for all quantum numbers of particle i except color, i.e., n i , m i , k i and σ i .