Deuteron-like heavy dibaryons from Lattice QCD

TL;DR

This lattice QCD study investigates heavy-flavor dibaryons, finding some stable bound states below two-baryon thresholds and predicting their masses, thus suggesting potential new exotic nuclei involving heavy quarks.

Contribution

First lattice QCD calculation of heavy-flavor deuteron-like dibaryons, identifying stable states and providing precise mass predictions, expanding understanding of exotic nuclear matter.

Findings

Stable bound states found for $ ext{Ω}_c ext{Ω}_{cc}$, $ ext{Ω}_b ext{Ω}_{bb}$, and $ ext{Ω}_{ccb} ext{Ω}_{cbb}$.

Masses of these stable dibaryons are predicted with high precision.

Binding strength increases with heavier dibaryons.

Abstract

We report the first lattice quantum chromodynamics (QCD) study of deuteron($np$)-like dibaryons with heavy quark flavours. These include particles with following dibaryon structures and valence quark contents: $\Sigma_c\Xi_{cc} (uucucc)$, $\Omega_c\Omega_{cc} (sscscc)$, $\Sigma_b\Xi_{bb} (uububb)$, $\Omega_b\Omega_{bb} (ssbsbb)$ and $\Omega_{ccb}\Omega_{cbb} (ccbcbb)$, and with spin ($J$)-parity ($P$), $J^{P} \equiv 1^{+}$. Using a state-of-the art lattice QCD calculation, after controlling relevant systematic errors, we unambiguously find that the ground state masses of dibaryons $\Omega_c\Omega_{cc} (sscscc)$, $\Omega_b\Omega_{bb} (ssbsbb)$ and $\Omega_{ccb}\Omega_{cbb} (ccbcbb)$ are below their respective two-baryon thresholds, suggesting the presence of bound states which are stable under strong and electromagnetic interactions. We also predict their masses precisely. For dibaryons $\Sigma_c\Xi_{cc} (uucucc)$, and $\Sigma_b\Xi_{bb} (uububb)$, we could not reach to a definitive conclusion about the presence of any bound state due to large systematics associated with these states. We also find that the binding of these dibaryons becomes stronger as they become heavier in mass. This study also opens up the possibility of the existence of many other exotic nuclei, which can be formed through the fusion of heavy baryons, similar to the formation of nuclei of elements in the Periodic Table.

Figures

Tables

Table 1. Table 1: Structures of spin-1 heavy dibaryons that we study in this work.

${𝒟}_{Qq}$	Interpolating fields
${𝒟}_{bc}$	$\frac{1}{\sqrt{2}}\left({\Omega }_{ccb}{\Omega }_{bbc}-{\Omega }_{bbc}{\Omega }_{ccb}\right)$
${𝒟}_{bs}$	$\frac{1}{\sqrt{2}}\left({\Omega }_b{\Omega }_{bb}-{\Omega }_{bb}{\Omega }_b\right)$
${𝒟}_{cs}$	$\frac{1}{\sqrt{2}}\left({\Omega }_c{\Omega }_{cc}-{\Omega }_{cc}{\Omega }_c\right)$
${𝒟}_{bu}$	$\frac{1}{\sqrt{2}}\left({\Sigma }_b{\Xi }_{bb}-{\Xi }_{bb}{\Sigma }_b\right)$
${𝒟}_{cu}$	$\frac{1}{\sqrt{2}}\left({\Sigma }_c{\Xi }_{cc}-{\Xi }_{cc}{\Sigma }_c\right)$

Table 2. Table 2: Energy differences between spin-1 heavy dibaryons and their two-baryon thresholds. The third column shows the predictions of masses for the stable dibaryons. The bottom part of the table is for dibaryons with unphysical quark, q 𝑞 q , with its mass ( m q subscript 𝑚 𝑞 m_{q} ) varying in between the charm quark mass and five times the bare bottom quark mass. The errors within brackets incorporate both statistical and systematic errors added in quadrature.

Dibaryon	Energy difference	Energy difference	Mass
	from spin-1/2	from spin-3/2
	threshold [MeV]	threshold [MeV]	[MeV]
