Fully-strange tetraquark $ss\bar{s}\bar{s}$ spectrum and possible experimental evidence

In this work we construct 36 tetraquark configurations for the $1S$-, $1P$-, and $2S$-wave states, and make a prediction of the mass spectrum for the tetraquark $ss\bar{s}\bar{s}$ system in the framework of a nonrelativistic potential quark model without the diquark-antidiquark approximation. The model parameters are well determined by our previous study of the strangeonium spectrum. We find that the resonances $f_0(2200)$ and $f_2(2340)$ may favor the assignments of ground states $T_{(ss\bar{s}\bar{s})0^{++}}(2218)$ and $T_{(ss\bar{s}\bar{s})2^{++}}(2378)$, respectively, and the newly observed $X(2500)$ at BESIII may be a candidate of the lowest mass $1P$-wave $0^{-+}$ state $T_{(ss\bar{s}\bar{s})0^{-+}}(2481)$. Signals for the other $0^{++}$ ground state $T_{(ss\bar{s}\bar{s})0^{++}}(2440)$ may also have been observed in the $\phi\phi$ invariant mass spectrum in $J/\psi\to\gamma\phi\phi$ at BESIII. The masses of the $J^{PC}=1^{--}$ $T_{ss\bar{s}\bar{s}}$ states are predicted to be in the range of $\sim 2.44-2.99$ GeV, which indicates that the $\phi(2170)$ resonance may not be a good candidate of the $T_{ss\bar{s}\bar{s}}$ state. This study may provide a useful guidance for searching for the $T_{ss\bar{s}\bar{s}}$ states in experiments.


I. INTRODUCTION
From the Review of Particle Physics (RPP) of Particle Data Group [1], above the mass range of 2.0 GeV one can see that there are several unflavored qq isoscaler states, such as f 0 (2200), f 2 (2150), f 2 (2300), f 2 (2340) etc., dominantly decaying into φφ, ηη, and/or KK final states. The decay modes indicate that these states might be good candidates for conventional ss meson resonances. Recently, we carried out a systematical study of the mass spectrum and strong decay properties of the ss system in Ref. [2]. It shows that these states cannot be easily accommodated by the conventional ss meson spectrum. While they may be candidates for tetraquark ssss (T ssss ) states, it is easy to understand that they can fall apart into φφ and ηη final states through quark rearrangements, or easily decay into KK final states through a pair of ss annihilations and then a pair of light quark creations. The mass analysis with the relativistic quark model in Ref. [3] supports the f 0 (2200), and f 2 (2340) to be assigned as the T ssss ground states with 0 ++ and 2 ++ , respectively. However, a relativized quark model calculation [4] only favors f 2 (2300) to be a T ssss state.
Some other candidates of the T ssss states from experiment are also suggested in the literature. For example the vector meson resonance φ(2170) listed in RPP [1] is suggested to be a 1 −− T ssss state based on the mass analysis of QCD sum rules [5][6][7][8][9], and flux-tube model [10]. The newly observed X(2239) resonance in the e + e − → K + K − process at BE-SIII [11] is suggested to be a candidate of the lowest 1 −− T ssss state in a relativized quark model [4]. Moreover, the newly observed resonances X(2500) observed in J/ψ → γφφ [12] and X(2060) observed in J/ψ → φηη ′ [13] at BESIII are suggested * E-mail: zhongxh@hunnu.edu.cn † E-mail: zhaoq@ihep.ac.cn to be 0 −+ and 1 +− T ssss states, respectively, according to the QCD sum rule studies [14,15]. The assignment of X(2500) is consistent with that in Ref. [4].
With the recent experimental progresses more quantitative studies on the T ssss states can be carried out and their evidences can also be searched for in experiments. Very recently, the LHCb Collaboration reported their results on the observations of tetraquark cccc (T cccc ) states [16]. A broad structure above the J/ψJ/ψ threshold ranging from 6.2 to 6.8 GeV and a narrower structure T cccc (6900) are observed with more than 5 σ of significance level. There are also some vague structure around 7.2 GeV to be confirmed. These observations could be evidences for genuine T cccc states [17][18][19][20][21].
The observations of the T cccc states above the J/ψJ/ψ threshold at LHCb may provide an important clue for the underlying dynamics for the cccc system. In particular, the narrowness of T cccc (6900) suggests that the there should be more profound mechanism that "slows down" the fall-apart decays of such a tetraquark system. Although this may be related to the properties of the static potential of heavy quark systems, more direct evidences are still needed to disentangle the dynamical features between the heavy and light flavor systems. As an analogy of the T cccc system, there might exist stable T ssss states above the φφ threshold, and can likely be observed in the di-φ mass spectrum. On the other hand, flavor mixings could be important for the light flavor systems and pure ssss states may not exist. To answer such questions, systematic calculations of the ssss system should be carried out. The BESIII experiments can provide a large data sample for the search of the T ssss states in J/ψ and ψ(2S ) decays. In theory, although there have been some predictions of the T ssss spectrum within the quark model [3,4,10] and QCD sum rules [17][18][19][20], most of the studies focus on some special states in a diquark-antidiquark picture. About the status of the tetraquark states, some recent review works can be referenced [22,23]. In this study we intend to provide a systematical calculation of the mass spectrum of the 1S , 1P and 2S -wave T ssss states without the diquark-antidiquark approximation in a nonrelativistic potential quark model (NRPQM).
The NRPQM is based on the Hamiltonian proposed by the Cornell model [24], which contains a linear confinement and a one-gluon-exchange (OGE) potential for quarkquark and quark-antiquark interactions. With the NRPQM, we have successfully described the ss, cc, and bb meson spectra [2,25,26], and sss, ccc and bbb baryon spectra [27,28]. Furthermore, we adopted the NRPQM for the study of both 1S and 1P-wave all-heavy tetraquark states with a Gaussian expansion method [21,29]. In this work we continue to extend this method to study the T ssss spectrum by constructing the full tetraquark configurations without the diquark-antidiquark approximation. With the parameters determined in our study of the ss spectrum [2], we obtain a relatively reliable prediction of the mass spectrum for 36 T ssss states, i.e., 4 1S -wave ground states, 20 1P-wave orbital excitations, and 12 2S -wave radial excitations.
The paper is organized as follows: a brief introduction to the tetraquark spectrum is given in Sec. II. In Sec. III, the numerical results and discussions are presented. A short summary is given in Sec. IV.

