Effect of relativistic kinematics on the stability of multiquarks

We discuss whether the bound nature of multiquark states in quark models could benefit from relativistic effects on the kinetic energy operator. For mesons and baryons, relativistic corrections to the kinetic energy lead to lower energies, and thus call for a retuning of the parameters of the model. For multiquark states, as well as their respective thresholds, a comparison is made of the results obtained with non-relativistic and relativistic kinetic energy. It is found that the binding energy is lower in the relativistic case. In particular, $QQ\bar q\bar q$ tetraquarks with double heavy flavor become stable for a larger ratio of the heavy to light quark masses; and the all-heavy tetraquarks $QQ\bar Q\bar Q$ that are not stable in standard non-relativistic quark models remain unstable when a relativistic form of kinetic energy is adopted


I. INTRODUCTION
In a rather celebrated paper [1], Isgur and his collaborators analyzed how to take out the "naive" and "nonrelativistic" out of the quark model. After this pioneering work, there have been several other studies on how to implement a minimal amount of relativity in the quark model, e.g., [2][3][4][5].
We are aware that this is just a small part of the problem. For instance, in the case of the positronium atom or ion, there are many effects, very often canceling each other, and the kinematics is just one of them. See, e.g., [6,7]. We nevertheless deemed appropriate to examine the role of relativistic kinematics in the quark model and to go beyond the case of mesons and baryons. The main motivation is the current interest about the stability of multiquark hadrons and the delicate interplay between the energy of collective configurations and their corresponding thresholds. This is precisely the aim of the present article, to study to which extent the stability of tetraquarks is influenced by relativistic kinematics. This implies that the mesons constituting the threshold and the tetraquarks are estimated consistently within the same framework, either non-relativistic (NR) or relativistic, or, more precisely, semi-relativistic (SR), to keep in mind that we take into account only a fraction of the relativistic effects. For a given potential, all hadron masses tend to decrease if one adopts a relativistic form of the kinetic energy. It * j-m.richard@ipnl.in2p3.fr † valcarce@usal.es ‡ javier.vijande@uv.es will be shown that the mesons are more affected than the tetraquarks, so that the binding energy with respect to the lowest threshold becomes smaller for bound tetraquarks.
For each hadron, we concentrate on its energy E, so that its mass is given by M = m i + E, where the m i are the constituent masses of its quarks. For the discussion about the stability of tetraquarks, we have the same cumulated constituent mass, m i , in the tetraquark and in the mesons entering its threshold.
The paper is organized as follows. In Sec. II, we review the formalism, focusing on how to estimate the matrix elements in a variational calculation using a basis of correlated Gaussians. Then, some applications are given for ordinary hadrons in Sec. III, and for tetraquarks in Sec. IV. In particular, we discuss how the relativistic effects influence the binding of doubly-heavy tetraquarks with respect to their dissociation into two flavored mesons and how the fully-heavy tetraquarks are affected by the choice of kinematics. Some conclusions are drawn in Sec. V.