${𝒟}_{bc}$	$-91\left(12\right)$	$-52\left(13\right)$	$19105\left(21\right)$
${𝒟}_{bs}$	$-287\left(45\right)$	$-29\left(13\right)$	$16004\left(17\right)$
${𝒟}_{cs}$	$-26\left(9\right)$	$-90\left(20\right)$	$6381\left(20\right)$
${𝒟}_{bu}$	$-350\left(110\right)$	$3\left(50\right)$
${𝒟}_{cu}$	$-8\left(17\right)$	$-75\left(46\right)$
${𝒟}_{bq(m_q=1.38m_c})$		$-60\left(10\right)$
${𝒟}_{bq\left(m_q=1.72m_c\right)}$		$-87\left(8\right)$
${𝒟}_{bq\left(m_q=2.07m_c\right)}$		$-101\left(8\right)$
${𝒟}_{bq\left(m_q=m_b\right)}$		$-109\left(5\right)$
${𝒟}_{bq\left(m_q=2m_b\right)}$		$-85\left(10\right)$
${𝒟}_{bq\left(m_q=5m_b\right)}$		$-15\left(8\right)$

Table 3. Table 3: Error budget in the calculation of energy splittings for 𝒟 b ​ c subscript 𝒟 𝑏 𝑐 \mathcal{D}_{bc} dibaryon.

$Source$	Error (MeV)
Statistical	10
Discretization	5
Scale setting	4
$m_b$ tuning	3
$m_c$ tuning	3
Fit window	4
Finite size	2
Electromagnetism	3
Total	10 (stat) & 9 (syst)

Equations

$${\mathcal{D}}_{qQ}=\frac{1}{\sqrt{2}}\bigg{(}\Omega_{qqQ}(C\gamma^{j})\Omega_{QQq}-\Omega_{QQq}(C\gamma^{j})\Omega_{qqQ}\bigg{)}$$

$$\Delta E=M_{\mathcal{D}}-M_{B_1{B}_2},$$

$${\mathcal{D}}^{\prime }(\tau )=\frac{\mathcal{D}(\tau )}{B_1(\tau )\times B_2(\tau )}\to \mathcal{A}{e}^{-\Delta E\tau }+...$$

$$\Delta M_{\mathcal{D}}=[M_{\mathcal{D}}^L-n_c{\overline{1S}}_c/2-n_b{\overline{1S}}_b/2]a^{-1}.$$

$$M_{\mathcal{D}}^c=\Delta M_{\mathcal{D}}^c+n_b({\overline{1S}}_c{)}_{phys}/2+n_c({\overline{1S}}_b{)}_{phys}/2.$$

Full text

A deuteron is a bound state of two baryons, a proton and a neutron, and is made of six light valence quarks. In the early Universe, deuterons were created and their stability is responsible for the creation of other elements. Interestingly, the strong interactions between quarks, which bring stability to deuterons, also allow various other six-quark combinations leading to the possible formation of many other dibaryons. However, no such strong-interaction-stable dibaryons, though speculated about many times Jaffe (1977); Mulders et al. (1980); Balachandran et al. (1984); Sakai et al. (2000); Ikeda and Sato (2007); Bashkanov et al. (2013); Shanahan et al. (2011); Clement (2017), have been observed yet Lees et al. (2019). Using a state-of-the-art first principles calculation of lattice quantum chromodynamics (QCD), here we report, for the first time, a definite prediction of the existence of other deuteron-like spin-1 dibaryons. We also predict their masses precisely. These new subatomic particles could either be made of six heavy quarks (charm and bottom) or heavy and strange quarks. Such dibaryons are stable against strong and electromagnetic decays, but, unlike the deuteron, they can decay through weak interactions. We also find that such dibaryons become more strong-interaction-stable as they become heavier. We expect that prediction from this calculation will initiate more theoretical works on heavy dibaryons, particularly to understand their binding mechanism, and may as well aid in discovering these new subatomic particles at future experimental facilities, such as at the upgraded Large Hadron Collider and future high energy heavy ion facilities. In fact this study opens up the possibility of the existence of many other exotic nuclei, which can be formed through the fusion of heavy baryons, similar to the formation of nuclei of elements in the Periodic Table. Formation of these hadrons perhaps also enhances the possibility of a quark-level analogue of nuclear fusion as discussed recently Karliner and Rosner (2017a). Further study of these exotic states can also provide information on the strong interaction dynamics at multiple scales.