A. Hamiltonian
We adopt a NRPQM to calculate the mass spectrum of the ssss system. In this model the Hamiltonian is given by where m i and T i stand for the constituent quark mass and kinetic energy of the ith quark, respectively; T G stands for the center-of-mass (c.m.) kinetic energy of the tetraquark system; r i j ≡ |r i − r j | is the distance between the ith and jth quark; and V i j (r i j ) stands for the effective potential between them. In this work the V i j (r i j ) adopts a widely used form [24][25][26][30][31][32][33][34][35][36]: where the confinement potential adopts the standard form of the Cornell potential [24], which includes the spinindependent linear confinement potential V Lin i j (r i j ) ∝ r i j and Coulomb-like potential V Coul i j (r i j ) ∝ 1/r i j : The constant C 0 stands for the zero point energy. While the spin-dependent potential V sd i j (r i j ) is the sum of the spin-spin contact hyperfine potential V S S i j , the spin-orbit potential V S S i j , and the tensor term V T i j : with In the above equations, S i stands for the spin of the ith quark, and L i j stands for the relative orbital angular momentum between the ith and jth quark. If the interaction occurs between two quarks or antiquarks, the λ i · λ j operator is defined as λ i · λ j ≡ 8 a=1 λ a i λ a j , while if the interaction occurs between a quark and an antiquark, the λ i · λ j operator is defined as where λ a * is the complex conjugate of the Gell-Mann matrix λ a . The parameters b i j and α i j denote the strength of the confinement and strong coupling of the one-gluon-exchange potential, respectively.
The five parameters m s , α ss , σ ss , b ss , and C 0 have been determined by fitting the mass spectrum of the strangeonium in our previous work [2]. The quark model parameters adopted in this work are collected in the Table I.

