Study of three-flavored heavy dibaryons using lattice QCD

We present results of the first lattice QCD calculation of three-flavored heavy dibaryons both in the flavor-symmetric and antisymmetric channels. These dibaryons have spin zero, and are constructed using various possible combinations of quark flavors with at least one of them as the charm or the bottom quark, i.e., namely, $H_c(cudcud), H_b(budbud), H_{bcs}(bcsbcs)$, $H_{csl}(cslcsl), H_{bsl}(bslbsl)$ and $H_{bcl}(bclbcl)$; $l\in u,d$. We compute the ground state masses of these dibaryons and the calculations are performed on three $N_f=2+1+1$ HISQ gauge ensembles of the MILC collaboration, with lattice spacings $a =$ 0.1207, 0.0888 and 0.0582 fm. A relativistic overlap action is employed for the valence light to charm quarks while a non-relativistic-QCD Hamiltonian with improved coefficients is used for the bottom quarks. Unlike the doubly heavy tetraquarks, one and two-flavored heavy dibaryons, for which lattice QCD calculations have predicted deeply bound strong-interactions-stable states, for these $H_c, H_b, H_{csl},H_{bsl}$ dibaryons we do not find any such deeply bound state. However, for $H_{bcs}$, our results indicate the presence of an energy level $29\pm 24$ MeV below the lowest two-baryon threshold, which could be relevant for its future experimental searches. Moreover, we find that the energy difference between the ground state of $H_{bcl}$ and its lowest threshold increases when $m_l>m_s$. Taken together, our findings indicate the possibility of the existence of the $H_{bcs}$ dibaryon while all other physical three-flavored dibaryons are much closer to their thresholds suggesting either they are weakly bound or unbound, resolving which requires further detail study. Our results also point that the binding of a dibaryon configuration becomes stronger with the increase of its valence quark masses which suggests an interesting aspect of strong interactions at multiple scales.