The particular dibaryons ($\mathcal{D}$)111We identify the two-flavour spin-1 dibaryons with the symbol $\mathcal{D}$ and name them as $\mathcal{D}_{q_{1}q_{2}}$, which are made of two baryons with valence quark contents $(q_{1}q_{1}q_{2})$ and $(q_{1}q_{2}q_{2})$. In this notation the deuteron is $\mathcal{D}_{ud}\equiv np(uududd)$. Such dibaryons may be called as $\mathcal{D}$-dibaryons. that we investigate are heavy quark analogues of deuteron ($np$). They have the spin-($J$)-parity ($P$) quantum numbers: $J^{P}=1^{+}$, with following dibaryon configurations: $\mathcal{D}_{cu}\equiv\Sigma_{c}\Xi_{cc}(uucucc)$, $\mathcal{D}_{cs}\equiv\Omega_{c}\Omega_{cc}(sscscc)$, $\mathcal{D}_{bu}\equiv\Sigma_{b}\Xi_{bb}(uububb)$, $\mathcal{D}_{bs}\equiv\Omega_{b}\Omega_{bb}(ssbsbb)$ and $\mathcal{D}_{bc}\equiv\Omega_{ccb}\Omega_{cbb}(ccbcbb)$. Here, $\Sigma_{q},\Xi_{qq},\Omega_{qq},\Omega_{q_{1}q_{2}q_{2}}$’s are heavy baryons with the usual nomenclature of the Particle Data Group Tanabashi et al. (2018), and $u,s,c,b$ inside brackets are various quark flavours. We find that $D_{cs}$, $D_{bs}$ and $D_{bc}$ are stable against strong and electromagnetic decays, and thus it will interesting to carry out further theoretical studies on their binding mechanisms as well as their production mechanisms. However, for $\mathcal{D}_{cu}$ and $\mathcal{D}_{bu}$, we find the ground state masses are consistent with their respective two-baryon thresholds, and therefore our results are not currently precise enough to reach a definitive conclusion on their stability. Incidentally, only recently tetra-Choi et al. (2008); Aaij et al. (2014); Bondar et al. (2012) and pentaquark Aaij et al. (2015, 2019) states have been discovered and those are made of heavy quarks. Recent model and lattice QCD studies also suggest the existence of other heavy tetraquark hadrons Karliner and Rosner (2017b); Eichten and Quigg (2017); Francis et al. (2017); Bicudo et al. (2017); Junnarkar et al. (2019); Francis et al. (2019a); Leskovec et al. (2019). It is thus natural to search for six-quark states made of heavy quarks and our predictions provide first indication for the existence of such heavy hadrons. Dibaryons have been studied through various models over the years Jaffe (1977); Balachandran et al. (1984); Mulders et al. (1980); Sakai et al. (2000); Ikeda and Sato (2007); Bashkanov et al. (2013); Shanahan et al. (2011); Clement (2017), However, it is quite crucial to have first principles lattice QCD based studies on these states to predict their masses and to understand their structures. In fact a few lattice QCD studies have already been carried out Inoue et al. (2011); Beane et al. (2011, 2012); Yamazaki et al. (2012); Gongyo et al. (2018); Francis et al. (2019b); Iritani et al. (2019). However, those are mainly focused on light quarks with spin-0 states Inoue et al. (2011); Beane et al. (2011); Gongyo et al. (2018); Francis et al. (2019b), studies of deuteron in Refs. Beane et al. (2012); Yamazaki et al. (2012) as well as studies of spin-2 states Iritani et al. (2019). This work is the first lattice study on heavy dibaryons. A lattice dibaryon calculation with heavy quarks is advantageous over the light counterparts in two ways. The two point correlators are less noisy in comparison with light dibaryons and secondly the signal-to-noise ratio of the two point correlators is far better in comparison to the light quark calculations. Both of these provide an added advantage in our calculation and help in reliable extraction of binding energies of these heavy dibaryons.

The lattice set up that we utilize here is similar to that was used in Refs. Basak et al. (2012); Mathur et al. (2018); Mathur and Padmanath (2019); Junnarkar et al. (2019). Below we elaborate it further for the sake of completeness.

A. Lattice ensembles: Three sets of dynamical 2+1+1 flavours ($u/d,s,c$) lattice ensembles, within a volume of about 3 fm, generated by MILC collaboration Bazavov et al. (2013) with HISQ fermion action Follana et al. (2007) are employed for this study. Lattice spacings, using $r_{1}$ parameter, for these ensembles are measured to be 0.1207(11) 0.0888(8) and 0.0582(5) fm, respectively Bazavov et al. (2013), which are also found to be consistent with scales obtained through Wilson flow Bazavov et al. (2016).

B. Quark actions: In the valence sector, from light to charm quarks, we utilize the overlap action which has exact chiral symmetry at finite lattice spacings Neuberger (1998a, b); Luscher (1998) and no $\mathcal{O}^{n}(ma),n=1,3,5,\cdots$ errors. A non-relativistic QCD (NRQCD) formulation Lepage et al. (1992) is adapted for the bottom quark with an $\mathcal{O}(\alpha_{s}v^{4})$ improved Hamiltonian with non-perturbatively tuned improvement coefficients Dowdall et al. (2012a). For reliable extraction of the ground states we employ a wall source and point sink.