B. Configurations classified in the quark model
To calculate the spectroscopy of a qqqq (q ∈ {s, c, b}) system, first we construct the configurations in the product space of flavor, color, spin, and spatial parts.
In the color space, there are two color-singlet bases |66 c and |33 c , their wave functions are given by |66 c = 1 2 √ 6 (rb + br)(br +rb) + (gr + rg)(ḡr +rḡ) In the spin space, there are six spin bases, which are denoted by χ S 12 S 34 S . Where S 12 stands for the spin quantum numbers for the diquark (q 1 q 2 ) (or antidiquark (q 1q2 )), while S 34 stands for the spin quantum number for the antidiquark (q 3q4 ) (or diquark (q 3 q 4 ) ). S is the total spin quantum number of the tetraquark qqqq system. The spin wave functions χ S 12 S 34 S S z with  0.135 C 0 (GeV) −0.519 a determined S z (S z stands for the third component of the total spin S) can be explicitly expressed as follows: In the spatial space, we define the relative Jacobi coordinates with the single-partial coordinates r i (i = 1, 2, 3, 4): R ≡ m 1 r 1 + m 2 r 2 + m 3 r 3 + m 4 r 4 Note that ξ 1 and ξ 2 stand for the relative Jacobi coordinates between two quarks q 1 and q 2 (or antiquarksq 1 andq 2 ), and two antiquarksq 3 andq 4 (or quarks q 3 and q 4 ), respectively. While ξ 3 stands for the relative Jacobi coordinate between diquark qq and anti-diquarkqq. Using the above Jacobi coordinates, it is easy to obtain basis functions that have well-defined symmetry under permutations of the pairs (12) and (34) [37].
In the Jacobi coordinate system, the spatial wave function Ψ NLM (ξ 1 , ξ 2 , ξ 3 , R) for a qqqq system with principal quantum number N and orbital angular momentum quantum numbers LM may be expressed as the linear combination of Φ(R)ψ α 1 (ξ 1 )ψ α 2 (ξ 2 )ψ α 3 (ξ 3 ): (20) where C α 1 ,α 2 ,α 3 stands for the combination coefficients, Φ(R) is the center-of-mass (c.m.) motion wave function. In the quantum number set α i ≡ {n ξ i , l ξ i , m ξ i }, n ξ i is the principal quantum number, l ξ i is the angular momentum, and m ξ i is its third component projection. The wave functions ψ α i (ξ i ), which account for the relative motions, can be written as where Y l ξ i m ξ i (ξ i ) is the spherical harmonic function, and R n ξ i l ξ i (ξ i ) is the radial part. It is seen that for an excited state, there are three spatial excitation modes corresponding to three independent internal wave functions ψ α i (ξ i ) (i = 1, 2, 3), which are denoted as ξ 1 , ξ 2 , and ξ 3 , respectively, in the present work. One point should be emphasized that considering the fact that the ssss system is composed of equal mass constituent quarks and antiquarks, we adopt a single set of Jacobi coordinates in this study as an approximation. In fact, the four-body wave function describing a scalar ssss state contains a small contribution of internal angular momentum. This contribution is neglected in our calculations. To precisely treat an N-body system, one can involve several different sets of Jacobi coordinates as those done in Refs. [38][39][40][41][42][43]; or adopt a single set of Jacobi coordinates X = (ξ i , ξ 2 , ..., ξ N−1 ) with non diagonal Gaussians e −XAX T as those done in Refs. [44][45][46][47], where A is a symmetric matrix. Taking into account the Pauli principle and color confinement for the four-quark system qqqq, we have 4 configurations for 1S -wave ground states, 20 configurations for the 1P-wave orbital excitations, and 12 configurations for the 2Swave radial excitations. The spin-parity quantum numbers, notations, and wave functions for these configurations are presented in Table II. With the wave functions for all the configurations, the mass matrix elements of the Hamiltonian can be worked out.
To work out the matrix elements in the coordinate space, we expand the radial part R n ξ i l ξ i (ξ i ) with a series of harmonic oscillator functions [21,27]: with It should be pointed out that if there are no radial excitations, the expansion method with harmonic oscillator wave functions are just the same as the Gaussian expansion method adopted in the literature [38,39].
For an ssss system, if we ensure that the spatial wave function with Jacobi coordinates can transform into the single particle coordinate system, the harmonic oscillator frequencies ω ξ i ℓ (i = 1, 2, 3) can be related to the harmonic oscillator stiffness factor K ℓ with ω ξ 1 ℓ = 2K ℓ /µ ξ 1 , ω ξ 2 ℓ = 2K ℓ /µ ξ 2 , and ω ξ 3 ℓ = 4K ℓ /µ ξ 3 . Considering the reduced masses µ ξ 1 = µ ξ 2 = m s /2, µ ξ 3 = m s for T (ssss) , one has ω ξ 1 ℓ = ω ξ 2 ℓ = ω ξ 3 ℓ = √ 4K ℓ /m s . It indicates that the harmonic oscillator frequencies ω ξ i ℓ for T (ssss) are not independent. According to the relation ω ξ 1 ℓ = ω ξ 2 ℓ = ω ξ 3 ℓ = ω ℓ , the expansion Then we introduce oscillator length parameters d ℓ that can be related to the harmonic oscillator frequencies ω ℓ with 1/d 2 ℓ = m s ω ℓ . Following the method of Refs. [38,39], we let the d ℓ parameters form a geometric progression where n is the number of harmonic oscillator functions, and a is the ratio coefficient. There are three parameters {d 1 , d n , n} to be determined through the variation method. It is found that with the parameter set {0.085 fm, 3.399 fm, 15} for the ssss system, we can obtain stable solutions. The numerical results should be independent of the parameter d 1 . To confirm this point, as done in the literature [40][41][42] we scale the parameter d 1 of the basis functions as d 1 → αd 1 . The mass of a T ssss state should be stable at a resonance energy insensitive to the scaling parameter α. As an example, we plot the masses of 12 2S -wave T ssss configurations as a function of the scaling factor α in Fig. 1. It is found that the numerical results are nearly independent of the scaling factor α. The stabilization of other states predicted in this work has also been examined by the same method.
With the mass matrix elements ready for each configuration, the mass of the tetraquark configuration and its spacial wave function can be determined by solving a generalized eigenvalue problem. The details can be found in our previous works [27,29]. Finally, the physical states can be obtained by diagonalizing the mass matrix of different configurations with the same J PC numbers.

