Phenomenological mass model for exotic hadrons and predictions for masses of non-strange dibaryons as hexaquarks

We investigate the mass spectra of exotic hadrons known as hexaquarks in the form of dibaryons. We use a phenomenological model based on an extended version of the G\"ursey-Radicati mass formula for hadrons to include non-charmed baryons, charmed baryons, and non-strange dibaryons to be able to predict masses of potential dibaryon states. We perform six numerical fits of this model to input data for three different sets of masses of baryons and dibaryons. We find that the model can fit some of the data sets well, especially the sets including charmed baryons and non-strange dibaryons, and observe that the predicted mass of one of the dibaryons is close to the measured mass of the observed hexaquark candidate $d^*(2380)$ reported by the WASA-at-COSY experiment. The predicted mass of the deuteron is slightly larger than its measured mass. Finally, for the data sets including charmed baryon and non-strange dibaryon masses, we find that the predicted masses of potential dibaryon states are all in the range from 1900 MeV to 3700 MeV.


I. INTRODUCTION
Quarks are subatomic particles that were proposed independently by Gell-Mann [1] and Zweig [2] in 1964, and verified by experiments a few years later. They are never found isolated in Nature, but rather exist in composite structures called hadrons. Observed hadrons normally consist of quark-antiquark pairs or collections of three quarks (or three antiquarks). However, apart from the constraint of hadrons having zero net color, the quark model does not impose any limitation on the number of quarks that can constitute a hadron. The existence of so-called exotic hadrons consisting of more than three quarks was also proposed by Gell-Mann [1] in 1964, and such hadrons have been observed by various particle accelerator experiments in recent years [3,4]. The most commonly discussed exotic hadrons are those composed of four, five or six quarks, known as tetraquarks, pentaquarks and hexaquarks, respectively (see e.g. Refs. [5,6]).
The proposition of six quarks combined into a single structure, although not necessarily called a hexaquark, was initially proposed by Dyson and Xuong [7] in 1964, as a so-called dibaryon state. They predicted non-strange S-wave dibaryon states D IJ , where I and J denote isospin and spin, respectively. The proposal was that the deuteron, D 01 (N N ), and a virtual state, D 10 (N N ), were contained in the 10 and 27 representations of the group SU(3) that is part of the decomposition SU(3) × SU(2) ⊂ SU (6). They also predicted four additional states, D 03 (∆∆), D 30 (∆∆), D 12 (N ∆), and D 21 (N ∆), based on group theoretical symmetry arguments. Yet another hexaquark was proposed by Jaffe in 1976 [8]. It was assumed to have the quark content uuddss and given the name dihyperon or the H-particle. Note that the name hyperon is nowadays normally used for baryons containing one or more strange quarks. The stability of the H-particle and many other dibaryon candidates was later studied further in several works by Leandri and Silvestre-Brac, see Ref. [9] and references therein. Also, in Refs. [10][11][12], the ideas of hiddencolor dibaryons and hidden-charm hexaquarks were discussed. For the pure hexaquark sector, simply meaning six quarks and no antiquarks, a large systematic study of dibaryon candidates was carried out [9] by means of a schematic chromomagnetic model (see e.g. Refs. [13][14][15]). Results pointed to dibaryon candidates most favorable for stability being composed of three different pairs of identical quarks and these states were also studied in Ref. [9]. Despite the possibility of such resonances not being excluded, the investigations entailed that stable dibaryon states are not favored -there always exists a two-baryon channel with lower energy. However, the search for hexaquarks revived in 2014, when the WASAat-COSY collaboration revealed observations of a narrow resonance-like structure around 2380 MeV by studying the double-pionic fusion channels pn → dπ 0 π 0 and pn → dπ + π − . This structure was named d * (2380) and data suggested it having quantum numbers I(J P ) = 0(3 + ) [16]. A serious candidate for the d * (2380) has later been thought to be the D 03 state, proposed fifty years earlier, having the proper quantum numbers and a predicted mass of about 2350 MeV. More recently, Huang et al. [17] investigated the six dibaryon candidates proposed by Dyson and Xuong, using dynamical calculations and a simple mass formula previously applied to baryons and pentaquarks. In particular, they investigated the possibility of the D 21 as a bound-state dibaryon candidate arXiv:2106.11978v2 [hep-ph] 22 Nov 2021 and concluded that models that obtain the experimental d * (2380) as the D 03 are not compatible with also obtaining the D 21 as a bound state. Hexaquarks have also been studied in chiral quark models [18,19], using QCD sum rules [20][21][22], and by constructing wave functions with effective potentials [23]. Finally, note that alternative interpretations of the observed d * (2380) structure have been proposed, see e.g. Refs. [24,25].
In this work, we will consider the simple mass formula for hadrons discussed in Ref. [17], which could predict the masses of hexaquarks using numerical fits to measured masses of various baryons. In general, we will use an extended version of this formula, fit it to larger data sets of baryons that have not been considered previously, and predict masses of the six non-strange S-wave baryon states D IJ that could constitute hexaquarks.
This work is organized as follows. In Sec. II, we present the phenomenological mass model based on the so-called extended Gürsey-Radicati mass formula for hadrons that can predict masses of hexaquarks, and describe the numerical fitting procedure. Then, in Sec. III, we introduce three different data sets for measured masses of baryons, perform six numerical fits of this model, and state and discuss the results of these fits in detail. The predicted masses of hexaquarks are also presented. Finally, in Sec. IV, we summarize our main results and conclude based on our stated results.