C. Quark mass tuning: We tune both the charm and bottom quark masses using the Fermilab prescription for heavy quarks El-Khadra et al. (1997) in which we equate the lattice-extracted spin-averaged kinetic masses of the $1S$ quarkonia states with their physical values Tanabashi et al. (2018). Following Ref. Chakraborty et al. (2015) we tune the strange quark mass to its physical value.

D. Dibaryon interpolators: The dibaryon interpolating operator with spin $(J=1)$ and antisymmetric in flavours $(q,Q)\in(s,c,b)$ is constructed as:

[TABLE]

where $\Omega_{qqQ}$ and $\Omega_{QQq}$ are spin-1/2 baryons defined as, $(\Omega_{qqQ})_{\alpha}=\epsilon^{abc}q^{a}_{\alpha}(x)q^{b}_{\mu}(x)(C\gamma_{5})_{\mu\nu}Q^{c}_{\nu}(x)$ and $(\Omega_{QQq})_{\alpha}=\epsilon^{abc}Q^{a}_{\alpha}(x)q^{b}_{\mu}(x)(C\gamma_{5})_{\mu\nu}Q^{c}_{\nu}(x)$. Here Latin letters indicate color, while Greek letters indicate the spinor degree of freedom. The various deuteron analogues with appropriate flavour antisymmetry are listed in Table 1.

In a lattice QCD formulation Wilson (1974), the mass ($m$) of a particle is extracted in two steps: first by calculating the two-point Euclidean time ($\tau$) correlator ($\langle C(\tau)\rangle$) of the interpolating source and sink operators, over many gauge configurations, and then extracting the exponent ($m$) by fitting the exponential decay ($\langle C(\tau)\rangle\,\sim e^{-m\tau}$) of that correlator at large Euclidean time. Following the dibaryon structure of deuteron ($d={1\over\sqrt{2}}(pn-np)$), the spin-1 flavour-antisymmetric dibaryon interpolating fields are constructed with the appropriate spin projection of two individual spin-1/2 baryons, as shown in Table  1. We then calculate the two-point correlators with those interpolating fields and from the exponential decay of these correlators we calculate the lowest energy states, i.e., the ground state masses of each dibaryons. Masses of individual baryons are also calculated similarly (See Supplemental Material [url] for more details on calculation procedure and error analysis, which includes Refs. Lepage et al. (2002); Mathur et al. (2018); Mathur and Padmanath (2019); Brown et al. (2014); Basak et al. (2014); Bazavov et al. (2016); El-Khadra et al. (1997); Luscher and Wolff (1990); McNeile et al. (2012); Dowdall et al. (2012b); Chakraborty et al. (2015); Basak et al. (2015); Borsanyi et al. (2015)).

The next step is to find the relative energy levels of the ground state dibaryons with respect to their two-baryon thresholds. Since baryon number is a conserved quantity in the Standard Model of particle physics, these spin-1 dibaryons can have only two strong decay channels: either to two spin-1/2 baryons or to two spin-3/2 baryons. Interestingly, we find that the ordering of masses for these two combinations is different with the charm and the bottom quarks. While the sum of masses of two spin-1/2 charmed baryons ($\Sigma_{c}$ and $\Xi_{cc}$ or $\Omega_{c}$ and $\Omega_{cc}$) are smaller than that of two spin-3/2 charmed baryons ($\Omega_{qqq,q=u,s}$ and $\Omega_{ccc}$), for the bottom quark this trend is opposite, i.e., $M_{\Sigma_{b}}+M_{\Xi_{bb}}>M_{\Delta_{uuu}}+M_{\Omega_{bbb}}$, $M_{\Omega_{b}}+M_{\Omega_{bb}}>M_{\Omega_{sss}}+M_{\Omega_{bbb}}$ and $M_{\Omega_{ccb}}+M_{\Omega_{cbb}}>M_{\Omega_{ccc}}+M_{\Omega_{bbb}}$. These observations are consistent with known experimental results and lattice determination of baryon masses Brown et al. (2014); Mathur et al. (2018); Mathur and Padmanath (2019). After computing the dibaryon masses we then calculate their mass differences from both the spin-1/2 and spin-3/2 two-baryon thresholds. When the mass difference of a dibaryon from its closest threshold is negative and the finite volume effects are convincingly small, that particular dibaryon is likely to be a bound state.