III. RESULTS AND DISCUSSIONS
Our predictions of the T ssss mass spectrum with the HOEM are given in Table III, where the components of different configurations for a physical state can be seen. For example, the two 0 ++ ground states are mixing states between two different configurations 1 S 0 ++ (66) c and 1 S 0 ++ (33) c due to a strong contribution of the confinement potential to the non-diagonal elements. To see the contributions from each part of the Hamiltonian to the mass of different configurations, we also present our results in Table IV. It is found that both the kinetic energy term T and the linear confinement potential term V Lin contribute a large positive value to the mass, while the Coulomb type potential V Coul has a large cancelation with these two terms. The spin-spin interaction term V S S , the tensor potential term V T , and/or the spin-orbit interaction term V LS have also sizeable contributions to some configurations. Thus, as a reliable calculation, both the spin-independent and spindependent potentials should be reasonably included for the ssss system. For clarity, our predicted T ssss spectrum is plotted in Fig. 2.

A. Discussions of the numerical method
Herein we discuss the differences of numerical results between the expansion method with the harmonic oscillator wave functions (HOEM) used in present work and the Gaussian expansion method (GEM) often adopted in the literature. For the 1S -, 1P-wave T (ssss) states, etc., there are no radial excitations. Thus, the GEM is the same as the HOEM. For the first radial excited 2S -wave T (ssss) states, the HOEM is different from the GEM because the trail harmonic oscillator wave functions are different from the Gaussian functions.
To see the differences between the two expansion methods we also give our predictions of the 2S -wave T (ssss) states based on the GEM. It should be mentioned that by fully ex-   Table II. Then we cannot numerically work out the masses for the following states of 0 +− (2 1 S 0 +− (66) c (ξ 1 ,ξ 2 ) and 2 1 S 0 +− (33) c (ξ 1 ,ξ 2 ) ), 1 ++ (2 3 S 1 ++ (33) c (ξ 1 ,ξ 2 ) ), and 2 +− (2 5 S 2 +− (33) c (ξ 1 ,ξ 2 ) ) listed in Table II. To overcome this problem, the spatial wave functions containing the radial excitations are expanded with the Gaussian functions, while the spatial wave functions containing no excitations are adopted 2628 95 25 95 2712 12 25 12 2633 the single Gaussian function as an approximation. We have tested the single Gaussian approximation in the calculations of the ground 1S T (ssss) states, the numerical values are reasonably consistent with those calculated with a series of Gaussian functions. The differences of the numerical results between these two methods are about 10 MeV.
Our numerical results for the 2S -wave T (ssss) states with the GEM are listed in Table IV and Table V. From Table IV, it is found that the numerical values for the 0 +− configuration 2 1 S 0 +− (66) c (ξ 1 ,ξ 2 ) and 0 ++ configuration 2 1 S 0 ++ (66) c (ξ 3 ) cal-culated with the HOEM are significantly different from those obtained with the GEM. For these two configurations, the predicted mass differences by the HOEM and GEM can reach up to ∼ 70 MeV. However, for the other 2S -wave T (ssss) configurations the numerical values of these two methods are comparable with each other. The differences of the predicted masses between these two methods are about 10 − 20 MeV. It should be mentioned that the Coulomb type potential V Coul for the 2S -wave states seems to be sensitive to the numerical methods as shown in Table IV. The average contributions of each part of the Hamiltonian to the ssss configurations with the HOEM. T stands for the contribution of the kinetic energy term. V Lin and V Coul stand for the contributions from the linear confinement potential and Coulomb type potential, respectively. V S S , V T , and V LS stand for the contributions from the spin-spin interaction term, the tensor potential term, and the spin-orbit interaction term, respectively. The second number in every column is calculated with the GEM.
In brief, most of the predictions are consistent with each other between the HOEM and GEM. The uncertainties from the numerical methods do not change our main predictions of the T (ssss) spectrum, Although some numerical results for the J PC = 0 +− and 0 ++ 2S states show a significant numerical method dependence (See Fig. 2), the GEM may give a slightly more accurate numerical result based on our tests of the charmonium spectrum. In the following, our discussions of the 2S states are based on the GEM calculations.