I. INTRODUCTION
Quantum chromodynamics (QCD), the theory of strong interactions of quarks and gluons, predicts a very rich energy spectra of hadronic states comprising various quark flavors from light to the bottom [1]. While most of the observed hadrons are classified as mesons and baryons within quark models, QCD also allows the existence of other bound state configurations of quarks (and antiquarks), which are generically known as exotic hadrons. Indeed the recent discovery of a large number of new subatomic particles, the so-called X,Y Z hadrons [2][3][4][5][6][7][8], including tetraquarks [9][10][11][12][13][14] and pentaquarks, [15,16] have confirmed the existence of a new class of subatomic particles in Nature. These discoveries naturally have created tremendous excitement to the field of hadron spectroscopy [1][2][3][4][5][6][7][8]. In terms of the number of valence quark content, these recent discoveries have so far been limited to four (tetra)-and five (penta)-quark states, and except the possible finding of a broad d * (2380) resonance [17] 1 no new six (hexa)-quark state has yet been discov-ered. It is therefore natural to investigate the existence of hadrons with six valence quark configurations within QCD, with the goal that QCD-predicted such states can guide in discovering them in future at high energy laboratories.
In our visible universe, so far deuteron is found to be the only stable six-quark bound state with a binding energy of about 2.2 MeV which has been modelled to be the result of a many body interactions of two nucleons [19,20]. The so-called H(udsuds)-dibaryon is another highly speculated compact six-quark bound state with strangeness S = −2, spin J = 0 and isospin I = 0. In the first calculation of H-dibaryon, a binding of 70 MeV was predicted using the MIT bag model [21]. Thereafter this state had seen a thorough investigation in the past four decades through various model studies [22][23][24][25]. However, till date, experimental searches have ruled out the existence of such a deeply bound state [26,27]. The H dibaryon has also been a subject of numerous lattice QCD calculations in recent years and the results have also ruled out a deeply bound state, and rather indicate an unbound or a very weakly bound state if at all it exists [28][29][30][31][32]. However, to confirm the existence of such a state through lattice QCD calculations, it is essential to perform chiral and continuum extrapolations to the physical limits, along with a detailed finite-volume amplitude analysis of the pole distribution in the scattering arXiv:2206.02942v2 [hep-lat] 4 Oct 2022 amplitude across the complex energy plane [33]. Light dibaryons have also been studied in other channels. For example, in the SU (3) f quark model, one expects attraction in dibaryon states without quark Pauli blocking and with attractive color-spin interaction from onegluon exchange, ΛΛ − ΣΣ − N Ξ (H-channel), ∆∆, and N Ω states [34,35]. Recent experimental and theoretical studies based on femtoscopy also suggest the possibility of a N Ω bound state in the J π = 2 + channel [36][37][38].
In the spirit of finding H-dibaryon, six-quark states in the isosinglet channel but with heavy quarks have also been explored in several model calculations with results not favoring a bound state [39][40][41][42]. On the contrary, recently in a first principles lattice QCD calculation it was found that strong-interactions-stable deuteron-like heavy dibaryons can exist if at least two of the quarks in a dibaryon have heavy flavors [43]. Lattice calculations have also been performed recently for single-flavored heavy dibaryons with charm as well as bottom flavors. While using HALQCD potential method the bindings for 1 S 0 Ω ccc -Ω ccc has been reported to be about −5 MeV [44], a direct calculation for Ω bbb -Ω bbb dibaryons finds a deep binding of about −90 MeV in the 1 S 0 channel [45].
Hence it is interesting to investigate the nature of quark mass dependence of H-like three-flavored dibaryons, particularly whether such a state is bound with one or more heavy quark content. Any clear indication on the existence of heavy H-like three-flavored dibaryon from lattice QCD calculations will be attractive for searching them at high energy laboratories. Moreover, a detailed study of quark mass dependence on the binding energy can further reveal the intriguing dynamics of heavy quarks in dibaryons which can further illuminate our knowledge of strong interactions at multiple scales.
Lattice QCD is an ideal tool for studying multi-hadron systems since in addition to being a first principles method, it is also possible to obtain quantitative results through lattice QCD at any quark masses, including at their unphysical values. As a result it provides a unique tool for systematic study of the quark mass dependence of binding energies, which otherwise is not possible to obtain. Moreover, with adequate computing resources it is possible to keep track of all the systematic uncertainties associated with such calculations. In recent years, beside the regular single-hadron energy spectra, lattice QCD methodology has been used successfully in studying multi-mesons as well as multi-baryons and nuclear systems [46][47][48][49][50][51][52][53]. These studies involve multifold challenges, namely evaluating a large number of Wick contractions, addressing the poor signal-to-noise ratio and the associated finite volume effects. Moreover, because of their exceedingly large numbers, the computational cost of Wick contractions for multi-baryon states can even exceed the cost of quark propagator computation and new algorithms are necessary to address this issue. New methods have indeed been developed, namely sink momentum projection [54], evaluating simultaneous contractions [55] and manipulation of permutation symmetry through ten-sor properties [56], which help to somewhat mitigate the computational cost. Besides Wick contractions, the issue of reliable ground state determination is also important due to the worsening of signal-to-noise ratio for multihadron systems [46]. The main challenge here is the need for a variational calculation with multi-baryon operators with good sources which can clearly separate the ground state from excited states. To overcome this issue a study has recently been performed employing both methods of point sources and distillation, and it reaches to the conclusion that more precise results for multi-hadron systems can indeed be obtained through operator smearing through distillation [31]. Another crucial challenge in searching for bound states is a detailed finite-volume amplitude analysis of the pole distribution in the scattering amplitude across the complex energy plane [33]. This requires calculations either at different physical volumes or in different momentum frames. Both of these require significant computational resources.
Working with heavy quarks can somewhat mitigate some of the aforementioned challenges associated with multi-baryon systems. For example, it is expected that the effect of chiral dynamics will be less severe in heavy multi-baryon systems and hence the signal-to-noise ratio in correlation functions may possibly be improved. Moreover, due to the presence of heavy hadrons at thresholds, a relatively large suppression of the finite volume effects on the extracted energy levels is expected for heavy dibaryons. Of course, light quarks in heavy dibaryons can still produce not-so-good signal-to-noise ratio and the presence of a light baryon at the threshold can also enhance the finite volume effects. While there have been numerous multi-baryon lattice QCD studies with light and strange quarks, there is almost no investigation on multi-baryons with heavy flavors until recently [43][44][45]. As mentioned before, in a first of its kind, deuteron-like dibaryons with heavy quarks were investigated recently and it was found that such states with charmed-bottom, strange-bottom and strange-charmed flavor combinations have large binding energies [43]. Large binding has also been found for the dibaryon which has only bottom quarks [45]. Interestingly, multiple lattice groups have predicted the large binding energies of tetraquarks with heavy quark contents [57][58][59][60]. Further, the exotics states that have been discovered in the last two decades are all have heavy flavors [1,[9][10][11][12][13][14][15][16].
Inspired by those theoretical and experimental studies, here we perform a pilot study of three-flavored heavy dibaryons using lattice QCD and report the findings. In particular, we investigate the spin-0 H-like dibaryons in the isosinglet (I = 0) and isotriplet (I = 1) channels by replacing the two strange quarks with the heavy flavors yielding dibaryons with the following quark contents: Qq 1 q 2 Qq 1 q 2 where Q ≡ c or b and q i ≡ u, d, s. These states are therefore charm and bottom quark analogues of H-dibaryon and also belong to SU (3) 27plet. Our results from this study do not support any physical bound state with deep bindings for any three-flavored H-like heavy dibaryons for isosinglet configurations. Energy levels for the isotriplet dibaryons are found to be even higher, suggesting possible scattering states. For example, for the cases of H c (cudcud), H b (budbud), and H csl (cslcsl), H bsl (bslbsl), H bsl (bclbsl); l ∈ u, d, we do not find any energy level far below their respective lowest threshold. However, for the physical isosinglet H bcs (bcsbcs) dibaryon we find an energy level consistently below its lowest elastic threshold in all lattice ensembles utilized in this work. When extrapolated to the continuum, this energy difference is found to be −29 (24) MeV. Although the large errobar prohibits us to reach to a conclusive evidence for the existence of H bcs (bcsbcs), this pilot study definitely indicates that if there is any deeply bound three-flavor dibaryon, then it has to be H bcs (bcsbcs). A further study with a large statistics and more control over systematics is now called-for to conclude about the binding nature of this state, and if a positive result is found from such a study there will be enough motivation for searching this dibaryon experimentally.
The paper is organized as below. In section II, we discuss the details of our lattice set up. Next we elaborate the relevant interpolating fields used in this study. In section III, we discuss the details of the analysis and present the results. Finally, in section IV, we provide a discussion of the results from this study and possible future outlooks.

II. LATTICE SET UP AND INTERPOLATING OPERATORS
The lattice set up that we employ for this work is similar to the one used in our previous works in Refs. [43,59,61,62], but we describe it here for completeness. The gauge ensembles utilized for this calculation are generated by the MILC collaboration with N f = 2+1+1 flavors of sea quarks using HISQ action [63]. We employ following three sets of lattice ensembles as listed in Table I. For these ensembles the strange and charm sea-quark masses were set to their physical values on these ensembles. The scale was set by the MILC collaboration using the r 1 parameter and the corresponding lattice spacings are listed in the last column of Table I. These values were also found to be consistent with the scales obtained through Wilson flow [64].
For the valence fermions, as in our earlier works [43,59,61,62], we employ a relativistic overlap action for light to charm quarks. A gauge fixed wall source is utilized to compute the valence propagators. The strange quark mass is tuned by setting the unphysical pseudoscalar mass ss to 688 MeV [65]. The charm quark mass is set to its physical value by equating the kinetic mass of the spin-averaged 1S-charmonia, 1 4 (3M J/ψ + M ηc ), to its experimental value. For the bottom quark, we use a non-relativistic QCD (NRQCD) formulation [66], where the Hamiltonian includes all the terms up to 1/(am b ) 2 and leading order term of 1/(am b ) 3 , with m b as the bare bottom quark mass. The NRQCD Hamiltonian is given by H = H 0 + ∆H, where the interaction term, ∆H, as used here, is given by, with c 1..6 as the tuned improvement coefficients [67]. The bottom quark mass is tuned by setting the kinetic mass of the spin-averaged 1S-bottomonia, 1 4 (3M Υ + M η b ), to its experimental value. The bottom quark propagators are computed following the usual NRQCD evolution of the above Hamiltonian.
As mentioned above, we use overlap fermions for the light valence quark propagators, and in Table II below, we list the range of valence quark masses and the corresponding pseudoscalar meson masses that we use for three different ensembles. We now elaborate the interpolating operators used in calculating the three-flavored heavy dibaryons in this work. These are local six-quark interpolating operators projected onto the antisymmetric and symmetric flavor representations and have already been explored in the literature in the context of searches for the H dibaryon in previous lattice calculations [59,68,69] as well as in model studies [22,23]. In the context of the present work, the charm and bottom quarks replace the strange quarks and accordingly the quantum numbers of the states. The interpolating operator is constructed as products of three diquarks as, where the alphabets with bold calligraphy, (a, b, c, d, e, f), indicate quark fields at the site ( x, t), C is the charge conjugation operator and P + ≡ (1 + γ 4 ) is used for the positive parity projection. With this notation, the antisymmetric combination is represented as [68], where the flavor Q ∈ (c, b) represents the heavy quark flavor of charm c or bottom b. Similarly the flavor q ∈ (c, s) when Q = b, and the flavor l is understood to be the light flavor. In addition to the antisymmetric channel, we also have computed the flavor symmetric channel, for which the interpolating operator is given by, One can easily notice that, by construction, these states are the heavy quark generalisation of the singlet and the 27-plet of the SU(3) flavor symmetry. By choosing the appropriate quark flavor for Q, q in Eq.(2) one obtains three possible flavor combinations namely H bcl , H bsl and H csl . We also consider the case where q = l, i.e, two degenerate light flavors with isospin symmetry. These dibaryons will be denoted as H b and H c corresponding to the case of Q = b and Q = c respectively. The non-interacting two-baryon thresholds for the above six-quark configurations related to these dibaryons involve both light and heavy single baryons, as will be discussed in the next section. Those single baryon correlators are computed using the standard interpolating operators for single baryons.

III. ANALYSIS AND RESULTS
In this section we discuss our calculations and present the results with relevant analysis. With the operators (O) so constructed, using the interpolating fields as given in Eqs. (2,3,4), we compute the single baryon and dibaryon two-point correlation functions between source (t i ) and sink (t f ) time-slices, For each case the ground state mass is obtained by fitting the respective average correlation function C O (τ ) with a single exponential at sufficiently large times (τ = t f −t i ). Coulomb gauge fixed wall sources are employed to obtain good overlap to the ground states. The single baryon correlators are utilized to evaluate the non-interacting two-baryon states. To evaluate the possible bindings of the dibaryon states it is foremost important to find first the threshold levels and particularly the lowest noninteracting two baryon energy level. Below we discuss that.

A. Threshold energy levels
We first discuss the relevant thresholds for the charm and bottom dibaryons, H c (cudcud) and H b (budbud). To represent the correlation functions C(τ ) showing their signal saturation, signal-to-noise ratio and possible fit ranges, we calculate the effective masses as defined below .
(6) Figure 1 shows the representative effective mass plots of the various possible non-interacting two-baryon correlators (C T (t)). These are obtained from the separate two baryon (B 1 and B 2 ) correlators as, The top left plot in Figure 1 shows the effective masses of possible non-interacting two-baryon threshold states corresponding to H c , namely Σ c Σ c , Λ c Λ c , and N Ξ cc , color coded in black (diamond), blue (square) and red (circle), respectively. Similarly, the top right plot shows the effective masses of Σ b Σ b , Λ b Λ b and N Ξ bb for the possible non-interacting two baryon states corresponding to the dibaryon H b . In both cases, the results are computed at the SU(3) symmetric point, which corresponds to am l = am s = 0.028 for the 48 3 × 144 lattice with the lattice spacing a = 0.0582 fm. In both cases, the lowest threshold can be seen to be that of N Ξ QQ . This is in contrast with the H dibaryon case, where the lowest threshold is that of the two non-interacting Λ baryons. In addition, the splitting between the Λ Q Λ Q and the N Ξ QQ increases as the heavy quark mass becomes heavier − from the charm quark to the bottom quark. This is also consistent with the known experimental results for heavy baryons [1], and lattice determination of single baryons [61,62,70] where experimental results are not available. Taken together experimental values of light and charmed baryons and lattice extracted values for bottom baryons, one arrives at the following numbers at the physical quark masses [1,61,62,70]: That is, whereas the lowest threshold state is ΛΛ for the light H-dibaryon, the lowest threshold states for H c and H b are N Ξ cc and N Ξ bb , respectively. We also further point out that in the literature for searches of bound heavy charm dibaryons [39], the ground state of the dibaryon is often compared incorrectly with the Λ c Λ c threshold instead of the correct threshold N Ξ cc . The Σ Q Σ Q thresholds in both cases turn out to be higher in energy, similar to the SU(3) case of the H dibaryon.
Similarly, in the bottom two plots of Fig 1 we show the representative effective masses of the various possible non-interacting thresholds for the dibaryons H bsl and H bcs . Here the possible elastic thresholds for H bsl are the two-baryons Ξ Ξ bb , ΣΩ bb , Ξ b Ξ b , Σ b Ω b , and for H bcs , those are Ω c Ω cbb , Ω cc Ω bb , Ω cb Ω cb , Ω b Ω ccb . From the ordering of states it is clear that the lowest elastic thresholds for these dibaryons are Ξ Ξ bb and Ω c Ω cbb , respectively.
We identify the lowest thresholds of other dibaryons using the experimental values of the single baryons [1] as well as lattice-determined values of them when experimental results are not available, particularly for the bottom baryons [61,62,70]. In Table III, we tabulate all the possible elastic threshold states for the dibaryons that we study in this work. The second column shows the possible non-interacting two-baryon states in ascending order of energy and the third column shows the lowest threshold state.
The non-interacting two-baryon threshold energy levels (E T ) corresponding to the two-baryon combina-  tions of Table III, for a given quark mass combinations (q 1 q 2 q 3 ), are calculated by adding the single baryon masses extracted at those quark masses, The extracted single baryon masses are found to be consistent with our previous calculations in Refs. [61,62]. As mentioned previously, the lowest two-baryon noninteracting states at light quarks and heavy quarks are different. This is illustrated in Figure 2 for the possible non-interacting two-baryon states for the dibaryon H bcl . We show the variation in terms of the ratio of pseu- doscalar meson mass at a quark mass to the η b mass. For comparison purpose we keep the lowest threshold for all cases at the same level but maintain the relative energy differences between various thresholds. While below the charm quark mass the lowest and highest threshold states are Σ c Ω cbb and Ξ cb Ξ cb , respectively, it is completely opposite at the bottom quark mass. Similar level crossings are also found for other dibaryon threshold energy levels. It is therefore crucial to identify the relative positions of threshold energy levels at a given quark mass for studying heavy dibaryons. It will be very interesting to find a phenomenological explanation behind the minimization of total energy of these threshold levels leading to the observed ordering as shown in the Figure 2 B. Calculation of energy differences The ground state energies (E H ) of the dibaryons are obtained by fitting the correlators constructed with the operators as mentioned in Eqs. (2,3,4), with a single exponential form at large times: C(τ ) ∼ e −E0τ . We then calculate the energy difference (∆E H ) between the ground state energy of a dibaryon (E 0 H ) and the elastic threshold energy level (E 0 T ) as, For each dibaryon, the whole process of calculating ∆E H is performed through a bootstrap method. We also calculate these energy differences, ∆E H , by taking the ratio of the dibaryon correlators to the twobaryon correlators, as A fitting to the above ratio-correlator can also give the energy difference ∆E H . While such a ratio-correlator offers the advantage of reducing the systematic errors, one must be careful in using it as it can possible produce a fake plateau in R(τ ) due to the saturation of different energy states at different time windows. We therefore mostly extract the ∆E H values through direct fitting of individual correlators Eq.(10) and the ratio correlators method Eq.(11) is used for consistency checks. We now present our results for these heavy dibaryons.

C. Flavor-antisymmetric Hc and Hb
We first present our results for the flavor antisymmetric H c (cudcud) and H b (budbud) dibaryons. In Figure 3 we show the representative effective masses of H c (top row) and H b (bottom row) at two light quark masses corresponding to the pseudoscalar masses 480 and 550 MeV. The lowest threshold energy levels extracted from the two baryon non-interacting states, N Ξ cc and N Ξ bb , at those quark masses, for H c and H b respectively, are shown by the red horizontal lines. It is evident that both the dibaryon energy levels, within the statistical error, overlap with their respective two-baryon noninteracting states at large times. Correspondingly ∆E, between the lowest energy state and the two-baryon elastic thresholds N Ξ QQ , is found to be consistent with zero for all the light quark masses considered here. At even lighter quark masses signal-to-noise for the dibaryon correlation functions found to be much poorer and with the large error they overlap further more with the two-baryon non-interacting energy levels. Because of the large statistical error we do not include the data below the 480 MeV pion mass in this pilot study, and hence are unable to conclude on the relative positions and nature of these dibaryon energy levels with respect to their respective thresholds at the physical quark masses. However, given the trend of the results that we find from higher to the lower quark mass it is highly unlikely that there are any deeper bound state at the lighter quark masses both for the three-flavored H c (cudcud) and H b (budbud) dibaryons.

D. Flavor-antisymmetric Hcsl, Hbsl, Hbcl and Hbcs
In Figure 4, we show the representative effective masses of H csl (cslcsl); l ⊂ u, d, where the strange and the charm quark masses are set to their physical values whereas the light quark mass (m l ) is varied. The pseudoscalar meson masses corresponding to these light quark masses are shown inside the figures. The non-interacting twobaryon state that is lowest in energy is Σ l Ω cc . At each light quark mass (m l ) we compute its energy by adding the baryon masses of Σ(lls) and Ω(ccs), and the thresholds thus obtained are represented by the horizontal lines. We find that for all ranges of m l the effective masses of H csl , at large times, are either consistent with the twobaryon non-interacting state or stay above that. The effective masses for the cases with m l < m s are found to be quite noisy with large errorbars and are not shown here. Within the large statistical error above time-slices 1.5 fm, they overlap with the thresholds indicating the absence of any energy level much below the thresholds.
Here again, it is very likely that the physical states H csu and H csd are either resonances or a scattering states or loosely bound states near the two-baryon thresholds. To identify that one needs much more statistics and scattering amplitudes analysis of these finite volume energy levels.
In Fig. 5, similarly we show the effective masses of the dibaryons H bsl where the strange and the bottom quark masses are set to their physical values and the third quark mass (m l ) is varied over a range. Here the elastic threshold state is Ξ(ssl)Ξ(bbl) and is shown by the horizontal line in each plot. At m l = m s we find the lowest energy for H bsl is consistent with the elastic threshold though its central value lies just below that. Here again, at the lighter quark masses we do not find any energy level much below the threshold that we can distinguish from the threshold within the statistics used in this pilot study. It is very likely that there is no deeply bound state of H bsl dibaryons at the physical quark masses. However, there are hints of an energy level below the threshold for each unphysically large values of m l (> m c ), and the energy splitting (|∆E|) between them increases as m l increases further.
Next we discuss the dibaryons H bcl (bclbcl); l ⊂ u, d. In Fig. 6 we show their representative effective masses, where the charm and the bottom quark masses are set to their physical values and the other quark mass (m l ) is varied over a range. As shown in Figure 2, the corresponding lowest thresholds are different at different quark masses and those are shown by the horizontal lines. Here the main observation is that unlike the previous cases, we find an energy level consistently below the threshold, particularly when m l is large. We fit the H bcl correlators with a single exponential and extract its ground state energy E H bcl at a large time. The fit ranges and fit values with one standard deviation (σ) are shown by the magenta band. We find it to be consistently lower than the threshold values for most quark masses for m s ≤ m l : E H bcl < E th bcl . The energy difference ∆E = E th bcl − E H bcl increases as the quark mass m l increases from light to the bottom quark masses.
At m l = m s , the relevant state is the physical threeflavored H bcs dibaryon which is particularly interesting. We find that at large times the lowest energy level is below but consistent with the threshold within 1.5-σ errorband. We extract the ∆E value of H bcs on three lattice spacings and show the results in Fig 7, and also tabulate that in Table IV. A continuum extrapolation with the form A + Ba 2 yields a value of ∆E bcs | cont = −29 (24)   MeV. This form of extrapolation is justified by the usage of overlap fermions which have no O(ma) error. One can also use the O(a 2 log(a)) term. However, with only three data points, inclusion of such terms is not possible for us. Though the continuum value is consistent with zero within 1.5-σ band, it is clearly noticeable that the ground state energy of H bcs is consistently below the non-interacting threshold state Ω c Ω bb on all three lattice ensembles. The extrapolated continuum result suggests that the H bcs dibaryon possibly has non-zero binding. With the given statistics and with fits without considering correlations between the thresholds and dibaryon correlators, it will not be possible to reach a definite conclusion about the nature of binding for H bcs . Moreover a detail finite volume amplitude analysis of the extracted energy levels is essential to reach a definite conclusion which is beyond the scope of this calculation. Nevertheless, findings from this pilot study are definitely encouraging for pursuing a more quantitative study in the future to achieve that goal.
At the lighter quark masses the signal-to-noise ratio found to be much poorer. Though the central values of the extracted energy levels for H bcl dibaryons are always found to be below the lowest threshold, with the given statistics they are consistent with the lowest threshold. One needs a detail finite volume amplitude analysis with more statistics to find if there is a loosely bound state, or a resonance at threshold or a scattering state for H bcu and H bcd dibaryons.  0.0582 −32 (18) In Table V we show ∆E values for H bcl at various pseudoscalar meson masses corresponding to the quark masses m s ≤ m l ≤ m b . In Fig. 8 we show the variation of ∆E bcl . Following HQET it is expected that the variation of ∆E will scale with heavy quark masses [71]. We thus plot ∆E bcl as a function of pseudoscalar masses which also scale with quark masses in the heavy quark limit. The x-axis is normalized by the mass of η b so that at the bottom quark its value becomes 1. The errorbands shown in the figure are obtained by fitting ∆E bcl as a function of x = m ps /η b with the forms ∆E(x) = A + Bx (cyan), ∆E(x) = A + Bx + Cx 2 (magenta) and ∆E(x) = A + B log(x) (grey). It is interesting to note that a logarithmic form also fits the data very well which could be phenomenologically interesting to consider for other splittings for their heavy quark mass dependence. From the above discussion it is quite apparent that when any of the quark mass in a three-flavored dibaryon H q1q2q3 becomes heavier its binding tends to increase. Therefore the strongest binding is expected for the case when m q1 = m q2 = m q3 = m b . In Fig. 9 we show the effective mass plot of this case which shows that the lowest energy state lies much below the corresponding elastic threshold. The extracted ∆E for this case is found to be −99 (8) The energy difference (∆E) between the ground state of the H bcq dibaryons and the lowest energy level of the non-interacting two-baryon states at various values of the quark mass mq in between the strange to the bottom quarks. This figure shows that the finite volume ∆E, which is related to the infinite volume binding energy of H bcq , increases with mq, leading to the result that while the physical dibaryon H bcs is more likely bound or weekly bound, other heavier unphysical dibaryons are more likely strongly bound as mq increases. Errorbands represent the fitting forms: ∆E(x) = A + Bx (cyan), ∆E(x) = A + Bx + Cx 2 (magenta) and ∆E(x) = A + Blog(x) (grey), where x is the ratio of the pseudoscalar meson mass at mq to the η b mass (x = mps/mη b ). of deuteron-like heavy dibaryons when all quark masses are set at the bottom quark mass [43], and also with the binding energy, −89 −17 +12 (12) MeV, of the single flavored heavy dibaryon at the bottom quark mass [45]. It indicates that at very heavy quark masses, bindings of the single, two and three-flavored dibaryons are similar.
The conclusions on H bcl dibaryons are the following: (I) We find a finite volume energy level below the threshold for H bcs for all the ensembles used here. A continuum extrapolation yields this energy difference from threshold ∆E bcs | cont = −29(24) MeV (Table IV and Figure 7).
(II) This energy difference increases as the quark mass m l increases and it becomes maximum at the bottom quark mass (m l = m b ) (Table V and Figure  8).
(III) There is no finite volume energy level much below the lowest thresholds for the physical H bcl ; l ⊂ u, d, dibaryons. However, there is an indication for a finite-volume energy level close to the threshold which we could not resolve and needs to be investigated further to find whether that energy level is associated with a closely bound state, or a resonance at threshold or a scattering state.

E. Flavor-symmetric three-flavored heavy dibaryons
We now discuss the results of the flavored-symmetric cases. In Fig. 10 we plot the effective masses of the flavored-symmetric H bcl dibaryon at m l = m s (top plot) and m l = m c (bottom plot), while keeping the m c and m b at their physical values, and compared them with that of the flavored-antisymmetric dibaryons. The data with black squares and blue circles represent the flavoredsymmetric and antisymmetric cases, respectively, while the red line represents the non-interacting lowest energy levels of the two baryons. A general feature that we find is that the extracted lowest energy levels for the flavorsymmetric configurations are always found to be higher than that of the antisymmetric cases. We observe that the symmetric H bcl states, at m l = m s and m l = m c , are above their respective elastic threshold energy. Most possibly they are scattering states or resonances above the threshold, and a detailed scattering amplitude analysis of the extracted energy levels is necessary to determine that. However, at the very heavy quark masses, particularly when m l = m c = m b , we observe large bindings, but always smaller than the corresponding symmetric cases. For other flavored-symmetric states, for example, for H bsl and H csl , we also observe that the lowest energy levels are always higher than their corresponding non-interacting threshold energy levels. Since we do not find any signature of a distinguishable extra energy level below the thresholds for any of the symmetric cases we will not discuss them further in this work.

F. Finite volume effects
We extracted dibaryon energy levels on Euclidean lattices at finite volume (3 fermi box extent). These cannot be directly associated with the physical states. In order  to do that one needs to perform a finite volume analysis through the scattering amplitude analysis of these finite volume energy levels [33]. However, for multihadron states with heavy quarks, it has been noted in Refs [43,59] that the finite volume corrections to infinite volume binding energy of the relevant hadronic state receives a non-trivial large suppression from the masses of the non-interacting heavy hadrons [72][73][74], as Here k ∞ is the binding momentum of the infinite volume state, E F V is the energy level computed on a cubic lattice, and (m 1 , m 2 ) are the masses of the two noninteracting hadrons with the threshold energy m 1 + m 2 .
Therefore, ∆ F V is expected to be smaller for larger values of m 1 and m 2 , that is when the threshold state consists of two heavy hadrons as in the dibaryons that we are studying. However, as shown in Table III, for most of these three-flavored dibaryons, one is a light baryon out of the two-baryons at the threshold. In particular, for H c and H b , because of the presence of nucleon at the threshold (N Ξ QQ ), the combination m 1 + m 2 may not provide a stronger suppression in comparison to the case of tetraquarks [59] and two-flavored heavy dibaryons [43]. Nevertheless, the volume suppression still is expected to be larger than that of light dibaryons. For the cases of H bcl and H bcs , volume suppression would be even larger. For the unphysical dibaryons as shown in Fig. 8, the presence of two heavy baryons will bring back larger volume suppression and one can argue that the energy levels mentioned in Table V are expected to be closer to their infinite volume limits. Nevertheless, it will be important to perform a finite volume analysis, in particular for H c and H b , as was performed in Ref. [31], where infinite volume binding energy was computed by locating the bound state pole in the scattering amplitude [33]. Such an analysis is not possible within the framework of our current set up and we would like to pursue that in future. Beside the statistical and finite volume effects, other systematic errors are also involved in this work, namely, mixed action partially quenching, discretization, scale setting, mass tuning, fit window and electromagnetism. In Ref. [43] we estimated such errors can be as large as 10 MeV. The parameters set used in this work are similar to that Ref. [43] and hence we expect similar systematic errors, particularly for the dibaryons H bcs , and those heavier in masses. For other dibaryons involving light quarks these systematics are expected to be larger and without addressing them properly it is not possible to reach a definitive conclusion for their natures of binding.

IV. SUMMARY AND DISCUSSION
In this work, we report the first lattice QCD study of three-flavored heavy dibaryons both in the flavored-symmetric and antisymmetric channels. These states are the heavy quark analogues of the much investigated H-dibaryon and are in the SU(3) 27plet of quark flavors.
From this pilot study we summarize our findings below.
First, in the flavor-symmetric channel, for the physical dibaryons H c (cudcud), H b (budbud), H csl (cslcsl), H bsl (bslbsl); l ∈ u, d, within our statistics we do not find any energy level much below their respective lowest elastic thresholds, which suggests that there is no deeply bound dibaryons in these channels. Most likely they are either loosely bound states near their respective thresholds or resonances just above the thresholds or scattering states. On the other hand, for H bcs (bcsbcs), we find an energy level below the corresponding lowest non-interacting thresh-old, Ω c Ω cbb . An extrapolation of the energy difference ∆E H bcs between the ground state of this dibaryon from the non-interacting Ω c Ω cbb , yields ∆E H bcs = −29 ± 24 MeV. Though this result on ∆E H bcs has a large error, and is consistent with zero within 1.5 standard deviation, there is a clear trend that it is consistently below the lowest threshold energy level in all the three lattice ensembles employed in this work. Since H bcs is a physical state and could be an attractive dibaryon candidate to be searched in future at high energy laboratories, this finding of the possibility of an energy level below the threshold from this pilot study is very interesting and calls for an extension with more statistics and better control over systematics. However, for these threeflavored dibaryons, when the light quark mass is set to an unphysically high value, for example for H bcl with m l > m c , while keeping the charm and bottom quark masses at their physical values, we always find an energy level much below the respective threshold energy level. That clearly indicates the possibility of strong binding of those unphysical dibaryons. Moreover, the energy difference from the respective elastic threshold becomes deeper as the quark mass m l increases, as shown in Fig. [8].
For the dibaryon H bsl , we also find the presence of an energy level below the elastic threshold at m l ∼ m b , though somewhat closer to the threshold than that of H bcl . For the dibaryon H csl , within the statistics employed in this study we do not find an energy level below its lowest threshold that can be distinguished from the lowest threshold for any value of quark masses m l employed in this work. For the flavor-symmetric channels the corresponding energy levels are observed to be always higher than those of flavor-symmetric cases, suggesting possible scattering states or resonances above the thresholds. Taken together all results, we can summarize that for the three-flavored dibaryons H q1q2q3 , there is no deeply bound state if any of the quark mass (m qi ) is below the charm quark mass. However, we find strong indications of a shallow level below the threshold for the physical H bcs state which needs to be probed further. Moreover, an energy level below the threshold always emerges when all the three quark masses become heavier than the charm quark mass, and the binding increases with the increase of quark masses.
We would also like to point out that there are different dynamics as far as binding is concerned for the threeflavored light and heavy dibaryons. That is reflected through the presence of different types of two baryons at their respective elastic thresholds. For the H-dibaryon, which is the lightest three-flavored dibaryon, the elastic threshold state is ΛΛ (with M (ΛΛ) − M (N Ξ) = −21.7 MeV). This also continues to be the case at the SU(3) point. However, for the heavy three-flavored dibaryons, H c and H b , the lowest thresholds are N Ξ cc and N Ξ bb , respectively (with M (Λ c Λ c ) − M (N Ξ cc ) = 13.05 MeV and M (Λ b Λ b ) − M (N Ξ bb ) = 158(30)) [1,61,62,75]. The presence of a doubly heavy baryon lowers the threshold for a heavy three-flavored dibaryon.
The results obtained in this work, when are taken together with the findings in doubly-heavy two-flavored deuteron-like dibaryons [43], all-heavy single-flavored dibaryons [45] and doubly heavy tetraquarks [57][58][59][60], point to an interesting dynamics of the heavy multihadron systems. A common pattern emerges that for a doubly heavy multiquark hadron, the heavier the two heavy quarks the stronger is the binding. However, the mass of other quarks (or antiquarks) towards the strong binding of these systems are quite intriguing. In the case of doubly heavy tetraquarks, the heavier the heavy quarks (or anti-quarks) and lighter the light antiquarks (or quarks), the stronger is the binding [31,59]. On the contrary, for the heavy dibaryons, binding increases when all the quarks are heavier, that is, in the presence of a light quark the binding decreases. In addition, while for various two-flavored dibaryons with two heavy quarks, the third quark can still be lighter to have an energy level below the elastic threshold, for the three-flavored case only H bcs shows such behaviour. All other physical three-flavored H-dibaryons are most likely either unbound or very weakly bound. That is, the two-flavored heavy dibaryons have stronger binding than that of threeflavored heavy dibaryons. However, when all the quarks become much heavier (m q1,q2,q3 ∼ m b ) the one-, two-and three-flavored dibaryons all exhibit similar strong binding.
The study pursued here is the first effort to investigate the three-flavored heavy dibaryons. Given the amount of theoretical and experimental efforts put into the exploration of the H dibaryon state, our motivation has been to elucidate the trend of the lattice ground state energy levels with respect to the elastic thresholds as the strange quark becomes heavier. In doing so, our hope has been to identify a possible favorable three-flavored channels in charm and/or bottom sectors which may exhibit a bound state, and guide in discovering them in future given the large experimental data being collected and to be collected for heavy hadron spectroscopy at various laboratories. This pilot study indicates that the dibaryon H b (bcsbcs) is possibly such a bound state. Considering the feasibility of discovering it in high energy experimental laboratories, it will thus be worthwhile to pursue a more detailed study in future to get a definite conclusion on the binding of this state. That can be accomplished with the variational method combined with the use of distillation method for dibaryon systems as in Ref. [31] or with the potential method [50]. Along with that, as mentioned earlier, a detailed finite-volume analysis is needed to discern the pole distribution in the scattering amplitude across the complex energy plane. We will pursue such a study in future.