II. FORMALISM
Let us start with the one-body Hamiltonian in three dimensions with K(p, m) being either K NR = p 2 /(2 m) or K SR = (p 2 + m 2 ) 1/2 − m. The ground state is sought at variationnally using a trial wave function Ψ(r) = Y 00 (r)u(r)/r arXiv:2103.02222v1 [hep-ph] 3 Mar 2021 The computation involves the matrix elements of normalization, potential and kinetic energy between two basis functions w(a, r) and w(b, r), which are noted as where c = (a + b)/2 andc = 2 a b/(a + b) are suitable averages. By straightforward calculation, one obtains g V (c) = 3/(2 c) for an harmonic potential, 2/ √ c π for a linear one, 2 c/π for a Coulomb one, etc. As for the kinetic energy, in the NR case, one gets f NR (c) = 3c/(4 m), while in the SR case, it is where some care is required in the computation of the Bessel function K 1 , at least with some of the available softwares.
The variational expansion (2) converges remarkably well. For instance, for m = 1 and V (r) = r, for which the ground state energy is exactly the negative of the first root of the Airy function, one can estimate whose vanishing ensures that the Temple-Kato lower bound coincides with the variational upper bound, as explained, e.g., in [8]. For a single Gaussian (N = 1), one gets δ 2 H = 0.027, and for N = 6, δ 2 = 0.00002. It has been checked that a similar convergence is obtained for other potentials and/or kinetic-energy operators. For a meson, in the center of mass, the NR Hamiltonian is just (1) with m in K NR (p, m) being the reduced mass. For the SR Hamiltonian, one simply adds up K SR (p, m 1 ) and K SR (p, m 2 ). Hence, the estimate of the energy and wave function is rather similar to that of the one-body case.
For a baryon, in the center of mass, the Hamiltonian reads where r ij = |r j − r i |, and p i = 0. In the NR case, it is customary, but not compulsory, to introduce Jacobi coordinates such that the total kinetic energy operator becomes diagonal. For our purpose, it is sufficient to introduce an universal set of variables in momentum space and their conjugate so that the individual momenta p i are linear combination of P , p x and p y , and the distances r ij linear combinations of x and y. The most general spatial wave function corresponding to an overall S-wave reads or, in abbreviated form, |ψ = γ i |A i , where the symmetric, definite-positive, matrices A i contain the range coefficients. As in the one-body case, the minimisation is reached in two steps: for a given set of A i , the coefficients γ i and the variational energy Ψ|H 3 |Ψ / Ψ|Ψ are given by an eigenvalue equation, and this energy is minimized by varying the A i .
The matrix elements of interest are obtained by standard techniques of Gaussian integration [9,10]. If a pair corresponds to a separation Similarly, for the kinetic energy, if, say The extension to tetraquarks is straightforward. One can introduce momenta such as and the conjugate variables in position space, and express the individual momenta and the relative distances in terms of them. Then, Eqs. (9)-(11) are easily extended from dimension 2 × 2 to 3 × 3.
The convergence of the variational expansion (9) is as good as for the one-body case. The Gaussian expansion method, indeed, is well documented, as discussed, e.g., in a paper written by some of the best experts of few-body physics, with benchmark calculations in various fields of physics [11]. In the case of baryons (or multiquarks), care should be taken that all degrees of freedom are incorporated, and this is achieved if one introduces 2 × 2 (3 × 3 for tetraquarks) matrices which are not restricted to a scalar or diagonal form. For instance, if a system of three equal masses is described with an expansion γ i exp(−a i (x 2 + y 2 )/2), one will never account for the possibility of internal orbital momenta x = y > 0.
Note that the formalism does not need to be modified to handle unequal masses. On can still use the Jacobi coordinates (7) and (12) in momentum space and their conjugates in position space. Let us insist on that the center-of-mass energy is exactly removed if one adopts a variational wave function that is invariant under translations. For instance, in atomic physics, many estimates of M + e − e − ions are made in a frame attached to the positive nucleus, and the corrections are taken into account by the "mass-polarization" term. The beautiful proof of the stability of the positronium molecule by two pioneers of quantum physics [12] was achieved with a trial wavefunction that depends only on some relative separations r j − r i . Nevertheless, this proof was unjustly criticized by arguing that the center of mass is not removed beforehand in the Hamiltonian [13].