0 ++ states
In the 1S -wave mutiplets, the two 0 ++ ground states include T (ssss)0 ++ (2218) and T (ssss)0 ++ (2440). Their mass splitting reaches up to about 200 MeV. These two states have a strong mixing between the two color structures |66 c and |33 c . Their masses are much larger than the mass threshold of φφ. Thus, they may easily decay into φφ pair through quark rearrangements. The mass of the lowest 0 ++ T ssss in our model is close to the prediction of 2203 MeV in the relativistic diquark-antidiquark model [3]. However, it turns out to be much higher than the predicted value 1716 MeV by the relativized quark model with a diquark-antidiquark approximation [4]. There might be some crucial dynamics missing in the diquark-antidiquark approximation. As a test of the diquark-antidiquark approximation we adopt the approximation as done in Ref. [4] and calculate the mass of the 0 ++ T ssss state with the same potential model parameters. We obtain a mass of 1758 MeV, which is comparable with the prediction of Ref. [4], but is obviously smaller than the results without the diquark-antidiquark approximation.
It should be mentioned that f 0 (2200) is listed in RPP [1] as a well-established state. It has been seen in the K 0 s K 0 s , K + K − and ηη, and may be assigned to T (ssss)0 ++ (2218). Some qualitative features can be expected: (i) The 0 ++ T ssss state can decay into ηη, η ′ η ′ , and ηη ′ through quark rearrangements via the ss component in the η and η ′ mesons. An approximate branching ratio fraction can be examined: BR(ηη) : BR(η ′ η ′ ) : BR(ηη ′ ) ≃ sin 4 α P : cos 4 α P : 2 sin 2 α P cos 2 α P ≃ 0.24 : 0.25 : 0.50, with α P ≡ arctan √ 2 + θ P ≃ 44.7 • and without including the phase space factors. (ii) The 0 ++ states may also easily decay into K 0 s K 0 s and K + K − final states through annihilating a pair of ss and creating a pair of light qq. (iii) It is interesting to note that no conventional 0 ++ ss states are predicted around 2.2 GeV in most literatures [2].