II. MODEL AND FITTING PROCEDURE
A. Extended Gürsey-Radicati model Already in 1964, Gürsey and Radicati proposed a simple mass formula for baryons, based on the breaking of SU(6) ⊃ SU(2) s × SU(3) f spin-flavor symmetry [26]. As pointed out in Ref. [27], the Gürsey-Radicati (GR) mass formula describes quite well the way in which symmetry breaking affects the mass spectrum of baryons, despite its simplicity. In the original work of Gürsey and Radicati, the mass formula is given by (up to differences in notation) where J, Y , and I denote spin, hypercharge, and isospin, respectively, and M 0 , A, B, and C are phenomenological model parameters. These parameters are obtained from fits to baryon data. An extension to Eq. (1) was introduced and studied in Ref. [28], and further investigated in Ref. [29], where applications to pentaquark systems are considered. The extended version of the GR mass formula reads [28] where C 2 is the eigenvalue of the SU(3) f Casimir operator, N c is the number of constituent charm quarks in the considered baryon, and D and E are two additional model parameters. In Refs. [28,29], the parameters are fixed using the baryon spectrum and pentaquark masses are then calculated under the assumption that these parameter values are universally valid beyond baryons. Furthermore, the extended GR formula (2) is used in the context of hexaquarks in Ref. [17], when studying the possible dibaryon candidates D IJ suggested by Dyson and Xuong in 1964. For these dibaryon states, a simple formula for the eigenvalue of the Casimir operator is presented as [7] C 2 = 12 + 2I(I + 1) .
B. Fitting procedure Now, we continue to the fitting procedure of baryon data. Following the approach of Huang et al. [17], we use the extended GR mass formula (2) to calculate the mass spectrum of dibaryons previously predicted by Dyson and Xuong [7], imposing universality of model parameters.
The six model parameters are fitted to experimental data on hadron masses, retrieved from the Particle Data Group (PDG) [30] of 2020. In Ref. [17], the data set used consists of eight non-charmed baryons and six non-strange dibaryons, where of the latter only two have experimentally measured and verified masses. To distinguish our study from theirs, we additionally include eight charmed baryons, and furthermore, we only use the two verified dibaryon candidates N N and d * (2380) for D 01 and D 03 , respectively. In Tab. I, the baryons along with their corresponding masses and quantum numbers are presented, and likewise in Tab. II, the corresponding properties of the dibaryon candidates are found.
These data are divided into three sets, for each of which separate fits and calculations will be carried out. The division of the data is as follows: Set I: Data for the eight non-charmed baryons in the upper half of Tab I.
Set II: Data for the eight non-charmed baryons and the eight charmed baryons in Tab I.
Set III: Data for the eight non-charmed baryons, the eight charmed baryons, and the two dibaryon candidates in Tabs. I and II.
The parameters for each data set are then determined by minimizing the χ 2 -function where M exp.
i are the tabulated experimental masses, σ exp.
i the corresponding experimental uncertainties, and M i the input masses for data sets I, II, and III, respectively. The minimization of the χ 2 -function (4) is performed using both the built-in function Minimize in Mathematica [31] and a numerical method called basin-hopping, which is an iterative technique for globally optimizing a scalar function of one or several variables. Basin-hopping is a two-phase method consisting of a stochastic global stepping algorithm and local minimization at each step [32], which exist as built-in functions in the SciPy library. For the minimization, several numerical methods can be used, but in this work the above-mentioned function in Mathematica and the so-called BFGS-algorithm are chosen.