A. Mesons
As a first illustration, we compare in Fig. 1 the energy calculated with a potential V (r) = r and constituent masses (m, m). The energy scale is given by the string tension set to σ = 1. In the NR case, all energies are proportional to m 1/3 , while in the SR case, they are recomputed for each m. Not surprisingly, the SR energy is lower than the NR one. Indeed, We now turn to the case of a pure Coulomb interaction. Both NR and SR energies are proportional to m, so we set m = 1 and deal with the well-studied "Herbst" Hamiltonian and its NR analog with p 2 /2. As discussed in the literature [4,[14][15][16][17], the Herbst Hamiltonian becomes delicate when α approaches 2/π. It is not our aim to enter the mathematical subtleties of (14), but to stress that • in the Coulomb regime, the relativistic effects are governed by the strength of the attractive 1/r ij terms, • when this strength increases, the convergence of the variational expansion becomes more delicate. In Fig. 2 is shown a comparison of the ground state of the Herbst Hamiltonian (14) and its NR version. Now for a realistic potential such as V (r) = −0.4/r + 0.2 r, where V is in GeV and r in GeV −1 , the SR and NR energies of (m, m) mesons significantly differ for light quarks, say m < 1 GeV, and become close in the heavy quark regime, say m ∼ 3 -5 GeV. But for very heavy constituents, the Coulomb interaction becomes dominant, and the ratio of energies increases again. This is shown in Fig. 3.
For all above potentials, the decrease of the meson energy when going from the NR to SR case is accompanied by a significant increase of the wave function at the origin,Ψ(0) = u (0)/ √ 4 π. For instance for a binary system with masses m 1 = m 2 = 1 bound by V (r) = r, one gets exactly |u (0)| 2 = 1, while the SR analog is |u (0)| 2 1.8. For a Coulomb interaction −g/r , the effect is even larger: the ratio of the |u (0)| 2 is nearly 3 for g = 1/2.
For quark models, this dramatic enhancement of the short-range correlations implies a drastic re-tuning of the spin-spin part of the interaction, whether or not it is treated as first-order perturbation, or included nonperturbatively in the model. For the latter option, an example is the AL1 potential [18], according to which the interaction of a quark of mass m i and an antiquark of

B. Baryons
For symmetric baryons (m, m, m), the same pattern is observed as for mesons. For instance, in a simple linear potential r ij /2, the ground state is found at E ∼ 3.863 in the NR case for m = 1, and is lowered to 3.522 in the SR case. This corresponds to a decrease by about 9%, very similar to the decrease by 8% observed for a two-body system with masses m 1 = m 2 = 1 and bound by V (r) = r.
Other mass combinations have been studied, with the results given in Table I. There is a smooth transition from the SR to the NR regime. We now consider a doubly-heavy baryon (M, M, m) with m = 0.5 and M = 5 whose quarks are bound by a pairwise interaction V (r ij )/2 and the same Coulombplus-linear potential V (r) = −0.4/r + 0.2 r as earlier for mesons. We obtain E 0.652 in the SR case, lower than the 0.780 in the NR case.
If one adopts a naive diquark model, i.e., a two-step method in which one first solves for QQ using the QQ potential alone, and then for QQ-q using V (r) where r is the distance from the light quark to the center of the diquark, 1 one obtains 0.746 in the NR case, while the SR scheme gives 0.473. This means that the distortion induced by the naive diquark model is seemingly amplified in semi-relativistic calculations.

IV. TETRAQUARKS
We now study tetraquarks, first in simple toy models, and then in more realistic potentials AL1 and AL1N tuned to reproduce the ordinary hadrons. We concentrate on two types of mass distribution for QQqq: M M mm for which binding is expected if the mass ratio is large enough, and the "fully-heavy" configurations, QQQQ, for which there are somewhat conflicting predictions in the literature.

A. Doubly-flavored tetraquarks
For QQqq, it is known for decades that in a static and flavor-independent potential, the system becomes bound if the mass ratio M/m is large enough. Moreover, in the case of QQūd with isospin I = 0 , the binding is favored by the chromomagnetic interaction. See, e.g., [19], and the recent discussion about various dynamical effects [20].
We start with a purely chromoelectric model, namely a potential and utilize the usual notation T =33 and M = 66 for the color states in the qq-qq basis. The energy scale is fixed by imposing a potential of unit strength for a quark-antiquark pair forming a meson. In Fig. 4 Fig. 5, the critical mass ratios of the SR and NR cases become closer, but still, stability occurs for a smaller value of the mass ratio M/m in the NR case than in the SR one.
We have checked that with T-M mixing, still with a purely linear potential, the same pattern is observed.
We have repeated the study with a purely Coulombic interaction, that would ideally describe a bound state of very heavy quarks, provided that none of these quarks decay weakly too fast. This corresponds to the Hamiltonian The study is restricted to a color T. The pattern is similar to the one observed for a linear interaction, but the critical mass ratio for tetraquark binding is significantly lowered. For a weak coupling g = 0.1, see Fig. 6, the NR and SR energies are very close, as expected. Their comparison is just a check of the consistency of our computation scheme.
The results corresponding to a somewhat larger coupling g = 0.5 are shown in Fig. 7. The tetraquark energy is moderately lowered by relativistic effects, much less than the threshold energy. Hence the mass ratio required for stability is higher, and when stability is reached, the binding of the tetraquark is smaller. We now do the calculation with the aforementioned AL1 model (15), including the mixing of the T and M color components. The doubly-heavy tetraquarks have been already studied by several authors with this model [21,22], with the good surprise that a tiny binding is obtained for ccūd, and, of course, a more pronounced binding for the heavier analogs where one or two c quarks are replaced by b. This is shown in Fig. 8: the curve labeled AL1-NR crosses the threshold th-NR before M = m c , if one fixes beforehand m = m q . Two remarks here are in order: • a value of the variational energy E > E th simply means that binding is not found. The lowest energy of the 4-body Hamiltonian is E th . This is confirmed by the observation that the color content of the variational wave function tends to 1/3 T and 2/3 M, corresponding to a meson-meson decomposition.
• the energy with a frozen M color wave function, not shown, is higher than for T, except near M = m, as discussed, e.g., in [20]. The SR analogs are also shown in Fig. 8, and the tetraquark is bound for any value of M/m. This is due to an unrealistic strength for the hyperfine component, which is attractive for theūd pair, as discussed in (15).
Using the more realistic AL1N potential, we obtain the pattern shown in Fig. 9 for the tetraquark vs. its threshold. It can be seen that once again, binding is obtained for a larger value of the mass ratio M/m than in the NR case.

B. All-heavy tetraquarks
We now consider the case of tetraquarks with two heavy quarks and two heavy antiquarks. There is a flurry of calculations, in particular following the announcement by LHCb of the peak in the J/ψ-J/ψ distribution, and its interpretation as a cccc resonance. One of the questions raised by the LHCb discovery is whether there exist bound states of cccc or bbbb. Most of calculations are done in a simplified scheme with a diquark and an antidiquark. As shown elsewhere [20], this is not a good approxima- tion to the standard quark model, 2 and hence we do not include the corresponding papers in our discussion. Using a standard (color-dependent, pairwise) potential and an expansion on an harmonic-oscillator basis, Llyod and Vary [23] found a bound cccc bound state, but their result was not confirmed by other authors. 3 In [20,24] and references therein, an explanation is given on why in the chromoelectric limit QQQQ is not bound, while the electric analog e + e + e − e − is stable against dissociation into two e + e − atoms.
For instance, if we adopt the above AL1 potential, the lowest bbbb state is estimated at 18.872 GeV, above the threshold 18.848 GeV, and this energy would decrease toward this threshold if the variational expansion were pushed further. The unbound character is reinforced by the observation that the color content is nearly exactly 1/3 for T and 2/3 for M, which corresponds to a singletsinglet content. The analog for SR form of kinetic energy is 18.792 GeV for the variational estimate, again with 33% T and 67% M, above the threshold at 18.772 GeV. So the relativistic form of kinetic energy does not rescue the binding of QQQQ.

V. OUTLOOK
Clearly, estimating the mass and properties of ordinary and exotic hadrons requires sophisticated tools where long-range and short-range aspects of the dynamics are well accounted for. Simple models are nevertheless useful to probe some mechanisms, provided these models are solved carefully.
In this article, we revisited the effect of relativistic kinematics in the quark model by studying how the results are modified by replacing the NR form of kinetic operator p 2 /(2 m) by p 2 + m 2 − m. In the meson sector, most relativistic effects can be absorbed by a tuning of the parameters. Such change is illustrated in Eq. (16). The readjustment of the constituent masses is not very dramatic, and could be supplemented by minor changes of the static potential. More important is the necessary modification of the chromomagnetic term, because the relativistic corrections influence the short-range correlations.
The same is true for baryons. The study becomes more delicate for tetraquarks. For a given potential, both the threshold energy and the multiquark energy are lowered by the change of kinetic energy. It is observed that the effect is more pronounced for the former, so that the binding energy decreases. In particular, for QQqq configurations in a given potential, the mass ratio M/m at which the system becomes stable is larger for the semi-relativistic models than for the non-relativistic ones. Besides, fullyheavy tetraquarks remain above threshold if relativistic kinematics is used.
In other words, the threshold "benefits" more from the relativistic corrections than the collective configuration. This is similar to what is observed for some symmetry breakings. For instance, if one starts from a QQqq system with masses {M, M, m, m}, and breaks the particle identity, then one usually observes that the system {M, M , m, m } with M > M and m > m is less bound with respect to the {M , m } + {M, m} threshold than the more symmetric system with respect to its threshold 2 {M, m}.