We compute the aforementioned mass differences at multiple lattice spacings within our lattice set-up. Our results are shown in Figure 1. Mass differences from spin-1/2 and spin-3/2 thresholds are shown by red and green colours respectively. The upper plot is for the dibaryon $\mathcal{D}_{bc}(\Omega_{ccb}\Omega_{cbb})$, middle one is for $\mathcal{D}_{bs}(\Omega_{ssb}\Omega_{sbb})$ and the bottom one is for $\mathcal{D}_{cs}(\Omega_{ssc}\Omega_{scc})$. It should be noted that these dibaryon energy levels are already computed at their physical quark masses and therefore only require a continuum extrapolation and finite volume corrections. For the continuum extrapolation, the energy levels are computed at three lattice spacings, indicated by different marker styles, and we use a linear in $a^{2}$ ansatz as well as with an $a^{2}\,\rm{ln}(a)$ term. To be noted that use of overlap action, which has exact chiral symmetry on lattice, ensures no $\mathcal{O}^{n}(ma),n=1,3,5,\cdots$ errors, which then assures that higher order discretization errors are smaller, particularly at the finest lattice. In addition, since we calculate mass splittings, rather than masses, errors due to finite lattices are in good control. The errors shown include both statistical as well as systematic errors which are then added in quadrature to get the final errors (see Supplemental Material for details on error analysis).

From these figures it is quite apparent that the ground state masses of $\mathcal{D}_{bc},\mathcal{D}_{bs}$ and $\mathcal{D}_{cs}$ lie below their respective closest two baryon thresholds by about 4, 2 and 3 sigma errors, respectively. Next, it is natural to ask if these mass differences obtained at finite volume lattices ($L^{3}$) are indeed the physical binding energies that hold these dibaryons together from decaying to individual baryons. To sort this out one needs to study finite volume corrections of these lattice computed energy levels. Fortunately for multi-hadron systems this has already been worked out Beane et al. (2004); Davoudi and Savage (2011); Briceño et al. (2013) and the finite volume corrections, $\Delta_{FV}$, for such systems was found to be $\sim{\mathcal{O}}(e^{-k_{\infty}L})/L$, where $k_{\infty}=\sqrt{(m_{1}+m_{2})B_{\infty}}$, $m_{1},m_{2}$ being the masses of threshold states and $B_{\infty}$ the infinite volume binding energy. Since here $m_{1},m_{2}$ are masses of two baryons with multiple heavy quarks, $\Delta_{FV}$ receives a large suppression even when $B_{\infty}$ is of a few MeV. This assures that the binding energies for these dibaryons, particularly for $\mathcal{D}_{bc}$, will be very close to the extracted mass differences that we showed in Figure 1. Since we use the strange, charm and bottom quark masses already at their physical values, the extracted mass differences, after assuming negligible finite volume corrections, are indeed the physical binding energies of these dibaryons. However, to confirm the bound state nature of these dibaryons, particularly for $\mathcal{D}_{bs}$ and $\mathcal{D}_{cs}$, one needs to carry out a detail finite volume study Luscher and Wolff (1990) and we will address that in future.

For dibaryons, $\mathcal{D}_{cu}$ and $\mathcal{D}_{bu}$ we first perform chiral extrapolations with a constant plus a term linear in $m_{\pi}^{2}$. Due to the presence of light quarks, the signal to noise ratios in the correlation functions of these states are rather poor. Moreover for $\mathcal{D}_{bu}$, one decay product is the resonance state $\Delta_{uuu}$, which needs to be addressed with adequate finite volume study Luscher and Wolff (1990). With the current lattice set up, it is therefore difficult for us to make a precise conclusive statement about the stability of these two dibaryons. In future such a study can be carried out with the availability of more computing resources. However, following the example of deuteron, if we assume they are also bound, our results suggest that their binding energies will be much smaller compared to other three dibaryons mentioned above.

The final values of dibaryon masses are calculated by adding the known two-baryon threshold masses to the continuum value of binding energies that we extracted. For masses of yet-to-be discovered baryons we use their lattice values as calculated in this work and in Refs. Brown et al. (2014); Mathur et al. (2018); Mathur and Padmanath (2019). We also use subtraction method Dowdall et al. (2012b); Brown et al. (2014); Mathur et al. (2018); Mathur and Padmanath (2019) utilizing spin-average value of $1S$ charmonia and bottomonia and find results consistent with above mass estimates. The final values of dibaryon masses are shown in the third column of Table 2. In Fig. 2 we show the relative energy levels of these dibaryons with respect to their spin-1/2 and spin-3/2 two-baryon non-interacting strong-decay thresholds. It is quite apparent that the dibaryons $\mathcal{D}_{bc},\mathcal{D}_{bs}$ and $\mathcal{D}_{cs}$ lie below their respective closest thresholds (horizontal zero line) by about 4, 2 and 3 sigma errors, respectively, while for the other two, $\mathcal{D}_{bu}$, and $\mathcal{D}_{bu}$, we could not reach to a definitive conclusion due to large errors. It may be noted here that we have obtained our results from simulations with a finite number of statistical ensembles, and on space-time grids with finite lattice spacings and finite volume, and hence these are associated with both statistical and systematic errors. Our final results are obtained after carefully addressing those errors and we elaborate that in SI.