A. Data sets
After having divided the data into three sets, we fit the parameters of the extended GR mass formula (2) with Eq. (3) in two ways for each set. Firstly, we use the experimental uncertainties of the hadron masses -these fits are called unprimed and denoted I, II, and III. Secondly, we do what is called a "1 % error fit", i.e. setting each uncertainty σ i to one percent of the corresponding experimentally determined mass. This can sometimes yield a better fit, although the precision of the results cannot be better than just 1 %. These fits are called primed and denoted I , II , and III .

B. Numerical fits and results
In Tab. III, resulting parameter values for each case are presented along with the value of the minimized χ 2function (4). We then use the GR formula (2) with these parameters together with quantum numbers from Tab. II, to predict the mass spectrum of the dibaryon candidates. The predictions for the dibaryon mass spectrum are presented in Tab. IV. For fits III and III , the experimental masses of D 01 and D 03 are used as inputs, hence they are not included as predictions, whereas the resulting predictions for the masses of the other four dibaryons, i.e. D 10 , D 30 , D 12 , and D 21 , are all found to be in the rather narrow spectrum of (1900, 2700) MeV. In fact, for fits II-III , it turns out that the predictions for the masses of all six dibaryons lay in the interval (1900, 2900) MeV, except for the predicted mass of D 30 for fits II and II that is slightly larger. We start by examining Tab. III, where the values of the minimized χ 2 -function and the values of the phenomenological model parameters for the six different fits are presented. One striking feature is the difference between the minimized values for the primed and the unprimed fits. The unprimed fits I, II, and III have χ 2 ∼ 10 3 − 10 5 , which is very large and indicative of a poor fit, while the primed fits I , II , and III have χ 2 ∼ 10 0 − 10 1 , suggesting that much better parameter values have been found. This is thought to arise mainly due to N and N N having significantly smaller experimental uncertainties (∼ 10 −6 MeV) than all other included baryons, making their contributions dominate the χ 2 -function and hence giving a larger discrepancy for the other data points. Behavior like this is precisely the motivation behind also performing a 1 % error fit in the first place. Concerning the parameter values, we note that results for the fits I and I deviate significantly from the rest of the fits. For instance, the value of the parameter M 0 is much smaller, about 280 MeV, compared to all other fits having M 0 > 950 MeV. This is unexpected, since M 0 in a sense should represent the total mass of the constituents. Moreover, the parameters A and D -coefficients for the spin and Casimir eigenvalue terms, respectively -are radically different for the first two fits. For fits I and I , A is negative, about −280 MeV, while for the rest of the fits, A is positive and of smaller magnitude, about 20 MeV for fits II and II and (40 − 50) MeV for fits III and III . The parameter D also differs widely, about 350 MeV for fits I and I , compared to (40 − 50) MeV for fits II and II and (10 − 20) MeV for fits III and III . Finally, the parameters B and C -coefficients for the hypercharge and isospin terms -obtain similar values for all six fits, i.e. B ∈ (−200, −190) MeV and C ∈ (33, 41) MeV. For reference, we note that the four fits that yield more comparable results (i.e. fits II, II , III, and III ) agree better with the results in Ref. [17] than the other two fits (i.e. fits I and I ), although their method for obtaining the parameter values is not mentioned.
We observe that the results of fits I and I , using only non-charmed baryons, deviate significantly from the re- maining four fits. Since the parameter values obtained for these fits were found to be far off, this is not entirely unexpected. Although most fits follow a general trend for the masses, with D 30 and D 21 being the two heaviest, both fits I and I show a much more extreme spread with practically all predicted values being rendered unrealistic. However, the remaining fits show a more reasonable behavior, regarding both the primed and the unprimed ones. One explanation for this could be that the heavier charmed baryons, contributing with larger weight to the χ 2 -function, result in more consistent parameter values. This might indicate that including charmed baryons is beneficial, provided that the result from data set II reproduces the known dibaryon candidates well. In particular, looking at fit II , we note that the predicted value for the mass of the D 03 is about 2324 MeV, which is, in fact, rather close to the experimental value of the mass for the d * (2380), and basically within our 1 % error margin for the primed fits. However, the prediction for the D 01 , around 2100 MeV, is too large compared to the experimental value of N N at 1876 MeV. Note that none of the fits were able to describe the deuteron to within even 200 MeV precision. One way to investigate the validity of the model is to look at how well it reproduces already known results and data. Since there are only two experimentally determined masses for dibaryon candidates to compare with, we instead look at the well-known baryon spectrum. Figures 1-3 show the pulls of the calculated baryon masses for the primed fits of the data sets I, II, and III, i.e. the deviation from the experimental values divided by the uncertainty. On a general level, we make one immediate observation. The pulls for all of the primed fits are of order unity for data sets I, II, and III, while the individual pulls of the unprimed fits can be as large as 10 3 (but not presented in Figs. 1-3). This agrees with the previously mentioned discrepancies for χ 2 , giving an overall better fit at the expense of losing precision for the baryons that have the smallest uncertainties. Another noticeable feature, perhaps not unexpected, is that the data sets  including charmed baryons and dibaryon candidates, in general, yield larger pulls for the non-charmed baryons.
One explanation for this is that the charmed baryons are heavier, due to the charm quark being more massive than the up, down, and strange quarks, and hence shift the spectra towards larger mass. The same logic applies to the dibaryon candidates, naturally being more massive than conventional baryons. One could therefore argue that the inclusion of charmed baryons contributes to making the parameter values more universal. Although it should be emphasized that the approach of using the GR mass formula with 1 % error is far from the most precise method conceivable, it is still evident that the basic features of the baryon spectrum are well described by this model. The primed fits give pulls in the range ±3 for all baryons as well as D 01 and D 03 , which is of the order one per mille.