2 ++ states
There is only one 2 ++ state T (ssss)2 ++ (2378) in the 1S -wave states. This state lies between the two 0 ++ ground states, and has a pure |33 c color structure. T (ssss)2 ++ (2378) may have large decay rates into the φφ, ηη and η ′ η ′ final states through quark rearrangements, and/or into K ( * )K( * ) final states through the annihilation of ss and creation of a pair of nonstrange qq. It should be mentioned that with the diquark-antidiquark approximation, the mass of the 2 ++ state is predicted to be 2192 MeV, which is about 200 MeV lower than the four-body calculation results.
The f 2 (2340) resonance listed in RPP [1] may be assigned to T (ssss)2 ++ (2378). Besides the measured mass 2345 +50 −40 MeV, the observed decay modes φφ and ηη are consistent with the expectation of the tetraquark scenario. On the other hand, as a conventional ss state the f 2 (2340) cannot be easily accommodated by the quark model expectation [2]. The relativistic quark model calculation of Ref. [3] also supports the f 2 (2340) to be assigned as the T ssss ground state with 2 ++ . To confirm this assignment, the other main decay modes of T (ssss)2 ++ (2378) such as ηη ′ , η ′ η ′ , K ( * )K( * ) should be investigated in experiment.

1 +− states
In the 1S -wave multiplets T (ssss)1 +− (2323) is the only state with C = −1, and has a pure |33 c color structure. Its mass is about 100 MeV larger than the lowest 1S -wave state T (ssss)0 ++ (2218). Its mass is about 200 ∼ 300 MeV larger than that predicted by the QCD sum rules [15] and the relativized quark model [4] in the diquark picture. The mass of the 1 +− state may be notably underestimated in the diquark picture. and BESIII [59], which could be signals of the 1 −− ssss tetraquark states [60]. For the heavier states T (ssss)1 −− (2766) and T (ssss)1 −− (2984), they can also decay into ΞΞ baryon pair through a qq pair production in vacuum. Thus, experimental search for these states in e + e − → ΞΞ should be very interesting.

3 −− state
There is only one 3 −− state T (ssss)3 −− (2719) predicted in the NRPQM. This state has a pure color structure |33 c , and also a pure orbital excitation between the diquark (ss) and antidiquark (ss). Our predicted mass is about 60 MeV larger than that predicted by the relativized quark model [4] with a diquark approximation. The 3 −− states may easily decay into φ f ′ 2 (1525) in an S wave by the quark rearrangements. Since it has a high spin, it may be produced relatively easier in pp or pp collisions.

IV. SUMMARY
In this work we calculate the mass spectra for the 1S , 1P and 2S -wave T ssss states in a nonrelativistic potential quark model without the often-adopted diquark-antidiquark approximation. The 1S -wave ground states lie in the mass range of ∼ 2.21 − 2.44 GeV, while the 1P-and 2S -wave states scatter in a rather wide mass range of ∼ 2.44 − 2.99 GeV. For the 2Swave states, except for the highest state T (ssss)0 ++ (3155) all the other states lie in a relatively narrow range of ∼ 2.78 − 2.98 GeV. We find that most of the physical states are mixed states with different configurations.
For the ssss system it shows that both the kinetic energy T and the linear confinement potential V Lin contribute a large positive value to the mass, while the Coulomb type potential V Coul has a large cancellation with the these two terms. The spin-spin interaction V S S , tensor potential V T , and/or the spin-orbit interaction term V LS also have sizeable contributions to some configurations. Some T ssss states may have shown hints in experiment. For instance, the observed decay modes and masses of f 0 (2200) and f 2 (2340) listed in RPP [1] could be good candidates for the ground states T (ssss)0 ++ (2218) and T (ssss)2 ++ (2378), respectively. The newly observed X(2500) at BESIII may be a candidate for the lowest mass 1P-wave 0 −+ state T (ssss)0 −+ (2481). Another 0 ++ ground state T (ssss)0 ++ (2440) may have shown signals in the φφ channel at BESIII [12,48]. Our calculation shows that φ(2170) may not favor a vector state of T ssss , because of the much higher mass obtained in our model.
It should be stressed that as a flavor partner of T cccc , the T ssss system may have very different dynamic features that need further studies. One crucial point is that the strange quark is rather light and the light flavor mixing effects could become non-negligible. It suggests that strong couplings between T ssss and open strangeness channels could be sizeable. As a consequence, mixings between T ssss and T sqsq would be inevitable. For an S -wave strong coupling, it may also lead to configuration mixings which can be interpreted as hadronic molecules for a near-threshold structure. In such a sense, this study can set up a reference on the basis of orthogonal states. More elaborated dynamics can be investigated by including the hadron interactions in the Hamiltonian. For states with exotic quantum numbers, experimental searches for their signals can be carried out at BESIII and Belle-II.