Interestingly, our results point out that the strong binding energies of these heavy dibaryons increase as they become heavier in mass. To confirm this pattern we also calculate similar dibaryons, $\mathcal{D}_{bq}$, with bottom quark at its physical value while varying the other quark mass ($m_{q}$) between the charm quark mass to five times the bare bottom quark mass (for $m_{q}>m_{b}$, we use the same tuning coefficients as for $m_{q}=m_{b}$). Of course, such quarks are unphysical and so are these dibaryons, but since they obey the same strong dynamics they can clarify the pattern of stability. In Table 2, at the bottom part, we tabulate the binding energies of these unphysical dibaryons (these results are obtained only at our finest lattice).

To depict it more clearly, in Fig. 3 we plot these binding energies for dibaryons $\mathcal{D}_{bq}(\Omega_{bbb}\Omega_{qqq})$ with $q$ varying from light to all the way to the bottom quark. It clearly shows that the spin-1 dibaryons become more stable when they are heavier. We even calculate such dibaryons with much higher quark masses (which are of course unphysical) and observe that their binding energies decrease when $m_{q}>m_{b}$ and they vanish at very large quark masses ($m_{q}\rightarrow\infty$). This suggests that the combination $\mathcal{D}_{bq}(\Omega_{bbb}\Omega_{qqq,q=b})$ has the maximum binding.

The strong-binding of these dibaryons and the pattern that it increases with the increase of quark mass is pointing towards an interesting strong dynamics in the heavy hadrons. It is to be noted that for the doubly heavy tetraquarks lattice studies Francis et al. (2017); Junnarkar et al. (2019); Francis et al. (2019a); Leskovec et al. (2019) found a complete opposite patterns where the strong-binding increases as the light quark becomes lighter. This suggests that the strong-binding mechanisms of spin-1 heavy tetraquarks and spin-1 heavy dibaryons are different and certainly this needs to be investigated further, perhaps similar to the way as investigated in Ref. Francis et al. (2019a) for the heavy tetraquarks.

The stability of these dibaryons against the strong and electromagnetic decays opens up the possibility of finding these and even more complicated higher nuclei with many heavy quarks, similar to nuclei of various elements in our periodic table. Similar to the role of deuteron in the nuclear fusion cycle for the creation of elements, these dibaryons can help the fusion processes of heavy baryons to produce nuclei with heavy quarks. Such nuclei can be studied theoretically in future with adequate computational resources and it may well be possible to discover them in future higher energy heavy ion facilities. It will also be interesting to see if the formation of such states has any implication in cosmology. Formation of such hadrons also enhances the speculation on the possibility of a quark-level analogue of nuclear fusion which was discussed recently in Ref. Karliner and Rosner (2017a). For example, formation of $\mathcal{D}_{bs}$ through fusion of $\Omega_{bb}$ and $\Omega_{b}$, as well as through fusion of $\Omega_{bbb}$ and $\Omega_{sss}$, are highly exothermic with the release of energy about 300 and 30 MeV, respectively. We also find that reactions such as $\mathcal{B}^{1/2}_{qbb}+\mathcal{B}^{1/2}_{qqb}\rightarrow\mathcal{B}^{3/2}_{bbb}+\mathcal{B}^{3/2}_{qqq},q\equiv c,s,u/d$, and $\mathcal{B}^{3/2}_{ccc}+\mathcal{B}^{3/2}_{qqq}\rightarrow\mathcal{B}^{1/2}_{qcc}+\mathcal{B}^{1/2}_{qqc},q\equiv s,u/d$, are highly exothermic (here we represent $\mathcal{B}^{J}q_{1}q_{2}q_{3}$ as a baryon with spin $J$ and quark contents $q_{1},q_{2}$ and $q_{3}$).

Since these dibaryons involve quarks with masses over a wide range, studies of their properties will be helpful to understand the strong dynamics at multiple scales. Due to the presence of multiple heavy quarks they will decay via various possible weak decay processes. For example, $\mathcal{D}_{cb}$ can decay through $b\rightarrow c$, $b\rightarrow s$ and $c\rightarrow s$ to various light baryons and multiple mesons which can interfere among themselves. Detailed analysis of these multiple ways of decay may also be an ideal place to study the hadronic-interference in weak decay processes. Furthermore, this work, particularly the strong-binding of heavier dibaryons, motivate us to study other heavy dibaryons, namely spin-0 states as well as heavy quark analogues of the H-dibaryon. Results on those dibaryons will be reported in future which, taking together with the finding here, may well be able to infer important information on the dynamics and binding mechanisms of heavy hadrons and perhaps also for the nuclei with heavy quarks.