D. Comparison with other works
Since the approach of this study is inspired by the work of Huang et al. in Ref. [17], it is reasonable to compare our results to theirs. When fixing the parameters of the GR mass formula, they use experimental masses of the same non-charmed baryons that we present in the upper half of Tab. I as well as the dibaryons listed in Tab. II. In addition to D 01 and D 03 with masses assumed to be 1876 MeV and 2380 MeV, respectively, they have also included masses of D 10 , D 12 , and D 21 , albeit followed by a question mark (presumably to signal that these masses are not experimentally verified). Since there is no clear reference to where the values originate, we have chosen to leave them out of our analysis. In any case, the parameter values obtained in their work are (in terms of our notation): M 0 1026.2 MeV, A 56.883 MeV, B −194.70 MeV, C 33.218 MeV, and D 9.4085 MeV [17]. We shall also note that they obtain a separate M 0 2091.9 MeV for dibaryons, whereas we have simply used 2M 0 instead. The values of these parameters are in decent agreement with our obtained values for fits II-III , i.e. in terms of sign and order of magnitude.
Regarding the mass spectra, our results are generally larger than those obtained in Ref. [17]. For the masses, they obtain 1883 MeV for D 10 , 2394 MeV for D 30 , 2168 MeV for D 12 , and 2182 MeV for D 21 . This is not too surprising, since, as we have mentioned, inclusion of heavier charmed baryons in the data sets naturally should shift the mass spectrum somewhat. However, our results do agree to the extent that D IJ is heavier than D JI for I > J in all cases, meaning that forming bound states of D 01 , D 03 , and D 12 would be preferred. They present a threshold for the D 21 being a bound state at about 2171 MeV, which is lower than the value they obtain using the GR approach and also lower than what we obtain. This would indicate that no such state is possible. Further, their more extensive dynamical calculations yield concordant results, after which it is concluded that a bound state D 21 cannot be obtained within these models, while at the same time accounting for the experimentally observed d * (2380).

IV. SUMMARY AND CONCLUSIONS
Inspired by other works, we have used an extension of the Gürsey-Radicati mass formula, presented in Ref. [28], to predict the mass spectrum for some nonstrange dibaryon candidates. Numerical fits of parameters have been performed based on three different data sets, using both experimental and 1 % uncertainties for each set, resulting in a total of six fits.
Using only non-charmed baryon masses as input, very large fluctuations in the dibaryon mass spectra have been observed. The other fits, including charmed baryon and non-strange dibaryon masses, have given a more coher-ent prediction with some results very close to the experimentally observed mass of the hexaquark candidate d * (2380). The best prediction is 2324 MeV for the mass of D 03 (∆∆), whereas 2100 MeV for the mass of D 01 (N N ). Hence, it is not possible to obtain the mass of the deuteron to a better precision than 200 MeV. In general, the fits that use experimental uncertainties have resulted in very large residual values, i.e. χ 2 ∼ 10 3 − 10 5 , while the 1 % error fits have given χ 2 ∼ 10 0 −10 1 . For all fits, including charmed baryon and non-strange dibaryon masses as input, all predicted masses of dibaryon states are in the region (1900, 3700) MeV.