Acknowledgements: Computations were carried out on the Cray-XC30 of ILGTI, TIFR, and on the Gaggle/Pride clusters of the Department of Theoretical Physics, TIFR. We would like to thank A. Dighe, M. Padmanath, S. Pavaskar and S. Raychaudhuri for many useful discussions. We are thankful to the MILC collaboration and in particular to S. Gottlieb for providing us with the HISQ lattices. N. M. would also like to thank A. Salve, K. Ghadiali and P. M. Kulkarni for computational support.

Supplementary Information

Below we detail the procedure that we utilize for computing energy differences between the ground state dibaryons from their non-interacting two-baryon thresholds. Afterwards we address statistical and systematic errors associated with this calculation.

Calculating mass differences and masses:    Euclidean two-point correlators for single and dibaryons are calculated using wall source smearing and from those we extract their ground state masses. We calculate the mass difference that we showed in Table 2 as,

[TABLE]

where $M_{\mathcal{D}}$ is the ground state mass obtained from the two-point dibaryon correlators, while $M_{B_{1}B_{2}}=M_{B_{1}}+M_{B_{2}}$ is the mass of the non-interacting two-baryon ($B_{1}$ and $B_{2}$) threshold. We calculate this mass difference by two methods: first by fitting the single baryons and dibaryon data sets separately and then computing their difference in Jackknife sample. In the second method we take the Jackknife ratio of the dibaryon correlator, $\mathcal{D}(\tau)$, to the two-baryon correlators, $B_{1}(\tau)\times B_{2}(\tau)$) as:

[TABLE]

A fit to the ratio correlator ($\mathcal{D}^{\prime}(\tau)$) then directly yields the mass splitting with respect to the relevant threshold ($M_{B_{1}B_{2}}$). However, one must be careful in using ratio method as it can produce spurious effects due to saturation of different particle at different time slices. Results from both methods are found to be consistent in our case.

Continuum extrapolation:    Due to exact chiral symmetry on lattice, overlap fermions have no $\mathcal{O}^{n}(ma),n=1,3,5,\cdots$ errors. Hence the first two terms that enter in chiral extrapolations are $a^{2}$ and $a^{2}\rm{ln}(a)$. Since we have only three lattice spacings it is impossible to fit them together and so we fit them separately. However, with the given precision and with data at the finest lattice spacing it is not possible to distinguish two fittings with those two forms. We also fit them in a constrained fitting similar to those in Ref. Lepage et al. (2002).

Final results on dibaryon masses:    The final values of dibaryon masses are calculated using two different subtraction methods. In the first method we add the continuum extrapolated mass differences with the two-baryon threshold masses. For example, for $\mathcal{D}_{bc}$ we add the mass difference of $52(13)$ MeV with the spin-3/2 $\Omega_{bbb}$ and $\Omega_{ccc}$ masses. For single baryons, which have not been discovered yet, we have calculated their masses and found those to be consistent with Refs. Mathur et al. (2018); Mathur and Padmanath (2019); Brown et al. (2014). In the second method, following Refs. Mathur et al. (2018); Mathur and Padmanath (2019), we first calculate subtracted masses on the lattice as

[TABLE]

Here $n_{c}$ and $n_{b}$ are the number of $c$ and $b$ valence quarks in dibaryons and $\overline{1S}_{c}$ and $\overline{1S}_{b}$ are the lattice calculated spin average $\overline{1S}$ charmonia and bottomonia masses respectively. These subtracted masses are then extrapolated to their continuum limit ($\Delta M^{c}_{\mathcal{D}}$) and finally the physical results are obtained by adding physical values of spin average masses to that as

[TABLE]

We find results from both methods are consistent with each other.

Error analysis:    Results obtained at finite volume lattices have both statistical as well as systematic errors. We address each of those below which are similar to the error analysis in Refs. Mathur et al. (2018); Mathur and Padmanath (2019); Brown et al. (2014).

Statistical error: We use wall source to obtain dibaryon correlators and it helps to obtain better ground state plateau. We find a statistical uncertainty of 10 MeV while calculating mass difference for $\mathcal{D}_{bc}$.

Fitting window error: With stable plateau we find uncertainty due to different fitting windows to be maximum of 4 MeV.

Discretization error: To obtain reliable results for these dibaryons with heavy quarks, a crucial issue is the control of discretization errors. The following three methods help us to reduce discretization errors: $i$) continuum extrapolation, as mentioned above, of the results obtained at three lattice spacings, the finest one being at 0.06 fm, $ii$) use of overlap action and $iii$) extraction of mass differences, which has less discretization error than masses. As mentioned previously, the first two terms that enters in the continuum extrapolation are $a^{2}$ and $a^{2}\,\rm{ln}(a)$. With three data points we fit these two forms separately and also together in a constrained fitting. Differences in central values at the continuum limit obtained with these different fittings are included in the discretization errors. After continuum extrapolation, we find the maximum discretization would be less than 4 MeV for the mass difference in $\mathcal{D}_{bc}$.

Scale setting error: On these set of ensembles we have also determined scales by measuring the $\Omega_{sss}$ baryon mass and those were found to be consistent with the determinations using $r_{1}$ parameter Basak et al. (2014). Measurement of scale with Wilson flow Bazavov et al. (2016) was also found to be consistent with the scale used here. The scale setting uncertainty in the mass difference (Eq. (2)) for $\mathcal{D}_{cs}$ is found to be $\sim$ 4 MeV.

Quark mass tuning error: We use the Fermilab method of heavy quarks El-Khadra et al. (1997) to tune the charm and bottom quark masses. In this method, we calculate the kinetic masses, instead of pole masses, of the spin average $\overline{1S}$ quarkonia and equate those with their experimental values. We perform this process corresponding to the central value of the scale and also with central $\pm$ error values. For each of the tuned masses, we calculate hadron masses and include the variation as the errors due to quark mass tuning. Our estimate for the mass tuning errors in the energy difference for $\mathcal{D}_{bc}$ due to the charm and bottom quarks are about 2-3 MeV.

Finite size effects: The finite volume corrections, $\Delta_{FV}$, for these multi-hadron systems were found to be $\sim{\mathcal{O}}(e^{-k_{\infty}L})/L$, where $k_{\infty}=\sqrt{(m_{1}+m_{2})B_{\infty}}$, $m_{1},m_{2}$ being the masses of threshold states and $B_{\infty}$ infinite volume binding energy. Since here $m_{1},m_{2}$ are masses of two heavy baryons, $\Delta_{FV}$ would be very small even if $B_{\infty}$ is of a few MeV. Following our works on heavy baryons in Ref. Mathur et al. (2018) and similar works in Ref. Brown et al. (2014), we include an uncertainty of 2 MeV from finite volume effects while calculating the mass difference for $\mathcal{D}_{bc}$. However, the fine volume effects for dibaryons with light-quarks would be bigger and one needs to perform a dedicated finite volume study for that Luscher and Wolff (1990).

Other sources: For dibaryons $\mathcal{D}_{cu}$ and $\mathcal{D}_{bu}$, we had to carry out chiral extrapolations and we perform that with constant plus $m_{\pi}^{2}$ terms. However, since the pion masses are not so light, particularly on the fine lattice, it is difficult to control the systematic associated with the chiral extrapolation. The dibaryon $\mathcal{D}_{bu}$ also involves its decay to $\Delta$-resonance which needs to be treated with finite volume study and that is beyond the scope of this work. We have thus included large errors associated with it coming from chiral and possible finite volume corrections which in turn limits the precision of our predictions for these baryons. In future with the availability of more computer resources this limitation could be addressed. No chiral extrapolation is necessary for $\mathcal{D}_{cs}$, $\mathcal{D}_{bs}$ and $\mathcal{D}_{bc}$ for which we have made conclusive prediction.

The unphysical sea quark mass effects are expected to be within a percent level for heavier dibaryons with no effective valence light quark content McNeile et al. (2012); Dowdall et al. (2012b); Chakraborty et al. (2015). This effect however would be larger for $\mathcal{D}_{cu}$ and $\mathcal{D}_{bu}$ dibaryons. Errors due to mixed action effects are found to be small within this lattice set up Basak et al. (2015) and are expected to vanish in the continuum limit. Errors from electromagnetism are expected to be within 3 MeV for baryons Borsanyi et al. (2015) and should be similar for these dibaryons.

Adding all these systematic errors in quadrature, we found that for the energy difference (Eq. 2), i.e., for the binding energy, in the case of $\mathcal{D}_{bc}$ is less than 8 MeV. For other dibaryons we also estimated those similarly and included those in table II. We find that systematic errors for dibaryons $\mathcal{D}_{cu}$ and $\mathcal{D}_{bu}$ are too big to reach a definitive conclusion about their bindings.

Below we summarize the error budget for $\mathcal{D}_{bc}$ dibaryon.